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