/[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 979 - (show annotations) (download) (as text)
Wed Dec 20 14:34:16 2006 UTC (13 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 120946 byte(s)
Added simplified ASC_PANIC call that uses var-args, added throughout relation_util.c.
Fixed var_filter_t stuff in djex and fvex.
More assertions in integrator.c
Added output of initial state from lsode.c (hoping that's a good idea?)
Fixed output code from relman_diff2.
Added asc_panic_nofunc for non var-arg CPPs.
Disabled -O3 flag in building C++ API
Added __getitem__ and __getattr__ methods in Simuluation for simplified python syntax (eg M.x instead M.sim.x)
Integrator::analyse throws exceptions on error now.

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