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