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

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