/[ascend]/trunk/base/generic/compiler/relation_util.c
ViewVC logotype

Annotation of /trunk/base/generic/compiler/relation_util.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 910 - (hide annotations) (download) (as text)
Thu Oct 26 13:35:25 2006 UTC (13 years, 9 months ago) by johnpye
File MIME type: text/x-csrc
File size: 121036 byte(s)
In instantiate.c, made new blackbox code tolerant of blackboxes that don't need initialisation.
Removed some debug output.
Expanded 'extfntest.py' a little bit, for ease of testing.
Converted 'blackbox is experimental' warnings to one-time-only.
Minor change to way that webbrowser is invoked under linux.
1 johnpye 709 /* ASCEND modelling environment
2     Copyright (C) 1997, 2006 Carnegie Mellon University
3     Copyright (C) 1993, 1994 Joseph James Zaher, Benjamin Andrew Allan
4     Copyright (C) 1993 Joseph James Zaher
5     Copyright (C) 1990 Thomas Guthrie Epperly, Karl Michael Westerberg
6    
7     This program is free software; you can redistribute it and/or modify
8     it under the terms of the GNU General Public License as published by
9     the Free Software Foundation; either version 2, or (at your option)
10     any later version.
11    
12     This program is distributed in the hope that it will be useful,
13     but WITHOUT ANY WARRANTY; without even the implied warranty of
14     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15     GNU General Public License for more details.
16    
17     You should have received a copy of the GNU General Public License
18     along with this program; if not, write to the Free Software
19     Foundation, Inc., 59 Temple Place - Suite 330,
20     Boston, MA 02111-1307, USA.
21     *//**
22     @file
23     Relation utility functions for Ascend
24    
25     This module defines the dimensionality checking and some other
26 johnpye 908 relation auxillaries for ASCEND.
27 johnpye 709 *//*
28     Last in CVS: $Revision: 1.44 $ $Date: 1998/04/23 23:51:09 $ $Author: ballan $
29     */
30    
31 aw0a 1 #include <math.h>
32     #include <errno.h>
33     #include <stdarg.h>
34 johnpye 399 #include <utilities/ascConfig.h>
35     #include <utilities/ascMalloc.h>
36     #include <utilities/ascPanic.h>
37 johnpye 908 #include <general/mathmacros.h>
38 johnpye 399 #include <general/list.h>
39     #include <general/dstring.h>
40     #include "compiler.h"
41     #include "symtab.h"
42     #include "functype.h"
43     #include "safe.h"
44     #include "fractions.h"
45     #include "dimen.h"
46 johnpye 669 #include "expr_types.h"
47 johnpye 399 #include "vlist.h"
48     #include "dimen_io.h"
49     #include "instance_enum.h"
50     #include "bintoken.h"
51     #include "find.h"
52     #include "atomvalue.h"
53     #include "instance_name.h"
54     #include "relation_type.h"
55 johnpye 908 #include "extfunc.h"
56     #include "rel_blackbox.h"
57 johnpye 399 #include "relation.h"
58     #include "relation_util.h"
59 johnpye 908 #include "rel_blackbox.h"
60 johnpye 399 #include "instance_io.h"
61     #include "instquery.h"
62     #include "visitinst.h"
63     #include "mathinst.h"
64     #include "rootfind.h"
65     #include "func.h"
66 aw0a 1
67 johnpye 709 #ifndef NDEBUG
68     #include "relation_io.h"
69     #endif
70    
71 johnpye 694 /*------------------------------------------------------------------------------
72     DATA TYPES AND GLOBAL VARS
73     */
74 aw0a 1
75     int g_check_dimensions_noisy = 1;
76     #define GCDN g_check_dimensions_noisy
77    
78 johnpye 694 /** global relation pointer to avoid passing a relation recursively */
79     static struct relation *glob_rel;
80 aw0a 1
81 johnpye 669 /**
82 johnpye 694 The following global variables are used thoughout the
83     functions called by RelationFindroot.
84 aw0a 1
85 johnpye 694 These should probably be located at the top of this
86     file alonge with glob_rel. [OK, let it be so then. -- JP]
87 johnpye 669 */
88 johnpye 846 static unsigned long glob_varnum;
89 johnpye 694 static int glob_done;
90 aw0a 1
91 johnpye 694 /* some data structurs...*/
92 aw0a 1
93     struct dimnode {
94     dim_type d;
95     enum Expr_enum type;
96     short int_const;
97     double real_const;
98     struct fraction power;
99     };
100    
101 johnpye 694 struct ds_soln_list {
102     int length,capacity;
103     double *soln;
104     };
105 aw0a 1
106 johnpye 669 /*
107 johnpye 694 Define the following if you want ASCEND to panic when it hits a
108     relation error in this file. This will help with debugging (GDB).
109    
110     Comment it out for ERROR_REPORTER to be used instead.
111     */
112     #define RELUTIL_CHECK_ABORT
113    
114     /*------------------------------------------------------------------------------
115     forward declarations
116     */
117    
118     static struct fraction real_to_frac(double real);
119     char *tmpalloc(int nbytes);
120     int ArgsForRealToken(enum Expr_enum type);
121     static int IsZero(struct dimnode *node);
122    
123     /* bunch of support functions for RelationFindRoots */
124     static double RootFind(struct relation *rel, double *lower_bound, double *upper_bound,
125     double *nominal,double *tolerance,unsigned long varnum,int *status);
126     static int CalcResidGivenValue(int *mode, int *m, unsigned long *varnum,double *val, double *u, double *f, double *g);
127     int RelationInvertTokenTop(struct ds_soln_list *soln_list);
128     int RelationInvertToken(struct relation_term **term,struct ds_soln_list *soln_list,enum safe_err *not_safe);
129     static void SetUpInvertTokenTop(struct relation_term **invert_side,double *value);
130     static int SetUpInvertToken(struct relation_term *term,struct relation_term **invert_side,double *value);
131     static int SearchEval_Branch(struct relation_term *term);
132     static void InsertBranchResult(struct relation_term *term, double value);
133     static void remove_soln( struct ds_soln_list *sl, int ndx);
134     static void append_soln( struct ds_soln_list *sl, double soln);
135     static struct relation *RelationTmpTokenCopy(CONST struct relation *src);
136     static int RelationTmpCopySide(union RelationTermUnion *old,unsigned long len,union RelationTermUnion *arr);
137     static struct relation *RelationCreateTmp(unsigned long lhslen, unsigned long rhslen, enum Expr_enum relop);
138    
139     /* the following appear only to be used locally, so I've made them static -- JP */
140    
141 johnpye 908 static int RelationCalcDerivative(struct Instance *i, unsigned long vindex, double *grad);
142 johnpye 694 /**<
143     * This calculates the derivative of the relation df/dx (f = lhs-rhs)
144 johnpye 908 * where x is the VINDEX-th entry in the relation's var list.
145 johnpye 694 * The var list is a gl_list_t indexed from 1 to length.
146     * Non-zero return value implies a problem.<br><br>
147     *
148     * Notes: This function is a possible source of floating point
149     * exceptions and should not be used during compilation.
150     */
151    
152     static enum safe_err
153 johnpye 908 RelationCalcDerivativeSafe(struct Instance *i, unsigned long vindex, double *grad);
154 johnpye 694 /**<
155     * Calculates the derivative safely.
156     * Non-zero return value implies a problem.
157     */
158    
159     #ifndef NDEBUG
160     static int relutil_check_inst_and_res(struct Instance *i, double *res);
161     #endif
162    
163     #ifndef NDEBUG
164     # define CHECK_INST_RES(i,res,retval) if(!relutil_check_inst_and_res(i,res)){return retval;}
165     #else
166     # define CHECK_INST_RES(i,res,retval) ((void)0)
167     #endif
168    
169     /*------------------------------------------------------------------------------
170     SOME STUFF TO DO WITH DIMENSIONS
171     */
172    
173     /*
174 johnpye 669 @TODO what this does needs to be documented here
175     */
176 aw0a 1 static void apply_term_dimensions(struct relation *rel,
177 johnpye 694 struct relation_term *rt,
178     struct dimnode *first,
179     struct dimnode *second,
180     int *con,
181     int *wild
182     ){
183 aw0a 1 enum Expr_enum type;
184    
185     switch(type=RelationTermType(rt)) {
186     case e_zero:
187     CopyDimensions(WildDimension(),&(first->d));
188     first->real_const = 0.0;
189     first->type = type;
190     break;
191    
192     case e_int:
193     CopyDimensions(Dimensionless(),&(first->d));
194     first->int_const = (short)TermInteger(rt);
195     first->type = type;
196     break;
197    
198     case e_real:
199     CopyDimensions(TermDimensions(rt),&(first->d));
200     first->real_const = TermReal(rt);
201     first->type = type;
202     break;
203    
204     case e_var: {
205     struct Instance *var = RelationVariable(rel,TermVarNumber(rt));
206     CopyDimensions(RealAtomDims(var),&(first->d));
207     first->type = type;
208     break;
209     }
210     case e_func: {
211     enum Func_enum id = FuncId(TermFunc(rt));
212     switch( id ) {
213     case F_ABS:
214     case F_HOLD:
215     /* no checking or scaling */
216     break;
217    
218     case F_SQR:
219     /* no checking, just simple scaling */
220     first->d = ScaleDimensions(&(first->d),CreateFraction(2,1));
221     break;
222    
223     case F_CUBE:
224     /* no checking, just simple scaling */
225     first->d = ScaleDimensions(&(first->d),CreateFraction(3,1));
226     break;
227    
228     case F_SQRT:
229     /* no checking, just simple scaling */
230     first->d = ScaleDimensions(&(first->d),CreateFraction(1,2));
231     break;
232    
233     case F_CBRT:
234     /* no checking, just simple scaling */
235     first->d = ScaleDimensions(&(first->d),CreateFraction(1,3));
236     break;
237    
238     case F_EXP:
239     case F_LN:
240     case F_LNM:
241 johnpye 123 case F_LOG10:
242 aw0a 1 #ifdef HAVE_ERF
243     case F_ERF:
244     #endif /* HAVE_ERF */
245     case F_SINH:
246     case F_COSH:
247     case F_TANH:
248     case F_ARCSINH:
249     case F_ARCCOSH:
250     case F_ARCTANH:
251     /**
252     *** first must now be dimensionless. It will
253     *** end up dimensionless as well.
254     **/
255     if( IsWild(&(first->d)) && !IsZero(first) ) {
256     if( !*wild ) *wild = TRUE;
257     if (GCDN) {
258     FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
259     FPRINTF(ASCERR," in function %s.\n",
260     FuncName(TermFunc(rt)));
261     }
262     } else if( !IsWild(&(first->d)) &&
263     CmpDimen(&(first->d),Dimensionless()) ) {
264     if( *con ) *con = FALSE;
265     if (GCDN) {
266     FPRINTF(ASCERR,"ERROR: Function %s called with\n",
267     FuncName(TermFunc(rt)));
268     FPRINTF(ASCERR," dimensions ");
269     WriteDimensions(ASCERR,&(first->d));
270     FPRINTF(ASCERR,".\n");
271     }
272     }
273     CopyDimensions(Dimensionless(),&(first->d));
274     break;
275    
276     case F_SIN:
277     case F_COS:
278     case F_TAN: {
279     /**
280     *** first must now be of dimension D_PLANE_ANGLE.
281     *** It will then be made dimensionless.
282     **/
283     if( IsWild(&(first->d)) && !IsZero(first) ) {
284     if( !*wild ) *wild = TRUE;
285     if (GCDN) {
286     FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
287     FPRINTF(ASCERR," in function %s.\n",
288     FuncName(TermFunc(rt)) );
289     }
290 johnpye 726 }else{
291 aw0a 1 if( !IsWild(&(first->d)) &&
292     CmpDimen(&(first->d),TrigDimension()) ) {
293     if( *con ) *con = FALSE;
294     if (GCDN) {
295     FPRINTF(ASCERR,"ERROR: Function %s called with\n",
296     FuncName(TermFunc(rt)));
297     FPRINTF(ASCERR," dimensions ");
298     WriteDimensions(ASCERR,&(first->d));
299     FPRINTF(ASCERR,".\n");
300     }
301     }
302     }
303     CopyDimensions(Dimensionless(),&(first->d));
304     break;
305     }
306    
307     case F_ARCSIN:
308     case F_ARCCOS:
309     case F_ARCTAN:
310     /**
311     *** first must now be dimensionless. It will
312     *** end up with dimension D_PLANE_ANGLE
313     **/
314     if( IsWild(&(first->d)) && !IsZero(first) ) {
315     if( !*wild ) *wild = TRUE;
316     if (GCDN) {
317     FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
318     FPRINTF(ASCERR," in function %s.\n",
319     FuncName(TermFunc(rt)));
320     }
321     } else if( !IsWild(&(first->d)) &&
322     CmpDimen(&(first->d),Dimensionless()) ) {
323     if( *con ) *con = FALSE;
324     if (GCDN) {
325     FPRINTF(ASCERR,"ERROR: Function %s called with\n",
326     FuncName(TermFunc(rt)));
327     FPRINTF(ASCERR," dimensions ");
328     WriteDimensions(ASCERR,&(first->d));
329     FPRINTF(ASCERR,".\n");
330     }
331     }
332     CopyDimensions(TrigDimension(),&(first->d));
333     break;
334     }
335     first->type = type;
336     break;
337     }
338    
339     case e_uminus:
340     first->type = type;
341     break;
342    
343     case e_times:
344     first->d = AddDimensions(&(first->d),&(second->d));
345     first->type = type;
346     break;
347    
348     case e_divide:
349     first->d = SubDimensions(&(first->d),&(second->d));
350     first->type = type;
351     break;
352    
353     case e_power: /* fix me and add ipower */
354     if( IsWild(&(second->d)) && !IsZero(second) ) {
355     if( !*wild ) *wild = TRUE;
356     if (GCDN) {
357     FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
358     FPRINTF(ASCERR," in exponent.\n");
359     }
360     } else if( !IsWild(&(second->d)) &&
361     CmpDimen(&(second->d),Dimensionless()) ) {
362     if( *con ) *con = FALSE;
363     if (GCDN) {
364     FPRINTF(ASCERR,"ERROR: Exponent has dimensions ");
365     WriteDimensions(ASCERR,&(second->d));
366     FPRINTF(ASCERR,".\n");
367     }
368     }
369     CopyDimensions(Dimensionless(),&(second->d));
370     switch( second->type ) {
371     case e_int:
372     if( !IsWild(&(first->d)) &&
373     CmpDimen(&(first->d),Dimensionless()) ) {
374     struct fraction power;
375     power = CreateFraction(second->int_const,1);
376     first->d = ScaleDimensions(&(first->d),power);
377     }
378     break;
379    
380     case e_real:
381     if( !IsWild(&(first->d)) &&
382     CmpDimen(&(first->d),Dimensionless()) ) {
383     struct fraction power;
384     power = real_to_frac(second->real_const);
385     first->d = ScaleDimensions(&(first->d),power);
386     }
387     break;
388    
389     /* what about e_zero? */
390     default:
391     if( IsWild(&(first->d)) && !IsZero(first) ) {
392     if( !*wild ) *wild = TRUE;
393     if (GCDN) {
394     FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
395     FPRINTF(ASCERR," raised to a non-constant power.\n");
396     }
397     } else if( !IsWild(&(first->d)) &&
398     CmpDimen(&(first->d),Dimensionless()) ) {
399     if( *con ) *con = FALSE;
400     if (GCDN) {
401     FPRINTF(ASCERR,"ERROR: Dimensions ");
402     WriteDimensions(ASCERR,&(first->d));
403     FPRINTF(ASCERR," are\n");
404     FPRINTF(ASCERR," raised to a non-constant power.\n");
405     }
406     }
407     CopyDimensions(Dimensionless(),&(first->d));
408     break;
409    
410     }
411     first->type = type;
412     break;
413    
414     case e_plus:
415     case e_minus:
416     if( IsWild(&(first->d)) && IsZero(first) ) {
417     /* first wild zero */
418     CopyDimensions(&(second->d),&(first->d));
419     first->type = second->type;
420     if( second->type==e_int )
421     first->int_const = second->int_const;
422     if( second->type==e_real )
423     first->real_const = second->real_const;
424     } else if( IsWild(&(first->d)) && !IsZero(first) ) {
425     /* first wild non-zero */
426     if( IsWild(&(second->d)) && !IsZero(second) ) {
427     /* second wild non-zero */
428     if( !*wild ) *wild = TRUE;
429     if (GCDN) {
430     FPRINTF(ASCERR,"ERROR: %s has wild dimensions on\n",
431     type==e_plus ? "Addition":"Subtraction");
432     FPRINTF(ASCERR," left and right hand sides.\n");
433     }
434     first->type = type;
435     } else if( !IsWild(&(second->d)) ) {
436     /* second not wild */
437     if( !*wild ) *wild = TRUE;
438     if (GCDN) {
439     FPRINTF(ASCERR,"ERROR: %s has wild dimensions on\n",
440     type==e_plus ? "Addition":"Subtraction");
441     FPRINTF(ASCERR," left hand side.\n");
442     }
443     CopyDimensions(&(second->d),&(first->d));
444     first->type = type;
445     }
446     } else if( !IsWild(&(first->d)) ) {
447     /* first not wild */
448     if( IsWild(&(second->d)) && !IsZero(second) ) {
449     /* second wild non-zero */
450     if( !*wild ) *wild = TRUE;
451     if (GCDN) {
452     FPRINTF(ASCERR,"ERROR: %s has wild dimensions on\n",
453     type==e_plus ? "Addition":"Subtraction");
454     FPRINTF(ASCERR," right hand side.\n");
455     }
456     first->type = type;
457     } else if ( !IsWild(&(second->d)) ) {
458     /* second not wild */
459     if( CmpDimen(&(first->d),&(second->d)) ) {
460     if( *con ) *con = FALSE;
461     if (GCDN) {
462     FPRINTF(ASCERR,"ERROR: %s has dimensions ",
463     type==e_plus ? "Addition":"Subtraction");
464     WriteDimensions(ASCERR,&(first->d));
465     FPRINTF(ASCERR," on left\n");
466     FPRINTF(ASCERR," and dimensions ");
467     WriteDimensions(ASCERR,&(second->d));
468     FPRINTF(ASCERR," on right.\n");
469     }
470     }
471     first->type = type;
472     }
473     }
474     break;
475    
476     default:
477     FPRINTF(ASCERR,"ERROR: Unknown relation term type.\n");
478     if( *con ) *con = FALSE;
479     first->type = type;
480     break;
481     }
482     }
483    
484 johnpye 694 /**
485     @TODO what this does needs to be documented here
486     */
487 aw0a 1 int RelationCheckDimensions(struct relation *rel, dim_type *dimens)
488     {
489     struct dimnode *stack, *sp;
490     int consistent = TRUE;
491     int wild = FALSE;
492     unsigned long c, len;
493    
494 johnpye 908
495     /* FIXME blackbox checkDimens: this is where we need to check on relation type
496     (which we can't without the Instance) and get the output var
497     of blackbox for dims. */
498 aw0a 1 if ( !IsWild(RelationDim(rel)) ) { /* don't do this twice */
499     CopyDimensions(RelationDim(rel),dimens);
500     return 2;
501     }
502 johnpye 709 sp = stack = ASC_NEW_ARRAY(struct dimnode,RelationDepth(rel));
503 aw0a 1 switch( RelationRelop(rel) ) {
504     case e_less:
505     case e_lesseq:
506     case e_greater:
507     case e_greatereq:
508     case e_equal:
509     case e_notequal:
510     /* Working on the left-hand-side */
511     len = RelationLength(rel,TRUE);
512     for( c = 1; c <= len; c++ ) {
513     struct relation_term *rt;
514     rt = (struct relation_term *)RelationTerm(rel,c,TRUE);
515     sp += 1-ArgsForRealToken(RelationTermType(rt));
516     apply_term_dimensions(rel,rt,sp-1,sp,&consistent,&wild);
517     } /* stack[0].d contains the dimensions of the lhs expression */
518    
519     /* Now working on the right-hand_side */
520     len = RelationLength(rel,FALSE);
521     for( c = 1; c <= len; c++ ) {
522     struct relation_term *rt;
523     rt = (struct relation_term *) RelationTerm(rel,c,FALSE);
524     sp += 1-ArgsForRealToken(RelationTermType(rt));
525     apply_term_dimensions(rel,rt,sp-1,sp,&consistent,&wild);
526     } /* stack[1].d contains the dimensions of the rhs expression */
527    
528     if( IsWild(&(stack[0].d)) || IsWild(&(stack[1].d)) ) {
529     if( IsWild(&(stack[0].d)) && !IsZero(&(stack[0])) ) {
530     if( !wild ) wild = TRUE;
531     if (GCDN) {
532     FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
533     FPRINTF(ASCERR," on left hand side.\n");
534     }
535     }
536     if( IsWild(&(stack[1].d)) && !IsZero(&(stack[1])) ) {
537     if( !wild ) wild = TRUE;
538     if (GCDN) {
539     FPRINTF(ASCERR,"ERROR: Relation has wild dimensions\n");
540     FPRINTF(ASCERR," on right hand side.\n");
541     }
542     }
543 johnpye 726 }else{
544 aw0a 1 if( CmpDimen(&(stack[0].d),&(stack[1].d)) ) {
545     if( consistent ) consistent = FALSE;
546     if (GCDN) {
547     FPRINTF(ASCERR,"ERROR: Relation has dimensions ");
548     WriteDimensions(ASCERR,&(stack[0].d));
549     FPRINTF(ASCERR," on left\n");
550     FPRINTF(ASCERR," and dimensions ");
551     WriteDimensions(ASCERR,&(stack[1].d));
552     FPRINTF(ASCERR," on right.\n");
553     }
554     }
555     }
556     break;
557     case e_maximize:
558     case e_minimize:
559     /* Working on the left-hand-side */
560     len = RelationLength(rel,TRUE);
561     for( c = 1; c <= len; c++ ) {
562     struct relation_term *rt;
563     rt = (struct relation_term *) RelationTerm(rel,c,TRUE);
564     sp += 1-ArgsForRealToken(RelationTermType(rt));
565     apply_term_dimensions(rel,rt,sp-1,sp,&consistent,&wild);
566     } /* stack[0].d contains the dimensions of the lhs expression */
567    
568     if( IsWild(&(stack[0].d)) && !IsZero(&(stack[0])) ) {
569     if( !wild ) wild = TRUE;
570     if (GCDN) {
571     FPRINTF(ASCERR,"ERROR: Objective has wild dimensions.\n");
572     }
573     }
574     break;
575    
576     default:
577     FPRINTF(ASCERR,"ERROR: Unknown relation type.\n");
578     if( consistent ) consistent = FALSE;
579     break;
580     }
581     CopyDimensions(&(stack[0].d),dimens);
582     ascfree(stack);
583     return( consistent && !wild );
584     }
585    
586 johnpye 669 /*------------------------------------------------------------------------------
587     CALCULATION FUNCTIONS
588     */
589    
590     /** @NOTE ANY function calling RelationBranchEvaluator should set
591 aw0a 1 glob_rel to point at the relation being evaluated. The calling
592 johnpye 669 function should also set glob_rel = NULL when it is done.
593     */
594 aw0a 1 static double RelationBranchEvaluator(struct relation_term *term)
595     {
596     assert(term != NULL);
597     switch(RelationTermType(term)) {
598     case e_func:
599 johnpye 846 /* CONSOLE_DEBUG("Evaluating term using FuncEval..."); */
600 aw0a 1 return FuncEval(TermFunc(term),
601     RelationBranchEvaluator(TermFuncLeft(term)) );
602     case e_var:
603     return TermVariable(glob_rel , term);
604     case e_int:
605     return (double)TermInteger(term);
606     case e_real:
607     return TermReal(term);
608     case e_zero:
609     return 0.0;
610 johnpye 669 case e_diff:
611     return 0.0;
612 aw0a 1
613     case e_plus:
614     return (RelationBranchEvaluator(TermBinLeft(term)) +
615     RelationBranchEvaluator(TermBinRight(term)));
616     case e_minus:
617     return (RelationBranchEvaluator(TermBinLeft(term)) -
618     RelationBranchEvaluator(TermBinRight(term)));
619     case e_times:
620     return (RelationBranchEvaluator(TermBinLeft(term)) *
621     RelationBranchEvaluator(TermBinRight(term)));
622     case e_divide:
623     return (RelationBranchEvaluator(TermBinLeft(term)) /
624     RelationBranchEvaluator(TermBinRight(term)));
625     case e_power:
626     case e_ipower:
627     return pow( RelationBranchEvaluator(TermBinLeft(term)) ,
628     RelationBranchEvaluator(TermBinRight(term)) );
629     case e_uminus:
630     return - RelationBranchEvaluator(TermBinLeft(term));
631     default:
632     FPRINTF(ASCERR, "error in RelationBranchEvaluator routine\n");
633     FPRINTF(ASCERR, "relation term type not recognized\n");
634     return 0.0;
635     }
636     }
637    
638 johnpye 669 /**
639     This function is passed a relation pointer, r, a pointer, pos, to a
640     position in the postfix version of the relation (0<=pos<length), and
641     a flag, lhs, telling whether we are interested in the left(=1) or
642     right(=0) side of the relation. This function will tranverse and
643     evaluate the subtree rooted at pos and will return the value as a
644     double. To do its evaluation, this function goes backwards through
645     the postfix representation of relation and calls itself at each
646     node--creating a stack of function calls.
647    
648     @NOTE: This function changes the value of pos---to the position of
649     the deepest leaf visited
650     */
651 aw0a 1 static double
652     RelationEvaluatePostfixBranch(CONST struct relation *r,
653     unsigned long *pos,
654     int lhs)
655     {
656     CONST struct relation_term *term; /* the current term */
657     CONST struct Func *funcptr; /* a pointer to a function */
658     double x, y; /* temporary values */
659    
660     term = NewRelationTerm(r, *pos, lhs);
661     assert(term != NULL);
662     switch( RelationTermType(term) ) {
663     case e_zero:
664     case e_real:
665     return TermReal(term);
666 johnpye 669 case e_diff:
667     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Can't evaluate DIFF(...)");
668     return 0;
669 aw0a 1 case e_int:
670     return TermInteger(term);
671     case e_var:
672     return TermVariable(r, term);
673     case e_plus:
674     (*pos)--;
675     y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y==right-side of '+' */
676     (*pos)--;
677     return RelationEvaluatePostfixBranch(r, pos, lhs) + y;
678     case e_minus:
679     (*pos)--;
680     y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y==right-side of '-' */
681     (*pos)--;
682     return RelationEvaluatePostfixBranch(r, pos, lhs) - y;
683     case e_times:
684     (*pos)--;
685     y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y==right-side of '*' */
686     (*pos)--;
687     return RelationEvaluatePostfixBranch(r, pos, lhs) * y;
688     case e_divide:
689     (*pos)--;
690     y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y is the divisor */
691     (*pos)--;
692     return RelationEvaluatePostfixBranch(r, pos, lhs) / y;
693     case e_power:
694     (*pos)--;
695     y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y is the power */
696     (*pos)--;
697     x = RelationEvaluatePostfixBranch(r, pos, lhs); /* x is the base */
698     return pow(x, y);
699     case e_ipower:
700     (*pos)--;
701     y = RelationEvaluatePostfixBranch(r, pos, lhs); /* y is the power */
702     (*pos)--;
703     x = RelationEvaluatePostfixBranch(r, pos, lhs); /* x is the base */
704     return asc_ipow(x, (int)y);
705     case e_uminus:
706     (*pos)--;
707     return -1.0 * RelationEvaluatePostfixBranch(r, pos, lhs);
708     case e_func:
709     funcptr = TermFunc(term);
710     (*pos)--;
711     return FuncEval(funcptr, RelationEvaluatePostfixBranch(r, pos, lhs));
712     default:
713     Asc_Panic(2, NULL,
714     "Don't know this type of relation type\n"
715     "in function RelationEvaluatePostfixBranch\n");
716 johnpye 846
717 aw0a 1 break;
718     }
719     }
720    
721 johnpye 726 #define RECURSE RelationEvaluatePostfixBranchSafe
722 aw0a 1 static double
723     RelationEvaluatePostfixBranchSafe(CONST struct relation *r,
724     unsigned long *pos,
725     int lhs,
726 johnpye 726 enum safe_err *serr
727     ){
728 aw0a 1 CONST struct relation_term *term; /* the current term */
729     double x, y; /* temporary values */
730    
731     term = NewRelationTerm(r, *pos, lhs);
732     assert(term != NULL);
733     switch( RelationTermType(term) ) {
734     case e_zero:
735     case e_real:
736     return TermReal(term);
737 johnpye 669 case e_diff:
738     return 0; /* for the moment, all derivatives evaluate to zero */
739 aw0a 1 case e_int:
740     return TermInteger(term);
741     case e_var:
742     return TermVariable(r, term);
743     case e_plus:
744 johnpye 726 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = RHS of '+' */
745     (*pos)--; return safe_add_D0(RECURSE(r,pos,lhs,serr), y, serr);
746 aw0a 1 case e_minus:
747 johnpye 726 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = RHS of '-' */
748     (*pos)--; return safe_sub_D0(RECURSE(r,pos,lhs,serr), y, serr);
749 aw0a 1 case e_times:
750 johnpye 726 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = RHS of '*' */
751     (*pos)--; return safe_mul_D0(RECURSE(r,pos,lhs,serr), y, serr);
752 aw0a 1 case e_divide:
753 johnpye 726 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = the divisor (RHS of '/') */
754     (*pos)--; return safe_div_D0(RECURSE(r,pos,lhs,serr), y, serr);
755 aw0a 1 case e_power:
756 johnpye 726 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* = the power (RHS of '^') */
757     (*pos)--; x = RECURSE(r, pos, lhs, serr); /* = the base */
758 aw0a 1 return safe_pow_D0(x, y, serr);
759     case e_ipower:
760 johnpye 726 (*pos)--; y = RECURSE(r, pos, lhs, serr); /* y is the power (RHS of '^') */
761     (*pos)--; x = RECURSE(r, pos, lhs, serr); /* x is the base */
762 aw0a 1 return safe_ipow_D0(x, (int)y, serr);
763     case e_uminus:
764 johnpye 726 (*pos)--; return -1.0 * RECURSE(r, pos, lhs, serr);
765 aw0a 1 case e_func:
766 johnpye 726 (*pos)--; return FuncEvalSafe(TermFunc(term),RECURSE(r,pos,lhs,serr),serr);
767 aw0a 1 default:
768 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
769 aw0a 1 }
770     }
771 johnpye 726 #undef RECURSE
772 aw0a 1
773 johnpye 669 /**
774     Yet another function for calculating the residual of a relation.
775     This function also uses the postfix version of the relations, but it
776     manages a stack(array) of doubles and calculates the residual in this
777     stack; therefore the function is not recursive. If the funtion
778     cannot allocate memory for its stack, it returns 0.0, so there is
779     currently no way of knowing if this function failed.
780     */
781 aw0a 1 static double
782     RelationEvaluateResidualPostfix(CONST struct relation *r)
783     {
784     unsigned long t; /* the current term in the relation r */
785     int lhs; /* looking at left(=1) or right(=0) hand side */
786     double *res_stack; /* the stack we use for evaluating the residual */
787     long s = -1; /* the top position in the stacks */
788     unsigned long length_lhs, length_rhs;
789     CONST struct relation_term *term;
790     CONST struct Func *funcptr;
791    
792    
793     length_lhs = RelationLength(r, 1);
794     length_rhs = RelationLength(r, 0);
795     if( (length_lhs+length_rhs) == 0 ) return 0.0;
796    
797     /* create the stacks */
798     res_stack = tmpalloc_array((1+MAX(length_lhs,length_rhs)),double);
799     if( res_stack == NULL ) return 0.0;
800    
801     lhs = 1;
802     t = 0;
803     while (1) {
804     if( lhs && (t >= length_lhs) ) {
805     /* finished processing left hand side, switch to right if it exists */
806     if( length_rhs ) {
807     lhs = t = 0;
808     }
809     else {
810     /* do not need to check for s>=0, since we know that
811     * (length_lhs+length_rhs>0) and that (length_rhs==0), the
812     * length_lhs must be > 0, thus s>=0
813     */
814     return (res_stack[s]);
815     }
816     }
817     else if( (!lhs) && (t >= length_rhs) ) {
818     /* finished processing right hand side */
819     if( length_lhs ) {
820     /* we know length_lhs and length_rhs are both > 0, since if
821     * length_rhs == 0, we would have exited above.
822     */
823     return (res_stack[s-1] - res_stack[s]);
824     }
825     else {
826     /* do not need to check for s>=0, since we know that
827     * (length_lhs+length_rhs>0) and that (length_lhs==0), the
828     * length_rhs must be > 0, thus s>=0
829     */
830     return (-1.0 * res_stack[s]);
831     }
832     }
833    
834     term = NewRelationTerm(r, t++, lhs);
835     switch( RelationTermType(term) ) {
836     case e_zero:
837 johnpye 737 s++; res_stack[s] = 0.0; break;
838 aw0a 1 case e_real:
839 johnpye 737 s++; res_stack[s] = TermReal(term); break;
840 aw0a 1 case e_int:
841 johnpye 737 s++; res_stack[s] = TermInteger(term); break;
842 aw0a 1 case e_var:
843 johnpye 737 s++; res_stack[s] = TermVariable(r, term); break;
844 aw0a 1 case e_plus:
845 johnpye 737 res_stack[s-1] += res_stack[s]; s--; break;
846 aw0a 1 case e_minus:
847 johnpye 737 res_stack[s-1] -= res_stack[s]; s--; break;
848 aw0a 1 case e_times:
849 johnpye 737 res_stack[s-1] *= res_stack[s]; s--; break;
850 aw0a 1 case e_divide:
851 johnpye 737 res_stack[s-1] /= res_stack[s]; s--; break;
852 aw0a 1 case e_uminus:
853 johnpye 737 res_stack[s] *= -1.0; break;
854 aw0a 1 case e_power:
855     case e_ipower:
856 johnpye 737 res_stack[s-1] = pow(res_stack[s-1], res_stack[s]); s--; break;
857 aw0a 1 case e_func:
858     funcptr = TermFunc(term);
859     res_stack[s] = FuncEval(funcptr, res_stack[s]);
860     break;
861     default:
862 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
863 aw0a 1 }
864     }
865     }
866    
867 johnpye 726 /*------------------------------------------------------------------------------
868     GRADIENT AND DERIVATIVE CALCULATIONS
869     */
870 aw0a 1
871 johnpye 669 /**
872     This function evaluates the residual and the gradient for the relation
873     r. The calling function must provide a pointer to a double for the
874     residual and an array of doubles for the gradients. This function
875     assumes r exists and that the pointers to the residual and gradients
876     are not NULL. This function returns 0 is everythings goes o.k., and
877     1 otherwise (out of memory). The function computes the gradients by
878     maintaining a n stacks, where n = (number-of-variables-in-r + 1)
879     The +1 is for the residual. The stacks come from a single array which
880     this function gets by calling tmpalloc_array. Two macros are defined
881     to make referencing this array easier.
882     */
883 aw0a 1 static int
884     RelationEvaluateResidualGradient(CONST struct relation *r,
885     double *residual,
886     double *gradient)
887     {
888     unsigned long t; /* the current term in the relation r */
889     unsigned long num_var; /* the number of variables in the relation r */
890     unsigned long v; /* the index of the variable we are looking at */
891     int lhs; /* looking at left(=1) or right(=0) hand side of r */
892     double *stacks; /* the memory for the stacks */
893     unsigned long stack_height; /* height of each stack */
894     long s = -1; /* the top position in the stacks */
895     double temp, temp2; /* temporary variables to speed gradient calculatns */
896     unsigned long length_lhs, length_rhs;
897     CONST struct relation_term *term;
898     CONST struct Func *fxnptr;
899    
900     num_var = NumberVariables(r);
901     length_lhs = RelationLength(r, 1);
902     length_rhs = RelationLength(r, 0);
903     if( (length_lhs + length_rhs) == 0 ) {
904     for( v = 0; v < num_var; v++ ) gradient[v] = 0.0;
905     *residual = 0.0;
906     return 0;
907     }
908     else {
909     stack_height = 1 + MAX(length_lhs,length_rhs);
910     }
911    
912     /* create the stacks */
913     stacks = tmpalloc_array(((num_var+1)*stack_height),double);
914     if( stacks == NULL ) return 1;
915    
916     #define res_stack(s) stacks[(s)]
917     #define grad_stack(v,s) stacks[((v)*stack_height)+(s)]
918    
919     lhs = 1;
920     t = 0;
921     while(1) {
922     if( lhs && (t >= length_lhs) ) {
923     /* need to switch to the right hand side--if it exists */
924     if( length_rhs ) {
925     lhs = t = 0;
926     }
927     else {
928     /* Set the pointers we were passed to the tops of the stacks.
929     * We do not need to check for s>=0, since we know that
930     * (length_lhs+length_rhs>0) and that (length_rhs==0), the
931     * length_lhs must be > 0, thus s>=0
932     */
933     for( v = 1; v <= num_var; v++ ) gradient[v-1] = grad_stack(v,s);
934     *residual = res_stack(s);
935     return 0;
936     }
937     }
938     else if( (!lhs) && (t >= length_rhs) ) {
939     /* we have processed both sides, quit */
940     if( length_lhs ) {
941     /* Set the pointers we were passed to lhs - rhs
942     * We know length_lhs and length_rhs are both > 0, since if
943     * length_rhs == 0, we would have exited above.
944     */
945     for( v = 1; v <= num_var; v++ ) {
946     gradient[v-1] = grad_stack(v,s-1) - grad_stack(v,s);
947     }
948     *residual = res_stack(s-1) - res_stack(s);
949     return 0;
950     }
951     else {
952     /* Set the pointers we were passed to -1.0 * top of stacks.
953     * We do not need to check for s>=0, since we know that
954     * (length_lhs+length_rhs>0) and that (length_lhs==0), the
955     * length_rhs must be > 0, thus s>=0
956     */
957     for( v = 1; v <= num_var; v++ ) {
958     gradient[v-1] = -grad_stack(v,s);
959     }
960     *residual = -res_stack(s);
961     return 0;
962     }
963     }
964    
965     term = NewRelationTerm(r, t++, lhs);
966     switch( RelationTermType(term) ) {
967     case e_zero:
968     s++;
969     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
970     res_stack(s) = 0.0;
971     break;
972     case e_real:
973     s++;
974     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
975     res_stack(s) = TermReal(term);
976     break;
977     case e_int:
978     s++;
979     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
980     res_stack(s) = TermInteger(term);
981     break;
982     case e_var:
983     s++;
984     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
985     grad_stack(TermVarNumber(term),s) = 1.0;
986     res_stack(s) = TermVariable(r, term);
987     break;
988     case e_plus:
989     /* d(u+v) = du + dv */
990     for( v = 1; v <= num_var; v++ ) grad_stack(v,s-1) += grad_stack(v,s);
991     res_stack(s-1) += res_stack(s);
992     s--;
993     break;
994     case e_minus:
995     /* d(u-v) = du - dv */
996     for( v = 1; v <= num_var; v++ ) grad_stack(v,s-1) -= grad_stack(v,s);
997     res_stack(s-1) -= res_stack(s);
998     s--;
999     break;
1000     case e_times:
1001     /* d(u*v) = u*dv + v*du */
1002     for( v = 1; v <= num_var; v++ ) {
1003     grad_stack(v,s-1) = ((res_stack(s-1) * grad_stack(v,s)) +
1004     (res_stack(s) * grad_stack(v,s-1)));
1005     }
1006     res_stack(s-1) *= res_stack(s);
1007     s--;
1008     break;
1009     case e_divide:
1010     /* d(u/v) = du/v - u*dv/(v^2) = (1/v) * [du - (u/v)*dv] */
1011     res_stack(s) = 1.0 / res_stack(s); /* 1/v */
1012     res_stack(s-1) *= res_stack(s); /* u/v */
1013     for( v = 1; v <= num_var; v++ ) {
1014     grad_stack(v,s-1) = (res_stack(s) *
1015     (grad_stack(v,s-1) -
1016     (res_stack(s-1) * grad_stack(v,s))));
1017     }
1018     s--;
1019     break;
1020     case e_uminus:
1021     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = -grad_stack(v,s);
1022     res_stack(s) = -res_stack(s);
1023     break;
1024     case e_power:
1025     /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1026     /* First compute: v*u^(v-1) */
1027     temp = res_stack(s) * pow( res_stack(s-1), (res_stack(s) - 1.0) );
1028     /* Now compute: ln(u) */
1029     temp2 = FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1030     /* Next compute: u^v */
1031     res_stack(s-1) = pow(res_stack(s-1), res_stack(s));
1032     /* Compute: [ln(u)] * [u^v] */
1033     temp2 *= res_stack(s-1);
1034     /* Finally, compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1035     for( v = 1; v <= num_var; v++ ) {
1036     grad_stack(v,s-1) = ((temp * grad_stack(v,s-1)) +
1037     (temp2 * grad_stack(v,s)));
1038     }
1039     s--;
1040     break;
1041     case e_ipower:
1042     /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1043     /* First compute: v*u^(v-1) */
1044     temp = asc_d1ipow( res_stack(s-1), ((int)res_stack(s)) );
1045     /* Now compute: ln(u) */
1046     temp2 = FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1047     /* Next compute: u^v */
1048     res_stack(s-1) = asc_ipow( res_stack(s-1), ((int)res_stack(s)) );
1049     /* Compute: [ln(u)] * [u^v] */
1050     temp2 *= res_stack(s-1);
1051     /* Finally, compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1052     for( v = 1; v <= num_var; v++ ) {
1053     grad_stack(v,s-1) = ((temp * grad_stack(v,s-1)) +
1054     (temp2 * grad_stack(v,s)));
1055     }
1056     s--;
1057     break;
1058     case e_func:
1059     /*
1060     funcptr = TermFunc(term);
1061     for (v = 0; v < num_var; v++) {
1062     grad_stack(v,s) = FuncDeriv(funcptr, grad_stack(v,s));
1063     }
1064     res_stack(s) = FuncEval(funcptr, res_stack(s)); */
1065     fxnptr = TermFunc(term);
1066     temp = FuncDeriv( fxnptr, res_stack(s) );
1067     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) *= temp;
1068     res_stack(s) = FuncEval( fxnptr, res_stack(s) );
1069     break;
1070     default:
1071 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
1072 aw0a 1 break;
1073     }
1074     }
1075     #undef grad_stack
1076     #undef res_stack
1077     }
1078    
1079     static int
1080     RelationEvaluateResidualGradientSafe(CONST struct relation *r,
1081     double *residual,
1082     double *gradient,
1083     enum safe_err *serr)
1084     {
1085     unsigned long t; /* the current term in the relation r */
1086     unsigned long num_var; /* the number of variables in the relation r */
1087     unsigned long v; /* the index of the variable we are looking at */
1088     int lhs; /* looking at left(=1) or right(=0) hand side of r */
1089     double *stacks; /* the memory for the stacks */
1090     unsigned long stack_height; /* height of each stack */
1091     long s = -1; /* the top position in the stacks */
1092     double temp, temp2; /* temporary variables to speed gradient calculatns */
1093     unsigned long length_lhs, length_rhs;
1094     CONST struct relation_term *term;
1095     CONST struct Func *fxnptr;
1096    
1097     num_var = NumberVariables(r);
1098     length_lhs = RelationLength(r, 1);
1099     length_rhs = RelationLength(r, 0);
1100     if( (length_lhs + length_rhs) == 0 ) {
1101     for( v = 0; v < num_var; v++ ) gradient[v] = 0.0;
1102     *residual = 0.0;
1103     return 0;
1104     }
1105     else {
1106     stack_height = 1 + MAX(length_lhs,length_rhs);
1107     }
1108    
1109     /* create the stacks */
1110     stacks = tmpalloc_array(((num_var+1)*stack_height),double);
1111     if( stacks == NULL ) return 1;
1112    
1113     #define res_stack(s) stacks[(s)]
1114     #define grad_stack(v,s) stacks[((v)*stack_height)+(s)]
1115    
1116     lhs = 1;
1117     t = 0;
1118     while(1) {
1119     if( lhs && (t >= length_lhs) ) {
1120     /* need to switch to the right hand side--if it exists */
1121     if( length_rhs ) {
1122     lhs = t = 0;
1123     }
1124     else {
1125     /* Set the pointers we were passed to the tops of the stacks.
1126     * We do not need to check for s>=0, since we know that
1127     * (length_lhs+length_rhs>0) and that (length_rhs==0), the
1128     * length_lhs must be > 0, thus s>=0
1129     */
1130     for( v = 1; v <= num_var; v++ ) gradient[v-1] = grad_stack(v,s);
1131     *residual = res_stack(s);
1132     return 0;
1133     }
1134     }
1135     else if( (!lhs) && (t >= length_rhs) ) {
1136     /* we have processed both sides, quit */
1137     if( length_lhs ) {
1138     /* Set the pointers we were passed to lhs - rhs
1139     * We know length_lhs and length_rhs are both > 0, since if
1140     * length_rhs == 0, we would have exited above.
1141     */
1142     for( v = 1; v <= num_var; v++ ) {
1143     gradient[v-1] = safe_sub_D0(grad_stack(v,s-1),grad_stack(v,s),serr);
1144     }
1145     *residual = safe_sub_D0(res_stack(s-1), res_stack(s), serr);
1146     return 0;
1147     }
1148     else {
1149     /* Set the pointers we were passed to -1.0 * top of stacks.
1150     * We do not need to check for s>=0, since we know that
1151     * (length_lhs+length_rhs>0) and that (length_lhs==0), the
1152     * length_rhs must be > 0, thus s>=0
1153     */
1154     for( v = 1; v <= num_var; v++ ) {
1155     gradient[v-1] = -grad_stack(v,s);
1156     }
1157     *residual = -res_stack(s);
1158     return 0;
1159     }
1160     }
1161    
1162     term = NewRelationTerm(r, t++, lhs);
1163     switch( RelationTermType(term) ) {
1164     case e_zero:
1165     s++;
1166     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1167     res_stack(s) = 0.0;
1168     break;
1169     case e_real:
1170     s++;
1171     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1172     res_stack(s) = TermReal(term);
1173     break;
1174     case e_int:
1175     s++;
1176     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1177     res_stack(s) = TermInteger(term);
1178     break;
1179     case e_var:
1180     s++;
1181     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = 0.0;
1182     grad_stack(TermVarNumber(term),s) = 1.0;
1183     res_stack(s) = TermVariable(r, term);
1184     break;
1185     case e_plus:
1186     /* d(u+v) = du + dv */
1187     for( v = 1; v <= num_var; v++ ) {
1188     grad_stack(v,s-1)=safe_add_D0(grad_stack(v,s-1),grad_stack(v,s),serr);
1189     }
1190     res_stack(s-1) = safe_add_D0(res_stack(s-1),res_stack(s),serr);
1191     s--;
1192     break;
1193     case e_minus:
1194     /* d(u-v) = du - dv */
1195     for( v = 1; v <= num_var; v++ ) {
1196     grad_stack(v,s-1)=safe_sub_D0(grad_stack(v,s-1),grad_stack(v,s),serr);
1197     }
1198     res_stack(s-1) = safe_sub_D0(res_stack(s-1),res_stack(s),serr);
1199     s--;
1200     break;
1201     case e_times:
1202     /* d(u*v) = u*dv + v*du */
1203     for( v = 1; v <= num_var; v++ ) {
1204     grad_stack(v,s-1) =
1205     safe_add_D0(safe_mul_D0(res_stack(s-1),grad_stack(v,s),serr),
1206     safe_mul_D0(res_stack(s),grad_stack(v,s-1),serr),
1207     serr);
1208     }
1209     res_stack(s-1) = safe_mul_D0(res_stack(s-1),res_stack(s),serr);
1210     s--;
1211     break;
1212     case e_divide:
1213     /* d(u/v) = du/v - u*dv/(v^2) = (1/v) * [du - (u/v)*dv] */
1214     res_stack(s) = safe_rec(res_stack(s),serr); /* 1/v */
1215     res_stack(s-1) = safe_mul_D0(res_stack(s-1),res_stack(s),serr); /* u/v */
1216     for( v = 1; v <= num_var; v++ ) {
1217     grad_stack(v,s-1) =
1218     safe_mul_D0(res_stack(s),
1219     safe_sub_D0(grad_stack(v,s-1),
1220     safe_mul_D0(res_stack(s-1),
1221     grad_stack(v,s),
1222     serr),serr),serr);
1223     }
1224     s--;
1225     break;
1226     case e_uminus:
1227     for( v = 1; v <= num_var; v++ ) grad_stack(v,s) = -grad_stack(v,s);
1228     res_stack(s) = -res_stack(s);
1229     break;
1230     case e_power:
1231     /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1232     /* v*u^(v-1) */
1233     temp = safe_pow_D1( res_stack(s-1), res_stack(s), 0, serr );
1234     /* ln(u)*u^v */
1235     temp2 = safe_pow_D1( res_stack(s-1), res_stack(s), 1, serr );
1236     /* Compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1237     for( v = 1; v <= num_var; v++ ) {
1238     grad_stack(v,s-1) =
1239     safe_add_D0(safe_mul_D0(temp, grad_stack(v,s-1), serr),
1240     safe_mul_D0(temp2, grad_stack(v,s), serr), serr);
1241     }
1242     /* u^v */
1243     res_stack(s-1) = safe_pow_D0( res_stack(s-1), res_stack(s), serr );
1244     s--;
1245     break;
1246     case e_ipower:
1247     /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1248     /* v*u^(v-1) */
1249     temp = safe_ipow_D1( res_stack(s-1), res_stack(s), 0, serr );
1250     /* ln(u)*u^v */
1251     temp2 = safe_ipow_D1( res_stack(s-1), res_stack(s), 1, serr );
1252     /* Compute: [v*u^(v-1)] * [du] + [ln(u)*u^v] * [dv] */
1253     for( v = 1; v <= num_var; v++ ) {
1254     grad_stack(v,s-1) =
1255     safe_add_D0(safe_mul_D0(temp, grad_stack(v,s-1), serr),
1256     safe_mul_D0(temp2, grad_stack(v,s), serr), serr);
1257     }
1258     /* Next compute: u^v */
1259     res_stack(s-1) = safe_ipow_D0( res_stack(s-1), res_stack(s), serr );
1260     s--;
1261     break;
1262     case e_func:
1263     fxnptr = TermFunc(term);
1264     temp = FuncDerivSafe( fxnptr, res_stack(s), serr);
1265     for( v = 1; v <= num_var; v++ ) {
1266     grad_stack(v,s) = safe_mul_D0( grad_stack(v,s), temp, serr );
1267     }
1268     res_stack(s) = FuncEvalSafe( fxnptr, res_stack(s), serr);
1269     break;
1270     default:
1271 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
1272 aw0a 1 }
1273     }
1274     #undef grad_stack
1275     #undef res_stack
1276     }
1277    
1278 johnpye 669 /**
1279     This function evaluates and returns the derivative of the
1280     relation r with respect to the variable whose index is pos.
1281     This function assumes r exists and that pos is within the proper range.
1282     The function computes the gradients by maintaining 2 stacks, one for
1283     the residual and one for the derivative. The stacks come from a
1284     single array which this gets by calling tmpalloc_array. Two macros
1285     are defined to make referencing this array easier. Of the malloc fails,
1286     this function returns 0.0, so there is currently no way to know if
1287     the function failed.
1288     */
1289 aw0a 1 static double
1290     RelationEvaluateDerivative(CONST struct relation *r,
1291     unsigned long pos)
1292     {
1293     unsigned long t; /* the current term in the relation r */
1294     int lhs; /* looking at left(=1) or right(=0) hand side of r */
1295     double *stacks; /* the memory for the stacks */
1296     unsigned long stack_height; /* height of each stack */
1297     long s = -1; /* the top position in the stacks */
1298     unsigned long length_lhs, length_rhs;
1299     CONST struct relation_term *term;
1300     CONST struct Func *fxnptr;
1301    
1302     length_lhs = RelationLength(r, 1);
1303     length_rhs = RelationLength(r, 0);
1304     if( (length_lhs + length_rhs) == 0 ) {
1305     return 0.0;
1306     }
1307     else {
1308     stack_height = 1 + MAX(length_lhs,length_rhs);
1309     }
1310    
1311     /* create the stacks */
1312     stacks = tmpalloc_array((2*stack_height),double);
1313     if( stacks == NULL ) return 0.0;
1314    
1315     #define res_stack(s) stacks[(s)]
1316     #define grad_stack(s) stacks[stack_height+(s)]
1317    
1318     lhs = 1;
1319     t = 0;
1320     while(1) {
1321     if( lhs && (t >= length_lhs) ) {
1322     /* need to switch to the right hand side--if it exists */
1323     if( length_rhs ) {
1324     lhs = t = 0;
1325     }
1326     else {
1327     /* do not need to check for s>=0, since we know that
1328     * (length_lhs+length_rhs>0) and that (length_rhs==0), the
1329     * length_lhs must be > 0, thus s>=0
1330     */
1331     return grad_stack(s);
1332     }
1333     }
1334     else if( (!lhs) && (t >= length_rhs) ) {
1335     /* we have processed both sides, quit */
1336     if( length_lhs ) {
1337     /* we know length_lhs and length_rhs are both > 0, since if
1338     * length_rhs == 0, we would have exited above.
1339     */
1340     return (grad_stack(s-1) - grad_stack(s));
1341     }
1342     else {
1343     /* do not need to check for s>=0, since we know that
1344     * (length_lhs+length_rhs>0) and that (length_lhs==0), the
1345     * length_rhs must be > 0, thus s>=0
1346     */
1347     return (-1.0 * grad_stack(s));
1348     }
1349     }
1350    
1351     term = NewRelationTerm(r, t++, lhs);
1352     switch( RelationTermType(term) ) {
1353     case e_zero:
1354     s++;
1355     grad_stack(s) = res_stack(s) = 0.0;
1356     break;
1357     case e_real:
1358     s++;
1359     grad_stack(s) = 0.0;
1360     res_stack(s) = TermReal(term);
1361     break;
1362     case e_int:
1363     s++;
1364     grad_stack(s) = 0.0;
1365     res_stack(s) = TermInteger(term);
1366     break;
1367     case e_var:
1368     s++;
1369     grad_stack(s) = ( (pos == TermVarNumber(term)) ? 1.0 : 0.0 );
1370     res_stack(s) = TermVariable(r, term);
1371     break;
1372     case e_plus:
1373     /* d(u+v) = du + dv */
1374     grad_stack(s-1) += grad_stack(s);
1375     res_stack(s-1) += res_stack(s);
1376     s--;
1377     break;
1378     case e_minus:
1379     /* d(u-v) = du - dv */
1380     grad_stack(s-1) -= grad_stack(s);
1381     res_stack(s-1) -= res_stack(s);
1382     s--;
1383     break;
1384     case e_times:
1385     /* d(u*v) = u*dv + v*du */
1386     grad_stack(s-1) = ((res_stack(s-1) * grad_stack(s)) +
1387     (res_stack(s) * grad_stack(s-1)));
1388     res_stack(s-1) *= res_stack(s);
1389     s--;
1390     break;
1391     case e_divide:
1392     /* d(u/v) = du/v - u*dv/(v^2) = [du - (u/v)*dv]/v */
1393     res_stack(s-1) = res_stack(s-1) / res_stack(s);
1394     grad_stack(s-1) = ((grad_stack(s-1) - (res_stack(s-1) * grad_stack(s))) /
1395     res_stack(s));
1396     s--;
1397     break;
1398     case e_uminus:
1399     grad_stack(s) = -grad_stack(s);
1400     res_stack(s) = -res_stack(s);
1401     break;
1402     case e_power:
1403     /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1404     /* First we compute: v*u^(v-1)*du */
1405     grad_stack(s-1) *= (pow(res_stack(s-1), (res_stack(s) - 1.0)) *
1406     res_stack(s));
1407     /* Now compute: ln(u)*dv */
1408     grad_stack(s) *= FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1409     /* Next compute: u^v */
1410     res_stack(s-1) = pow( res_stack(s-1), res_stack(s) );
1411     /* Finally, compute: [v*u^(v-1)*du] + [u^v] * [ln(u)*dv] */
1412     grad_stack(s-1) += (res_stack(s-1) * grad_stack(s));
1413     s--;
1414     break;
1415     case e_ipower:
1416     /* d(x^y) = y * dx * x^(y-1) + ln(x) * dy * x^y */
1417     /* First we compute: v*u^(v-1)*du */
1418     grad_stack(s-1) *= asc_d1ipow( res_stack(s-1), ((int)res_stack(s)) );
1419     /* Now compute: ln(u)*dv */
1420     grad_stack(s) *= FuncEval( LookupFuncById(F_LN), res_stack(s-1) );
1421     /* Next compute: u^v */
1422     res_stack(s-1) = asc_ipow( res_stack(s-1), ((int)res_stack(s)) );
1423     /* Finally, compute: [v*u^(v-1)*du] + [u^v] * [ln(u)*dv] */
1424     grad_stack(s-1) += (res_stack(s-1) * grad_stack(s));
1425     s--;
1426     break;
1427     case e_func:
1428     fxnptr = TermFunc(term);
1429     grad_stack(s) *= FuncDeriv( fxnptr, res_stack(s) );
1430     res_stack(s) = FuncEval( fxnptr, res_stack(s) );
1431     break;
1432     default:
1433 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
1434 aw0a 1 break;
1435     }
1436     }
1437     #undef grad_stack
1438     #undef res_stack
1439     }
1440    
1441     static double
1442     RelationEvaluateDerivativeSafe(CONST struct relation *r,
1443     unsigned long pos,
1444     enum safe_err *serr)
1445     {
1446     unsigned long t; /* the current term in the relation r */
1447     int lhs; /* looking at left(=1) or right(=0) hand side of r */
1448     double *stacks; /* the memory for the stacks */
1449     unsigned long stack_height; /* height of each stack */
1450     long s = -1; /* the top position in the stacks */
1451     unsigned long length_lhs, length_rhs;
1452     CONST struct relation_term *term;
1453     CONST struct Func *fxnptr;
1454    
1455     length_lhs = RelationLength(r, 1);
1456     length_rhs = RelationLength(r, 0);
1457     if( (length_lhs + length_rhs) == 0 ) {
1458     return 0.0;
1459     }
1460     else {
1461     stack_height = 1 + MAX(length_lhs,length_rhs);
1462     }
1463    
1464     /* create the stacks */
1465     stacks = tmpalloc_array((2*stack_height),double);
1466     if( stacks == NULL ) return 0.0;
1467    
1468     #define res_stack(s) stacks[(s)]
1469     #define grad_stack(s) stacks[stack_height+(s)]
1470    
1471     lhs = 1;
1472     t = 0;
1473     while(1) {
1474     if( lhs && (t >= length_lhs) ) {
1475     /* need to switch to the right hand side--if it exists */
1476     if( length_rhs ) {
1477     lhs = t = 0;
1478     }
1479     else {
1480     /* do not need to check for s>=0, since we know that
1481     * (length_lhs+length_rhs>0) and that (length_rhs==0), the
1482     * length_lhs must be > 0, thus s>=0
1483     */
1484     return grad_stack(s);
1485     }
1486     }
1487     else if( (!lhs) && (t >= length_rhs) ) {
1488     /* we have processed both sides, quit */
1489     if( length_lhs ) {
1490     /* we know length_lhs and length_rhs are both > 0, since if
1491     * length_rhs == 0, we would have exited above.
1492     */
1493     return safe_sub_D0(grad_stack(s-1), grad_stack(s), serr);
1494     }
1495     else {
1496     /* do not need to check for s>=0, since we know that
1497     * (length_lhs+length_rhs>0) and that (length_lhs==0), the
1498     * length_rhs must be > 0, thus s>=0
1499     */
1500     return (-grad_stack(s));
1501     }
1502     }
1503    
1504     term = NewRelationTerm(r, t++, lhs);
1505     switch( RelationTermType(term) ) {
1506     case e_zero:
1507     s++;
1508     grad_stack(s) = res_stack(s) = 0.0;
1509     break;
1510     case e_real:
1511     s++;
1512     grad_stack(s) = 0.0;
1513     res_stack(s) = TermReal(term);
1514     break;
1515     case e_int:
1516     s++;
1517     grad_stack(s) = 0.0;
1518     res_stack(s) = TermInteger(term);
1519     break;
1520     case e_var:
1521     s++;
1522     grad_stack(s) = ( (pos == TermVarNumber(term)) ? 1.0 : 0.0 );
1523     res_stack(s) = TermVariable(r, term);
1524     break;
1525     case e_plus:
1526     /* d(u+v) = du + dv */
1527     grad_stack(s-1) = safe_add_D0( grad_stack(s-1), grad_stack(s), serr );
1528     res_stack(s-1) = safe_add_D0( res_stack(s-1), res_stack(s), serr );
1529     s--;
1530     break;
1531     case e_minus:
1532     /* d(u-v) = du - dv */
1533     grad_stack(s-1) = safe_sub_D0( grad_stack(s-1), grad_stack(s), serr );
1534     res_stack(s-1) = safe_sub_D0( res_stack(s-1), res_stack(s), serr );
1535     s--;
1536     break;
1537     case e_times:
1538     /* d(u*v) = u*dv + v*du */
1539     grad_stack(s-1) =
1540     safe_add_D0(safe_mul_D0( res_stack(s-1), grad_stack(s), serr),
1541     safe_mul_D0( res_stack(s), grad_stack(s-1), serr), serr);
1542     res_stack(s-1) = safe_mul_D0( res_stack(s-1), res_stack(s), serr );
1543     s--;
1544     break;
1545     case e_divide:
1546     /* d(u/v) = du/v - u*dv/(v^2) = [du - (u/v)*dv]/v */
1547     res_stack(s-1) = safe_div_D0( res_stack(s-1), res_stack(s), serr);
1548     grad_stack(s-1) =
1549     safe_div_D0(safe_sub_D0(grad_stack(s-1),
1550     safe_mul_D0(res_stack(s-1),grad_stack(s),serr),
1551     serr),
1552     res_stack(s),serr);
1553     s--;
1554     break;
1555     case e_uminus:
1556     grad_stack(s) = -grad_stack(s);
1557     res_stack(s) = -res_stack(s);
1558     break;
1559     case e_power:
1560     /* d(u^v) = v * u^(v-1) * du + ln(u) * u^v * dv */
1561     grad_stack(s-1) =
1562     safe_add_D0(safe_mul_D0(safe_pow_D1(res_stack(s-1),
1563     res_stack(s),0,serr),
1564     grad_stack(s-1), serr),
1565     safe_mul_D0(safe_pow_D1(res_stack(s-1),
1566     res_stack(s),1,serr),
1567     grad_stack(s),serr), serr);
1568     /* u^v */
1569     res_stack(s-1) = safe_pow_D0( res_stack(s-1), res_stack(s), serr);
1570     s--;
1571     break;
1572     case e_ipower:
1573     /* d(x^y) = y * dx * x^(y-1) + ln(x) * dy * x^y */
1574     grad_stack(s-1) =
1575     safe_add_D0(safe_mul_D0(safe_ipow_D1(res_stack(s-1),
1576     res_stack(s),0,serr),
1577     grad_stack(s-1), serr),
1578     safe_mul_D0(safe_ipow_D1(res_stack(s-1),
1579     res_stack(s),1,serr),
1580     grad_stack(s),serr), serr);
1581     /* u^v */
1582     res_stack(s-1) = safe_ipow_D0( res_stack(s-1), res_stack(s), serr);
1583     s--;
1584     break;
1585     case e_func:
1586     fxnptr = TermFunc(term);
1587     grad_stack(s) = safe_mul_D0(FuncDerivSafe( fxnptr, res_stack(s), serr ),
1588     grad_stack(s), serr);
1589     res_stack(s) = FuncEvalSafe( fxnptr, res_stack(s), serr);
1590     break;
1591     default:
1592 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
1593 aw0a 1 }
1594     }
1595     #undef grad_stack
1596     #undef res_stack
1597     }
1598    
1599 johnpye 669 /*------------------------------------------------------------------------------
1600     RELATION & RELATION TERM QUERIES FUNCTIONS (FOR USE BY EXTERNAL CODE)
1601     */
1602 aw0a 1
1603 johnpye 726 /**
1604 johnpye 846 @return =, <, >, etc, etc. not e_token, e_glassbox, etc
1605 johnpye 726 */
1606 aw0a 1 enum Expr_enum RelationRelop(CONST struct relation *rel)
1607     {
1608     AssertAllocatedMemory(rel,sizeof(struct relation));
1609     return rel->share->s.relop;
1610     }
1611    
1612 johnpye 726 /**
1613     This query only applies to TokenRelations and OpcodeRelations.
1614     */
1615 aw0a 1 unsigned long RelationLength(CONST struct relation *rel, int lhs)
1616     {
1617     assert(rel!=NULL);
1618     AssertAllocatedMemory(rel,sizeof(struct relation));
1619     if (lhs){
1620     if (RTOKEN(rel).lhs) return (RTOKEN(rel).lhs_len);
1621     else return 0;
1622     }
1623     if (RTOKEN(rel).rhs) return (RTOKEN(rel).rhs_len);
1624     else return 0;
1625     }
1626    
1627 johnpye 908 struct gl_list_t *RelationBlackBoxFormalArgs(CONST struct relation *rel)
1628     {
1629     assert(rel!=NULL);
1630     return RelationBlackBoxCache(rel)->formalArgList;
1631     }
1632    
1633     struct BlackBoxCache * RelationBlackBoxCache(CONST struct relation *rel)
1634     {
1635     assert(rel!=NULL);
1636     AssertAllocatedMemory(rel,sizeof(struct relation));
1637     return ((struct BlackBoxData *) (rel->externalData))->common;
1638     }
1639    
1640     struct BlackBoxData * RelationBlackBoxData(CONST struct relation *rel)
1641     {
1642     assert(rel!=NULL);
1643     AssertAllocatedMemory(rel,sizeof(struct relation));
1644     return (struct BlackBoxData *) rel->externalData;
1645     }
1646    
1647 aw0a 1 /*
1648 johnpye 726 This query only applies to TokenRelations. It assumes the
1649     user still thinks tokens number from [1..len].
1650     */
1651 aw0a 1 CONST struct relation_term *RelationTerm(CONST struct relation *rel,
1652     unsigned long int pos, int lhs)
1653     {
1654     assert(rel!=NULL);
1655     AssertAllocatedMemory(rel,sizeof(struct relation));
1656     if (lhs){
1657     if (RTOKEN(rel).lhs)
1658     return A_TERM(&(RTOKEN(rel).lhs[pos-1]));
1659     else return NULL;
1660     }
1661     else{
1662     if (RTOKEN(rel).rhs)
1663     return A_TERM(&(RTOKEN(rel).rhs[pos-1]));
1664     else return NULL;
1665     }
1666     }
1667    
1668 johnpye 726 /**
1669     This query only applies to TokenRelations. It assumes the
1670     clued user thinks tokens number from [0..len-1], which they do.
1671     */
1672 aw0a 1 CONST struct relation_term
1673     *NewRelationTermF(CONST struct relation *rel, unsigned long pos, int lhs)
1674     {
1675     assert(rel!=NULL);
1676     AssertAllocatedMemory(rel,sizeof(struct relation));
1677     if (lhs){
1678     if (RTOKEN(rel).lhs != NULL)
1679     return A_TERM(&(RTOKEN(rel).lhs[pos]));
1680     else return NULL;
1681 johnpye 726 }else{
1682 aw0a 1 if (RTOKEN(rel).rhs != NULL)
1683     return A_TERM(&(RTOKEN(rel).rhs[pos]));
1684     else return NULL;
1685     }
1686     }
1687    
1688 johnpye 726 /**
1689     This query only applies to sides from TokenRelations. It assumes the
1690     clued user thinks tokens number from [0..len-1], which they do,
1691     and that the side came from a token relation instance.
1692     */
1693 aw0a 1 CONST struct relation_term
1694     *RelationSideTermF(CONST union RelationTermUnion *side, unsigned long pos)
1695     {
1696     assert(side!=NULL);
1697     return A_TERM(&(side[pos]));
1698     }
1699    
1700 johnpye 726 /**
1701     This query only applies to TokenRelations. It assumes the
1702     clued user thinks tokens number from [0..len-1], which they do.
1703     */
1704 aw0a 1 enum Expr_enum RelationTermTypeF(CONST struct relation_term *term)
1705     {
1706     AssertMemory(term);
1707     return term->t;
1708     }
1709    
1710     unsigned long TermVarNumber(CONST struct relation_term *term)
1711     {
1712     assert(term&&term->t == e_var);
1713     AssertMemory(term);
1714     return V_TERM(term)->varnum;
1715     }
1716    
1717 johnpye 726 long TermInteger(CONST struct relation_term *term){
1718 aw0a 1 assert(term&&(term->t==e_int));
1719     AssertMemory(term);
1720     return I_TERM(term)->ivalue;
1721     }
1722    
1723 johnpye 726 double TermReal(CONST struct relation_term *term){
1724 aw0a 1 assert(term&&( term->t==e_real || term->t==e_zero));
1725     AssertMemory(term);
1726     return R_TERM(term)->value;
1727     }
1728    
1729     double
1730 johnpye 726 TermVariable(CONST struct relation *rel, CONST struct relation_term *term){
1731 aw0a 1 return
1732     RealAtomValue((struct Instance*)RelationVariable(rel,TermVarNumber(term)));
1733     }
1734    
1735 johnpye 726 CONST dim_type *TermDimensions(CONST struct relation_term *term){
1736 aw0a 1 assert( term && (term->t==e_real || term->t == e_int || term->t == e_zero) );
1737     AssertMemory(term);
1738     if (term->t==e_real) return R_TERM(term)->dimensions;
1739     if (term->t==e_int) return Dimensionless();
1740     if (term->t==e_zero) return WildDimension();
1741     return NULL;
1742     }
1743    
1744 johnpye 726 CONST struct Func *TermFunc(CONST struct relation_term *term){
1745 aw0a 1 assert(term&&(term->t == e_func));
1746     AssertMemory(term);
1747     return F_TERM(term)->fptr;
1748     }
1749    
1750 johnpye 726 struct relation_term *RelationINF_Lhs(CONST struct relation *rel){
1751 aw0a 1 return RTOKEN(rel).lhs_term;
1752     }
1753    
1754 johnpye 726 struct relation_term *RelationINF_Rhs(CONST struct relation *rel){
1755 aw0a 1 return RTOKEN(rel).rhs_term;
1756     }
1757    
1758 johnpye 908 struct ExternalFunc *RelationBlackBoxExtFunc(CONST struct relation *rel)
1759     {
1760 aw0a 1 assert(rel!=NULL);
1761 johnpye 908 return RBBOX(rel).efunc;
1762 aw0a 1 }
1763    
1764 johnpye 726 /*------------------------------------------------------------------------------
1765     GLASSBOX RELATION QUERIES
1766    
1767     For picking apart a GlassBox relation.
1768     */
1769    
1770 johnpye 908 struct ExternalFunc *RelationGlassBoxExtFunc(CONST struct relation *rel){
1771 aw0a 1 assert(rel!=NULL);
1772     return RGBOX(rel).efunc;
1773     }
1774    
1775 johnpye 726 int GlassBoxRelIndex(CONST struct relation *rel){
1776 aw0a 1 assert(rel!=NULL);
1777     return RGBOX(rel).index;
1778     }
1779    
1780 johnpye 726 int *GlassBoxArgs(CONST struct relation *rel){
1781 aw0a 1 assert(rel!=NULL);
1782     return RGBOX(rel).args;
1783     }
1784    
1785    
1786 johnpye 726 /*------------------------------------------------------------------------------
1787     GENERAL RELATION QUERIES
1788 aw0a 1
1789 johnpye 726 Not all of these queries may be applied to all relation types.
1790     Those that cannot be are so marked.
1791     */
1792    
1793     CONST struct gl_list_t *RelationVarList(CONST struct relation *rel){
1794 aw0a 1 return (CONST struct gl_list_t *)rel->vars;
1795     }
1796    
1797 johnpye 726 dim_type *RelationDim(CONST struct relation *rel){
1798 aw0a 1 assert(rel!=NULL);
1799     return rel->d;
1800     }
1801    
1802 johnpye 908 int SetRelationDim(struct relation *rel, CONST dim_type *d)
1803     {
1804 aw0a 1 if (!rel) return 1;
1805 johnpye 908 rel->d = (dim_type*)d;
1806 aw0a 1 return 0;
1807     }
1808    
1809 johnpye 726 double RelationResidual(CONST struct relation *rel){
1810 aw0a 1 assert(rel!=NULL);
1811     return rel->residual;
1812     }
1813    
1814 johnpye 726 void SetRelationResidual(struct relation *rel, double value){
1815 aw0a 1 assert(rel!=NULL);
1816     rel->residual = value;
1817     }
1818    
1819 johnpye 726 double RelationMultiplier(CONST struct relation *rel){
1820 aw0a 1 assert(rel!=NULL);
1821     return rel->multiplier;
1822     }
1823    
1824 johnpye 726 void SetRelationMultiplier(struct relation *rel, double value){
1825 aw0a 1 assert(rel!=NULL);
1826     rel->multiplier = value;
1827     }
1828    
1829 johnpye 726 double RelationNominal(CONST struct relation *rel){
1830 aw0a 1 assert(rel!=NULL);
1831     return rel->nominal;
1832     }
1833    
1834 johnpye 726 void SetRelationNominal(struct relation *rel, double value){
1835 aw0a 1 assert(rel!=NULL);
1836     rel->nominal = (fabs(value) > 0.0) ? fabs(value) : rel->nominal;
1837     }
1838    
1839    
1840 johnpye 726 int RelationIsCond(CONST struct relation *rel){
1841 aw0a 1 if ( rel != NULL) {
1842     return rel->iscond;
1843     }
1844     return 0;
1845     }
1846    
1847 johnpye 726 void SetRelationIsCond(struct relation *rel){
1848 aw0a 1 if ( rel != NULL) {
1849     rel->iscond = 1;
1850 johnpye 726 }else{
1851     ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation");
1852 aw0a 1 }
1853     }
1854    
1855     unsigned long NumberVariables(CONST struct relation *rel)
1856     {
1857     unsigned long n;
1858     assert(rel!=NULL);
1859     n = (rel->vars!=NULL) ? gl_length(rel->vars) : 0;
1860     return n;
1861     }
1862    
1863     struct Instance *RelationVariable(CONST struct relation *rel,
1864 johnpye 726 unsigned long int varnum
1865     ){
1866 aw0a 1 assert(rel!=NULL);
1867     return (struct Instance *)gl_fetch(rel->vars,varnum);
1868     }
1869    
1870     static void CalcDepth(CONST struct relation *rel,
1871 johnpye 726 int lhs,
1872     unsigned long int *depth,
1873     unsigned long int *maxdepth
1874     ){
1875 aw0a 1 unsigned long c,length;
1876     CONST struct relation_term *term;
1877     length = RelationLength(rel,lhs);
1878     for(c=0;c<length;c++){
1879     term = NewRelationTerm(rel,c,lhs);
1880     switch(RelationTermType(term)){
1881     case e_zero:
1882     case e_int:
1883     case e_real:
1884     case e_var:
1885     if (++(*depth) > *maxdepth) *maxdepth = *depth;
1886     break;
1887     case e_func:
1888     case e_uminus:
1889     break;
1890     case e_plus:
1891     case e_minus:
1892     case e_times:
1893     case e_divide:
1894     case e_power:
1895     case e_ipower:
1896     (*depth)--;
1897     break;
1898     default:
1899 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
1900 aw0a 1 }
1901     }
1902     }
1903    
1904 johnpye 726 unsigned long RelationDepth(CONST struct relation *rel){
1905 aw0a 1 unsigned long depth=0,maxdepth=0;
1906     switch(RelationRelop(rel)){
1907     case e_equal:
1908     case e_notequal:
1909     case e_less:
1910     case e_greater:
1911     case e_lesseq:
1912     case e_greatereq:
1913     CalcDepth(rel,1,&depth,&maxdepth);
1914     CalcDepth(rel,0,&depth,&maxdepth);
1915     assert(depth == 2);
1916     break;
1917     case e_maximize:
1918     case e_minimize:
1919     CalcDepth(rel,1,&depth,&maxdepth);
1920     assert(depth == 1);
1921     break;
1922     default:
1923 johnpye 726 Asc_Panic(2, __FUNCTION__,"Unknown relation term type");
1924 aw0a 1 }
1925     return maxdepth;
1926     }
1927    
1928 johnpye 669 /*------------------------------------------------------------------------------
1929     EQUATION SCALING
1930 aw0a 1
1931 johnpye 669 "Documentation will be added at a later date" -- someone last century
1932     */
1933    
1934 johnpye 726 static double FindMaxAdditiveTerm(struct relation_term *s){
1935 aw0a 1 enum safe_err serr;
1936 ben.allan 118 double lhs, rhs;
1937 johnpye 85
1938 aw0a 1 switch (RelationTermType(s)) {
1939     case e_plus:
1940     case e_minus:
1941 ben.allan 118 /** note these used to be inlined with max, but a bug in gcc323 caused it to be split out. */
1942     lhs = FindMaxAdditiveTerm(TermBinLeft(s));
1943     rhs = FindMaxAdditiveTerm(TermBinRight(s));
1944     return MAX(fabs(lhs), fabs(rhs));
1945 aw0a 1 case e_uminus:
1946     return (FindMaxAdditiveTerm(TermUniLeft(s)));
1947     case e_times:
1948     return (FindMaxAdditiveTerm(TermBinLeft(s))*
1949     FindMaxAdditiveTerm(TermBinRight(s)));
1950     case e_divide:
1951     /* bug patch / 0 */
1952     return safe_div_D0(FindMaxAdditiveTerm(TermBinLeft(s)) ,
1953     RelationBranchEvaluator(TermBinRight(s)),&serr);
1954     default:
1955     return RelationBranchEvaluator(s);
1956     }
1957     }
1958    
1959 johnpye 726 static double FindMaxFromTop(struct relation *s){
1960 ben.allan 118 double lhs;
1961     double rhs;
1962 aw0a 1 if (s == NULL) {
1963     return 0;
1964     }
1965 ben.allan 118 /** note these used to be inlined with max, but a bug in gcc323 caused it to be split out. */
1966     lhs = FindMaxAdditiveTerm(Infix_LhsSide(s));
1967     rhs = FindMaxAdditiveTerm(Infix_RhsSide(s));
1968     return MAX(fabs(lhs), fabs(rhs));
1969 aw0a 1 }
1970    
1971 johnpye 726 /* send in relation */
1972     double CalcRelationNominal(struct Instance *i){
1973 aw0a 1 enum Expr_enum reltype;
1974 johnpye 815 static int msg=0;
1975 johnpye 85
1976     char *iname;
1977     iname = WriteInstanceNameString(i,NULL);
1978     ascfree(iname);
1979    
1980 jds 97 glob_rel = NULL;
1981 aw0a 1 if (i == NULL){
1982     FPRINTF(ASCERR, "error in CalcRelationNominal routine\n");
1983     return (double)0;
1984     }
1985     if (InstanceKind(i) != REL_INST) {
1986     FPRINTF(ASCERR, "error in CalcRelationNominal routine\n");
1987     return (double)0;
1988     }
1989     glob_rel = (struct relation *)GetInstanceRelation(i,&reltype);
1990     if (glob_rel == NULL) {
1991     FPRINTF(ASCERR, "error in CalcRelationNominal routine\n");
1992     return (double)0;
1993     }
1994    
1995     if (reltype == e_token) {
1996     double temp;
1997     temp = FindMaxFromTop(glob_rel);
1998     if (isnan(temp) || !finite(temp)) {
1999     glob_rel = NULL;
2000     return (double)1;
2001     }
2002     if ( temp > 0) { /* this could return some really small numbers */
2003     glob_rel = NULL;
2004     return temp;
2005     }
2006     }
2007 johnpye 694 if (reltype == e_blackbox){
2008 johnpye 908 #ifdef BBOXWHINE
2009     ERROR_REPORTER_HERE(ASC_PROG_ERR,"blackbox not implemented yet (assuming 1.0) (%s)",__FUNCTION__); /* FIXME nominal blackbox */
2010     #endif
2011     /* should pull the nominal from the assigned lhs variable. */
2012     }
2013     if (reltype == e_glassbox){
2014 johnpye 815 if(!msg){
2015     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"glassbox not implemented yet (assuming 1.0) (%s)",__FUNCTION__);
2016     msg=2;
2017     }
2018 johnpye 908 }
2019     if (reltype == e_opcode){
2020 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
2021 ben.allan 444 }
2022 aw0a 1 glob_rel = NULL;
2023     return (double)1;
2024     }
2025    
2026 johnpye 726 void PrintScale(struct Instance *i){
2027 aw0a 1 if (InstanceKind(i) == REL_INST) {
2028     double j;
2029     j = CalcRelationNominal(i);
2030     PRINTF(" scale constant = %g\n", j);
2031     }
2032     }
2033    
2034 johnpye 726 void PrintRelationNominals(struct Instance *i){
2035 aw0a 1 VisitInstanceTree(i,PrintScale, 0, 0);
2036     }
2037    
2038 johnpye 669 /*------------------------------------------------------------------------------
2039     CALCULATION ROUTINES
2040     */
2041    
2042 aw0a 1 /**
2043     * Load ATOM values into an array of doubles.
2044     * The array of doubles is indexed from 0 while the
2045     * var list is indexed from 1. The ultimate client of
2046     * the array calling this function thinks vars index from 0.
2047     */
2048 johnpye 669 static
2049 johnpye 726 void RelationLoadDoubles(struct gl_list_t *varlist, double *vars){
2050 aw0a 1 unsigned long c;
2051     vars--; /* back up the pointer so indexing from 1 puts data right */
2052     for (c= gl_length(varlist); c > 0; c--) {
2053     vars[c] = RealAtomValue((struct Instance *)gl_fetch(varlist,c));
2054     }
2055     }
2056    
2057 johnpye 669 /**
2058     only called on token relations.
2059     */
2060 johnpye 726 int RelationCalcResidualBinary(CONST struct relation *r, double *res){
2061 aw0a 1 double *vars;
2062     double tres;
2063     int old_errno;
2064    
2065     if (r == NULL || res == NULL) {
2066     return 1;
2067     }
2068     vars = tmpalloc_array(gl_length(r->vars),double);
2069     if (vars == NULL) {
2070     return 1;
2071     }
2072     RelationLoadDoubles(r->vars,vars);
2073     old_errno = errno; /* push C global errno */
2074     errno = 0;
2075     if (BinTokenCalcResidual(RTOKEN(r).btable,RTOKEN(r).bindex,vars,&tres)) {
2076     if (errno == 0) { /* pop if unchanged */
2077     errno = old_errno;
2078     }
2079     return 1;
2080     }
2081     if (!finite(tres) || errno == EDOM || errno == ERANGE ) {
2082     if (errno == 0) { /* pop if unchanged */
2083     errno = old_errno;
2084     }
2085     return 1;
2086     }
2087     if (errno == 0) { /* pop if unchanged */
2088     errno = old_errno;
2089     }
2090     *res = tres;
2091     return 0;
2092     }
2093    
2094     enum safe_err
2095 johnpye 669 RelationCalcResidualPostfixSafe(struct Instance *i, double *res){
2096     struct relation *r;
2097     enum Expr_enum reltype;
2098     enum safe_err status = safe_ok;
2099     unsigned long length_lhs, length_rhs;
2100 aw0a 1
2101 johnpye 669 /* CONSOLE_DEBUG("..."); */
2102    
2103 johnpye 694 CHECK_INST_RES(i,res,1);
2104    
2105 johnpye 669 r = (struct relation *)GetInstanceRelation(i, &reltype);
2106 aw0a 1
2107 johnpye 669 /* CONSOLE_DEBUG("..."); */
2108    
2109     if( r == NULL ) {
2110     ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation");
2111     return safe_problem;
2112     }
2113    
2114     /* CONSOLE_DEBUG("..."); */
2115    
2116     switch(reltype){
2117     case e_token:
2118     length_lhs = RelationLength(r, 1);
2119     length_rhs = RelationLength(r, 0);
2120    
2121     if(length_lhs > 0){
2122     length_lhs--;
2123     *res = RelationEvaluatePostfixBranchSafe(r, &length_lhs, 1,&status);
2124     /* CONSOLE_DEBUG("res = %g",res); */
2125     }else{
2126     *res = 0.0;
2127     }
2128    
2129     if(length_rhs > 0){
2130     length_rhs--;
2131     *res -= RelationEvaluatePostfixBranchSafe(r, &length_rhs, 0,&status);
2132     /* CONSOLE_DEBUG("res = %g",res); */
2133     }
2134    
2135     safe_error_to_stderr(&status);
2136     break;
2137     case e_blackbox:
2138 johnpye 910 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Attempting to evaluate blackbox..."); */
2139 johnpye 669 if ( RelationCalcResidualPostfix(i,res) != 0) {
2140     status = safe_problem;
2141     safe_error_to_stderr(&status);
2142     }
2143     break;
2144     case e_glassbox:
2145     ERROR_REPORTER_HERE(ASC_PROG_ERR,"'e_glassbox' relation not supported");
2146     status = safe_problem;
2147     break;
2148     case e_opcode:
2149     ERROR_REPORTER_HERE(ASC_PROG_ERR,"'e_opcode' relation not supported");
2150     status = safe_problem;
2151     break;
2152     default:
2153     if(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2154     status = safe_problem;
2155     }else{
2156 johnpye 726 Asc_Panic(2, __FUNCTION__, "reached end of routine!");
2157 johnpye 669 }
2158     }
2159     return status;
2160 aw0a 1 }
2161    
2162    
2163     int
2164 johnpye 726 RelationCalcResidualPostfix(struct Instance *i, double *res){
2165 aw0a 1 struct relation *r;
2166     enum Expr_enum reltype;
2167 johnpye 697 unsigned long length_lhs, length_rhs;
2168 aw0a 1
2169 johnpye 694 CHECK_INST_RES(i,res,1);
2170    
2171 aw0a 1 r = (struct relation *)GetInstanceRelation(i, &reltype);
2172 johnpye 709 if(r == NULL){
2173 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2174 aw0a 1 return 1;
2175     }
2176 johnpye 709
2177     /*
2178     struct Instance *p;
2179     p = InstanceParent(i,1);
2180     char *tmp;
2181     tmp = WriteRelationString(i,p,NULL,NULL,relio_ascend,NULL);
2182     CONSOLE_DEBUG("Evaluating residual for '%s'",tmp);
2183     ASC_FREE(tmp);
2184     */
2185    
2186 aw0a 1 if( reltype == e_token ) {
2187     length_lhs = RelationLength(r, 1);
2188     length_rhs = RelationLength(r, 0);
2189     if( length_lhs > 0 ) {
2190     length_lhs--;
2191     *res = RelationEvaluatePostfixBranch(r, &length_lhs, 1);
2192     }
2193     else {
2194     *res = 0.0;
2195     }
2196     if( length_rhs > 0 ) {
2197     length_rhs--;
2198     *res -= RelationEvaluatePostfixBranch(r, &length_rhs, 0);
2199     }
2200     return 0;
2201 johnpye 709 }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2202 johnpye 694
2203     if (reltype == e_blackbox){
2204 johnpye 908 return BlackBoxCalcResidual(i, res, r);
2205     }
2206     if (reltype == e_glassbox){
2207 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"glassbox not implemented yet (%s)",__FUNCTION__);
2208 johnpye 908 }
2209     if (reltype == e_opcode) {
2210 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"opcode not supported (%s)",__FUNCTION__);
2211     }
2212 ben.allan 624
2213 aw0a 1 return 1;
2214 johnpye 709 }else{
2215 johnpye 694 Asc_Panic(2, __FUNCTION__,"reached end of routine");
2216 aw0a 1 }
2217     }
2218    
2219 johnpye 726 int RelationCalcExceptionsInfix(struct Instance *i){
2220 aw0a 1 enum Expr_enum reltype;
2221     double res;
2222     int result = 0;
2223     int old_errno;
2224    
2225     glob_rel = NULL;
2226 johnpye 694
2227     CHECK_INST_RES(i,&res,-1);
2228    
2229 aw0a 1 glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);
2230     if( glob_rel == NULL ) {
2231 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation");
2232 aw0a 1 return -1;
2233     }
2234     if( reltype == e_token ) {
2235     if (Infix_LhsSide(glob_rel) != NULL) {
2236     old_errno = errno;
2237     errno = 0; /* save the last errno, because we don't know why */
2238     res = RelationBranchEvaluator(Infix_LhsSide(glob_rel));
2239     if (!finite(res) || errno == EDOM || errno == ERANGE) {
2240     result |= RCE_ERR_LHS;
2241     if (isnan(res)) {
2242     result |= RCE_ERR_LHSNAN;
2243 johnpye 726 }else{
2244 aw0a 1 if (!finite(res)) {
2245     result |= RCE_ERR_LHSINF;
2246     }
2247     }
2248     }
2249     if (errno == 0) {
2250     errno = old_errno;
2251     } /* else something odd happened in evaluation */
2252     }
2253     if(Infix_RhsSide(glob_rel) != NULL) {
2254     res = RelationBranchEvaluator(Infix_RhsSide(glob_rel));
2255     if (!finite(res)) {
2256     result |= RCE_ERR_RHS;
2257     if (isnan(res)) {
2258     result |= RCE_ERR_RHSNAN;
2259 johnpye 726 }else{
2260 aw0a 1 if (!finite(res)) {
2261     result |= RCE_ERR_LHSINF;
2262     }
2263     }
2264     }
2265     }
2266     glob_rel = NULL;
2267     return result;
2268 johnpye 694 }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2269     ERROR_REPORTER_HERE(ASC_PROG_ERR,"relation type not implemented (%s)",__FUNCTION__);
2270 aw0a 1 glob_rel = NULL;
2271     return -1;
2272 johnpye 694 }else{
2273 johnpye 726 Asc_Panic(2, __FUNCTION__,"reached end of routine");
2274 aw0a 1 }
2275     }
2276    
2277 johnpye 694
2278 johnpye 726 int RelationCalcResidualInfix(struct Instance *i, double *res){
2279 aw0a 1 enum Expr_enum reltype;
2280     glob_rel = NULL;
2281    
2282 johnpye 694 CHECK_INST_RES(i,res,1);
2283    
2284 aw0a 1 glob_rel = (struct relation *)GetInstanceRelation(i, &reltype);
2285     if( glob_rel == NULL ) {
2286 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL relation\n");
2287 aw0a 1 return 1;
2288     }
2289     if( reltype == e_token ) {
2290     if(Infix_LhsSide(glob_rel) != NULL) {
2291     *res = RelationBranchEvaluator(Infix_LhsSide(glob_rel));
2292 johnpye 726 }else{
2293 aw0a 1 *res = 0.0;
2294     }
2295     if(Infix_RhsSide(glob_rel) != NULL) {
2296     *res -= RelationBranchEvaluator(Infix_RhsSide(glob_rel));
2297     }
2298     glob_rel = NULL;
2299     return 0;
2300 johnpye 694 }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH) {
2301     ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
2302 aw0a 1 glob_rel = NULL;
2303     return 1;
2304 johnpye 694 }else{
2305 johnpye 726 Asc_Panic(2, __FUNCTION__,"reached end of routine");
2306 aw0a 1 }
2307     }
2308    
2309    
2310 johnpye 669 /*
2311     There used to be a stoopid comment here so I removed it.
2312     */
2313 aw0a 1 int
2314 johnpye 726 RelationCalcResidualPostfix2(struct Instance *i, double *res){
2315 aw0a 1 struct relation *r;
2316     enum Expr_enum reltype;
2317    
2318 johnpye 694 CHECK_INST_RES(i,res,1);
2319    
2320 aw0a 1 r = (struct relation *)GetInstanceRelation(i, &reltype);
2321     if( r == NULL ) {
2322 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2323 aw0a 1 return 1;
2324     }
2325    
2326 johnpye 694 if( reltype == e_token ){
2327 aw0a 1 *res = RelationEvaluateResidualPostfix(r);
2328     return 0;
2329 johnpye 694 }else if (reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2330     ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
2331 aw0a 1 return 1;
2332 johnpye 694 }else{
2333 johnpye 726 Asc_Panic(2, __FUNCTION__,"reached end of routine");
2334 aw0a 1 }
2335     }
2336    
2337    
2338 johnpye 669 /*
2339     simply call the version that calculates the gradient and the residual,
2340     then ignore the residual
2341     */
2342 aw0a 1 int
2343 johnpye 726 RelationCalcGradient(struct Instance *r, double *grad){
2344 aw0a 1 double residual;
2345     return RelationCalcResidGrad(r, &residual, grad);
2346     }
2347    
2348 johnpye 669 /*
2349     simply call the version that calculates the gradient and the residual,
2350     then ignore the residual
2351     */
2352 aw0a 1 enum safe_err
2353 johnpye 726 RelationCalcGradientSafe(struct Instance *r, double *grad){
2354 aw0a 1 double residual;
2355    
2356     return RelationCalcResidGradSafe(r, &residual, grad);
2357     }
2358    
2359    
2360     int
2361 johnpye 726 RelationCalcResidGrad(struct Instance *i, double *residual, double *gradient){
2362 aw0a 1 struct relation *r;
2363     enum Expr_enum reltype;
2364    
2365 johnpye 694 CHECK_INST_RES(i,residual,1);
2366     CHECK_INST_RES(i,residual,1);
2367    
2368 aw0a 1 r = (struct relation *)GetInstanceRelation(i, &reltype);
2369     if( r == NULL ) {
2370 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null relation\n");
2371 aw0a 1 return 1;
2372     }
2373    
2374 johnpye 694 if(reltype == e_token ){
2375 aw0a 1 return RelationEvaluateResidualGradient(r, residual, gradient);
2376 johnpye 694
2377     }else if(reltype >= TOK_REL_TYPE_LOW && reltype <= TOK_REL_TYPE_HIGH){
2378     ERROR_REPORTER_HERE(ASC_PROG_ERR,"reltype not implemented (%s)",__FUNCTION__);
2379 aw0a 1 return 1;
2380 johnpye 694
2381     }else{
2382     Asc_Panic(2, __FUNCTION__, "reached end of routine");
2383 aw0a 1 }
2384     }
2385    
2386     enum safe_err
2387 johnpye 726 RelationCalcResidGradSafe(struct Instance *i
2388     , double *residual, double *gradient
2389     ){
2390 aw0a 1 struct relation *r;
2391     enum Expr_enum reltype;
2392     enum safe_err not_safe = safe_ok;
2393     int dummy_int;
2394    
2395     #ifndef NDEBUG
2396     if( i == NULL ) {
2397 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null instance\n");
2398 aw0a 1 not_safe = safe_problem;
2399     return not_safe;
2400     }
2401     if( residual == NULL || gradient == NULL ) {
2402 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"null pointer\n");
2403 aw0a 1 not_safe = safe_problem;
2404     return not_safe;
2405     }
2406     if( InstanceKind(i) != REL_INST ) {
2407 johnpye 694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"not relation\n");
2408 aw0a 1 not_safe = safe_problem;
2409     return not_safe;
2410     }
2411     #endif
2412     r = (struct relation *)GetInstanceRelation(i, &reltype);
2413     if( r == NULL