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