/[ascend]/trunk/ascend/compiler/test/test_autodiff.c
ViewVC logotype

Contents of /trunk/ascend/compiler/test/test_autodiff.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2336 - (show annotations) (download) (as text)
Sun Dec 26 03:44:24 2010 UTC (11 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 25236 byte(s)
Suppress some debug output.
In mtx_vector, initialise some members of mtx_vector class, because otherwise they trigger valgrind errors.
1 /* ASCEND modelling environment
2 Copyright (C) 2007 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 */
19 /**
20 @file
21 Unit test functions reverse ad & second derivatives
22 Created by: Mahesh Narayanamurthi
23 Revised by: Ben Allan - Rev.1
24 Mahesh Narayanamurthi - Rev. 2
25 */
26 #include <CUnit/CUnit.h>
27
28 #include <string.h>
29 #include <stdlib.h>
30 #include <stdio.h>
31
32 #include <ascend/general/env.h>
33 #include <ascend/general/ospath.h>
34 #include <ascend/general/list.h>
35 #include <ascend/general/ltmatrix.h>
36
37 #include <ascend/general/platform.h>
38 #include <ascend/utilities/ascEnvVar.h>
39 #include <ascend/utilities/error.h>
40
41 #include <ascend/compiler/ascCompiler.h>
42 #include <ascend/compiler/module.h>
43 #include <ascend/compiler/parser.h>
44 #include <ascend/compiler/library.h>
45 #include <ascend/compiler/symtab.h>
46 #include <ascend/compiler/simlist.h>
47 #include <ascend/compiler/instquery.h>
48 #include <ascend/compiler/parentchild.h>
49 #include <ascend/compiler/atomvalue.h>
50 #include <ascend/compiler/relation_io.h>
51 #include <ascend/compiler/reverse_ad.h>
52 #include <ascend/compiler/relation_util.h>
53 #include <ascend/compiler/mathinst.h>
54 #include <ascend/compiler/watchpt.h>
55 #include <ascend/compiler/initialize.h>
56 #include <ascend/compiler/name.h>
57 #include <ascend/compiler/visitinst.h>
58 #include <ascend/compiler/functype.h>
59 #include <ascend/compiler/safe.h>
60 #include <ascend/compiler/qlfdid.h>
61 #include <ascend/compiler/instance_io.h>
62 // FIXME include anything else?
63
64 #include <test/assertimpl.h>
65 #include <test/test_globals.h>
66
67 #define RAD_TOL 1e-05
68
69 extern struct FilePath* ASC_TEST_DIR;
70
71 struct FileList
72 {
73 FILE * yacas; /* Meant for YACAS */
74 FILE * safeder; /* Meant for Safe Der */
75 FILE * nonsafeder; /* Meant for Non-Safe Der */
76 };
77
78 struct TestFileList
79 {
80 FILE * FirstDerYacas;
81 FILE * SecondDerYacas;
82 };
83
84 struct DiffTestData
85 {
86 FILE * outfile; /* Meant for Logging Debug Info */
87 FILE * varfile; /* Meant for recording the Variable Values */
88 struct FileList FirstDer;
89 struct FileList SecondDer;
90 int use_yacas;
91 FILE * first_yacas;
92 FILE * second_yacas;
93 struct Instance *root;
94 const char *infile;
95 int d0errors;
96 int d1errors;
97 int d1errors_yacas;
98 int d2errors_yacas;
99 };
100
101
102
103 #define INITDTD(d, log, varfile, FirstDer, SecondDer, use_yacas, YacasInFirst, YacasInSecond ,instroot, inputname) \
104 do { \
105 d.outfile = log; \
106 d.varfile = varfile; \
107 d.FirstDer = FirstDer;\
108 d.SecondDer = SecondDer;\
109 d.use_yacas = use_yacas; \
110 d.first_yacas = YacasInFirst;\
111 d.second_yacas = YacasInSecond;\
112 d.root = instroot; \
113 d.infile = inputname; \
114 d.d0errors= 0; \
115 d.d1errors= 0; \
116 d.d1errors_yacas= 0; \
117 d.d2errors_yacas= 0; \
118 } while (0)
119
120 #define LOG(d, ...) do { \
121 if ( d != NULL && d-> outfile != NULL) { \
122 fprintf(d->outfile, __VA_ARGS__); \
123 } \
124 } while ( 0 )
125
126 /** Temporary Declarations */
127 static void AutomateDiffTest(struct Instance *inst, VOIDPTR ptr);
128
129
130 static void test_autodiff(void){
131 #define OUTENV "/../ascend/compiler/test/LOG.html"
132 #define VARFILE "/../ascend/compiler/test/Vars.txt"
133
134 #define SAFEDER_2ND "/../ascend/compiler/test/Safes2nd.txt"
135 #define NONSAFEDER_2ND "/../ascend/compiler/test/Nonsafes2nd.txt"
136 #define YACAS_2ND "/../ascend/compiler/test/Yacas2nd.txt"
137
138 #define SAFEDER_1ST "/../ascend/compiler/test/Safes1st.txt"
139 #define NONSAFEDER_1ST "/../ascend/compiler/test/Nonsafes1st.txt"
140 #define YACAS_1ST "/../ascend/compiler/test/Yacas1st.txt"
141
142 #define YACAS_IN_1ST "/../ascend/compiler/test/FirstDeriv.txt"
143 #define YACAS_IN_2ND "/../ascend/compiler/test/SecondDeriv.txt"
144
145 #define CASEFILE "test/reverse_ad/allmodels.a4c"
146
147 #define USE_YACAS_ENV "ASC_YACAS_GEN"
148
149 int status;
150 int use_yacas = 0;
151
152 struct module_t *m;
153
154 struct Name *name;
155 enum Proc_enum pe;
156 struct Instance *root;
157
158 struct DiffTestData data;
159
160 char FILE_PATH[PATH_MAX];
161
162 struct FilePath* out_osp;
163 struct FilePath* varfile_osp;
164
165 struct FilePath* safe_osp_2nd;
166 struct FilePath* yacas_osp_2nd;
167 struct FilePath* nonsafe_osp_2nd;
168
169 struct FilePath* safe_osp_1st;
170 struct FilePath* yacas_osp_1st;
171 struct FilePath* nonsafe_osp_1st;
172
173 struct FilePath* first_yacas_osp;
174 struct FilePath* second_yacas_osp;
175
176
177 FILE * out = NULL;
178 FILE * varfile = NULL;
179
180 FILE * first_yacas = NULL;
181 FILE * second_yacas = NULL;
182
183 struct FileList FirstDer;
184 struct FileList SecondDer;
185
186 FirstDer.yacas = FirstDer.safeder = FirstDer.nonsafeder = NULL;
187 SecondDer.yacas = SecondDer.safeder = SecondDer.nonsafeder = NULL;
188
189 Asc_CompilerInit(1);
190 Asc_PutEnv(ASC_ENV_LIBRARY "=models");
191
192
193 // FIXME Use the environment variables here
194
195
196 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
197 strncat(FILE_PATH,OUTENV,PATH_MAX-strlen(FILE_PATH));
198 out_osp = ospath_new(FILE_PATH);
199 out = ospath_fopen(out_osp,"w");
200 CU_ASSERT_PTR_NOT_NULL_FATAL(out);
201
202 /** @TODO Open the following streams only if Environment Variable is set */
203 if(getenv(USE_YACAS_ENV) != NULL ){
204 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
205 strncat(FILE_PATH,VARFILE,PATH_MAX-strlen(FILE_PATH));
206 varfile_osp = ospath_new(FILE_PATH);
207 varfile = ospath_fopen(varfile_osp,"w");
208 CU_ASSERT_PTR_NOT_NULL_FATAL(varfile);
209
210 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
211 strncat(FILE_PATH,YACAS_2ND,PATH_MAX-strlen(FILE_PATH));
212 yacas_osp_2nd = ospath_new(FILE_PATH);
213 SecondDer.yacas = ospath_fopen(yacas_osp_2nd,"w");
214 CU_ASSERT_PTR_NOT_NULL_FATAL(SecondDer.yacas);
215
216 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
217 strncat(FILE_PATH,SAFEDER_2ND,PATH_MAX-strlen(FILE_PATH));
218 safe_osp_2nd = ospath_new(FILE_PATH);
219 SecondDer.safeder = ospath_fopen(safe_osp_2nd,"w");
220 CU_ASSERT_PTR_NOT_NULL_FATAL(SecondDer.safeder);
221
222
223 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
224 strncat(FILE_PATH,NONSAFEDER_2ND,PATH_MAX-strlen(FILE_PATH));
225 nonsafe_osp_2nd = ospath_new(FILE_PATH);
226 SecondDer.nonsafeder = ospath_fopen(nonsafe_osp_2nd,"w");
227 CU_ASSERT_PTR_NOT_NULL_FATAL(SecondDer.nonsafeder);
228
229 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
230 strncat(FILE_PATH,YACAS_1ST,PATH_MAX-strlen(FILE_PATH));
231 yacas_osp_1st = ospath_new(FILE_PATH);
232 FirstDer.yacas = ospath_fopen(yacas_osp_1st,"w");
233 CU_ASSERT_PTR_NOT_NULL_FATAL(FirstDer.yacas);
234
235 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
236 strncat(FILE_PATH,SAFEDER_1ST,PATH_MAX-strlen(FILE_PATH));
237 safe_osp_1st = ospath_new(FILE_PATH);
238 FirstDer.safeder = ospath_fopen(safe_osp_1st,"w");
239 CU_ASSERT_PTR_NOT_NULL_FATAL(FirstDer.safeder);
240
241
242 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
243 strncat(FILE_PATH,NONSAFEDER_1ST,PATH_MAX-strlen(FILE_PATH));
244 nonsafe_osp_1st = ospath_new(FILE_PATH);
245 FirstDer.nonsafeder = ospath_fopen(nonsafe_osp_1st,"w");
246 CU_ASSERT_PTR_NOT_NULL_FATAL(FirstDer.nonsafeder);
247
248 use_yacas=1;
249 }
250 else{
251 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
252 strncat(FILE_PATH,YACAS_IN_1ST,PATH_MAX-strlen(FILE_PATH));
253 first_yacas_osp = ospath_new(FILE_PATH);
254 first_yacas = ospath_fopen(first_yacas_osp,"r");
255 CU_ASSERT_PTR_NOT_NULL_FATAL(first_yacas);
256
257 ospath_strncpy(ASC_TEST_DIR,FILE_PATH,PATH_MAX);
258 strncat(FILE_PATH,YACAS_IN_2ND,PATH_MAX-strlen(FILE_PATH));
259 second_yacas_osp = ospath_new(FILE_PATH);
260 second_yacas = ospath_fopen(second_yacas_osp,"r");
261 CU_ASSERT_PTR_NOT_NULL_FATAL(second_yacas);
262 }
263
264 /* load the file */
265 m = Asc_OpenModule(CASEFILE,&status);
266 CU_ASSERT(status == 0);
267
268 /* parse it */
269 CU_ASSERT(0 == zz_parse());
270
271 /* find the model */
272 symchar * type = AddSymbol("allmodels");
273 CU_ASSERT(FindType(type)!=NULL);
274
275 /* instantiate it */
276 struct Instance *sim = SimsCreateInstance(type, AddSymbol("sim1"), e_normal, NULL);
277 CU_ASSERT_FATAL(sim!=NULL);
278
279 /*root of Simulation object*/
280 root = GetSimulationRoot(sim);
281
282 /** Call on_load */
283 error_reporter_tree_start();
284
285 name = CreateIdName(AddSymbol("on_load"));
286 pe = Initialize(root,name,"sim1",ASCERR,0, NULL, NULL);
287 int haserror=0;
288 if(error_reporter_tree_has_error()){
289 haserror=1;
290 }
291
292 error_reporter_tree_end();
293
294 if(pe == Proc_all_ok){
295 if(haserror){
296 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Method 'on_load' returned all_ok status but has error(s).");
297 }else{
298 ERROR_REPORTER_NOLINE(ASC_USER_SUCCESS,"Method 'on_load' returned 'all_ok' and output no errors.\n");
299 }
300 CONSOLE_DEBUG("Method 'on_load' : COMPLETED OK");
301 }else{
302 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Method 'on_load' returned the following error: %d",pe);
303 }
304
305 /* check for vars and rels */
306
307 INITDTD(data, out, varfile, FirstDer, SecondDer,use_yacas, first_yacas, second_yacas, root, CASEFILE);
308 if(data.outfile!=NULL){
309 fprintf(data.outfile,"<html><head><title>LOG File of Test-Suite for Automatic Differentiation</title></head><body style='font-size: 9pt;font-family: serif;'>\n");
310 fprintf(data.outfile,"<center><h3>LOG File of Test-Suite for Automatic Differentiation</h3></center></br></br>\n");
311 fprintf(data.outfile,"<center><h4><u>Legend</u></h4><b>RAD_TOL = %21.17g</br>\n",RAD_TOL);
312 fprintf(data.outfile,"Difference(a,b) = a - b</br>\n");
313 fprintf(data.outfile,"Error(a,b) = (a - b)/a </br>\n");
314 fprintf(data.outfile,"<font color='gold'>gold</font> = Residual Error </br>\n");
315 fprintf(data.outfile,"<font color='indigo'>indigo</font>= Gradient Error </br>\n");
316 fprintf(data.outfile,"<font color='purple'>purple</font>= Gradient Mismatch </br>\n");
317 fprintf(data.outfile,"<font color='red'>red</font>= Second Partials Mismatch</b> </center></br></br>\n");
318 }
319
320 VisitInstanceTreeTwo(root, AutomateDiffTest, 0, 0, &data);
321 sim_destroy(sim);
322 Asc_CompilerDestroy();
323
324 if(data.outfile!=NULL){
325 fprintf(data.outfile,"</br><center><u><h4>SUMMARY:</h4></u></br>\n");
326 fprintf(data.outfile,"<b>No. of Residual Errors</b> =<b> %d </b></br>\n",data.d0errors);
327 fprintf(data.outfile,"<b>No. of First RAD Errors</b> =<b> %d </b></br>\n",data.d1errors);
328 fprintf(data.outfile,"<b>No. of Total Errors</b> =<b> %d</b></br></br>\n",data.d0errors+data.d1errors);
329 fprintf(data.outfile,"<b>No. of First YACAS Mismatches</b> =<b> %d</b></br>\n",data.d1errors_yacas);
330 fprintf(data.outfile,"<b>No. of Second YACAS Mismatches</b> =<b> %d</b></br>\n",data.d2errors_yacas);
331 fprintf(data.outfile,"<b>No. of Total Mismatches</b> =<b> %d</b></center></br>\n</body>\n</html>"
332 ,data.d1errors_yacas + data.d2errors_yacas
333 );
334 }
335
336
337 if (out_osp!=NULL){
338 ospath_free(out_osp);
339
340 if (out != NULL)
341 {
342 fclose(out);
343 data.outfile = NULL;
344 }
345 }
346
347 if (use_yacas){
348 if (varfile_osp!=NULL){
349 ospath_free(varfile_osp);
350
351 if (varfile !=NULL)
352 {
353 fclose(varfile);
354 data.varfile = NULL;
355
356 }
357 }
358
359 if (yacas_osp_2nd!=NULL){
360 ospath_free(yacas_osp_2nd);
361
362 if (SecondDer.yacas != NULL)
363 {
364 fclose(SecondDer.yacas);
365 data.SecondDer.yacas = NULL;
366 }
367 }
368
369 if (safe_osp_2nd!=NULL){
370 ospath_free(safe_osp_2nd);
371
372 if (SecondDer.safeder !=NULL)
373 {
374 fclose(SecondDer.safeder);
375 data.SecondDer.safeder = NULL;
376 }
377 }
378
379 if (nonsafe_osp_2nd!=NULL){
380 ospath_free(nonsafe_osp_2nd);
381
382 if (SecondDer.nonsafeder !=NULL)
383 {
384 fclose(SecondDer.nonsafeder);
385 data.SecondDer.nonsafeder = NULL;
386 }
387 }
388
389
390 if (yacas_osp_1st!=NULL){
391 ospath_free(yacas_osp_1st);
392
393 if (FirstDer.yacas != NULL)
394 {
395 fclose(FirstDer.yacas);
396 data.FirstDer.yacas = NULL;
397 }
398 }
399
400 if (safe_osp_1st!=NULL){
401 ospath_free(safe_osp_1st);
402
403
404 if (FirstDer.safeder !=NULL)
405 {
406 fclose(FirstDer.safeder);
407 data.FirstDer.safeder = NULL;
408 }
409 }
410
411 if (nonsafe_osp_1st!=NULL){
412 ospath_free(nonsafe_osp_1st);
413
414 if (FirstDer.nonsafeder !=NULL)
415 {
416 fclose(FirstDer.nonsafeder);
417 data.FirstDer.nonsafeder = NULL;
418
419 }
420 }
421 }
422 else{
423 if (first_yacas_osp!=NULL){
424 ospath_free(first_yacas_osp);
425
426
427 if (first_yacas !=NULL)
428 {
429 fclose(first_yacas);
430 data.first_yacas = NULL;
431 }
432 }
433
434 if (second_yacas_osp!=NULL){
435 ospath_free(second_yacas_osp);
436
437 if (second_yacas !=NULL)
438 {
439 fclose(second_yacas);
440 data.second_yacas = NULL;
441
442 }
443 }
444 }
445 CU_ASSERT( 0 == (data.d0errors + data.d1errors) );
446 CONSOLE_DEBUG("For non-fatal errors refer ascend/compiler/test/LOG.html");
447 }
448
449
450 /**
451 Files into which information is extracted
452 outfile - comment
453 safeder - all safe 2nd derivative values
454 nonsafeder - all non-safe 2nd derivative values
455 outfile - End of a Relation
456 yacas - O/P to YACAS
457
458 */
459 static void AutomateDiffTest(struct Instance *inst, VOIDPTR ptr){
460 #define RETURN if (rname != NULL) { ASC_FREE(rname); } return
461 double residual_rev,residual_fwd;
462 double *gradients_rev,*gradients_fwd,*deriv_2nd;
463
464 float yacas_first_der = 0.0; /* compiler complains type mismatch when using */
465 float yacas_second_der = 0.0; /*double to read using fscanf(File*,"%21.17g"...) */
466
467
468 double err;
469
470
471 unsigned long num_var;
472 unsigned long i,j;
473
474 enum Expr_enum reltype;
475
476 struct relation *r;
477
478 int32 status;
479
480 char *rname = NULL;
481 char *infix_rel = NULL;
482 char *varname;
483 char buf[20];
484
485 struct RXNameData myrd = {"x",NULL,""};
486 struct DiffTestData *data = (struct DiffTestData*) ptr;
487 struct Instance *var_inst;
488
489 if (inst==NULL || InstanceKind(inst)!=REL_INST){
490 return;
491 }
492 if (data != NULL && data->outfile != NULL) {
493 rname = WriteInstanceNameString(inst, data->root);
494 }
495
496 CONSOLE_DEBUG("<<<<<<<<<<<<<<<<< The relation that follows: %s >>>>>>>>>>>>>>>>>>>>",rname);
497
498 r = (struct relation *)GetInstanceRelation(inst, &reltype);
499
500 if (r == NULL) {
501 LOG(data, "<font color='orange'><b>! skipping instance with null struct relation:</b></font> %s</br>\n", rname);
502 RETURN;
503 }
504
505 if (reltype != e_token) {
506 LOG(data, "<font color='orange'><b>! skipping non-token relation</b></font> %s</br>\n", rname);
507 RETURN;
508 }
509
510 LOG(data, "</br><font color='green'><b> Evaluating token relation </b></font> <b> %s </b> </br>\n", rname);
511
512 infix_rel = WriteRelationString(inst,
513 data->root,
514 (WRSNameFunc)RelationVarXName,
515 NULL,
516 relio_yacas,
517 NULL);
518
519 LOG(data,"</br></br><b>Relation: %s </br>\n",infix_rel); // Do I need to escape this??
520
521 num_var = NumberVariables(r); // FIXME or should we use rel_n_incidences
522
523 for(i=0; i<num_var; i++) {
524 varname = RelationVarXName(r,i+1,&myrd);
525 if (varname!=NULL){
526 LOG(data,"%s:=",varname);
527 var_inst = RelationVariable(r,i+1);
528 LOG(data,"%21.17g</br>\n",RealAtomValue(var_inst));
529 }
530 }
531
532 LOG(data,"</b>\n");
533
534 gradients_rev = ASC_NEW_ARRAY(double,num_var); // or rel_n_incidences
535 CU_ASSERT_PTR_NOT_NULL_FATAL(gradients_rev);
536
537 gradients_fwd = ASC_NEW_ARRAY(double,num_var); // or rel_n_incidences
538 CU_ASSERT_PTR_NOT_NULL_FATAL(gradients_fwd);
539
540 deriv_2nd = ASC_NEW_ARRAY(double,num_var); // or rel_n_incidences
541 CU_ASSERT_PTR_NOT_NULL_FATAL(deriv_2nd);
542
543 /** @todo log the infix form of the relation and associated variables values */ /*FIXME*/
544 if(data->use_yacas){
545 if(data->FirstDer.yacas!=NULL && data->SecondDer.yacas!=NULL && data->varfile!=NULL){
546 /** Print Variables Values First */
547 for(i=0; i<num_var; i++) {
548 varname = RelationVarXName(r,i+1,&myrd);
549 if (varname!=NULL){
550 fprintf(data->varfile,"%s:=",varname);
551 var_inst = RelationVariable(r,i+1);
552 fprintf(data->varfile,"%21.17g\n",RealAtomValue(var_inst));
553 }
554 }
555 fprintf(data->varfile,"@ Relation: %s follows\n",rname);
556
557 if (infix_rel!=NULL){
558 for(i=0; i<num_var;i++){
559 varname = RelationVarXName(r,i+1,&myrd);
560 strncpy(buf,varname,20);
561
562 // Generating Output for First Derivative Yacas Input file
563 fprintf(data->FirstDer.yacas,"ToStdout() [Echo({D(%s) ",buf);
564 fprintf(data->FirstDer.yacas,"%s});];\n",infix_rel);
565
566
567 // Generating Output for Second Derivative Yacas Input file
568 for (j=0;j<num_var;j++){
569 varname = RelationVarXName(r,j+1,&myrd);
570 fprintf(data->SecondDer.yacas,"ToStdout() [Echo({D(%s) ",buf);
571 fprintf(data->SecondDer.yacas,"D(%s) ",varname);
572 fprintf(data->SecondDer.yacas,"%s});];\n",infix_rel);
573 }
574 }
575 }
576 else{
577 ERROR_REPORTER_HERE(ASC_PROG_FATAL,"Infix Relation is NULL");
578 }
579 }
580 }
581 /** Testing non-safe routines */
582
583 /* we need to sigfpe trap this code or use the safe versions. */
584 RelationCalcResidGradRev(inst,&residual_rev,gradients_rev);
585 RelationCalcResidGrad(inst,&residual_fwd,gradients_fwd);
586
587 LOG(data,"</br> <b> Table of Values for Residuals </b> </br>\n");
588 LOG(data,"\n<table BORDER>\n");
589 LOG(data,"<tr><td>ASCEND(NONSAFE,REV)</td><td>ASCEND(NONSAFE,FWD)</td><td>Error</td></tr>\n");
590 CU_ASSERT(fabs(residual_rev - residual_fwd) <= RAD_TOL);
591 if ( fabs(residual_rev - residual_fwd) > RAD_TOL) {
592 data->d0errors ++;
593 LOG(data,"<tr bgcolor='yellow'><td><font color='red'>%21.17g</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%.4g</font></td></tr>\n", residual_rev,residual_fwd, fabs(residual_rev - residual_fwd));
594 }
595 else{
596 LOG(data,"<tr><td>%21.17g</td><td>%21.17g</td><td>%.4g</td></tr>\n",residual_rev,residual_fwd,0.0);
597 }
598 LOG(data,"</table>\n");
599
600
601
602 LOG(data,"</br> <b> Table of Values for First Derivatives </b> </br>\n");
603 LOG(data,"\n<table BORDER>\n");
604 for(i=0; i<num_var; i++) {
605
606 if(data->use_yacas && data->FirstDer.nonsafeder!=NULL){
607 /* Recording Reverse AD Non-Safe Derivatives to file */
608 fprintf(data->FirstDer.nonsafeder,"%21.17g\n",gradients_rev[i]);
609 }
610 else if(data->first_yacas!=NULL){
611 /* Benchmarking Non-Safe Gradient Errors against Yacas */
612 if(!feof(data->first_yacas)){
613 fscanf(data->first_yacas,"%g\n",&yacas_first_der);
614
615 if(yacas_first_der!=0.0){
616 err = fabs((double)yacas_first_der - gradients_rev[i]) / (double)yacas_first_der;
617 }else{
618 err = gradients_rev[i];
619 }
620 err = fabs(err);
621 LOG(data,"<tr><td>Column</td><td>ASCEND(NONSAFE,REV)</td><td>YACAS</td><td>Percentage Mismatch</td></tr>\n");
622 CU_ASSERT(err <= RAD_TOL);
623 if (err > RAD_TOL) {
624 data->d1errors_yacas ++;
625 LOG(data,"<tr bgcolor='yellow'><td><font color='red'>%lu</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%.4g</font></td></tr>\n", i,gradients_rev[i],yacas_first_der, err*100);
626 }
627 else{
628 LOG(data,"<tr><td>%lu</td><td>%21.17g</td><td>%21.17g</td><td>%.4g</td></tr>\n", i,gradients_rev[i],yacas_first_der,0.0);
629 }
630 }
631 }
632
633 if (gradients_fwd[i] != 0.0) {
634 err = (gradients_fwd[i] - gradients_rev[i]) / gradients_fwd[i];
635 } else {
636 err = gradients_rev[i]; // These are totally different quantities should I make err as 1?
637 }
638 err = fabs(err);
639 LOG(data,"<tr><td>Column</td><td>ASCEND(NONSAFE,REV)</td><td>ASCEND(NONSAFE,FWD)</td><td>Percentage Mismatch</td></tr>\n");
640 CU_ASSERT(err <= RAD_TOL);
641 if (err > RAD_TOL) {
642 data->d1errors ++;
643 LOG(data,"<tr bgcolor='yellow'><td><font color='red'>%lu</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%.4g</font></td></tr>\n", i,gradients_rev[i],gradients_fwd[i], err*100);
644 }
645 else {
646 LOG(data,"<tr><td>%lu</td><td>%21.17g</td><td>%21.17g</td><td>%21.17g</td></tr>\n", i,gradients_rev[i],gradients_fwd[i],0.0);
647 }
648 }
649 LOG(data,"</table>\n");
650
651 /*Non Safe Second Derivative Calculations*/
652 if(data->use_yacas && data->SecondDer.nonsafeder!=NULL){
653 fprintf(data->SecondDer.nonsafeder,"@ Relation: %s Follows\n",rname);
654 }
655
656
657
658 LOG(data,"</br> <b> Table of Values for Second Derivatives </b> </br>\n");
659 LOG(data,"\n<table BORDER>\n");
660 LOG(data,"<tr><td>Row</td><td>Column</td><td>ASCEND (NON-SAFE)</td><td>Yacas</td><td>Percentage Mismatch</td></tr>\n");
661 for(i=0; i<num_var; i++){
662 RelationCalcSecondDeriv(inst,deriv_2nd,i);
663 if(data->use_yacas && data->SecondDer.nonsafeder!=NULL){
664 /** @todo log calculated values and indiceshere.*/ /*FIXME*/
665 for(j=0; j<num_var; j++){
666 fprintf(data->SecondDer.nonsafeder,"%21.17g\n",deriv_2nd[j]); /*FIXME*/
667 }
668 }
669 else if (data->second_yacas!=NULL){
670 /* Benchmarking Non-Safe Second Derivative Errors against Yacas */
671 for(j=0; j<num_var; j++){
672 if(!feof(data->second_yacas)){
673 fscanf(data->second_yacas,"%g\n",&yacas_second_der);
674
675 if(yacas_second_der!=0.0){
676 err = ((double)yacas_second_der - deriv_2nd[j]) / (double)yacas_second_der; /* todo err scaling */
677 }else{
678 err = deriv_2nd[j];
679 }
680 err = fabs(err);
681 CU_ASSERT(err <= RAD_TOL);
682 if (err > RAD_TOL) {
683 data->d2errors_yacas ++;
684 LOG(data,"<tr bgcolor='yellow'><td><font color='red'>%lu</font></td><td><font color='red'>%lu</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%.4g</font></td></tr>\n", i,j,deriv_2nd[j],yacas_second_der,err*100);
685 }
686 else{
687 LOG(data,"<tr><td>%lu</td><td>%lu</td><td>%21.17g</td><td>%21.17g</td><td>%21.17g</td></tr>\n", i,j,deriv_2nd[j],yacas_second_der,0.0);
688 }
689 }
690 }
691 }
692 }
693
694 LOG(data,"</table>\n");
695
696 /** Testing safe routines */
697
698 status = (int32) RelationCalcResidGradRevSafe(inst,&residual_rev,gradients_rev);
699 safe_error_to_stderr( (enum safe_err *)&status );
700
701 status = RelationCalcResidGradSafe(inst,&residual_fwd,gradients_fwd);
702 safe_error_to_stderr( (enum safe_err *)&status );
703
704
705
706 LOG(data,"</br> <b> Table of Values for Residuals </b> </br>\n");
707 LOG(data,"\n<table BORDER>\n");
708 LOG(data,"<tr><td>ASCEND(SAFE,REV)</td><td>ASCEND(SAFE,FWD)</td><td>Error</td></tr>\n");
709 CU_ASSERT(fabs(residual_rev - residual_fwd) <= RAD_TOL);
710 if ( fabs(residual_rev - residual_fwd) > RAD_TOL) {
711 data->d0errors ++;
712 LOG(data,"<tr bgcolor='yellow'><td><font color='red'>%21.17g</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%.4g</font></td></tr>\n", residual_rev,residual_fwd, fabs(residual_rev - residual_fwd));
713 }
714 else{
715 LOG(data,"<tr><td>%21.17g</td><td>%21.17g</td><td>%.4g</td></tr>\n", residual_rev,residual_fwd,0.0);
716 }
717 LOG(data,"</table>\n");
718
719
720
721 LOG(data,"</br> <b> Table of Values for First Derivatives </b> </br>\n");
722 LOG(data,"\n<table BORDER>\n");
723 LOG(data,"<tr><td>Column</td><td>ASCEND(SAFE,REV)</td><td>ASCEND(SAFE,FWD)</td><td>Percentage Mismatch</td></tr>\n");
724 for(i=0;i<num_var;i++) {
725
726 if(data->use_yacas && data->FirstDer.safeder!=NULL){
727 /* Recording Reverse AD Safe Derivatives to file */
728 fprintf(data->FirstDer.safeder,"%21.17g\n",gradients_rev[i]);
729 }
730
731
732 if (gradients_fwd[i] != 0.0) {
733 err = ( gradients_fwd[i] - gradients_rev[i] ) / gradients_fwd[i];
734 } else {
735 err = gradients_rev[i];
736 }
737 err = fabs(err);
738 CU_ASSERT(err <= RAD_TOL);
739 if (err > RAD_TOL) {
740 data->d1errors ++;
741 LOG(data,"<tr bgcolor='yellow'><td><font color='red'>%lu</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%21.17g</font></td><td><font color='red'>%.4g</font></td></tr>\n", i,gradients_rev[i],gradients_fwd[i], err*100);
742 }
743 else{
744 LOG(data,"<tr><td>%lu</td><td>%21.17g</td><td>%21.17g</td><td>%21.17g</td></tr>\n", i,gradients_rev[i],gradients_fwd[i],0.0);
745 }
746 }
747 LOG(data,"</table>\n");
748
749
750 /*Safe Second Derivative Calculations*/
751 if(data->use_yacas && data->SecondDer.safeder!=NULL){
752 fprintf(data->SecondDer.safeder,"@ Relation: %s Follows\n",rname);
753 }
754 for(i=0; i<num_var; i++){
755 status = (int32) RelationCalcSecondDerivSafe(inst,deriv_2nd,i);
756 safe_error_to_stderr( (enum safe_err *)&status );
757
758 if(data->use_yacas && data->SecondDer.safeder!=NULL){
759 /** @todo log calculated values here.*/ /*FIXME*/
760 for(j=0; j<num_var; j++){
761 fprintf(data->SecondDer.safeder,"%21.17g\n",deriv_2nd[j]); /*FIXME*/
762 }
763 }
764 }
765
766 /** End of Relation */
767 LOG(data,"</br><b>@ End of Relation </b>%p{%s}</br><hr></hr>",r,rname);
768
769
770 /** Freeing used memeory */
771 if(gradients_rev) {
772 ASC_FREE(gradients_rev);
773 }
774
775 if(gradients_fwd) {
776 ASC_FREE(gradients_fwd);
777 }
778
779 if(deriv_2nd) {
780 ASC_FREE(deriv_2nd);
781 }
782
783
784 RETURN;
785
786 #undef RETURN
787
788 }
789 #undef OUTENV
790 #undef VARFILE
791 #undef CASEFILE
792
793 #undef SAFEDER_2ND
794 #undef NONSAFEDER_2ND
795 #undef YACAS_2ND
796
797 #undef SAFEDER_1ST
798 #undef NONSAFEDER_1ST
799 #undef YACAS_1ST
800
801 #undef YACAS_IN_1ST
802 #undef YACAS_IN_2ND
803
804
805 /*===========================================================================*/
806 /* Registration information */
807
808 /* the list of tests */
809
810 #define TESTS(T,X) \
811 T(autodiff)
812
813 /* you shouldn't need to change the following */
814
815 #define TESTDECL(TESTFN) {#TESTFN,test_##TESTFN}
816
817 #define X ,
818
819 static CU_TestInfo autodiff_test_list[] = {
820 TESTS(TESTDECL,X)
821 X CU_TEST_INFO_NULL
822 };
823
824 static CU_SuiteInfo suites[] = {
825 {"compiler_autodiff", NULL, NULL, autodiff_test_list},
826 CU_SUITE_INFO_NULL
827 };
828
829 /*-------------------------------------------------------------------*/
830 CU_ErrorCode test_register_compiler_autodiff(void){
831 return CU_register_suites(suites);
832 }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22