/[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 3321 - (show annotations) (download) (as text)
Thu Jan 11 04:04:39 2018 UTC (5 months, 1 week ago) by jpye
File MIME type: text/x-csrc
File size: 23280 byte(s)
yacas test mostly working again, trouble with substitution before derivative determined.

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, see <http://www.gnu.org/licenses/>.
16 */
17 /**
18 @file
19 Unit test functions reverse ad & second derivatives
20 Created by: Mahesh Narayanamurthi
21 Revised by: Ben Allan - Rev.1
22 Mahesh Narayanamurthi - Rev. 2
23
24 TODO Note that this test case plays a DOUBLE ROLE of both preparing
25 the test files (yacas-output-*.txt) and then reading them and comparing
26 values with those returned by the computer algebra system YACAS. If you
27 want to change and of the model files being tested here, you will also
28 need to update the yacas output files against which the derivatives will
29 be compared.
30
31 To update the YACAS output files, use the following commands:
32
33 pushd ~/ascend/trunk \
34 && scons -j4 test ascend models solvers \
35 && LD_LIBRARY_PATH=. ASC_YACAS_GEN=1 test/test compiler_autodiff \
36 && cd ascend/compiler/test \
37 && python yacasgen.py \
38 && dirn=`pwd` \
39 && cd /usr/share/yacas \
40 && yacas $dirn/yacas-input-1st.txt > $dirn/yacas-output-1st.txt \
41 && yacas $dirn/yacas-input-2nd.txt > $dirn/yacas-output-2nd.txt \
42 && cd $dirn/../../.. && LD_LIBRARY_PATH=. test/test compiler_autodiff;
43 popd
44 */
45 #include <string.h>
46 #include <stdlib.h>
47 #include <stdio.h>
48
49 #include <ascend/general/env.h>
50 #include <ascend/general/ospath.h>
51 #include <ascend/general/list.h>
52 #include <ascend/general/ltmatrix.h>
53
54 #include <ascend/general/platform.h>
55 #include <ascend/utilities/ascEnvVar.h>
56 #include <ascend/utilities/error.h>
57
58 #include <ascend/compiler/ascCompiler.h>
59 #include <ascend/compiler/module.h>
60 #include <ascend/compiler/parser.h>
61 #include <ascend/compiler/library.h>
62 #include <ascend/compiler/symtab.h>
63 #include <ascend/compiler/simlist.h>
64 #include <ascend/compiler/instquery.h>
65 #include <ascend/compiler/parentchild.h>
66 #include <ascend/compiler/atomvalue.h>
67 #include <ascend/compiler/relation_io.h>
68 #include <ascend/compiler/reverse_ad.h>
69 #include <ascend/compiler/relation_util.h>
70 #include <ascend/compiler/mathinst.h>
71 #include <ascend/compiler/watchpt.h>
72 #include <ascend/compiler/initialize.h>
73 #include <ascend/compiler/name.h>
74 #include <ascend/compiler/visitinst.h>
75 #include <ascend/compiler/functype.h>
76 #include <ascend/compiler/safe.h>
77 #include <ascend/compiler/qlfdid.h>
78 #include <ascend/compiler/instance_io.h>
79 // FIXME include anything else?
80
81 #include <test/common.h>
82 #include <test/assertimpl.h>
83 #include <test/test_globals.h>
84
85 #define RAD_TOL 1e-05
86
87 #define AUTODIFF_DEBUG
88 #ifdef AUTODIFF_DEBUG
89 # define MSG CONSOLE_DEBUG
90 #else
91 # define MSG(ARGS...) ((void)0)
92 #endif
93
94 extern char ASC_TEST_PATH[PATH_MAX];
95
96 struct FileList
97 {
98 FILE * yacas; /* Meant for YACAS */
99 FILE * safeder; /* Meant for Safe Der */
100 FILE * nonsafeder; /* Meant for Non-Safe Der */
101 };
102
103 struct TestFileList
104 {
105 FILE * FirstDerYacas;
106 FILE * SecondDerYacas;
107 };
108
109 struct DiffTestData
110 {
111 FILE * outfile; /* Meant for Logging Debug Info */
112 FILE * varfile; /* Meant for recording the Variable Values */
113 struct FileList FirstDer;
114 struct FileList SecondDer;
115 int use_yacas;
116 FILE * first_yacas;
117 FILE * second_yacas;
118 struct Instance *root;
119 const char *infile;
120 int d0errors;
121 int d1errors;
122 int d1errors_yacas;
123 int d2errors_yacas;
124 int numrels;
125 };
126
127 #define INITDTD(d, log, varfile, FirstDer, SecondDer, use_yacas, YacasInFirst, YacasInSecond ,instroot, inputname) \
128 if(1){\
129 d.outfile = log; \
130 d.varfile = varfile; \
131 d.FirstDer = FirstDer;\
132 d.SecondDer = SecondDer;\
133 d.use_yacas = use_yacas; \
134 d.first_yacas = YacasInFirst;\
135 d.second_yacas = YacasInSecond;\
136 d.root = instroot; \
137 d.infile = inputname; \
138 d.d0errors= 0; \
139 d.d1errors= 0; \
140 d.d1errors_yacas= 0; \
141 d.d2errors_yacas= 0; \
142 d.numrels = 0;\
143 }
144
145 #define LOG(d, ...)\
146 if(d != NULL && d->outfile != NULL){ \
147 fprintf(d->outfile,__VA_ARGS__); \
148 }
149
150 /** Temporary Declarations */
151 static void AutomateDiffTest(struct Instance *inst, VOIDPTR ptr);
152
153 static void test_autodiff(void){
154 #define OUTENV "../ascend/compiler/test/LOG.html"
155 #define VARFILE "../ascend/compiler/test/varnames.txt"
156
157 #define SAFEDER_1ST "../ascend/compiler/test/derivs-safe-1st.txt"
158 #define SAFEDER_2ND "../ascend/compiler/test/derivs-safe-2nd.txt"
159 #define NONSAFEDER_1ST "../ascend/compiler/test/derivs-nonsafe-1st.txt"
160 #define NONSAFEDER_2ND "../ascend/compiler/test/derivs-nonsafe-2nd.txt"
161
162 #define YACAS_PREP_1ST "../ascend/compiler/test/yacas-prep-1st.txt"
163 #define YACAS_PREP_2ND "../ascend/compiler/test/yacas-prep-2nd.txt"
164 #define YACAS_OUT_1ST "../ascend/compiler/test/yacas-output-1st.txt"
165 #define YACAS_OUT_2ND "../ascend/compiler/test/yacas-output-2nd.txt"
166
167 #if 0
168 # define CASEFILE "test/reverse_ad/allmodels.a4c"
169 # define TYPENAME "allmodels"
170 #else
171 # define CASEFILE "test/ipopt/dummy.a4c"
172 # define TYPENAME "dummy"
173 #endif
174
175 #define USE_YACAS_ENV "ASC_YACAS_GEN"
176
177 int status;
178 int use_yacas = 0;
179
180 //struct module_t *m;
181
182 struct Name *name;
183 enum Proc_enum pe;
184 struct Instance *root;
185
186 struct DiffTestData data;
187
188 FILE *outfile = NULL;
189 FILE *varfile = NULL;
190
191 FILE *first_yacas = NULL;
192 FILE *second_yacas = NULL;
193
194 struct FileList FirstDer;
195 struct FileList SecondDer;
196
197 FirstDer.yacas = FirstDer.safeder = FirstDer.nonsafeder = NULL;
198 SecondDer.yacas = SecondDer.safeder = SecondDer.nonsafeder = NULL;
199
200 char env1[2*PATH_MAX];
201 Asc_CompilerInit(1);
202 /* set the needed environment variables so that models, solvers can be found */
203 snprintf(env1,2*PATH_MAX,ASC_ENV_LIBRARY "=%s","models");
204 CU_TEST(0 == Asc_PutEnv(env1));
205 CU_TEST(0 == Asc_PutEnv(ASC_ENV_SOLVERS "=solvers/ipopt"));
206
207 Asc_PutEnv(ASC_ENV_LIBRARY "=models");
208
209
210 // FIXME Use the environment variables here
211 MSG("ASC_TEST_PATH = '%s'",ASC_TEST_PATH);
212
213 struct FilePath *rootfp;
214 {
215 struct FilePath *tmp = ospath_new(ASC_TEST_PATH);
216 rootfp = ospath_getabs(tmp);
217 ospath_free(tmp);
218 }
219
220 #define OPENTESTFILE(FNAME,VAR,MODE) {\
221 struct FilePath *tmp = ospath_new_noclean(FNAME);\
222 struct FilePath *fp2 = ospath_concat(rootfp,tmp);\
223 ospath_cleanup(fp2);\
224 VAR = ospath_fopen(fp2,MODE);\
225 char *s=ospath_str(fp2);\
226 CU_ASSERT_PTR_NOT_NULL_FATAL(VAR);\
227 MSG("Opened test file '%s'",s);\
228 ASC_FREE(s);\
229 ospath_free(tmp);\
230 ospath_free(fp2);\
231 }
232
233 OPENTESTFILE(OUTENV,outfile,"w");
234
235 /** @TODO Open the following streams only if Environment Variable is set */
236 if(getenv(USE_YACAS_ENV) && atol(getenv(USE_YACAS_ENV))){
237 MSG("Generating YACAS input files for verification of ASCEND's values");
238 OPENTESTFILE(VARFILE,varfile,"w");
239 OPENTESTFILE(YACAS_PREP_2ND,SecondDer.yacas,"w");
240 OPENTESTFILE(SAFEDER_2ND,SecondDer.safeder,"w");
241 OPENTESTFILE(NONSAFEDER_2ND,SecondDer.nonsafeder,"w");
242 OPENTESTFILE(YACAS_PREP_1ST,FirstDer.yacas,"w");
243 OPENTESTFILE(SAFEDER_1ST,FirstDer.safeder,"w");
244 OPENTESTFILE(NONSAFEDER_1ST,FirstDer.nonsafeder,"w");
245 use_yacas=1;
246 }else{
247 MSG("Testing against pre-calculated derivatives from YACAS.");
248 OPENTESTFILE(YACAS_OUT_1ST,first_yacas,"r");
249 OPENTESTFILE(YACAS_OUT_2ND,second_yacas,"r");
250 }
251
252 ospath_free(rootfp);
253
254 /* load the file */
255 Asc_OpenModule(CASEFILE,&status);
256 CU_ASSERT_FATAL(status == 0);
257
258 /* parse it */
259 CU_ASSERT_FATAL(0 == zz_parse());
260
261 /* find the model */
262 symchar * type = AddSymbol(TYPENAME);
263 CU_ASSERT_FATAL(FindType(type)!=NULL);
264
265 /* instantiate it */
266 struct Instance *sim = SimsCreateInstance(type, AddSymbol("sim1"), e_normal, NULL);
267 CU_ASSERT_FATAL(sim!=NULL);
268
269 /*root of Simulation object*/
270 root = GetSimulationRoot(sim);
271
272 /** Call on_load */
273 error_reporter_tree_t *tree1 = error_reporter_tree_start(0);
274
275 name = CreateIdName(AddSymbol("on_load"));
276 pe = Initialize(root,name,"sim1",ASCERR,0, NULL, NULL);
277 CU_TEST_FATAL(0==error_reporter_tree_has_error(tree1));
278 error_reporter_tree_end(tree1);
279
280 CU_TEST_FATAL(pe == Proc_all_ok);
281
282 /* check for vars and rels */
283
284 INITDTD(data, outfile, varfile, FirstDer, SecondDer,use_yacas, first_yacas, second_yacas, root, CASEFILE);
285 if(data.outfile!=NULL){
286 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");
287 fprintf(data.outfile,"<center><h3>LOG File of Test-Suite for Automatic Differentiation</h3></center></br></br>\n");
288 fprintf(data.outfile,"<center><h4><u>Legend</u></h4><b>RAD_TOL = %21.17g</br>\n",RAD_TOL);
289 fprintf(data.outfile,"Difference(a,b) = a - b</br>\n");
290 fprintf(data.outfile,"Error(a,b) = (a - b)/a </br>\n");
291 fprintf(data.outfile,"<font color='gold'>gold</font> = Residual Error </br>\n");
292 fprintf(data.outfile,"<font color='indigo'>indigo</font>= Gradient Error </br>\n");
293 fprintf(data.outfile,"<font color='purple'>purple</font>= Gradient Mismatch </br>\n");
294 fprintf(data.outfile,"<font color='red'>red</font>= Second Partials Mismatch</b> </center></br></br>\n");
295 }
296
297 VisitInstanceTreeTwo(root, AutomateDiffTest, 0, 0, &data);
298 sim_destroy(sim);
299 Asc_CompilerDestroy();
300
301 if(data.outfile!=NULL){
302 fprintf(data.outfile,"</br><center><u><h4>SUMMARY:</h4></u></br>\n");
303 fprintf(data.outfile,"<b>No. of Residual Errors</b> =<b> %d </b></br>\n",data.d0errors);
304 fprintf(data.outfile,"<b>No. of First RAD Errors</b> =<b> %d </b></br>\n",data.d1errors);
305 fprintf(data.outfile,"<b>No. of Total Errors</b> =<b> %d</b></br></br>\n",data.d0errors+data.d1errors);
306 fprintf(data.outfile,"<b>No. of First YACAS Mismatches</b> =<b> %d</b></br>\n",data.d1errors_yacas);
307 fprintf(data.outfile,"<b>No. of Second YACAS Mismatches</b> =<b> %d</b></br>\n",data.d2errors_yacas);
308 fprintf(data.outfile,"<b>No. of Total Mismatches</b> =<b> %d</b></center></br>\n</body>\n</html>"
309 ,data.d1errors_yacas + data.d2errors_yacas
310 );
311 }
312 MSG("TOTAL OF %d RELATIONS TESTED:",data.numrels);
313 MSG(" Residual errors: %d",data.d0errors);
314 MSG(" First deriv errors: %d",data.d1errors);
315 MSG(" First YACAS mismatches: %d",data.d1errors_yacas);
316 MSG(" Second YACAS mismatches: %d",data.d2errors_yacas);
317
318 #define OSPCLEAN(FPTR) if(FPTR!=NULL){fclose(FPTR);data.FPTR = NULL;}
319 OSPCLEAN(outfile);
320 if(use_yacas){
321 OSPCLEAN(varfile);
322 OSPCLEAN(SecondDer.yacas);
323 OSPCLEAN(SecondDer.safeder);
324 OSPCLEAN(SecondDer.nonsafeder);
325 OSPCLEAN(FirstDer.yacas);
326 OSPCLEAN(FirstDer.safeder);
327 OSPCLEAN(FirstDer.nonsafeder);
328 }else{
329 OSPCLEAN(first_yacas);
330 OSPCLEAN(second_yacas);
331 }
332
333 CU_ASSERT( 0 == (data.d0errors + data.d1errors) );
334 MSG("For non-fatal errors refer ascend/compiler/test/LOG.html");
335 }
336
337
338 /*------------------------------------------------------------------------------
339 Routine to conduct derivative tests on a single relation
340 */
341 static void AutomateDiffTest(struct Instance *inst, VOIDPTR ptr){
342 double err,residual_rev,residual_fwd;
343 double *gradients_rev,*gradients_fwd,*deriv_2nd;
344 float yacas_first_der = 0.0; /* compiler complains type mismatch when using */
345 float yacas_second_der = 0.0; /*double to read using fscanf(File*,"%21.17g"...) */
346 unsigned long num_var,i,j;
347 enum Expr_enum reltype;
348 struct relation *r;
349 int32 status;
350 char *rname = NULL, *infix_rel = NULL, *varname;
351 char buf[20];
352 struct RXNameData myrd = {"x",NULL,""};
353 struct DiffTestData *data = (struct DiffTestData*) ptr;
354 struct Instance *var_inst;
355 #define RETURN if(rname != NULL)ASC_FREE(rname); return
356
357 if(inst==NULL || InstanceKind(inst)!=REL_INST)return;
358
359 if (data != NULL && data->outfile != NULL) {
360 rname = WriteInstanceNameString(inst, data->root);
361 }
362
363 MSG("Relation %d: '%s'",data->numrels,rname);
364 data->numrels++;
365
366 /* get the relation, check type */
367 r = (struct relation *)GetInstanceRelation(inst, &reltype);
368 if(r == NULL){
369 LOG(data, "<font color='orange'><b>! skipping instance with null struct relation:</b></font> %s</br>\n", rname);
370 RETURN;
371 }
372 if(reltype != e_token){
373 LOG(data, "<font color='orange'><b>! skipping non-token relation</b></font> %s</br>\n", rname);
374 RETURN;
375 }
376 LOG(data, "</br><font color='green'><b> Evaluating token relation </b></font> <b> %s </b> </br>\n", rname);
377
378 /* convert relation to YACAS code string */
379 infix_rel = WriteRelationString(inst, data->root
380 , (WRSNameFunc)RelationVarXName, NULL, relio_yacas, NULL
381 );
382 num_var = NumberVariables(r); // FIXME or should we use rel_n_incidences
383
384 /* write out the values of each variable; log rel + var vals */
385 LOG(data,"</br></br><b>Relation: %s </br>\n",infix_rel); // Do I need to escape this??
386 for(i=0; i<num_var; i++) {
387 varname = RelationVarXName(r,i+1,&myrd);
388 if (varname!=NULL){
389 LOG(data,"%s:=",varname);
390 var_inst = RelationVariable(r,i+1);
391 LOG(data,"%21.17g</br>\n",RealAtomValue(var_inst));
392 }
393 }
394 LOG(data,"</b>\n");
395
396 gradients_rev = ASC_NEW_ARRAY(double,num_var); // or rel_n_incidences
397 CU_TEST_FATAL(NULL!=gradients_rev);
398 gradients_fwd = ASC_NEW_ARRAY(double,num_var); // or rel_n_incidences
399 CU_TEST_FATAL(NULL!=gradients_fwd);
400 deriv_2nd = ASC_NEW_ARRAY(double,num_var); // or rel_n_incidences
401 CU_TEST_FATAL(NULL!=deriv_2nd);
402
403 /* write out data files for use by yacas */
404 if(data->use_yacas){
405 if(data->FirstDer.yacas!=NULL && data->SecondDer.yacas!=NULL && data->varfile!=NULL){
406 /** Print Variables Values First */
407 for(i=0; i<num_var; i++) {
408 varname = RelationVarXName(r,i+1,&myrd);
409 if (varname!=NULL){
410 fprintf(data->varfile,"%s==",varname);
411 var_inst = RelationVariable(r,i+1);
412 fprintf(data->varfile,"%-21.17g\n",RealAtomValue(var_inst));
413 }
414 }
415 fprintf(data->varfile,"@ Relation: %s\n",rname);
416
417 if (infix_rel!=NULL){
418 for(i=0; i<num_var;i++){
419 varname = RelationVarXName(r,i+1,&myrd);
420 strncpy(buf,varname,20);
421
422 // Generating Output for First Derivative Yacas Input file
423 fprintf(data->FirstDer.yacas
424 ,"D(%s) (%s)\n"
425 ,buf,infix_rel
426 );
427
428 // Generating Output for Second Derivative Yacas Input file
429 for (j=0;j<num_var;j++){
430 varname = RelationVarXName(r,j+1,&myrd);
431 fprintf(data->SecondDer.yacas
432 ,"D(%s) D(%s) (%s)\n"
433 ,buf,varname,infix_rel
434 );
435 }
436 }
437 }
438 else{
439 ERROR_REPORTER_HERE(ASC_PROG_FATAL,"Infix Relation is NULL");
440 }
441 }
442 }
443 ASC_FREE(infix_rel);
444
445 /*--- non-safe evaluation routines ---*/
446
447 /* TODO we need to sigfpe trap this code or use the safe versions. */
448
449 RelationCalcResidGradRev(inst,&residual_rev,gradients_rev);
450 RelationCalcResidGrad(inst,&residual_fwd,gradients_fwd);
451
452 LOG(data,"</br> <b> Table of Values for Residuals </b> </br>\n");
453 LOG(data,"\n<table BORDER>\n");
454 LOG(data,"<tr><td>ASCEND(NONSAFE,REV)</td><td>ASCEND(NONSAFE,FWD)</td><td>Error</td></tr>\n");
455
456 CU_ASSERT(fabs(residual_rev - residual_fwd) <= RAD_TOL);
457
458 if(fabs(residual_rev - residual_fwd) > RAD_TOL) {
459 data->d0errors ++;
460 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));
461 }else{
462 LOG(data,"<tr><td>%21.17g</td><td>%21.17g</td><td>%.4g</td></tr>\n",residual_rev,residual_fwd,0.0);
463 }
464 LOG(data,"</table>\n");
465
466 if(data->first_yacas!=NULL){
467 char s[PATH_MAX+1];
468 if(NULL==fgets(s,PATH_MAX,data->first_yacas)){
469 MSG("line: %s",s);
470 CU_FAIL_FATAL("Expected line from " YACAS_OUT_1ST);
471 }
472 }
473
474 LOG(data,"</br> <b> Table of Values for First Derivatives </b> </br>\n");
475 LOG(data,"\n<table BORDER>\n");
476 for(i=0; i<num_var; i++) {
477 if(data->use_yacas && data->FirstDer.nonsafeder!=NULL){
478 /* Recording Reverse AD Non-Safe Derivatives to file */
479 fprintf(data->FirstDer.nonsafeder,"%21.17g\n",gradients_rev[i]);
480
481 }else if(data->first_yacas!=NULL){
482 /* Benchmarking Non-Safe Gradient Errors against Yacas */
483 if(!feof(data->first_yacas)){
484 fscanf(data->first_yacas,"%g\n",&yacas_first_der);
485
486 if(yacas_first_der!=0.0){
487 err = fabs((double)yacas_first_der - gradients_rev[i]) / (double)yacas_first_der;
488 }else{
489 err = gradients_rev[i];
490 }
491 err = fabs(err);
492 LOG(data,"<tr><td>Column</td><td>ASCEND(nonsafe,rev)</td><td>YACAS</td><td>Percentage Mismatch</td></tr>\n");
493 CU_TEST(err < RAD_TOL);
494 if (err > RAD_TOL) {
495 MSG("dR/dx%lu_yacas = %g",i,yacas_first_der);
496 MSG("dR/dx%lu_rev = %g",i,gradients_rev[i]);
497 MSG("error = %e, tolerance = %e",err,RAD_TOL);
498 data->d1errors_yacas ++;
499 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);
500 }else{
501 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);
502 }
503 }
504 }
505
506 if(gradients_fwd[i] != 0.0){
507 err = (gradients_fwd[i] - gradients_rev[i]) / gradients_fwd[i];
508 }else{
509 err = gradients_rev[i]; // These are totally different quantities should I make err as 1?
510 }
511 err = fabs(err);
512 LOG(data,"<tr><td>Column</td><td>ASCEND(nonsafe,rev)</td><td>ASCEND(nonsafe,fwd)</td><td>Percentage Mismatch</td></tr>\n");
513 //CU_ASSERT(err <= RAD_TOL);
514 if(err > RAD_TOL){
515 MSG("Failed tolerance in first deriv #%lu",i);
516 CU_FAIL("Error exceeded tolerance");
517 data->d1errors ++;
518 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);
519 }else{
520 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);
521 }
522 }
523 LOG(data,"</table>\n");
524
525 /*Non Safe Second Derivative Calculations*/
526 if(data->use_yacas && data->SecondDer.nonsafeder!=NULL){
527 fprintf(data->SecondDer.nonsafeder,"@ Relation: %s\n",rname);
528 }
529
530 if(data->second_yacas!=NULL){
531 char s[PATH_MAX+1];
532 fgets(s,PATH_MAX,data->second_yacas);
533 char s2[PATH_MAX];
534 snprintf(s2,PATH_MAX,"@ Relation: %s\n",rname);
535 if(0!=strcmp(s,s2)){
536 MSG("line: <%s>",s);
537 MSG("expected: <%s>",s2);
538 CU_FAIL_FATAL("wrong relation read from second-derivs yacas file");
539 }
540 }
541
542 LOG(data,"</br> <b> Table of Values for Second Derivatives </b> </br>\n");
543 LOG(data,"\n<table BORDER>\n");
544 LOG(data,"<tr><td>Row</td><td>Column</td><td>ASCEND (nonsafe)</td><td>YACAS</td><td>Percentage Mismatch</td></tr>\n");
545 for(i=0; i<num_var; i++){
546 RelationCalcSecondDeriv(inst,deriv_2nd,i);
547 if(data->use_yacas && data->SecondDer.nonsafeder!=NULL){
548 /** @todo log calculated values and indiceshere.*/ /*FIXME*/
549 for(j=0; j<num_var; j++){
550 fprintf(data->SecondDer.nonsafeder,"%21.17g\n",deriv_2nd[j]); /*FIXME*/
551 }
552 }
553 else if (data->second_yacas!=NULL){
554 /* Benchmarking Non-Safe Second Derivative Errors against Yacas */
555 for(j=0; j<num_var; j++){
556 if(!feof(data->second_yacas)){
557 fscanf(data->second_yacas,"%g\n",&yacas_second_der);
558
559 if(yacas_second_der!=0.0){
560 err = ((double)yacas_second_der - deriv_2nd[j]) / (double)yacas_second_der; /* todo err scaling */
561 }else{
562 err = deriv_2nd[j];
563 }
564 err = fabs(err);
565 if(err > RAD_TOL) {
566 MSG("d2R/dx%ludx%lu_yacas = %g",i,j,yacas_second_der);
567 MSG("d2R/dx%ludx%lu_rev = %g",i,j,deriv_2nd[j]);
568 CU_FAIL(err <= RAD_TOL);
569 data->d2errors_yacas ++;
570 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);
571 }
572 else{
573 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);
574 }
575 }
576 }
577 }
578 }
579
580 LOG(data,"</table>\n");
581
582 /** Testing safe evaluation routines */
583
584 status = (int32) RelationCalcResidGradRevSafe(inst,&residual_rev,gradients_rev);
585 safe_error_to_stderr( (enum safe_err *)&status );
586 status = RelationCalcResidGradSafe(inst,&residual_fwd,gradients_fwd);
587 safe_error_to_stderr( (enum safe_err *)&status );
588
589 LOG(data,"</br> <b> Table of Values for Residuals </b> </br>\n");
590 LOG(data,"\n<table BORDER>\n");
591 LOG(data,"<tr><td>ASCEND(SAFE,REV)</td><td>ASCEND(SAFE,FWD)</td><td>Error</td></tr>\n");
592 CU_ASSERT(fabs(residual_rev - residual_fwd) <= RAD_TOL);
593 if ( fabs(residual_rev - residual_fwd) > RAD_TOL) {
594 data->d0errors ++;
595 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));
596 }
597 else{
598 LOG(data,"<tr><td>%21.17g</td><td>%21.17g</td><td>%.4g</td></tr>\n", residual_rev,residual_fwd,0.0);
599 }
600 LOG(data,"</table>\n");
601
602 LOG(data,"</br> <b> Table of Values for First Derivatives </b> </br>\n");
603 LOG(data,"\n<table BORDER>\n");
604 LOG(data,"<tr><td>Column</td><td>ASCEND(SAFE,REV)</td><td>ASCEND(SAFE,FWD)</td><td>Percentage Mismatch</td></tr>\n");
605 for(i=0;i<num_var;i++) {
606
607 if(data->use_yacas && data->FirstDer.safeder!=NULL){
608 /* Recording Reverse AD Safe Derivatives to file */
609 fprintf(data->FirstDer.safeder,"%21.17g\n",gradients_rev[i]);
610 }
611
612 if (gradients_fwd[i] != 0.0) {
613 err = ( gradients_fwd[i] - gradients_rev[i] ) / gradients_fwd[i];
614 } else {
615 err = gradients_rev[i];
616 }
617 err = fabs(err);
618 CU_ASSERT(err <= RAD_TOL);
619 if (err > RAD_TOL) {
620 data->d1errors ++;
621 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);
622 }
623 else{
624 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);
625 }
626 }
627 LOG(data,"</table>\n");
628
629 /*Safe Second Derivative Calculations*/
630 if(data->use_yacas && data->SecondDer.safeder!=NULL){
631 fprintf(data->SecondDer.safeder,"@ Relation: %s\n",rname);
632 }
633 for(i=0; i<num_var; i++){
634 status = (int32) RelationCalcSecondDerivSafe(inst,deriv_2nd,i);
635 safe_error_to_stderr( (enum safe_err *)&status );
636
637 if(data->use_yacas && data->SecondDer.safeder!=NULL){
638 /** @todo log calculated values here.*/ /*FIXME*/
639 for(j=0; j<num_var; j++){
640 fprintf(data->SecondDer.safeder,"%21.17g\n",deriv_2nd[j]); /*FIXME*/
641 }
642 }
643 }
644
645 /** End of Relation */
646 LOG(data,"</br><b>@ End of Relation </b>%p{%s}</br><hr></hr>",r,rname);
647
648 /** Freeing used memeory */
649 if(gradients_rev)ASC_FREE(gradients_rev);
650 if(gradients_fwd)ASC_FREE(gradients_fwd);
651 if(deriv_2nd)ASC_FREE(deriv_2nd);
652 RETURN;
653 #undef RETURN
654 }
655 /*===========================================================================*/
656 /* Registration information */
657
658 /* the list of tests */
659
660 #define TESTS(T) \
661 T(autodiff)
662
663 REGISTER_TESTS_SIMPLE(compiler_autodiff, TESTS)
664

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