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