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