/[ascend]/branches/relerrorlist/ascend/compiler/test/test_autodiff.c
ViewVC logotype

Contents of /branches/relerrorlist/ascend/compiler/test/test_autodiff.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3197 - (show annotations) (download) (as text)
Sun Apr 23 03:30:25 2017 UTC (17 months, 3 weeks ago) by jpye
File MIME type: text/x-csrc
File size: 22278 byte(s)
test suite almost done, one memory leak still to be fixed.

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

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