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