/[ascend]/trunk/ascend/compiler/relation_util.c
ViewVC logotype

Contents of /trunk/ascend/compiler/relation_util.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2130 - (show annotations) (download) (as text)
Fri Dec 18 01:32:28 2009 UTC (10 years, 6 months ago) by mahesh
File MIME type: text/x-csrc
File size: 130556 byte(s)
Merging all changes from Mahesh's GSOC project.
svn merge -r2273:2671 svn://ascend.cheme.cmu.edu/ascend/code/branches/mahesh .
See http://ascendwiki.cheme.cmu.edu/User:Mnm87 for details.
1 /* ASCEND modelling environment
2 Copyright (C) 1997, 2006 Carnegie Mellon University
3 Copyright (C) 1993, 1994 Joseph James Zaher, Benjamin Andrew Allan
4 Copyright (C) 1993 Joseph James Zaher
5 Copyright (C) 1990 Thomas Guthrie Epperly, Karl Michael Westerberg
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2, or (at your option)
10 any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 59 Temple Place - Suite 330,
20 Boston, MA 02111-1307, USA.
21 *//**
22 @file
23 Relation utility functions for Ascend
24
25 This module defines the dimensionality checking and some other
26 relation auxillaries for ASCEND.
27 *//*
28 Last in CVS: $Revision: 1.44 $ $Date: 1998/04/23 23:51:09 $ $Author: ballan $
29 */
30
31 #include <math.h>
32 #include <errno.h>
33 #include <stdarg.h>
34
35 #include "relation_util.h"
36
37
38 #include <ascend/utilities/ascMalloc.h>
39 #include <ascend/utilities/ascPanic.h>
40 #include <ascend/general/mathmacros.h>
41 #include <ascend/general/list.h>
42 #include <ascend/general/ltmatrix.h>
43 #include <ascend/general/dstring.h>
44 #include "symtab.h"
45 #include "vlist.h"
46 #include "dimen_io.h"
47 #include "instance_enum.h"
48 #include "bintoken.h"
49 #include "find.h"
50 #include "atomvalue.h"
51 #include "instance_name.h"
52 #include "rel_blackbox.h"
53 #include "relation.h"
54 #include "instance_io.h"
55 #include "instquery.h"
56 #include "visitinst.h"
57 #include "mathinst.h"
58 #include "rootfind.h"
59 #include "func.h"
60 #include "parentchild.h"
61
62
63 #ifndef NDEBUG
64 #include "relation_io.h"
65 #endif
66
67 /*------------------------------------------------------------------------------
68 DATA TYPES AND GLOBAL VARS
69 */
70
71 int g_check_dimensions_noisy = 1;
72 #define GCDN g_check_dimensions_noisy
73
74 /** global relation pointer to avoid passing a relation recursively */
75 static struct relation *glob_rel;
76
77 /**
78 The following global variables are used thoughout the
79 functions called by RelationFindroot.
80
81 These should probably be located at the top of this
82 file alonge with glob_rel. [OK, let it be so then. -- JP]
83 */
84 static unsigned long glob_varnum;
85 static int glob_done;
86
87 /* some data structurs...*/
88
89 struct dimnode {
90 dim_type d;
91 enum Expr_enum type;
92 short int_const;
93 double real_const;
94 struct fraction power;
95 };
96
97 struct ds_soln_list {
98 int length,capacity;
99 double *soln;
100 };
101
102 /*
103 Define the following if you want ASCEND to panic when it hits a
104 relation error in this file. This will help with debugging (GDB).
105
106 Comment it out for ERROR_REPORTER to be used instead.
107 */
108 #define RELUTIL_CHECK_ABORT
109
110 /*------------------------------------------------------------------------------
111 forward declarations
112 */
113
114 static struct fraction real_to_frac(double real);
115 int ArgsForRealToken(enum Expr_enum type);
116 static int IsZero(struct dimnode *node);
117
118 /* bunch of support functions for RelationFindRoots */
119 static double RootFind(struct relation *rel, double *lower_bound, double *upper_bound,
120 double *nominal,double *tolerance,unsigned long varnum,int *status);
121 static int CalcResidGivenValue(int *mode, int *m, unsigned long *varnum,double *val, double *u, double *f, double *g);
122 int RelationInvertTokenTop(struct ds_soln_list *soln_list);
123 int RelationInvertToken(struct relation_term **term,struct ds_soln_list *soln_list,enum safe_err *not_safe);
124 static void SetUpInvertTokenTop(struct relation_term **invert_side,double *value);
125 static int SetUpInvertToken(struct relation_term *term,struct relation_term **invert_side,double *value);
126 static int SearchEval_Branch(struct relation_term *term);
127 static void InsertBranchResult(struct relation_term *term, double value);
128 static void remove_soln( struct ds_soln_list *sl, int ndx);
129 static void append_soln( struct ds_soln_list *sl, double soln);
130 static struct relation *RelationTmpTokenCopy(CONST struct relation *src);
131 static int RelationTmpCopySide(union RelationTermUnion *old,unsigned long len,union RelationTermUnion *arr);
132 static struct relation *RelationCreateTmp(unsigned long lhslen, unsigned long rhslen, enum Expr_enum relop);
133
134 /* the following appear only to be used locally, so I've made them static -- JP */
135
136 static int RelationCalcDerivative(struct Instance *i, unsigned long vindex, double *grad);
137 /**<
138 * This calculates the derivative of the relation df/dx (f = lhs-rhs)
139 * where x is the VINDEX-th entry in the relation's var list.
140 * The var list is a gl_list_t indexed from 1 to length.
141 * Non-zero return value implies a problem.<br><br>
142 *
143 * Notes: This function is a possible source of floating point
144 * exceptions and should not be used during compilation.
145 */
146
147 static enum safe_err
148 RelationCalcDerivativeSafe(struct Instance *i, unsigned long vindex, double *grad);
149 /**<
150 * Calculates the derivative safely.
151 * Non-zero return value implies a problem.
152 */
153
154 #ifndef NDEBUG
155 static int relutil_check_inst_and_res(struct Instance *i, double *res);
156 #endif
157
158 #ifndef NDEBUG
159 # define CHECK_INST_RES(i,res,retval) if(!relutil_check_inst_and_res(i,res)){return retval;}
160 #else
161 # define CHECK_INST_RES(i,res,retval) ((void)0)
162 #endif
163
164 /*------------------------------------------------------------------------------
165 SOME STUFF TO DO WITH DIMENSIONS
166 */
167
168 /**
169 @TODO what this does needs to be documented here
170
171 UNUSED FUNCTION??
172 */
173 static void apply_term_dimensions(CONST struct relation *rel,
174 struct relation_term *rt,
175 struct dimnode *first,
176 struct dimnode *second,
177 int *con,
178 int *wild
179 ){
180 enum Expr_enum type;
181
182 switch(type=RelationTermType(rt)) {
183 case e_zero:
184 CopyDimensions(WildDimension(),&(first->d));
185 first->real_const = 0.0;
186 first->type = type;
187 break;
188
189 case e_int:
190 CopyDimensions(Dimensionless(),&(first->d));
191 first->int_const = (short)TermInteger(rt);
192 first->type = type;
193 break;
194
195 case e_real:
196 CopyDimensions(TermDimensions(rt),&(first->d));
197 first->real_const = TermReal(rt);
198 first->type = type;
199 break;
200
201 case e_var: {
202 struct Instance *var = RelationVariable(rel,TermVarNumber(rt));
203 CopyDimensions(RealAtomDims(var),&(first->d));
204 first->type = type;
205 break;
206 }
207 case e_func: {
208 enum Func_enum id = FuncId(TermFunc(rt));
209 switch( id ) {
210 case F_ABS:
211 case F_HOLD:
212 /* no checking or scaling */
213 break;
214
215 case F_SQR:
216 /* no checking, just simple scaling */
217 first->d = ScaleDimensions(&(first->d),CreateFraction(2,1));
218 break;
219
220 case F_CUBE:
221 /* no checking, just simple scaling */
222 first->d = ScaleDimensions(&(first->d),CreateFraction(3,1));
223 break;
224
225 case F_SQRT:
226 /* no checking, just simple scaling */
227 first->d = ScaleDimensions(&(first->d),CreateFraction(1,2));
228 break;
229
230 case F_CBRT:
231 /* no checking, just simple scaling */
232 first->d = ScaleDimensions(&(first->d),CreateFraction(1,3));
233 break;
234
235 case F_EXP:
236 case F_LN:
237 case F_LNM:
238 case F_LOG10:
239 #ifdef HAVE_ERF
240 case F_ERF:
241 #endif /* HAVE_ERF */
242 case F_SINH:
243 case F_COSH:
244 case F_TANH:
245 case F_ARCSINH:
246 case F_ARCCOSH:
247 case F_ARCTANH:
248 /**
249 *** first must now be dimensionless. It will
250 *** end up dimensionless as well.
251 **/
252 if( IsWild(&(first->d)) && !IsZero(first) ) {
253 if( !*wild ) *wild = TRUE;
254 if (GCDN) {
255 FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
256 FPRINTF(ASCERR," in function %s.\n",
257 FuncName(TermFunc(rt)));
258 }
259 } else if( !IsWild(&(first->d)) &&
260 CmpDimen(&(first->d),Dimensionless()) ) {
261 if( *con ) *con = FALSE;
262 if (GCDN) {
263 FPRINTF(ASCERR,"ERROR: Function %s called with\n",
264 FuncName(TermFunc(rt)));
265 FPRINTF(ASCERR," dimensions ");
266 WriteDimensions(ASCERR,&(first->d));
267 FPRINTF(ASCERR,".\n");
268 }
269 }
270 CopyDimensions(Dimensionless(),&(first->d));
271 break;
272
273 case F_SIN:
274 case F_COS:
275 case F_TAN: {
276 /**
277 *** first must now be of dimension D_PLANE_ANGLE.
278 *** It will then be made dimensionless.
279 **/
280 if( IsWild(&(first->d)) && !IsZero(first) ) {
281 if( !*wild ) *wild = TRUE;
282 if (GCDN) {
283 FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
284 FPRINTF(ASCERR," in function %s.\n",
285 FuncName(TermFunc(rt)) );
286 }
287 }else{
288 if( !IsWild(&(first->d)) &&
289 CmpDimen(&(first->d),TrigDimension()) ) {
290 if( *con ) *con = FALSE;
291 if (GCDN) {
292 FPRINTF(ASCERR,"ERROR: Function %s called with\n",
293 FuncName(TermFunc(rt)));
294 FPRINTF(ASCERR," dimensions ");
295 WriteDimensions(ASCERR,&(first->d));
296 FPRINTF(ASCERR,".\n");
297 }
298 }
299 }
300 CopyDimensions(Dimensionless(),&(first->d));
301 break;
302 }
303
304 case F_ARCSIN:
305 case F_ARCCOS:
306 case F_ARCTAN:
307 /**
308 *** first must now be dimensionless. It will
309 *** end up with dimension D_PLANE_ANGLE
310 **/
311 if( IsWild(&(first->d)) && !IsZero(first) ) {
312 if( !*wild ) *wild = TRUE;
313 if (GCDN) {
314 FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
315 FPRINTF(ASCERR," in function %s.\n",
316 FuncName(TermFunc(rt)));
317 }
318 } else if( !IsWild(&(first->d)) &&
319 CmpDimen(&(first->d),Dimensionless()) ) {
320 if( *con ) *con = FALSE;
321 if (GCDN) {
322 FPRINTF(ASCERR,"ERROR: Function %s called with\n",
323 FuncName(TermFunc(rt)));
324 FPRINTF(ASCERR," dimensions ");
325 WriteDimensions(ASCERR,&(first->d));
326 FPRINTF(ASCERR,".\n");
327 }
328 }
329 CopyDimensions(TrigDimension(),&(first->d));
330 break;
331 }
332 first->type = type;
333 break;
334 }
335
336 case e_uminus:
337 first->type = type;
338 break;
339
340 case e_times:
341 first->d = AddDimensions(&(first->d),&(second->d));
342 first->type = type;
343 break;
344
345 case e_divide:
346 first->d = SubDimensions(&(first->d),&(second->d));
347 first->type = type;
348 break;
349
350 case e_power: /* fix me and add ipower */
351 if( IsWild(&(second->d)) && !IsZero(second) ) {
352 if( !*wild ) *wild = TRUE;
353 if (GCDN) {
354 FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
355 FPRINTF(ASCERR," in exponent.\n");
356 }
357 } else if( !IsWild(&(second->d)) &&
358 CmpDimen(&(second->d),Dimensionless()) ) {
359 if( *con ) *con = FALSE;
360 if (GCDN) {
361 FPRINTF(ASCERR,"ERROR: Exponent has dimensions ");
362 WriteDimensions(ASCERR,&(second->d));
363 FPRINTF(ASCERR,".\n");
364 }
365 }
366 CopyDimensions(Dimensionless(),&(second->d));
367 switch( second->type ) {
368 case e_int:
369 if( !IsWild(&(first->d)) &&
370 CmpDimen(&(first->d),Dimensionless()) ) {
371 struct fraction power;
372 power = CreateFraction(second->int_const,1);
373 first->d = ScaleDimensions(&(first->d),power);
374 }
375 break;
376
377 case e_real:
378 if( !IsWild(&(first->d)) &&
379 CmpDimen(&(first->d),Dimensionless()) ) {
380 struct fraction power;
381 power = real_to_frac(second->real_const);
382 first->d = ScaleDimensions(&(first->d),power);
383 }
384 break;
385
386 /* what about e_zero? */
387 default:
388 if( IsWild(&(first->d)) && !IsZero(first) ) {
389 if( !*wild ) *wild = TRUE;
390 if (GCDN) {
391 FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
392 FPRINTF(ASCERR," raised to a non-constant power.\n");
393 }
394 } else if( !IsWild(&(first->d)) &&
395 CmpDimen(&(first->d),Dimensionless()) ) {
396 if( *con ) *con = FALSE;
397 if (GCDN) {
398 FPRINTF(ASCERR,"ERROR: Dimensions ");
399 WriteDimensions(ASCERR,&(first->d));
400 FPRINTF(ASCERR," are\n");
401 FPRINTF(ASCERR," raised to a non-constant power.\n");
402 }
403 }
404 CopyDimensions(Dimensionless(),&(first->d));
405 break;
406
407 }
408 first->type = type;
409 break;
410
411 case e_plus:
412 case e_minus:
413 if( IsWild(&(first->d)) && IsZero(first) ) {
414 /* first wild zero */
415 CopyDimensions(&(second->d),&(first->d));
416 first->type = second->type;
417 if( second->type==e_int )
418 first->int_const = second->int_const;
419 if( second->type==e_real )
420 first->real_const = second->real_const;
421 } else if( IsWild(&(first->d)) && !IsZero(first) ) {
422 /* first wild non-zero */
423 if( IsWild(&(second->d)) && !IsZero(second) ) {
424 /* second wild non-zero */
425 if( !*wild ) *wild = TRUE;
426 if (GCDN) {
427 FPRINTF(ASCERR,"ERROR: %s has wild dimensions on\n",
428 type==e_plus ? "Addition":"Subtraction");
429 FPRINTF(ASCERR," left and right hand sides.\n");
430 }
431 first->type = type;
432 } else if( !IsWild(&(second->d)) ) {
433 /* second not wild */
434 if( !*wild ) *wild = TRUE;
435 if (GCDN) {
436 FPRINTF(ASCERR,"ERROR: %s has wild dimensions on\n",
437 type==e_plus ? "Addition":"Subtraction");
438 FPRINTF(ASCERR," left hand side.\n");
439 }
440 CopyDimensions(&(second->d),&(first->d));
441 first->type = type;
442 }
443 } else if( !IsWild(&(first->d)) ) {
444 /* first not wild */
445 if( IsWild(&(second->d)) && !IsZero(second) ) {
446 /* second wild non-zero */
447 if( !*wild ) *wild = TRUE;
448 if (GCDN) {
449 FPRINTF(ASCERR,"ERROR: %s has wild dimensions on\n",
450 type==e_plus ? "Addition":"Subtraction");
451 FPRINTF(ASCERR," right hand side.\n");
452 }
453 first->type = type;
454 } else if ( !IsWild(&(second->d)) ) {
455 /* second not wild */
456 if( CmpDimen(&(first->d),&(second->d)) ) {
457 if( *con ) *con = FALSE;
458 if (GCDN) {
459 FPRINTF(ASCERR,"ERROR: %s has dimensions ",
460 type==e_plus ? "Addition":"Subtraction");
461 WriteDimensions(ASCERR,&(first->d));
462 FPRINTF(ASCERR," on left\n");
463 FPRINTF(ASCERR," and dimensions ");
464 WriteDimensions(ASCERR,&(second->d));
465 FPRINTF(ASCERR," on right.\n");
466 }
467 }
468 first->type = type;
469 }
470 }
471 break;
472
473 default:
474 FPRINTF(ASCERR,"ERROR: Unknown relation term type.\n");
475 if( *con ) *con = FALSE;
476 first->type = type;
477 break;
478 }
479 }
480
481 /**
482 @TODO what this does needs to be documented here
483
484 UNUSED FUNCTION??
485 */
486 int RelationCheckDimensions(struct Instance *relinst, dim_type *dimens)
487 {
488 CONST struct relation *rel;
489 enum Expr_enum reltype= 0;
490 struct dimnode *stack, *sp;
491 CONST dim_type* lhsdim;
492 struct Instance * lhs;
493 int consistent = TRUE;
494 int wild = FALSE;
495 unsigned long c, len;
496
497 rel = GetInstanceRelation(relinst, &reltype);
498 if ( !IsWild(RelationDim(rel)) ) { /* don't do this twice */
499 CopyDimensions(RelationDim(rel),dimens);
500 return 2;
501 }
502 if (reltype == e_glassbox || reltype == e_opcode) {
503 return 0;
504 }
505 if (reltype == e_blackbox) {
506 lhs = BlackBoxGetOutputVar(rel);
507 lhsdim = RealAtomDims(lhs);
508 CopyDimensions(lhsdim,dimens);
509 return( !IsWild(dimens) );
510 }
511 /* else token relation */
512 sp = stack = (struct dimnode *)
513 ascmalloc(RelationDepth(rel)*sizeof(struct dimnode));
514 switch( RelationRelop(rel) ) {
515 case e_less:
516 case e_lesseq:
517 case e_greater:
518 case e_greatereq:
519 case e_equal:
520 case e_notequal:
521 /* Working on the left-hand-side */
522 len = RelationLength(rel,TRUE);
523 for( c = 1; c <= len; c++ ) {
524 struct relation_term *rt;
525 rt = (struct relation_term *)RelationTerm(rel,c,TRUE);
526 sp += 1-ArgsForRealToken(RelationTermType(rt));
527 apply_term_dimensions(rel,rt,sp-1,sp,&consistent,&wild);
528 } /* stack[0].d contains the dimensions of the lhs expression */
529
530 /* Now working on the right-hand_side */
531 len = RelationLength(rel,FALSE);
532 for( c = 1; c <= len; c++ ) {
533 struct relation_term *rt;
534 rt = (struct relation_term *) RelationTerm(rel,c,FALSE);
535 sp += 1-ArgsForRealToken(RelationTermType(rt));
536 apply_term_dimensions(rel,rt,sp-1,sp,&consistent,&wild);
537 } /* stack[1].d contains the dimensions of the rhs expression */
538
539 if( IsWild(&(stack[0].d)) || IsWild(&(stack[1].d)) ) {
540 if( IsWild(&(stack[0].d)) && !IsZero(&(stack[0])) ) {
541 if( !wild ) wild = TRUE;
542 if (GCDN) {
543 FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
544 FPRINTF(ASCERR," on left hand side.\n");
545 }
546 }
547 if( IsWild(&(stack[1].d)) && !IsZero(&(stack[1])) ) {
548 if( !wild ) wild = TRUE;
549 if (GCDN) {
550 FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
551 FPRINTF(ASCERR," on right hand side.\n");
552 }
553 }
554 }else{
555 if( CmpDimen(&(stack[0].d),&(stack[1].d)) ) {
556 if( consistent ) consistent = FALSE;
557 if (GCDN) {
558 FPRINTF(ASCERR,"ERROR: Relation has dimensions ");
559 WriteDimensions(ASCERR,&(stack[0].d));
560 FPRINTF(ASCERR," on left\n");
561 FPRINTF(ASCERR," and dimensions ");
562 WriteDimensions(ASCERR,&(stack[1].d));
563 FPRINTF(ASCERR," on right.\n");
564 }
565 }
566 }
567 break;
568 case e_maximize:
569 case e_minimize:
570 /* Working on the left-hand-side */
571 len = RelationLength(rel,TRUE);
572 for( c = 1; c <= len; c++ ) {
573 struct relation_term *rt;
574 rt = (struct relation_term *) RelationTerm(rel,c,TRUE);
575 sp += 1-ArgsForRealToken(RelationTermType(rt));
576 apply_term_dimensions(rel,rt,sp-1,sp,&consistent,&wild);
577 } /* stack[0].d contains the dimensions of the lhs expression */
578
579 if( IsWild(&(stack[0].d)) && !IsZero(&(stack[0])) ) {
580 if( !wild ) wild = TRUE;
581 if (GCDN) {
582 FPRINTF(ASCERR,"ERROR: Objective has wild dimensions.\n");
583 }
584 }
585 break;
586
587 default:
588 FPRINTF(ASCERR,"ERROR: Unknown relation type.\n");
589 if( consistent ) consistent = FALSE;
590 break;
591 }
592 CopyDimensions(&(stack[0].d),dimens);
593 ascfree(stack);
594 return( consistent && !wild );
595 }
596
597 /*------------------------------------------------------------------------------
598 CALCULATION FUNCTIONS
599 */
600
601 /** @NOTE ANY function calling RelationBranchEvaluator should set
602 glob_rel to point at the relation being evaluated. The calling
603 function should also set glob_rel = NULL when it is done.
604 */
605 static double RelationBranchEvaluator(struct relation_term *term)
606 {
607 assert(term != NULL);
608 switch(RelationTermType(term)) {
609 case e_func:
610 /* CONSOLE_DEBUG("Evaluating term using FuncEval..."); */
611 return FuncEval(TermFunc(term),
612 RelationBranchEvaluator(TermFuncLeft(term)) );
613 case e_var:
614 return TermVariable(glob_rel , term);
615 case e_int:
616 return (double)TermInteger(term);
617 case e_real:
618 return TermReal(term);
619 case e_zero:
620 return 0.0;
621 case e_plus:
622 return (RelationBranchEvaluator(TermBinLeft(term)) +
623 RelationBranchEvaluator(TermBinRight(term)));
624 case e_minus:
625 return (RelationBranchEvaluator(TermBinLeft(term)) -
626 RelationBranchEvaluator(TermBinRight(term)));
627 case e_times:
628 return (RelationBranchEvaluator(TermBinLeft(term)) *
629 RelationBranchEvaluator(TermBinRight(term)));
630 case e_divide:
631 return (RelationBranchEvaluator(TermBinLeft(term)) /
632 RelationBranchEvaluator(TermBinRight(term)));
633 case e_power:
634 case e_ipower:
635 return pow( RelationBranchEvaluator(TermBinLeft(term)) ,
636 RelationBranchEvaluator(TermBinRight(term)) );
637 case e_uminus:
638 return - RelationBranchEvaluator(TermBinLeft(term));
639 default:
640 FPRINTF(ASCERR, "error in RelationBranchEvaluator routine\n");
641 FPRINTF(ASCERR, "relation term type not recognized\n");
642 return 0.0;
643 }
644 }
645
646 /**
647 This function is passed a relation pointer, r, a pointer, pos, to a
648 position in the postfix version of the relation (0<=pos<length), and
649 a flag, lhs, telling whether we are interested in the left(=1) or
650 right(=0) side of the relation. This function will tranverse and
651 evaluate the subtree rooted at pos and will return the value as a
652 double. To do its evaluation, this function goes backwards through
653 the postfix representation of relation and calls itself at each
654 node--creating a stack of function calls.
655
656 @NOTE: This function changes the value of pos---to the position of
657 the deepest leaf visited
658 */
659 static double
660 RelationEvaluatePostfixBranch(CONST struct relation *r,
661 unsigned long *pos,
662 int lhs)
663 {
664 CONST struct relation_term *term; /* the current term */
665 CONST struct Func *funcptr; /* a pointer to a function */
666 double x, y; /* temporary values */
667
668 term = NewRelationTerm(r, *pos, lhs);
669 assert(term != NULL);
670 switch( RelationTermType(term) ) {
671 case e_zero:
672 case e_real:
673 return TermReal(term);
674 case e_int:
675 return TermInteger(term);
676 case e_var:
677 return TermVariable(r, term);
678 case e_plus:
679 (*pos)--;
680 y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y==right-side of '+' */
681 (*pos)--;
682 return RelationEvaluatePostfixBranch(r, pos, lhs) + y;
683 case e_minus:
684 (*pos)--;
685 y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y==right-side of '-' */
686 (*pos)--;
687 return RelationEvaluatePostfixBranch(r, pos, lhs) - y;
688 case e_times:
689 (*pos)--;
690 y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y==right-side of '*' */
691 (*pos)--;
692 return RelationEvaluatePostfixBranch(r, pos, lhs) * y;
693 case e_divide:
694 (*pos)--;
695 y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y is the divisor */
696 (*pos)--;
697 return RelationEvaluatePostfixBranch(r, pos, lhs) / y;
698 case e_power:
699 (*pos)--;
700 y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y is the power */
701 (*pos)--;
702 x = RelationEvaluatePostfixBranch(r, pos, lhs); /* x is the base */
703 return pow(x, y);
704 case e_ipower:
705 (*pos)--;
706 y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y is the power */
707 (*pos)--;
708 x = RelationEvaluatePostfixBranch(r, pos, lhs); /* x is the base */
709 return asc_ipow(x, (int)y);
710 case e_uminus:
711 (*pos)--;
712 return -1.0 * RelationEvaluatePostfixBranch(r, pos, lhs);
713 case e_func:
714 funcptr = TermFunc(term);
715 (*pos)--;
716 return FuncEval(funcptr, RelationEvaluatePostfixBranch(r, pos, lhs));
717 default:
718 ASC_PANIC("unrecognised relation type");
719 break;
720 }
721 }
722
723 #define RECURSE RelationEvaluatePostfixBranchSafe
724 static double
725 RelationEvaluatePostfixBranchSafe(CONST struct relation *r,
726 unsigned long *pos,
727 int lhs,
728 enum safe_err *serr
729 ){
730 CONST struct relation_term *term; /* the current term */
731 double x, y; /* temporary values */
732
733 term = NewRelationTerm(r, *pos, lhs);
734 assert(term != NULL);
735 switch( RelationTermType(term) ) {
736 case e_zero:
737 case e_real:
738 return TermReal(term);
739 case e_int:
740 return TermInteger(term);
741 case e_var:
742 return TermVariable(r, term);
743 case e_plus:
744 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = RHS of '+' */
745 (*pos)--; return safe_add_D0(RECURSE(r,pos,lhs,serr), y, serr);
746 case e_minus:
747 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = RHS of '-' */
748 (*pos)--; return safe_sub_D0(RECURSE(r,pos,lhs,serr), y, serr);
749 case e_times:
750 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = RHS of '*' */
751 (*pos)--; return safe_mul_D0(RECURSE(r,pos,lhs,serr), y, serr);
752 case e_divide:
753 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = the divisor (RHS of '/') */
754 (*pos)--; return safe_div_D0(RECURSE(r,pos,lhs,serr), y, serr);
755 case e_power:
756 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = the power (RHS of '^') */
757 (*pos)--; x = RECURSE(r, pos, lhs, serr); /* = the base */
758 return safe_pow_D0(x, y, serr);
759 case e_ipower:
760 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* y is the power (RHS of '^') */
761 (*pos)--; x = RECURSE(r, pos, lhs, serr); /* x is the base */
762 return safe_ipow_D0(x, (int)y, serr);
763 case e_uminus:
764 (*pos)--; return -1.0 * RECURSE(r, pos, lhs, serr);
765 case e_func:
766 (*pos)--; return FuncEvalSafe(TermFunc(term),RECURSE(r,pos,lhs,serr),serr);
767 default:
768 ASC_PANIC("Unknown relation term type");
769 }
770 }
771 #undef RECURSE
772
773 /**
774 Yet another function for calculating the residual of a relation.
775 This function also uses the postfix version of the relations, but it
776 manages a stack(array) of doubles and calculates the residual in this
777 stack; therefore the function is not recursive. If the funtion
778 cannot allocate memory for its stack, it returns 0.0, so there is
779 currently no way of knowing if this function failed.
780 */
781 static double
782 RelationEvaluateResidualPostfix(CONST struct relation *r)
783 {
784 unsigned long t; /* the current term in the relation r */
785 int lhs; /* looking at left(=1) or right(=0) hand side */
786 double *res_stack; /* the stack we use for evaluating the residual */
787 long s = -1; /* the top position in the stacks */
788 unsigned long length_lhs, length_rhs;
789 CONST struct relation_term *term;
790 CONST struct Func *funcptr;
791
792
793 length_lhs = RelationLength(r, 1);
794 length_rhs = RelationLength(r, 0);
795 if( (length_lhs+length_rhs) == 0 ) return 0.0;
796
797 /* create the stacks */
798 res_stack = tmpalloc_array((1+MAX(length_lhs,length_rhs)),double);
799 if( res_stack == NULL ) return 0.0;
800
801 lhs = 1;
802 t = 0;
803 while (1) {
804 if( lhs && (t >= length_lhs) ) {
805 /* finished processing left hand side, switch to right if it exists */
806 if( length_rhs ) {
807 lhs = t = 0;
808 }
809 else {
810 /* do not need to check for s>=0, since we know that
811 * (length_lhs+length_rhs>0) and that (length_rhs==0), the
812 * length_lhs must be > 0, thus s>=0
813 */
814 return (res_stack[s]);
815 }
816 }
817 else if( (!lhs) && (t >= length_rhs) ) {
818 /* finished processing right hand side */
819 if( length_lhs ) {
820 /* we know length_lhs and length_rhs are both > 0, since if
821 * length_rhs == 0, we would have exited above.
822 */
823 return (res_stack[s-1] - res_stack[s]);
824 }
825 else {
826 /* do not need to check for s>=0, since we know that
827 * (length_lhs+length_rhs>0) and that (length_lhs==0), the
828 * length_rhs must be > 0, thus s>=0
829 */
830 return (-1.0 * res_stack[s]);
831 }
832 }
833
834 term = NewRelationTerm(r, t++, lhs);
835 switch( RelationTermType(term) ) {
836 case e_zero:
837 s++; res_stack[s] = 0.0; break;
838 case e_real:
839 s++; res_stack[s] = TermReal(term); break;
840 case e_int:
841 s++; res_stack[s] = TermInteger(term); break;
842 case e_var:
843 s++; res_stack[s] = TermVariable(r, term); break;
844 case e_plus:
845 res_stack[s-1] += res_stack[s]; s--; break;
846 case e_minus:
847 res_stack[s-1] -= res_stack[s]; s--; break;
848 case e_times:
849 res_stack[s-1] *= res_stack[s]; s--; break;
850 case e_divide:
851 res_stack[s-1] /= res_stack[s]; s--; break;
852 case e_uminus:
853 res_stack[s] *= -1.0; break;
854 case e_power:
855 case e_ipower:
856 res_stack[s-1] = pow(res_stack[s-1], res_stack[s]); s--; break;
857 case e_func:
858 funcptr = TermFunc(term);
859 res_stack[s] = FuncEval(funcptr, res_stack[s]);
860 break;
861 default:
862 ASC_PANIC("Unknown relation term type");
863 }
864 }
865 }
866
867 /*------------------------------------------------------------------------------
868 GRADIENT AND DERIVATIVE CALCULATIONS
869 */
870
871 /**
872 Compute residuals and gradients for a relation (compiler-side routine)
873
874 @param r relation for which residual and gradients are to be calculated.
875 @param residual pointer to a double in which the residual shall be stored (must have already been allocated)
876 @param gradient pointer to an array of doubles where the gradients can be stored (must already have been allocated)
877
878 @return 0 on success, 1 on out-of-memeory.
879
880 Computes the gradients by maintaining n stacks, where
881 n = (number-of-variables-in-r + 1)
882 The +1 is for the residual. The stacks come from a single array which
883 this function gets by calling tmpalloc_array. Two macros are defined
884 to make referencing this array easier.
885 */
886 static int
887 RelationEvaluateResidualGradient(CONST struct relation *r,
888 double *residual,
889 double *gradient)
890 {
891 unsigned long t; /* the current term in the relation r */
892 unsigned long num_var; /* the number of variables in the relation r */
893 unsigned long v; /* the index of the variable we are looking at */
894 int lhs; /* looking at left(=1) or right(=0) hand side of r */
895 double *stacks; /* the memory for the stacks */
896 unsigned long stack_height; /* height of each stack */
897 long s = -1; /* the top position in the stacks */
898 double temp, temp2; /* temporary variables to speed gradient calculatns */
899 unsigned long length_lhs, length_rhs;
900 CONST struct relation_term *term;
901 CONST struct Func *fxnptr;
902
903 num_var = NumberVariables(r);
904 length_lhs = RelationLength(r, 1);
905 length_rhs = RelationLength(r, 0);
906 if( (length_lhs + length_rhs) == 0 ) {
907 for( v = 0; v < num_var; v++ ) gradient[v] = 0.0;
908 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Relation with no LHS and no RHS: returning residual 0");
909 *residual = 0.0;
910 return 0;
911 }
912 else {
913 stack_height = 1 + MAX(length_lhs,length_rhs);
914 }
915
916 /* create the stacks */
917 stacks = tmpalloc_array(((num_var+1)*stack_height),double);
918 if( stacks == NULL ) return 1;
919
920 #define res_stack(s) stacks[(s)]
921 #define grad_stack(v,s) stacks[((v)*stack_height)+(s)]
922
923 lhs = 1;
924 t = 0;
925 while(1) {
926 if( lhs && (t >= length_lhs) ) {
927 /* need to switch to the right hand side--if it exists */
928 if( length_rhs ) {
929 lhs = t = 0;
930 }
931 else {
932 /* Set the pointers we were passed to the tops of the stacks.
933 * We do not need to check for s>=0, since we know that
934 * (length_lhs+length_rhs>0) and that (length_rhs==0), the
935 * length_lhs must be > 0, thus s>=0
936 */
937 for( v = 1; v <= num_var; v++ ) gradient[v-1] = grad_stack(v,s);
938 *residual = res_stack(s);
939 return 0;
940 }
941 }
942 else if( (!lhs) && (t >= length_rhs) ) {
943 /* we have processed both sides, quit */
944 if( length_lhs ) {
945 /* Set the pointers we were passed to lhs - rhs
946 * We know length_lhs and length_rhs are both > 0, since if
947 * length_rhs == 0, we would have exited above.
948 */
949 for( v = 1; v <= num_var; v++ ) {
950 gradient[v-1] = grad_stack(v,s-1) - grad_stack(v,s);
951 }
952 *residual = res_stack(s-1) - res_stack(s);
953 return 0;
954 }
955 else {
956 /* Set the pointers we were passed to -1.0 * top of stacks.
957 * We do not need to check for s>=0, since we know that
958 * (length_lhs+length_rhs>0) and that (length_lhs==0), the
959 * length_rhs must be > 0, thus s>=0
960 */
961 for( v = 1; v <= num_var; v++ ) {
962 gradient[v-1] = -grad_stack(v,s);
963 }
964 *residual = -res_stack(s);
965 return 0;
966 }
967 }
968
969 term = NewRelationTerm(r, t++, lhs);
970 switch( RelationTermType(term) ) {
971 case e_zero:
972 s++;
973 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
974 res_stack(s) = 0.0;
975 break;
976 case e_real:
977 s++;
978 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
979 res_stack(s) = TermReal(term);
980 break;
981 case e_int:
982 s++;
983 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
984 res_stack(s) = TermInteger(term);
985 break;
986 case e_var:
987 s++;
988 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
989 grad_stack(TermVarNumber(term),s) = 1.0;
990 res_stack(s) = TermVariable(r, term);
991 break;
992 case e_plus:
993 /* d(u+v) = du + dv */
994 for( v = 1; v <= num_var; v++ ) grad_stack(v,s-1) += grad_stack(v,s);
995 res_stack(s-1) += res_stack(s);
996 s--;
997 break;
998 case e_minus:
999 /* d(u-v) = du - dv */
1000 for( v = 1; v <= num_var; v++ ) grad_stack(v,s-1) -= grad_stack(v,s);
1001 res_stack(s-1) -= res_stack(s);
1002 s--;
1003 break;
1004 case e_times:
1005 /* d(u*v) = u*dv + v*du */
1006 for( v = 1; v <= num_var; v++ ) {
1007 grad_stack(v,s-1) = ((res_stack(s-1) * grad_stack(v,s)) +
1008 (res_stack(s) * grad_stack(v,s-1)));
1009 }
1010 res_stack(s-1) *= res_stack(s);
1011 s--;
1012 break;
1013 case e_divide:
1014 /* d(u/v) = du/v - u*dv/(v^2) = (1/v) * [du - (u/v)*dv] */
1015 res_stack(s) = 1.0 / res_stack(s); /* 1/v */
1016 res_stack(s-1) *= res_stack(s); /* u/v */
1017 for( v = 1; v <= num_var; v++ ) {
1018 grad_stack(v,s-1) = (res_stack(s) *
1019 (grad_stack(v,s-1) -
1020 (res_stack(s-1) * grad_stack(v,s))));
1021 }
1022 s--;
1023 break;
1024 case e_uminus:
1025 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = -grad_stack(v,s);
1026 res_stack(s) = -res_stack(s);
1027 break;
1028 case e_power:
1029 /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1030 /* First compute: v*u^(v-1) */
1031 temp = res_stack(s) * pow( res_stack(s-1), (res_stack(s) - 1.0) );
1032 /* Now compute: ln(u) */
1033 temp2 = FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1034 /* Next compute: u^v */
1035 res_stack(s-1) = pow(res_stack(s-1), res_stack(s));
1036 /* Compute: [ln(u)] * [u^v] */
1037 temp2 *= res_stack(s-1);
1038 /* Finally, compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1039 for( v = 1; v <= num_var; v++ ) {
1040 grad_stack(v,s-1) = ((temp * grad_stack(v,s-1)) +
1041 (temp2 * grad_stack(v,s)));
1042 }
1043 s--;
1044 break;
1045 case e_ipower:
1046 /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1047 /* First compute: v*u^(v-1) */
1048 temp = asc_d1ipow( res_stack(s-1), ((int)res_stack(s)) );
1049 /* Now compute: ln(u) */
1050 temp2 = FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1051 /* Next compute: u^v */
1052 res_stack(s-1) = asc_ipow( res_stack(s-1), ((int)res_stack(s)) );
1053 /* Compute: [ln(u)] * [u^v] */
1054 temp2 *= res_stack(s-1);
1055 /* Finally, compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1056 for( v = 1; v <= num_var; v++ ) {
1057 grad_stack(v,s-1) = ((temp * grad_stack(v,s-1)) +
1058 (temp2 * grad_stack(v,s)));
1059 }
1060 s--;
1061 break;
1062 case e_func:
1063 /*
1064 funcptr = TermFunc(term);
1065 for (v = 0; v < num_var; v++) {
1066 grad_stack(v,s) = FuncDeriv(funcptr, grad_stack(v,s));
1067 }
1068 res_stack(s) = FuncEval(funcptr, res_stack(s)); */
1069 fxnptr = TermFunc(term);
1070 temp = FuncDeriv( fxnptr, res_stack(s) );
1071 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) *= temp;
1072 res_stack(s) = FuncEval( fxnptr, res_stack(s) );
1073 break;
1074 default:
1075 ASC_PANIC("Unknown relation term type");
1076 break;
1077 }
1078 }
1079 #undef grad_stack
1080 #undef res_stack
1081 }
1082
1083 static int
1084 RelationEvaluateResidualGradientSafe(CONST struct relation *r,
1085 double *residual,
1086 double *gradient,
1087 enum safe_err *serr)
1088 {
1089 unsigned long t; /* the current term in the relation r */
1090 unsigned long num_var; /* the number of variables in the relation r */
1091 unsigned long v; /* the index of the variable we are looking at */
1092 int lhs; /* looking at left(=1) or right(=0) hand side of r */
1093 double *stacks; /* the memory for the stacks */
1094 unsigned long stack_height; /* height of each stack */
1095 long s = -1; /* the top position in the stacks */
1096 double temp, temp2; /* temporary variables to speed gradient calculatns */
1097 unsigned long length_lhs, length_rhs;
1098 CONST struct relation_term *term;
1099 CONST struct Func *fxnptr;
1100
1101 num_var = NumberVariables(r);
1102 length_lhs = RelationLength(r, 1);
1103 length_rhs = RelationLength(r, 0);
1104 if( (length_lhs + length_rhs) == 0 ) {
1105 for( v = 0; v < num_var; v++ ) gradient[v] = 0.0;
1106 *residual = 0.0;
1107 return 0;
1108 }
1109 else {
1110 stack_height = 1 + MAX(length_lhs,length_rhs);
1111 }
1112
1113 /* create the stacks */
1114 stacks = tmpalloc_array(((num_var+1)*stack_height),double);
1115 if( stacks == NULL ) return 1;
1116
1117 #define res_stack(s) stacks[(s)]
1118 #define grad_stack(v,s) stacks[((v)*stack_height)+(s)]
1119
1120 lhs = 1;
1121 t = 0;
1122 while(1) {
1123 if( lhs && (t >= length_lhs) ) {
1124 /* need to switch to the right hand side--if it exists */
1125 if( length_rhs ) {
1126 lhs = t = 0;
1127 }
1128 else {
1129 /* Set the pointers we were passed to the tops of the stacks.
1130 * We do not need to check for s>=0, since we know that
1131 * (length_lhs+length_rhs>0) and that (length_rhs==0), the
1132 * length_lhs must be > 0, thus s>=0
1133 */
1134 for( v = 1; v <= num_var; v++ ) gradient[v-1] = grad_stack(v,s);
1135 *residual = res_stack(s);
1136 return 0;
1137 }
1138 }
1139 else if( (!lhs) && (t >= length_rhs) ) {
1140 /* we have processed both sides, quit */
1141 if( length_lhs ) {
1142 /* Set the pointers we were passed to lhs - rhs
1143 * We know length_lhs and length_rhs are both > 0, since if
1144 * length_rhs == 0, we would have exited above.
1145 */
1146 for( v = 1; v <= num_var; v++ ) {
1147 gradient[v-1] = safe_sub_D0(grad_stack(v,s-1),grad_stack(v,s),serr);
1148 }
1149 *residual = safe_sub_D0(res_stack(s-1), res_stack(s), serr);
1150 return 0;
1151 }
1152 else {
1153 /* Set the pointers we were passed to -1.0 * top of stacks.
1154 * We do not need to check for s>=0, since we know that
1155 * (length_lhs+length_rhs>0) and that (length_lhs==0), the
1156 * length_rhs must be > 0, thus s>=0
1157 */
1158 for( v = 1; v <= num_var; v++ ) {
1159 gradient[v-1] = -grad_stack(v,s);
1160 }
1161 *residual = -res_stack(s);
1162 return 0;
1163 }
1164 }
1165
1166 term = NewRelationTerm(r, t++, lhs);
1167 switch( RelationTermType(term) ) {
1168 case e_zero:
1169 s++;
1170 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1171 res_stack(s) = 0.0;
1172 break;
1173 case e_real:
1174 s++;
1175 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1176 res_stack(s) = TermReal(term);
1177 break;
1178 case e_int:
1179 s++;
1180 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1181 res_stack(s) = TermInteger(term);
1182 break;
1183 case e_var:
1184 s++;
1185 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1186 grad_stack(TermVarNumber(term),s) = 1.0;
1187 res_stack(s) = TermVariable(r, term);
1188 break;
1189 case e_plus:
1190 /* d(u+v) = du + dv */
1191 for( v = 1; v <= num_var; v++ ) {
1192 grad_stack(v,s-1)=safe_add_D0(grad_stack(v,s-1),grad_stack(v,s),serr);
1193 }
1194 res_stack(s-1) = safe_add_D0(res_stack(s-1),res_stack(s),serr);
1195 s--;
1196 break;
1197 case e_minus:
1198 /* d(u-v) = du - dv */
1199 for( v = 1; v <= num_var; v++ ) {
1200 grad_stack(v,s-1)=safe_sub_D0(grad_stack(v,s-1),grad_stack(v,s),serr);
1201 }
1202 res_stack(s-1) = safe_sub_D0(res_stack(s-1),res_stack(s),serr);
1203 s--;
1204 break;
1205 case e_times:
1206 /* d(u*v) = u*dv + v*du */
1207 for( v = 1; v <= num_var; v++ ) {
1208 grad_stack(v,s-1) =
1209 safe_add_D0(safe_mul_D0(res_stack(s-1),grad_stack(v,s),serr),
1210 safe_mul_D0(res_stack(s),grad_stack(v,s-1),serr),
1211 serr);
1212 }
1213 res_stack(s-1) = safe_mul_D0(res_stack(s-1),res_stack(s),serr);
1214 s--;
1215 break;
1216 case e_divide:
1217 /* d(u/v) = du/v - u*dv/(v^2) = (1/v) * [du - (u/v)*dv] */
1218 res_stack(s) = safe_rec(res_stack(s),serr); /* 1/v */
1219 res_stack(s-1) = safe_mul_D0(res_stack(s-1),res_stack(s),serr); /* u/v */
1220 for( v = 1; v <= num_var; v++ ) {
1221 grad_stack(v,s-1) =
1222 safe_mul_D0(res_stack(s),
1223 safe_sub_D0(grad_stack(v,s-1),
1224 safe_mul_D0(res_stack(s-1),
1225 grad_stack(v,s),
1226 serr),serr),serr);
1227 }
1228 s--;
1229 break;
1230 case e_uminus:
1231 for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = -grad_stack(v,s);
1232 res_stack(s) = -res_stack(s);
1233 break;
1234 case e_power:
1235 /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1236 /* v*u^(v-1) */
1237 temp = safe_pow_D1( res_stack(s-1), res_stack(s), 0, serr );
1238 /* ln(u)*u^v */
1239 temp2 = safe_pow_D1( res_stack(s-1), res_stack(s), 1, serr );
1240 /* Compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1241 for( v = 1; v <= num_var; v++ ) {
1242 grad_stack(v,s-1) =
1243 safe_add_D0(safe_mul_D0(temp, grad_stack(v,s-1), serr),
1244 safe_mul_D0(temp2, grad_stack(v,s), serr), serr);
1245 }
1246 /* u^v */
1247 res_stack(s-1) = safe_pow_D0( res_stack(s-1), res_stack(s), serr );
1248 s--;
1249 break;
1250 case e_ipower:
1251 /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1252 /* v*u^(v-1) */
1253 temp = safe_ipow_D1( res_stack(s-1), res_stack(s), 0, serr );
1254 /* ln(u)*u^v */
1255 temp2 = safe_ipow_D1( res_stack(s-1), res_stack(s), 1, serr );
1256 /* Compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1257 for( v = 1; v <= num_var; v++ ) {
1258 grad_stack(v,s-1) =
1259 safe_add_D0(safe_mul_D0(temp, grad_stack(v,s-1), serr),
1260 safe_mul_D0(temp2, grad_stack(v,s), serr), serr);
1261 }
1262 /* Next compute: u^v */
1263 res_stack(s-1) = safe_ipow_D0( res_stack(s-1), res_stack(s), serr );
1264 s--;
1265 break;
1266 case e_func:
1267 fxnptr = TermFunc(term);
1268 temp = FuncDerivSafe( fxnptr, res_stack(s), serr);
1269 for( v = 1; v <= num_var; v++ ) {
1270 grad_stack(v,s) = safe_mul_D0( grad_stack(v,s), temp, serr );
1271 }
1272 res_stack(s) = FuncEvalSafe( fxnptr, res_stack(s), serr);
1273 break;
1274 default:
1275 ASC_PANIC("Unknown relation term type");
1276 }
1277 }
1278 #undef grad_stack
1279 #undef res_stack
1280 }
1281
1282 /**
1283 This function evaluates and returns the derivative of the
1284 relation r with respect to the variable whose index is pos.
1285 This function assumes r exists and that pos is within the proper range.
1286 The function computes the gradients by maintaining 2 stacks, one for
1287 the residual and one for the derivative. The stacks come from a
1288 single array which this gets by calling tmpalloc_array. Two macros
1289 are defined to make referencing this array easier. Of the malloc fails,
1290 this function returns 0.0, so there is currently no way to know if
1291 the function failed.
1292 */
1293 static double
1294 RelationEvaluateDerivative(CONST struct relation *r,
1295 unsigned long pos)
1296 {
1297 unsigned long t; /* the current term in the relation r */
1298 int lhs; /* looking at left(=1) or right(=0) hand side of r */
1299 double *stacks; /* the memory for the stacks */
1300 unsigned long stack_height; /* height of each stack */
1301 long s = -1; /* the top position in the stacks */
1302 unsigned long length_lhs, length_rhs;
1303 CONST struct relation_term *term;
1304 CONST struct Func *fxnptr;
1305
1306 length_lhs = RelationLength(r, 1);
1307 length_rhs = RelationLength(r, 0);
1308 if( (length_lhs + length_rhs) == 0 ) {
1309 return 0.0;
1310 }
1311 else {
1312 stack_height = 1 + MAX(length_lhs,length_rhs);
1313 }
1314
1315 /* create the stacks */
1316 stacks = tmpalloc_array((2*stack_height),double);
1317 if( stacks == NULL ) return 0.0;
1318
1319 #define res_stack(s) stacks[(s)]
1320 #define grad_stack(s) stacks[stack_height+(s)]
1321
1322 lhs = 1;
1323 t = 0;
1324 while(1) {
1325 if( lhs && (t >= length_lhs) ) {
1326 /* need to switch to the right hand side--if it exists */
1327 if( length_rhs ) {
1328 lhs = t = 0;
1329 }
1330 else {
1331 /* do not need to check for s>=0, since we know that
1332 * (length_lhs+length_rhs>0) and that (length_rhs==0), the
1333 * length_lhs must be > 0, thus s>=0
1334 */
1335 return grad_stack(s);
1336 }
1337 }
1338 else if( (!lhs) && (t >= length_rhs) ) {
1339 /* we have processed both sides, quit */
1340 if( length_lhs ) {
1341 /* we know length_lhs and length_rhs are both > 0, since if
1342 * length_rhs == 0, we would have exited above.
1343 */
1344 return (grad_stack(s-1) - grad_stack(s));
1345 }
1346 else {
1347 /* do not need to check for s>=0, since we know that
1348 * (length_lhs+length_rhs>0) and that (length_lhs==0), the
1349 * length_rhs must be > 0, thus s>=0
1350 */
1351 return (-1.0 * grad_stack(s));
1352 }
1353 }
1354
1355 term = NewRelationTerm(r, t++, lhs);
1356 switch( RelationTermType(term) ) {
1357 case e_zero:
1358 s++;
1359 grad_stack(s) = res_stack(s) = 0.0;
1360 break;
1361 case e_real:
1362 s++;
1363 grad_stack(s) = 0.0;
1364 res_stack(s) = TermReal(term);
1365 break;
1366 case e_int:
1367 s++;
1368 grad_stack(s) = 0.0;
1369 res_stack(s) = TermInteger(term);
1370 break;
1371 case e_var:
1372 s++;
1373 grad_stack(s) = ( (pos == TermVarNumber(term)) ? 1.0 : 0.0 );
1374 res_stack(s) = TermVariable(r, term);
1375 break;
1376 case e_plus:
1377 /* d(u+v) = du + dv */
1378 grad_stack(s-1) += grad_stack(s);
1379 res_stack(s-1) += res_stack(s);
1380 s--;
1381 break;
1382 case e_minus:
1383 /* d(u-v) = du - dv */
1384 grad_stack(s-1) -= grad_stack(s);
1385 res_stack(s-1) -= res_stack(s);
1386 s--;
1387 break;
1388 case e_times:
1389 /* d(u*v) = u*dv + v*du */
1390 grad_stack(s-1) = ((res_stack(s-1) * grad_stack(s)) +
1391 (res_stack(s) * grad_stack(s-1)));
1392 res_stack(s-1) *= res_stack(s);
1393 s--;
1394 break;
1395 case e_divide:
1396 /* d(u/v) = du/v - u*dv/(v^2) = [du - (u/v)*dv]/v */
1397 res_stack(s-1) = res_stack(s-1) / res_stack(s);
1398 grad_stack(s-1) = ((grad_stack(s-1) - (res_stack(s-1) * grad_stack(s))) /
1399 res_stack(s));
1400 s--;
1401 break;
1402 case e_uminus:
1403 grad_stack(s) = -grad_stack(s);
1404 res_stack(s) = -res_stack(s);
1405 break;
1406 case e_power:
1407 /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1408 /* First we compute: v*u^(v-1)*du */
1409 grad_stack(s-1) *= (pow(res_stack(s-1), (res_stack(s) - 1.0)) *
1410 res_stack(s));
1411 /* Now compute: ln(u)*dv */
1412 grad_stack(s) *= FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1413 /* Next compute: u^v */
1414 res_stack(s-1) = pow( res_stack(s-1), res_stack(s) );
1415 /* Finally, compute: [v*u^(v-1)*du] + [u^v] * [ln(u)*dv] */
1416 grad_stack(s-1) += (res_stack(s-1) * grad_stack(s));
1417 s--;
1418 break;
1419 case e_ipower:
1420 /* d(x^y) = y * dx * x^(y-1) + ln(x) * dy * x^y */
1421 /* First we compute: v*u^(v-1)*du */
1422 grad_stack(s-1) *= asc_d1ipow( res_stack(s-1), ((int)res_stack(s)) );
1423 /* Now compute: ln(u)*dv */
1424 grad_stack(s) *= FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1425 /* Next compute: u^v */
1426 res_stack(s-1) = asc_ipow( res_stack(s-1), ((int)res_stack(s)) );
1427 /* Finally, compute: [v*u^(v-1)*du] + [u^v] * [ln(u)*dv] */
1428 grad_stack(s-1) += (res_stack(s-1) * grad_stack(s));
1429 s--;
1430 break;
1431 case e_func:
1432 fxnptr = TermFunc(term);
1433 grad_stack(s) *= FuncDeriv( fxnptr, res_stack(s) );
1434 res_stack(s) = FuncEval( fxnptr, res_stack(s) );
1435 break;
1436 default:
1437 ASC_PANIC("Unknown relation term type");
1438 break;
1439 }
1440 }
1441 #undef grad_stack
1442 #undef res_stack
1443 }
1444
1445 static double
1446 RelationEvaluateDerivativeSafe(CONST struct relation *r,
1447 unsigned long pos,
1448 enum safe_err *serr)
1449 {
1450 unsigned long t; /* the current term in the relation r */
1451 int lhs; /* looking at left(=1) or right(=0) hand side of r */
1452 double *stacks; /* the memory for the stacks */
1453 unsigned long stack_height; /* height of each stack */
1454 long s = -1; /* the top position in the stacks */
1455 unsigned long length_lhs, length_rhs;
1456 CONST struct relation_term *term;
1457 CONST struct Func *fxnptr;
1458
1459 length_lhs = RelationLength(r, 1);
1460 length_rhs = RelationLength(r, 0);
1461 if( (length_lhs + length_rhs) == 0 ) {
1462 return 0.0;
1463 }
1464 else {
1465 stack_height = 1 + MAX(length_lhs,length_rhs);
1466 }
1467
1468 /* create the stacks */
1469 stacks = tmpalloc_array((2*stack_height),double);
1470 if( stacks == NULL ) return 0.0;
1471
1472 #define res_stack(s) stacks[(s)]
1473 #define grad_stack(s) stacks[stack_height+(s)]
1474
1475 lhs = 1;
1476 t = 0;
1477 while(1) {
1478 if( lhs && (t >= length_lhs) ) {
1479 /* need to switch to the right hand side--if it exists */
1480 if( length_rhs ) {
1481 lhs = t = 0;
1482 }
1483 else {
1484 /* do not need to check for s>=0, since we know that
1485 * (length_lhs+length_rhs>0) and that (length_rhs==0), the
1486 * length_lhs must be > 0, thus s>=0
1487 */
1488 return grad_stack(s);
1489 }
1490 }
1491 else if( (!lhs) && (t >= length_rhs) ) {
1492 /* we have processed both sides, quit */
1493 if( length_lhs ) {
1494 /* we know length_lhs and length_rhs are both > 0, since if
1495 * length_rhs == 0, we would have exited above.
1496 */
1497 return safe_sub_D0(grad_stack(s-1), grad_stack(s), serr);
1498 }
1499 else {
1500 /* do not need to check for s>=0, since we know that
1501 * (length_lhs+length_rhs>0) and that (length_lhs==0), the
1502 * length_rhs must be > 0, thus s>=0
1503 */
1504 return (-grad_stack(s));
1505 }
1506 }
1507
1508 term = NewRelationTerm(r, t++, lhs);
1509 switch( RelationTermType(term) ) {
1510 case e_zero:
1511 s++;
1512 grad_stack(s) = res_stack(s) = 0.0;
1513 break;
1514 case e_real:
1515 s++;
1516 grad_stack(s) = 0.0;
1517 res_stack(s) = TermReal(term);
1518 break;
1519 case e_int:
1520 s++;
1521 grad_stack(s) = 0.0;
1522 res_stack(s) = TermInteger(term);
1523 break;
1524 case e_var:
1525 s++;
1526 grad_stack(s) = ( (pos == TermVarNumber(term)) ? 1.0 : 0.0 );
1527 res_stack(s) = TermVariable(r, term);
1528 break;
1529 case e_plus:
1530 /* d(u+v) = du + dv */
1531 grad_stack(s-1) = safe_add_D0( grad_stack(s-1), grad_stack(s), serr );
1532 res_stack(s-1) = safe_add_D0( res_stack(s-1), res_stack(s), serr );
1533 s--;
1534 break;
1535 case e_minus:
1536 /* d(u-v) = du - dv */
1537 grad_stack(s-1) = safe_sub_D0( grad_stack(s-1), grad_stack(s), serr );
1538 res_stack(s-1) = safe_sub_D0( res_stack(s-1), res_stack(s), serr );
1539 s--;
1540 break;
1541 case e_times:
1542 /* d(u*v) = u*dv + v*du */
1543 grad_stack(s-1) =
1544 safe_add_D0(safe_mul_D0( res_stack(s-1), grad_stack(s), serr),
1545 safe_mul_D0( res_stack(s), grad_stack(s-1), serr), serr);
1546 res_stack(s-1) = safe_mul_D0( res_stack(s-1), res_stack(s), serr );
1547 s--;
1548 break;
1549 case e_divide:
1550 /* d(u/v) = du/v - u*dv/(v^2) = [du - (u/v)*dv]/v */
1551 res_stack(s-1) = safe_div_D0( res_stack(s-1), res_stack(s), serr);
1552 grad_stack(s-1) =
1553 safe_div_D0(safe_sub_D0(grad_stack(s-1),
1554 safe_mul_D0(res_stack(s-1),grad_stack(s),serr),
1555 serr),
1556 res_stack(s),serr);
1557 s--;
1558 break;
1559 case e_uminus:
1560 grad_stack(s) = -grad_stack(s);
1561 res_stack(s) = -res_stack(s);
1562 break;
1563 case e_power:
1564 /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1565 grad_stack(s-1) =
1566 safe_add_D0(safe_mul_D0(safe_pow_D1(res_stack(s-1),
1567 res_stack(s),0,serr),
1568 grad_stack(s-1), serr),
1569 safe_mul_D0(safe_pow_D1(res_stack(s-1),
1570 res_stack(s),1,serr),
1571 grad_stack(s),serr), serr);
1572 /* u^v */
1573 res_stack(s-1) = safe_pow_D0( res_stack(s-1), res_stack(s), serr);
1574 s--;
1575 break;
1576 case e_ipower:
1577 /* d(x^y) = y * dx * x^(y-1) + ln(x) * dy * x^y */
1578 grad_stack(s-1) =
1579 safe_add_D0(safe_mul_D0(safe_ipow_D1(res_stack(s-1),
1580 res_stack(s),0,serr),
1581 grad_stack(s-1), serr),
1582 safe_mul_D0(safe_ipow_D1(res_stack(s-1),
1583 res_stack(s),1,serr),
1584 grad_stack(s),serr), serr);
1585 /* u^v */
1586 res_stack(s-1) = safe_ipow_D0( res_stack(s-1), res_stack(s), serr);
1587 s--;
1588 break;
1589 case e_func:
1590 fxnptr = TermFunc(term);
1591 grad_stack(s) = safe_mul_D0(FuncDerivSafe( fxnptr, res_stack(s), serr ),
1592 grad_stack(s), serr);
1593 res_stack(s) = FuncEvalSafe( fxnptr, res_stack(s), serr);
1594 break;
1595 default:
1596 ASC_PANIC("Unknown relation term type");
1597 }
1598 }
1599 #undef grad_stack
1600 #undef res_stack
1601 }
1602
1603 /*------------------------------------------------------------------------------
1604 RELATION & RELATION TERM QUERIES FUNCTIONS (FOR USE BY EXTERNAL CODE)
1605 */
1606
1607 /**
1608 @return =, <, >, etc, etc. not e_token, e_glassbox, etc
1609 */
1610 enum Expr_enum RelationRelop(CONST struct relation *rel)
1611 {
1612 AssertAllocatedMemory(rel,sizeof(struct relation));
1613 return rel->share->s.relop;
1614 }
1615
1616 /**
1617 This query only applies to TokenRelations and OpcodeRelations.
1618 */
1619 unsigned long RelationLength(CONST struct relation *rel, int lhs){
1620 assert(rel!=NULL);
1621 AssertAllocatedMemory(rel,sizeof(struct relation));
1622 if(lhs){
1623 if(RTOKEN(rel).lhs){
1624 asc_assert(RTOKEN(rel).lhs_len >= 0);
1625 return RTOKEN(rel).lhs_len;
1626 }
1627 return 0;
1628 }
1629 if(RTOKEN(rel).rhs){
1630 asc_assert(RTOKEN(rel).rhs_len >= 0);
1631 return RTOKEN(rel).rhs_len;
1632 }
1633 return 0;
1634 }
1635
1636 struct gl_list_t *RelationBlackBoxArgNames(CONST struct relation *rel)
1637 {
1638 assert(rel!=NULL);
1639 return RelationBlackBoxCache(rel)->argListNames;
1640 }
1641
1642 struct Name *RelationBlackBoxDataName(CONST struct relation *rel)
1643 {
1644 assert(rel!=NULL);
1645 return RelationBlackBoxCache(rel)->dataName;
1646 }
1647
1648 struct BlackBoxCache * RelationBlackBoxCache(CONST struct relation *rel)
1649 {
1650 assert(rel!=NULL);
1651 AssertAllocatedMemory(rel,sizeof(struct relation));
1652 return ((struct BlackBoxData *) (rel->externalData))->common;
1653 }
1654
1655 struct BlackBoxData * RelationBlackBoxData(CONST struct relation *rel)
1656 {
1657 assert(rel!=NULL);
1658 AssertAllocatedMemory(rel,sizeof(struct relation));
1659 return (struct BlackBoxData *) rel->externalData;
1660 }
1661
1662 /*
1663 This query only applies to TokenRelations. It assumes the
1664 user still thinks tokens number from [1..len].
1665 */
1666 CONST struct relation_term *RelationTerm(CONST struct relation *rel,
1667 unsigned long int pos, int lhs)
1668 {
1669 assert(rel!=NULL);
1670 AssertAllocatedMemory(rel,sizeof(struct relation));
1671 if (lhs){
1672 if (RTOKEN(rel).lhs)
1673 return A_TERM(&(RTOKEN(rel).lhs[pos-1]));
1674 else return NULL;
1675 }
1676 else{
1677 if (RTOKEN(rel).rhs)
1678 return A_TERM(&(RTOKEN(rel).rhs[pos-1]));
1679 else return NULL;
1680 }
1681 }
1682
1683 /**
1684 This query only applies to TokenRelations. It assumes the
1685 clued user thinks tokens number from [0..len-1], which they do.
1686 */
1687 CONST struct relation_term
1688 *NewRelationTermF(CONST struct relation *rel, unsigned long pos, int lhs)
1689 {
1690 assert(rel!=NULL);
1691 AssertAllocatedMemory(rel,sizeof(struct relation));
1692 if (lhs){
1693 if (RTOKEN(rel).lhs != NULL)
1694 return A_TERM(&(RTOKEN(rel).lhs[pos]));
1695 else return NULL;
1696 }else{
1697 if (RTOKEN(rel).rhs != NULL)
1698 return A_TERM(&(RTOKEN(rel).rhs[pos]));
1699 else return NULL;
1700 }
1701 }
1702
1703 /**
1704 This query only applies to sides from TokenRelations. It assumes the
1705 clued user thinks tokens number from [0..len-1], which they do,
1706 and that the side came from a token relation instance.
1707 */
1708 CONST struct relation_term
1709 *RelationSideTermF(CONST union RelationTermUnion *side, unsigned long pos)
1710 {
1711 assert(side!=NULL);
1712 return A_TERM(&(side[pos]));
1713 }
1714
1715 /**
1716 This query only applies to TokenRelations. It assumes the
1717 clued user thinks tokens number from [0..len-1], which they do.
1718 */
1719 enum Expr_enum RelationTermTypeF(CONST struct relation_term *term)
1720 {
1721 AssertMemory(term);
1722 return term->t;
1723 }
1724
1725 unsigned long TermVarNumber(CONST struct relation_term *term)
1726 {
1727 assert(term&&term->t == e_var);
1728 AssertMemory(term);
1729 return V_TERM(term)->varnum;
1730 }
1731
1732 long TermInteger(CONST struct relation_term *term){
1733 assert(term&&(term->t==e_int));
1734 AssertMemory(term);
1735 return I_TERM(term)->ivalue;
1736 }
1737
1738 double TermReal(CONST struct relation_term *term){
1739 assert(term&&( term->t==e_real || term->t==e_zero));
1740 AssertMemory(term);
1741 return R_TERM(term)->value;
1742 }
1743
1744 double
1745 TermVariable(CONST struct relation *rel, CONST struct relation_term *term){
1746 return
1747 RealAtomValue((struct Instance*)RelationVariable(rel,TermVarNumber(term)));
1748 }
1749
1750 CONST dim_type *TermDimensions(CONST struct relation_term *term){
1751 assert( term && (term->t==e_real || term->t == e_int || term->t == e_zero) );
1752 AssertMemory(term);
1753 if (term->t==e_real) return R_TERM(term)->dimensions;
1754 if (term->t==e_int) return Dimensionless();
1755 if (term->t==e_zero) return WildDimension();
1756 return NULL;
1757 }
1758
1759 CONST struct Func *TermFunc(CONST struct relation_term *term){
1760 assert(term&&(term->t == e_func));
1761 AssertMemory(term);
1762 return F_TERM(term)->fptr;
1763 }
1764
1765 struct relation_term *RelationINF_Lhs(CONST struct relation *rel){
1766 return RTOKEN(rel).lhs_term;
1767 }
1768
1769 struct relation_term *RelationINF_Rhs(CONST struct relation *rel){
1770 return RTOKEN(rel).rhs_term;
1771 }
1772
1773 struct ExternalFunc *RelationBlackBoxExtFunc(CONST struct relation *rel)
1774 {
1775 assert(rel!=NULL);
1776 return RelationBlackBoxCache(rel)->efunc;
1777 }
1778
1779 /*------------------------------------------------------------------------------
1780 GLASSBOX RELATION QUERIES
1781
1782 For picking apart a GlassBox relation.
1783 */
1784
1785 struct ExternalFunc *RelationGlassBoxExtFunc(CONST struct relation *rel){
1786 assert(rel!=NULL);
1787 return RGBOX(rel).efunc;
1788 }
1789
1790 int GlassBoxRelIndex(CONST struct relation *rel){
1791 assert(rel!=NULL);
1792 return RGBOX(rel).index;
1793 }
1794
1795 int *GlassBoxArgs(CONST struct relation *rel){
1796 assert(rel!=NULL);
1797 return RGBOX(rel).args;
1798 }
1799
1800
1801 /*------------------------------------------------------------------------------
1802 GENERAL RELATION QUERIES
1803
1804 Not all of these queries may be applied to all relation types.
1805 Those that cannot be are so marked.
1806 */
1807
1808 CONST struct gl_list_t *RelationVarList(CONST struct relation *rel){
1809 return (CONST struct gl_list_t *)rel->vars;
1810 }
1811
1812 dim_type *RelationDim(CONST struct relation *rel){
1813 assert(rel!=NULL);
1814 return rel->d;
1815 }
1816
1817 int SetRelationDim(struct relation *rel, CONST dim_type *d)
1818 {
1819 if (!rel) return 1;
1820 rel->d = (dim_type*)d;
1821 return 0;
1822 }
1823
1824 double RelationResidual(CONST struct relation *rel){
1825 assert(rel!=NULL);
1826 return rel->residual;
1827 }
1828
1829 void SetRelationResidual(struct relation *rel, double value){
1830 assert(rel!=NULL);
1831 rel->residual = value;
1832 }
1833
1834 double RelationMultiplier(CONST struct relation *rel){
1835 assert(rel!=NULL);
1836 return rel->multiplier;
1837 }
1838
1839 void SetRelationMultiplier(struct relation *rel, double value){
1840 assert(rel!=NULL);
1841 rel->multiplier = value;
1842 }
1843
1844 double RelationNominal(CONST struct relation *rel){
1845 assert(rel!=NULL);
1846 return rel->nominal;
1847 }
1848
1849 void SetRelationNominal(struct relation *rel, double value){
1850 assert(rel!=NULL);
1851 rel->nominal = (fabs(value) > 0.0) ? fabs(value) : rel->nominal;
1852 }
1853
1854
1855 int RelationIsCond(CONST struct relation *rel){
1856 if ( rel != NULL) {
1857 return rel->iscond;
1858 }
1859 return 0;
1860 }
1861
1862 void SetRelationIsCond(struct relation *rel){
1863 if ( rel != NULL) {
1864 rel->iscond = 1;
1865 }else{
1866 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation");
1867 }
1868 }
1869
1870 unsigned long NumberVariables(CONST struct relation *rel)
1871 {
1872 unsigned long n;
1873 assert(rel!=NULL);
1874 n = (rel->vars!=NULL) ? gl_length(rel->vars) : 0;
1875 return n;
1876 }
1877
1878 struct Instance *RelationVariable(CONST struct relation *rel,
1879 unsigned long int varnum
1880 ){
1881 assert(rel!=NULL);
1882 return (struct Instance *)gl_fetch(rel->vars,varnum);
1883 }
1884
1885 static void CalcDepth(CONST struct relation *rel,
1886 int lhs,
1887 unsigned long int *depth,
1888 unsigned long int *maxdepth
1889 ){
1890 unsigned long c,length;
1891 CONST struct relation_term *term;
1892 length = RelationLength(rel,lhs);
1893 for(c=0;c<length;c++){
1894 term = NewRelationTerm(rel,c,lhs);
1895 switch(RelationTermType(term)){
1896 case e_zero:
1897 case e_int:
1898 case e_real:
1899 case e_var:
1900 if (++(*depth) > *maxdepth) *maxdepth = *depth;
1901 break;
1902 case e_func:
1903 case e_uminus:
1904 break;
1905 case e_plus:
1906 case e_minus:
1907 case e_times:
1908 case e_divide:
1909 case e_power:
1910 case e_ipower:
1911 (*depth)--;
1912 break;
1913 default:
1914 ASC_PANIC("Unknown relation term type");
1915 }
1916 }
1917 }
1918
1919 unsigned long RelationDepth(CONST struct relation *rel){
1920 unsigned long depth=0,maxdepth=0;
1921 switch(RelationRelop(rel)){
1922 case e_equal:
1923 case e_notequal:
1924 case e_less:
1925 case e_greater:
1926 case e_lesseq:
1927 case e_greatereq:
1928 CalcDepth(rel,1,&depth,&maxdepth);
1929 CalcDepth(rel,0,&depth,&maxdepth);
1930 assert(depth == 2);
1931 break;
1932 case e_maximize:
1933 case e_minimize:
1934 CalcDepth(rel,1,&depth,&maxdepth);
1935 assert(depth == 1);
1936 break;
1937 default:
1938 ASC_PANIC("Unknown relation term type");
1939 }
1940 return maxdepth;
1941 }
1942
1943 /*------------------------------------------------------------------------------
1944 EQUATION SCALING
1945
1946 "Documentation will be added at a later date" -- someone last century
1947 */
1948
1949 static double FindMaxAdditiveTerm(struct relation_term *s){
1950 enum safe_err serr;
1951 double lhs, rhs;
1952
1953 switch (RelationTermType(s)) {
1954 case e_plus:
1955 case e_minus:
1956 /** note these used to be inlined with max, but a bug in gcc323 caused it to be split out. */
1957 lhs = FindMaxAdditiveTerm(TermBinLeft(s));
1958 rhs = FindMaxAdditiveTerm(TermBinRight(s));
1959 return MAX(fabs(lhs), fabs(rhs));
1960 case e_uminus:
1961 return (FindMaxAdditiveTerm(TermUniLeft(s)));
1962 case e_times:
1963 return (FindMaxAdditiveTerm(TermBinLeft(s))*
1964 FindMaxAdditiveTerm(TermBinRight(s)));
1965 case e_divide:
1966 /* bug patch / 0 */
1967 return safe_div_D0(FindMaxAdditiveTerm(TermBinLeft(s)) ,
1968 RelationBranchEvaluator(TermBinRight(s)),&serr);
1969 default:
1970 return RelationBranchEvaluator(s);
1971 }
1972 }
1973
1974 static double FindMaxFromTop(struct relation *s){
1975 double lhs;
1976 double rhs;
1977 if (s == NULL) {
1978 return 0;
1979 }
1980 /** note these used to be inlined with max, but a bug in gcc323 caused it to be split out. */
1981 lhs = FindMaxAdditiveTerm(Infix_LhsSide(s));
1982 rhs = FindMaxAdditiveTerm(Infix_RhsSide(s));
1983 return MAX(fabs(lhs), fabs(rhs));
1984 }
1985
1986 /* send in relation */
1987 double CalcRelationNominal(struct Instance *i){
1988 enum Expr_enum reltype;
1989 struct Instance *c , *p;
1990 symchar *nomname;
1991
1992 char *iname;
1993 iname = WriteInstanceNameString(i,NULL);
1994 ascfree(iname);
1995
1996 glob_rel = NULL;
1997 if (i == NULL){
1998 FPRINTF(ASCERR, "error in CalcRelationNominal routine\n");
1999 return (double)0;
2000 }
2001 if (InstanceKind(i) != REL_INST) {
2002 FPRINTF(ASCERR, "error in CalcRelationNominal routine\n");
2003 return (double)0;
2004 }
2005 glob_rel = (struct relation *)GetInstanceRelation(i,&reltype);
2006 if (glob_rel == NULL) {
2007 FPRINTF(ASCERR, "error in CalcRelationNominal routine\n");
2008 return (double)0;
2009 }
2010
2011 if (reltype == e_token) {
2012 double temp;
2013 temp = FindMaxFromTop(glob_rel);
2014 if (asc_isnan(temp) || !asc_finite(temp)) {
2015 glob_rel = NULL;
2016 return (double)1;
2017 }
2018 if ( temp > 0) { /* this could return some really small numbers */
2019 glob_rel = NULL;
2020 return temp;
2021 }
2022 }
2023 if (reltype == e_blackbox){
2024 p = BlackBoxGetOutputVar(glob_rel);
2025 nomname = AddSymbol("nominal");
2026 glob_rel = NULL;
2027 c = ChildByChar(p,nomname);
2028 if (c == NULL) {
2029 ERROR_REPORTER_HERE(ASC_PROG_ERR,"nominal missing from standard var definition (assuming 1.0) (%s)",__FUNCTION__);
2030 return 1.0;
2031 } else {
2032 return( RealAtomValue(c) );
2033 }
2034 }
2035 if (reltype == e_glassbox){
2036 p = BlackBoxGetOutputVar(glob_rel);
2037 nomname = AddSymbol("nominal");
2038 glob_rel = NULL;
2039 c = ChildByChar(p,nomname);
2040 if (c == NULL) {
2041 ERROR_REPORTER_HERE(ASC_PROG_ERR,"nominal missing from standard var definition (assuming 1.0) (%s)",__FUNCTION__);
2042 return 1.0;
2043 } else {
2044 return( RealAtomValue(c) );
2045 }
2046 }
2047 if (reltype == e_opcode){
2048 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
2049 }
2050 glob_rel = NULL;
2051 return (double)1;
2052 }
2053
2054 void PrintScale(struct Instance *i){
2055 if (InstanceKind(i) == REL_INST) {
2056 double j;
2057 j = CalcRelationNominal(i);
2058 PRINTF(" scale constant = %g\n", j);
2059 }
2060 }
2061
2062 void PrintRelationNominals(struct Instance *i){
2063 VisitInstanceTree(i,PrintScale, 0, 0);
2064 }
2065
2066 /*------------------------------------------------------------------------------
2067 CALCULATION ROUTINES
2068 */
2069
2070 /**
2071 * Load ATOM values into an array of doubles.
2072 * The array of doubles is indexed from 0 while the
2073 * var list is indexed from 1. The ultimate client of
2074 * the array calling this function thinks vars index from 0.
2075 */
2076 static
2077 void RelationLoadDoubles(struct gl_list_t *varlist, double *vars){
2078 unsigned long c;
2079 vars--; /* back up the pointer so indexing from 1 puts data right */
2080 for (c= gl_length(varlist); c > 0; c--) {
2081 vars[c] = RealAtomValue((struct Instance *)gl_fetch(varlist,c));
2082 }
2083 }
2084
2085 /**
2086 only called on token relations.
2087 */
2088 int RelationCalcResidualBinary(CONST struct relation *r, double *res){
2089 double *vars;
2090 double tres;
2091 int old_errno;
2092
2093 if (r == NULL || res == NULL) {
2094 return 1;
2095 }
2096 vars = tmpalloc_array(gl_length(r->vars),double);
2097 if (vars == NULL) {
2098 return 1;
2099 }
2100 RelationLoadDoubles(r->vars,vars);
2101 old_errno = errno; /* push C global errno */
2102 errno = 0;
2103 if (BinTokenCalcResidual(RTOKEN(r).btable,RTOKEN(r).bindex,vars,&tres)) {
2104 if (errno == 0) { /* pop if unchanged */
2105 errno = old_errno;
2106 }
2107 return 1;
2108 }
2109 if (!asc_finite(tres) || errno == EDOM || errno == ERANGE ) {
2110 if (errno == 0) { /* pop if unchanged */
2111 errno = old_errno;
2112 }
2113 return 1;
2114 }
2115 if (errno == 0) { /* pop if unchanged */
2116 errno = old_errno;
2117 }
2118 *res = tres;
2119 return 0;
2120 }
2121
2122 enum safe_err
2123 RelationCalcResidualPostfixSafe(struct Instance *i, double *res){
2124 struct relation *r;
2125 enum Expr_enum reltype;
2126 enum safe_err status = safe_ok;
2127 unsigned long length_lhs, length_rhs;
2128
2129 CHECK_INST_RES(i,res,1);
2130
2131 r = (struct relation *)GetInstanceRelation(i, &reltype);
2132
2133 if( r == NULL ) {
2134 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation");
2135 return safe_problem;
2136 }
2137
2138 switch(reltype){
2139 case e_token:
2140 length_lhs = RelationLength(r, 1);
2141 length_rhs = RelationLength(r, 0);
2142
2143 if(length_lhs){
2144 length_lhs--;
2145 *res = RelationEvaluatePostfixBranchSafe(r, &length_lhs, 1,&status);
2146 }else{
2147 *res = 0.0;
2148 }
2149
2150 if(length_rhs){
2151 length_rhs--;
2152 *res -= RelationEvaluatePostfixBranchSafe(r, &length_rhs, 0,&status);
2153 }
2154
2155 safe_error_to_stderr(&status);
2156 break;
2157 case e_blackbox:
2158 if ( RelationCalcResidualPostfix(i,res) != 0) {
2159 CONSOLE_DEBUG("Problem evalutating Blackbox residual");
2160 status = safe_problem;
2161 safe_error_to_stderr(&status);
2162 }
2163 break;
2164 case e_glassbox:
2165 ERROR_REPORTER_HERE(ASC_PROG_ERR,"'e_glassbox' relation not supported");
2166 status = safe_problem;
2167 break;
2168 case e_opcode:
2169 ERROR_REPORTER_HERE(ASC_PROG_ERR,"'e_opcode' relation not supported");
2170 status = safe_problem;
2171 break;
2172 default:
2173 if(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2174 status = safe_problem;
2175 }else{
2176 ASC_PANIC("Reached end of routine!");
2177 }
2178 }
2179 return status;
2180 }
2181
2182 /* return 0 on success */
2183 int
2184 RelationCalcResidualPostfix(struct Instance *i, double *res){
2185 struct relation *r;
2186 enum Expr_enum reltype;
2187 unsigned long length_lhs, length_rhs;
2188
2189 CHECK_INST_RES(i,res,1);
2190
2191 r = (struct relation *)GetInstanceRelation(i, &reltype);
2192 if(r == NULL){
2193 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2194 return 1;
2195 }
2196
2197 /*
2198 struct Instance *p;
2199 p = InstanceParent(i,1);
2200 char *tmp;
2201 tmp = WriteRelationString(i,p,NULL,NULL,relio_ascend,NULL);
2202 CONSOLE_DEBUG("Evaluating residual for '%s'",tmp);
2203 ASC_FREE(tmp);
2204 */
2205
2206 switch(reltype){
2207 case e_token:
2208 length_lhs = RelationLength(r, 1);
2209 length_rhs = RelationLength(r, 0);
2210 if(length_lhs > 0){
2211 length_lhs--;
2212 *res = RelationEvaluatePostfixBranch(r, &length_lhs, 1);
2213 }else{
2214 *res = 0.0;
2215 }
2216 if(length_rhs > 0){
2217 length_rhs--;
2218 *res -= RelationEvaluatePostfixBranch(r, &length_rhs, 0);
2219 }
2220 return 0;
2221
2222 case e_blackbox:
2223 return BlackBoxCalcResidual(i, res, r);
2224
2225 case e_glassbox:
2226 case e_opcode:
2227 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode/glassbox not supported");
2228 return 1;
2229
2230 default: break;
2231 }
2232
2233 ERROR_REPORTER_HERE(ASC_PROG_ERR,"invalid relation type");
2234 return 1;
2235 }
2236
2237 int RelationCalcExceptionsInfix(struct Instance *i){
2238 enum Expr_enum reltype;
2239 double res;
2240 int result = 0;
2241 int old_errno;
2242
2243 glob_rel = NULL;
2244
2245 CHECK_INST_RES(i,&res,-1);
2246
2247 glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);
2248 if( glob_rel == NULL ) {
2249 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation");
2250 return -1;
2251 }
2252 if( reltype == e_token ) {
2253 if (Infix_LhsSide(glob_rel) != NULL) {
2254 old_errno = errno;
2255 errno = 0; /* save the last errno, because we don't know why */
2256 res = RelationBranchEvaluator(Infix_LhsSide(glob_rel));
2257 if (!asc_finite(res) || errno == EDOM || errno == ERANGE) {
2258 result |= RCE_ERR_LHS;
2259 if (asc_isnan(res)) {
2260 result |= RCE_ERR_LHSNAN;
2261 }else{
2262 if (!asc_finite(res)) {
2263 result |= RCE_ERR_LHSINF;
2264 }
2265 }
2266 }
2267 if (errno == 0) {
2268 errno = old_errno;
2269 } /* else something odd happened in evaluation */
2270 }
2271 if(Infix_RhsSide(glob_rel) != NULL) {
2272 res = RelationBranchEvaluator(Infix_RhsSide(glob_rel));
2273 if (!asc_finite(res)) {
2274 result |= RCE_ERR_RHS;
2275 if (asc_isnan(res)) {
2276 result |= RCE_ERR_RHSNAN;
2277 }else{
2278 if (!asc_finite(res)) {
2279 result |= RCE_ERR_LHSINF;
2280 }
2281 }
2282 }
2283 }
2284 glob_rel = NULL;
2285 return result;
2286 }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2287 ERROR_REPORTER_HERE(ASC_PROG_ERR,"relation type not implemented (%s)",__FUNCTION__);
2288 glob_rel = NULL;
2289 return -1;
2290 }
2291
2292 ASC_PANIC("reached end of routine");
2293 }
2294
2295
2296 int RelationCalcResidualInfix(struct Instance *i, double *res){
2297 enum Expr_enum reltype;
2298 glob_rel = NULL;
2299
2300 CHECK_INST_RES(i,res,1);
2301
2302 glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);
2303 if( glob_rel == NULL ) {
2304 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation\n");
2305 return 1;
2306 }
2307 if( reltype == e_token ) {
2308 if(Infix_LhsSide(glob_rel) != NULL) {
2309 *res = RelationBranchEvaluator(Infix_LhsSide(glob_rel));
2310 }else{
2311 *res = 0.0;
2312 }
2313 if(Infix_RhsSide(glob_rel) != NULL) {
2314 *res -= RelationBranchEvaluator(Infix_RhsSide(glob_rel));
2315 }
2316 glob_rel = NULL;
2317 return 0;
2318 }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2319 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
2320 glob_rel = NULL;
2321 return 1;
2322 }
2323
2324 ASC_PANIC("reached end of routine");
2325 }
2326
2327
2328 /*
2329 There used to be a stoopid comment here so I removed it.
2330 */
2331 int
2332 RelationCalcResidualPostfix2(struct Instance *i, double *res){
2333 struct relation *r;
2334 enum Expr_enum reltype;
2335
2336 CHECK_INST_RES(i,res,1);
2337
2338 r = (struct relation *)GetInstanceRelation(i, &reltype);
2339 if( r == NULL ) {
2340 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2341 return 1;
2342 }
2343
2344 if( reltype == e_token ){
2345 *res = RelationEvaluateResidualPostfix(r);
2346 return 0;
2347 }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2348 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
2349 return 1;
2350 }
2351
2352 ASC_PANIC("reached end of routine");
2353 }
2354
2355
2356 /*
2357 simply call the version that calculates the gradient and the residual,
2358 then ignore the residual
2359 */
2360 int
2361 RelationCalcGradient(struct Instance *r, double *grad){
2362 double residual;
2363 return RelationCalcResidGrad(r, &residual, grad);
2364 }
2365
2366 /*
2367 simply call the version that calculates the gradient and the residual,
2368 then ignore the residual
2369
2370 return 0 on success (as 'safe_ok' enum)
2371 */
2372 enum safe_err
2373 RelationCalcGradientSafe(struct Instance *r, double *grad){
2374 double residual;
2375 //CONSOLE_DEBUG("Gradient Evaluation Type: SAFE");
2376 return RelationCalcResidGradSafe(r, &residual, grad);
2377 }
2378
2379 /* return 0 on success, 1 on error */
2380 int RelationCalcResidGrad(struct Instance *i, double *residual, double *gradient){
2381 struct relation *r;
2382 enum Expr_enum reltype;
2383
2384 CHECK_INST_RES(i,residual,1);
2385 CHECK_INST_RES(i,residual,1);
2386
2387 r = (struct relation *)GetInstanceRelation(i, &reltype);
2388 if( r == NULL ) {
2389 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation");
2390 return 1;
2391 }
2392
2393 if(reltype == e_token ){
2394 return RelationEvaluateResidualGradient(r, residual, gradient);
2395 }
2396
2397 if(reltype == e_blackbox){
2398 return BlackBoxCalcResidGrad(i, residual, gradient, r);
2399 }
2400
2401 assert(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH);
2402 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype %d not implemented",reltype);
2403 return 1;
2404 }
2405
2406 enum safe_err RelationCalcResidGradSafe(struct Instance *i
2407 , double *residual, double *gradient
2408 ){
2409 struct relation *r;
2410 enum Expr_enum reltype;
2411 enum safe_err not_safe = safe_ok;
2412 int dummy_int;
2413
2414 #ifndef NDEBUG
2415 if( i == NULL ) {
2416 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null instance\n");
2417 not_safe = safe_problem;
2418 return not_safe;
2419 }
2420 if( residual == NULL || gradient == NULL ) {
2421 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null pointer\n");
2422 not_safe = safe_problem;
2423 return not_safe;
2424 }
2425 if( InstanceKind(i) != REL_INST ) {
2426 ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation\n");
2427 not_safe = safe_problem;
2428 return not_safe;
2429 }
2430 #endif
2431 r = (struct relation *)GetInstanceRelation(i, &reltype);
2432 if( r == NULL ) {
2433 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2434 not_safe = safe_problem;
2435 return not_safe;
2436 }
2437
2438
2439 if( reltype == e_token ) {
2440 dummy_int =
2441 RelationEvaluateResidualGradientSafe(r, residual, gradient, &not_safe);
2442 //CONSOLE_DEBUG("Relation Type: e_token");
2443 return not_safe;
2444 }
2445 if (reltype == e_blackbox){
2446 if (BlackBoxCalcResidGrad(i, residual, gradient, r) ) {
2447 not_safe = safe_problem;
2448 }
2449 CONSOLE_DEBUG("Relation Type: e_blackbox");
2450 return not_safe;
2451 }
2452 if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2453 if (reltype == e_glassbox){
2454 CONSOLE_DEBUG("Relation Type: e_glassbox");
2455 ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (%s)",__FUNCTION__);
2456 }
2457 if (reltype == e_opcode){
2458 CONSOLE_DEBUG("Relation Type: e_opcode");
2459 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
2460 }
2461 CONSOLE_DEBUG("Relation Type: other");
2462 not_safe = safe_problem;
2463 return not_safe;
2464 }
2465
2466 ASC_PANIC( "reached end of routine");
2467 }
2468
2469 /**-------------Reverse Automatic Differentiation---------------*/
2470
2471 /*
2472 simply call the version that calculates the gradient and the residual,
2473 then ignore the residual
2474 */
2475 int RelationCalcGradientRev(struct Instance *r, double *grad){
2476 double residual;
2477 return RelationCalcResidGradRev(r, &residual, grad);
2478 }
2479
2480 /*
2481 simply call the version that calculates the gradient and the residual,
2482 then ignore the residual
2483
2484 return 0 on success (as 'safe_ok' enum)
2485 */
2486 enum safe_err RelationCalcGradientRevSafe(struct Instance *r, double *grad){
2487 double residual;
2488 //CONSOLE_DEBUG("Gradient Evaluation Type: SAFE");
2489 return RelationCalcResidGradRevSafe(r, &residual, grad);
2490 }
2491
2492
2493 /* return 0 on success, 1 on error */
2494 int RelationCalcResidGradRev(struct Instance *i, double *residual, double *gradient){
2495 struct relation *r;
2496 enum Expr_enum reltype;
2497
2498 CHECK_INST_RES(i,residual,1);
2499 CHECK_INST_RES(i,residual,1);
2500
2501 r = (struct relation *)GetInstanceRelation(i, &reltype);
2502 if( r == NULL ) {
2503 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation");
2504 return 1;
2505 }
2506
2507 if(reltype == e_token ){
2508 RelationEvaluateResidualGradientRev(r, residual, gradient,0);
2509 return 0;
2510 }
2511
2512 if(reltype == e_blackbox){
2513 //return BlackBoxCalcResidGrad(i, residual, gradient, r);
2514 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Black Box Relation not implemented");
2515 }
2516
2517 assert(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH);
2518 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype %d not implemented",reltype);
2519 return 1;
2520 }
2521
2522 enum safe_err RelationCalcResidGradRevSafe(struct Instance *i
2523 , double *residual, double *gradient)
2524 {
2525 struct relation *r;
2526 enum Expr_enum reltype;
2527 enum safe_err not_safe = safe_ok;
2528
2529
2530 #ifndef NDEBUG
2531 if( i == NULL ) {
2532 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null instance\n");
2533 not_safe = safe_problem;
2534 return not_safe;
2535 }
2536 if( residual == NULL || gradient == NULL ) {
2537 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null pointer\n");
2538 not_safe = safe_problem;
2539 return not_safe;
2540 }
2541 if( InstanceKind(i) != REL_INST ) {
2542 ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation\n");
2543 not_safe = safe_problem;
2544 return not_safe;
2545 }
2546 #endif
2547 r = (struct relation *)GetInstanceRelation(i, &reltype);
2548 if( r == NULL ) {
2549 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2550 not_safe = safe_problem;
2551 return not_safe;
2552 }
2553
2554
2555 if( reltype == e_token ) {
2556 RelationEvaluateResidualGradientRevSafe(r, residual, gradient,0, &not_safe);
2557 //CONSOLE_DEBUG("Relation Type: e_token");
2558 return not_safe;
2559 }
2560 if (reltype == e_blackbox){
2561 /*if (BlackBoxCalcResidGrad(i, residual, gradient, r) ) {
2562 not_safe = safe_problem;
2563 }
2564 CONSOLE_DEBUG("Relation Type: e_blackbox");
2565 return not_safe;*/
2566 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Black Box Relation not implemented");
2567 }
2568 if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2569 if (reltype == e_glassbox){
2570 CONSOLE_DEBUG("Relation Type: e_glassbox");
2571 ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (%s)",__FUNCTION__);
2572 }
2573 if (reltype == e_opcode){
2574 CONSOLE_DEBUG("Relation Type: e_opcode");
2575 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
2576 }
2577 CONSOLE_DEBUG("Relation Type: other");
2578 not_safe = safe_problem;
2579 return not_safe;
2580 }
2581
2582 ASC_PANIC( "reached end of routine");
2583 }
2584
2585 /**----------------------------Second Derivative Calculations ------------------------------- */
2586
2587 /* return 0 on success, 1 on error */
2588 int RelationCalcSecondDeriv(struct Instance *i, double *deriv2nd, unsigned long var_index){
2589 struct relation *r;
2590 enum Expr_enum reltype;
2591
2592
2593 r = (struct relation *)GetInstanceRelation(i, &reltype);
2594 if( r == NULL ) {
2595 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation");
2596 return 1;
2597 }
2598
2599 if(reltype == e_token ){
2600 RelationEvaluateSecondDeriv(r,deriv2nd,var_index,0,NULL); //FIXME Check this
2601 return 0;
2602 }
2603
2604 if(reltype == e_blackbox){
2605 //return BlackBoxCalcResidGrad(i, residual, gradient, r);
2606 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Black Box Relation not implemented");
2607 }
2608
2609 assert(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH);
2610 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype %d not implemented",reltype);
2611 return 1;
2612 }
2613
2614 enum safe_err RelationCalcSecondDerivSafe(struct Instance *i, double *deriv2nd,unsigned long var_index)
2615 {
2616 struct relation *r;
2617 enum Expr_enum reltype;
2618 enum safe_err not_safe = safe_ok;
2619
2620
2621 #ifndef NDEBUG
2622 if( i == NULL ) {
2623 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null instance\n");
2624 not_safe = safe_problem;
2625 return not_safe;
2626 }
2627 if(deriv2nd == NULL ) {
2628 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null pointer\n");
2629 not_safe = safe_problem;
2630 return not_safe;
2631 }
2632 if( InstanceKind(i) != REL_INST ) {
2633 ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation\n");
2634 not_safe = safe_problem;
2635 return not_safe;
2636 }
2637 #endif
2638 r = (struct relation *)GetInstanceRelation(i, &reltype);
2639 if( r == NULL ) {
2640 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2641 not_safe = safe_problem;
2642 return not_safe;
2643 }
2644
2645
2646 if( reltype == e_token ) {
2647 RelationEvaluateSecondDerivSafe(r,
2648 deriv2nd,
2649 var_index,
2650 0,
2651 NULL,
2652 &not_safe);
2653 return not_safe;
2654 }
2655 if (reltype == e_blackbox){
2656 /*if (BlackBoxCalcResidGrad(i, residual, gradient, r) ) {
2657 not_safe = safe_problem;
2658 }
2659 CONSOLE_DEBUG("Relation Type: e_blackbox");
2660 return not_safe;*/
2661 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Black Box Relation not implemented");
2662 }
2663 if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2664 if (reltype == e_glassbox){
2665 CONSOLE_DEBUG("Relation Type: e_glassbox");
2666 ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (%s)",__FUNCTION__);
2667 }
2668 if (reltype == e_opcode){
2669 CONSOLE_DEBUG("Relation Type: e_opcode");
2670 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
2671 }
2672 CONSOLE_DEBUG("Relation Type: other");
2673 not_safe = safe_problem;
2674 return not_safe;
2675 }
2676
2677 ASC_PANIC( "reached end of routine");
2678 }
2679 /**------------------------------------------------------------------------------------------ */
2680
2681 /**----------------------------Hessian Calculations Routines------------------------------- */
2682
2683 /* return 0 on success, 1 on error */
2684 int RelationCalcHessianMtx(struct Instance *i, hessian_mtx *hess_mtx, unsigned long dimension){
2685 struct relation *r;
2686 enum Expr_enum reltype;
2687
2688 // CONSOLE_DEBUG("IN FUNCTION RelationCalcHessianMtx");
2689
2690 r = (struct relation *)GetInstanceRelation(i, &reltype);
2691 if( r == NULL ) {
2692 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation");
2693 return 1;
2694 }
2695
2696 if(reltype == e_token ){
2697 RelationEvaluateHessianMtx(r,hess_mtx,dimension);
2698 return 0;
2699 }
2700
2701 if(reltype == e_blackbox){
2702 //return BlackBoxCalcResidGrad(i, residual, gradient, r);
2703 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Black Box Relation not implemented");
2704 }
2705
2706 assert(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH);
2707 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype %d not implemented",reltype);
2708 return 1;
2709 }
2710
2711 enum safe_err RelationCalcHessianMtxSafe(struct Instance *i, hessian_mtx *hess_mtx,unsigned long dimension)
2712 {
2713 struct relation *r;
2714 enum Expr_enum reltype;
2715 enum safe_err not_safe = safe_ok;
2716
2717 // CONSOLE_DEBUG("IN FUNCTION RelationCalcHessianMtxSafe");
2718
2719 #ifndef NDEBUG
2720 if( i == NULL ) {
2721 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null instance\n");
2722 not_safe = safe_problem;
2723 return not_safe;
2724 }
2725 if( hess_mtx == NULL ) {
2726 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null pointer\n");
2727 not_safe = safe_problem;
2728 return not_safe;
2729 }
2730 if( InstanceKind(i) != REL_INST ) {
2731 ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation\n");
2732 not_safe = safe_problem;
2733 return not_safe;
2734 }
2735 #endif
2736 r = (struct relation *)GetInstanceRelation(i, &reltype);
2737 if( r == NULL ) {
2738 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2739 not_safe = safe_problem;
2740 return not_safe;
2741 }
2742
2743
2744 if( reltype == e_token ) {
2745 RelationEvaluateHessianMtxSafe(r,hess_mtx,dimension,&not_safe);
2746 return not_safe;
2747 }
2748 if (reltype == e_blackbox){
2749 /*if (BlackBoxCalcResidGrad(i, residual, gradient, r) ) {
2750 not_safe = safe_problem;
2751 }
2752 CONSOLE_DEBUG("Relation Type: e_blackbox");
2753 return not_safe;*/
2754 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Black Box Relation not implemented");
2755 }
2756 if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2757 if (reltype == e_glassbox){
2758 CONSOLE_DEBUG("Relation Type: e_glassbox");
2759 ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (%s)",__FUNCTION__);
2760 }
2761 if (reltype == e_opcode){
2762 CONSOLE_DEBUG("Relation Type: e_opcode");
2763 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
2764 }
2765 CONSOLE_DEBUG("Relation Type: other");
2766 not_safe = safe_problem;
2767 return not_safe;
2768 }
2769
2770 ASC_PANIC( "reached end of routine");
2771 }
2772
2773
2774 /**------------------------------------------------------------------------------------------ */
2775
2776
2777 /*
2778 calculate the derivative with respect to a single variable
2779 whose index is index, where 1<=index<=NumberVariables(r)
2780
2781 @TODO this appears only to be used in PrintGradients
2782 */
2783 int
2784 RelationCalcDerivative(struct Instance *i,
2785 unsigned long vindex,
2786 double *gradient)
2787 {
2788 struct relation *r;
2789 enum Expr_enum reltype;
2790
2791 CHECK_INST_RES(i,gradient,1);
2792
2793 r = (struct relation *)GetInstanceRelation(i, &reltype);
2794 if( r == NULL ) {
2795 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2796 return 1;
2797 }
2798 if( (vindex < 1) || (vindex > NumberVariables(r)) ) {
2799 ERROR_REPORTER_HERE(ASC_PROG_ERR,"index out of bounds\n");
2800 return 1;
2801 }
2802
2803 if( reltype == e_token ) {
2804 *gradient = RelationEvaluateDerivative(r, vindex);
2805 return 0;
2806 }
2807 else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2808 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not supported (%s)",__FUNCTION__);
2809 return 1;
2810 }
2811 ASC_PANIC( "reached end of routine");
2812 }
2813
2814 enum safe_err
2815 RelationCalcDerivativeSafe(struct Instance *i,
2816 unsigned long vindex,
2817 double *gradient)
2818 {
2819 struct relation *r;
2820 enum Expr_enum reltype;
2821 enum safe_err not_safe = safe_ok;
2822
2823 #ifndef NDEBUG
2824 if(!relutil_check_inst_and_res(i,gradient)){
2825 not_safe = safe_problem;
2826 return not_safe;
2827 }
2828 #endif
2829 r = (struct relation *)GetInstanceRelation(i, &reltype);
2830 if( r == NULL ) {
2831 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2832 not_safe = safe_problem;
2833 return not_safe;
2834 }
2835 if( (vindex < 1) || (vindex > NumberVariables(r)) ) {
2836 ERROR_REPORTER_HERE(ASC_PROG_ERR,"index out of bounds\n");
2837 not_safe = safe_problem;
2838 return not_safe;
2839 }
2840
2841 if( reltype == e_token ) {
2842 *gradient = RelationEvaluateDerivativeSafe(r, vindex, &not_safe);
2843 return not_safe;
2844 }
2845 else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2846 ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not supported (%s)",__FUNCTION__);
2847 not_safe = safe_problem;
2848 return not_safe;
2849 }
2850
2851 ASC_PANIC( "reached end of routine");
2852 }
2853
2854 /**
2855 Function for testing residual and gradient calulations
2856 */
2857 void PrintGradients(struct Instance *i){
2858 if (InstanceKind(i) == REL_INST) {
2859 double res, grads[1000];
2860 unsigned long vars, v;
2861 enum Expr_enum type;
2862 enum safe_err safe;
2863 struct relation* rel;
2864
2865 rel = (struct relation *)GetInstanceRelation(i,&type);
2866 vars = NumberVariables(rel);
2867
2868 /***** use the non safe versions *****/
2869 for( v = 0; v < vars; v++ ) {
2870 if( ! RelationCalcDerivative(i, v+1, &res) ) {
2871 PRINTF("derivative in%5ld =\t%g\n", v+1, res);
2872 }
2873 else {
2874 PRINTF("**** RelationCalcDerivative returned nonzero status\n");
2875 }
2876 }
2877
2878 if( ! RelationCalcResidGrad(i,&res,grads) ) {
2879 for (v = 0; v < vars; v++) {
2880 PRINTF("gradient in %6ld =\t%g\n", v+1, grads[v]);
2881 }
2882 PRINTF("residual from grad =\t%g\n", res);
2883 }
2884 else {
2885 PRINTF("**** RelationCalcResidGrad returned nonzero status\n");
2886 }
2887
2888 if( !RelationCalcResidualInfix(i,&res) ) {
2889 PRINTF(" infix residual =\t%g\n", res);
2890 }
2891 else {
2892 PRINTF("**** RelationCalcResidualInfix returned nonzero status\n");
2893 }
2894
2895 if( !RelationCalcResidualPostfix(i,&res) ) {
2896 PRINTF(" postfix residual =\t%g\n", res);
2897 }
2898 else {
2899 PRINTF("**** RelationCalcResidualPostfix returned nonzero status\n");
2900 }
2901
2902 if( !RelationCalcResidualPostfix2(i,&res) ) {
2903 PRINTF(" postfix2 residual =\t%g\n", res);
2904 }
2905 else {
2906 PRINTF("**** RelationCalcResidualPostfix2 returned nonzero status\n");
2907 }
2908
2909 /***** use the safe versions *****/
2910 for( v = 0; v < vars; v++ ) {
2911 if(safe_ok == (safe = RelationCalcDerivativeSafe(i, v+1, &res)) ) {
2912 PRINTF("safe deriv in%5ld =\t%g\n", v+1, res);
2913 }
2914 else {
2915 PRINTF("**** RelationCalcDerivativeSafe returned nonzero: %d\n", safe);
2916 }
2917 }
2918
2919 if(safe_ok == (safe = RelationCalcResidGradSafe(i,&res,grads)) ) {
2920 for (v = 0; v < vars; v++) {
2921 PRINTF("safe grad in%6ld =\t%g\n", v+1, grads[v]);
2922 }
2923 PRINTF("safe resid ala grad=\t%g\n", res);
2924 }
2925 else {
2926 PRINTF("**** RelationCalcResidGradSafe returned nonzero: %d\n", safe);
2927 }
2928
2929 if(safe_ok==(safe=RelationCalcResidGradRevSafe(i,&res,grads))){
2930 for (v = 0; v < vars; v++) {
2931 PRINTF("reverse safe grad in%6ld =\t%g\n", v+1, grads[v]);
2932 }
2933 PRINTF("reverse safe resid ala grad=\t%g\n", res);
2934 }
2935 else {
2936 PRINTF("**** RelationCalcResidGradRevSafe returned nonzero: %d\n", safe);
2937 }
2938
2939 if( ! RelationCalcResidGradRev(i,&res,grads) ) {
2940 for (v = 0; v < vars; v++) {
2941 PRINTF("reverse non-safe gradient in %6ld =\t%g\n", v+1, grads[v]);
2942 }
2943 PRINTF("residual from reverse non-safe grad =\t%g\n", res);
2944 }
2945 else {
2946 PRINTF("**** RelationCalcResidGradRev returned nonzero status\n");
2947 }
2948 /***** not implemented
2949 if( ! (safe = RelationCalcResidualInfixSafe(i,&res)) ) {
2950 PRINTF("safe infix residual=\t%g\n", res);
2951 }
2952 else {
2953 PRINTF("**** RelationCalcResidualInfixSafe returned nonzero: %d\n",
2954 safe);
2955 }
2956 *****/
2957
2958 if(safe_ok == (safe = RelationCalcResidualPostfixSafe(i,&res)) ) {
2959 PRINTF("safe postfix resid =\t%g\n", res);
2960 }
2961 else {
2962 PRINTF("**** RelationCalcResidualPostfixSafe returned nonzero: %d\n",
2963 safe);
2964 }
2965
2966 /***** not implemented
2967 if( ! (safe = RelationCalcResidualPostfix2Safe(i,&res)) ) {
2968 PRINTF("safe postfix2 resd =\t%g\n", res);
2969 }
2970 else {
2971 PRINTF("**** RelationCalcResidualPostfix2Safe returned nonzero: %d\n",
2972 safe);
2973 }
2974 *****/
2975
2976 PRINTF("\n");
2977 }
2978 }
2979 void PrintRelationGradients(struct Instance *i)
2980 {
2981 VisitInstanceTree(i,PrintGradients, 0, 0);
2982 }
2983
2984 /* this function may make an fpe for method 2 or 3.
2985 * list must be of nonnull struct relation * for
2986 * meth = m_BIN and struct Instance * for 1-3.
2987 */
2988 #define m_BIN 0
2989 #define m_PFS 1
2990 #define m_PF 2
2991 #define m_IF 3
2992 void TimeCalcResidual(struct gl_list_t *rlist,int method){
2993 unsigned long c,len;
2994 double res;
2995
2996 if (rlist==NULL) return;
2997 switch (method) {
2998 case m_BIN:
2999 for (c=1,len=gl_length(rlist); c <= len; c++) {
3000 RelationCalcResidualBinary(gl_fetch(rlist,c),&res);
3001 }
3002 break;
3003 case m_PFS:
3004 for (c=1,len=gl_length(rlist); c <= len; c++) {
3005 RelationCalcResidualPostfixSafe(gl_fetch(rlist,c),&res);
3006 }
3007 break;
3008 case m_PF:
3009 for (c=1,len=gl_length(rlist); c <= len; c++) {
3010 RelationCalcResidualPostfix(gl_fetch(rlist,c),&res);
3011 }
3012 break;
3013 case m_IF:
3014 for (c=1,len=gl_length(rlist); c <= len; c++) {
3015 RelationCalcResidualInfix(gl_fetch(rlist,c),&res);
3016 }
3017 break;
3018 default:
3019 break;
3020 }
3021 return;
3022 }
3023
3024 void PrintResidual(struct Instance *i){
3025 enum safe_err se;
3026 struct relation *rel;
3027 enum Expr_enum reltype;
3028 int errb;
3029 #ifndef M_PI
3030 #define M_PIE 3.141590271828
3031 #else
3032 #define M_PIE M_PI
3033 #endif
3034 double post=M_PIE,in=M_PIE,postsafe=M_PIE,binary=M_PIE;
3035
3036 if (InstanceKind(i) == REL_INST) {
3037 rel = (struct relation *)GetInstanceRelation(i,&reltype);
3038 if (reltype == e_token) {
3039 errb = RelationCalcResidualBinary(rel,&(binary));
3040 }else{
3041 errb = 1;
3042 }
3043 se = RelationCalcResidualPostfixSafe(i,&(postsafe));
3044 if (errb || se != safe_ok) {
3045 FPRINTF(ASCERR,"Skipping Postfix,Infix\n");
3046 }else{
3047 RelationCalcResidualPostfix(i,&(post));
3048 RelationCalcResidualInfix(i,&(in));
3049 }
3050 PRINTF("binary residual = %.18g\n",binary);
3051 PRINTF("postfix safe res = %.18g\n",postsafe);
3052 if (errb||se!= safe_ok) {
3053 PRINTF("postfix residual = %.18g\n",post);
3054 PRINTF(" infix residual = %.18g\n",in);
3055 }
3056 if(binary != postsafe) {
3057 PRINTF("!!!!!!!ERROR!!!!!!! %g \n", binary-post);
3058 }
3059 PRINTF("(Unchanged residuals = %.18g\n\n",M_PIE);
3060 }
3061 }
3062
3063 void PrintRelationResiduals(struct Instance *i){
3064 VisitInstanceTree(i,PrintResidual, 0, 0);
3065 }
3066
3067 /*==============================================================================
3068 'RELATIONFINDROOTS' AND SUPPORT FUNCTIONS
3069
3070 The following functions support RelationFindRoots which
3071 is the compiler-side implementation of our old solver-side DirectSolve
3072 function.
3073 */
3074
3075 double *RelationFindRoots(struct Instance *i,
3076 double lower_bound, double upper_bound,
3077 double nominal,
3078 double tolerance,
3079 unsigned long *varnum,
3080 int *able,
3081 int *nsolns
3082 ){
3083 struct relation *rel;
3084 double sideval;
3085 enum Expr_enum reltype;
3086 static struct ds_soln_list soln_list = {0,0,NULL};
3087 CONST struct gl_list_t *list;
3088
3089 /* check for recycle shutdown */
3090 if (i==NULL && varnum == NULL && able == NULL && nsolns == NULL) {
3091 if (soln_list.soln != NULL) {
3092 ascfree(soln_list.soln);
3093 soln_list.soln = NULL;
3094 soln_list.length = soln_list.capacity = 0;
3095 }
3096 RootFind(NULL,NULL,NULL,NULL,NULL,0L,NULL); /*clear brent recycle */
3097 RelationCreateTmp(0,0,e_nop); /* clear tmprelation recycle */
3098 return NULL;
3099 }
3100 /* check assertions */
3101 #ifndef NDEBUG
3102 if( i == NULL ) {
3103 FPRINTF(ASCERR, "error in RelationFindRoot: NULL instance\n");
3104 glob_rel = NULL;
3105 return NULL;
3106 }
3107 if (able == NULL){
3108 FPRINTF(ASCERR,"error in RelationFindRoot: NULL able ptr\n");
3109 glob_rel = NULL;
3110 return NULL;
3111 }
3112 if (varnum == NULL){
3113 FPRINTF(ASCERR,"error in RelationFindRoot: NULL varnum\n");
3114 glob_rel = NULL;
3115 return NULL;
3116 }
3117 if( InstanceKind(i) != REL_INST ) {
3118 FPRINTF(ASCERR, "error in RelationFindRoot: not relation\n");
3119 glob_rel = NULL;
3120 return NULL;
3121 }
3122 #endif
3123
3124 *able = FALSE;
3125 *nsolns = -1; /* nsolns will be -1 for a very unhappy root-finder */
3126 glob_rel = NULL;
3127 glob_done = 0;
3128 soln_list.length = 0; /* reset len to 0. if NULL to start, append mallocs */
3129 append_soln(&soln_list,0.0);
3130 rel = (struct relation *)GetInstanceRelation(i, &reltype);
3131 if( rel == NULL ) {
3132 FPRINTF(ASCERR, "error in RelationFindRoot: NULL relation\n");
3133 glob_rel = NULL; return NULL;
3134 }
3135 /* here we should switch and handle all types. at present we don't
3136 * handle anything except e_token
3137 */
3138 if( reltype != e_token ) {
3139 FPRINTF(ASCERR, "error in RelationFindRoot: non-token relation\n");
3140 glob_rel = NULL;
3141 return NULL;
3142 }
3143
3144 if (RelationRelop(rel) == e_equal){
3145 glob_rel = RelationTmpTokenCopy(rel);
3146 assert(glob_rel!=NULL);
3147 glob_done = 0;
3148 list = RelationVarList(glob_rel);
3149 if( *varnum >= 1 && *varnum <= gl_length(list)){
3150 glob_done = 1;
3151 }
3152 if (!glob_done) {
3153 FPRINTF(ASCERR, "error in FindRoot: var not found\n");
3154 glob_rel = NULL;
3155 return NULL;
3156 }
3157
3158 glob_varnum = *varnum;