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

Contents of /trunk/base/generic/compiler/relation_util.c

Parent Directory Parent Directory | Revision Log Revision Log


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