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