/[ascend]/trunk/ascend4/solver/extrel.c
ViewVC logotype

Annotation of /trunk/ascend4/solver/extrel.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 16242 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 /*
2     * External Relations Cache for solvers.
3     * by Kirk Andre Abbott
4     * Created: 08/10/94
5     * Version: $Revision: 1.10 $
6     * Version control file: $RCSfile: extrel.c,v $
7     * Date last modified: $Date: 1997/07/18 12:14:14 $
8     * Last modified by: $Author: mthomas $
9     *
10     * This file is part of the SLV solver.
11     *
12     * Copyright (C) 1994 Kirk Abbott
13     *
14     * The SLV solver is free software; you can redistribute
15     * it and/or modify it under the terms of the GNU General Public License as
16     * published by the Free Software Foundation; either version 2 of the
17     * License, or (at your option) any later version.
18     *
19     * The SLV solver is distributed in hope that it will be
20     * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22     * General Public License for more details.
23     *
24     * You should have received a copy of the GNU General Public License along with
25     * the program; if not, write to the Free Software Foundation, Inc., 675
26     * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
27     * COPYING is found in ../compiler.
28     */
29     #include "utilities/ascConfig.h"
30     #include "compiler/packages.h"
31     /* #include "solver/exprman.h" */
32     #include "solver/rel.h"
33     #include "solver/extrel.h"
34     #include "compiler/extcall.h"
35    
36     double g_external_tolerance = 1.0e-12;
37    
38     struct ExtRelCache *CreateExtRelCache(struct ExtCallNode *ext)
39     {
40     struct ExtRelCache *result;
41     unsigned long n_input_args, n_output_args;
42     int ninputs, noutputs;
43    
44     assert(ext!=NULL);
45     result = (struct ExtRelCache *)ascmalloc(sizeof(struct ExtRelCache));
46     result->nodestamp = ExternalCallNodeStamp(ext);
47     result->efunc = ExternalCallExtFunc(ext);
48     result->data = ExternalCallDataInstance(ext);
49     result->arglist = ExternalCallArgList(ext);
50     result->user_data = NULL; /* absolutely important to be
51     * initialized to NULL ! */
52    
53     n_input_args = NumberInputArgs(result->efunc);
54     n_output_args = NumberOutputArgs(result->efunc);
55    
56     /*
57     * We own the result of the LinearizeArgList call.
58     */
59     result->inputlist = LinearizeArgList(result->arglist,1,n_input_args);
60     ninputs = (int)gl_length(result->inputlist);
61     noutputs = (int)CountNumberOfArgs(result->arglist,n_input_args+1,
62     n_input_args+n_output_args);
63     result->ninputs = ninputs;
64     result->noutputs = noutputs;
65    
66     result->inputs = (double *)asccalloc(ninputs,sizeof(double));
67     result->outputs = (double *)asccalloc(noutputs,sizeof(double));
68     result->jacobian = (double *)asccalloc(ninputs*noutputs,sizeof(double));
69     /*
70     * Setup default flags for controlling calculations.
71     */
72     result->newcalc_done = (unsigned)1;
73     return result;
74     }
75    
76    
77     struct ExtRelCache *CreateCacheFromInstance(struct Instance *relinst)
78     {
79     struct ExtCallNode *ext;
80     struct ExtRelCache *cache;
81     CONST struct relation *reln;
82     enum Expr_enum type;
83    
84     reln = GetInstanceRelation(relinst,&type);
85     if (type!=e_blackbox) {
86     FPRINTF(stderr,"Invalid relation type in (CreateCacheFromInstance)\n");
87     return NULL;
88     }
89     ext = BlackBoxExtCall(reln);
90     cache = CreateExtRelCache(ext);
91     return cache;
92     }
93    
94     void ExtRel_DestroyCache(struct ExtRelCache *cache)
95     {
96     if (cache) {
97     cache->nodestamp=-1;
98     cache->efunc = NULL;
99     cache->data = NULL;
100     gl_destroy(cache->inputlist); /* we own the list */
101     ascfree(cache->inputs);
102     ascfree(cache->outputs);
103     ascfree(cache->jacobian);
104     cache->inputlist = NULL;
105     cache->inputs = NULL;
106     cache->outputs = NULL;
107     cache->jacobian = NULL;
108     ascfree(cache);
109     }
110     }
111    
112    
113     /*
114     *********************************************************************
115     * ExtRel_PreSolve:
116     *
117     * To deal with the first time we also want to do arguement
118     * checking, and then turn off the first_func_eval flag.
119     * Turn on the newcalc_done flag. The rationale behind this is
120     * as follows:
121     * The solver at the moment does not treat an external relation
122     * specially, i.e., as a block. It also calls for its functions
123     * a relation at a time. However the external relations compute
124     * all their outputs at once. So as not to do unnecessary
125     * recalculations we use these flag bits. We set newcalc_done
126     * initially to true, so as to force *at least* one calculation
127     * of the external relations. By similar reasoning first_func_eval (done)
128     * is set to false.
129     *********************************************************************
130     */
131     int ExtRel_PreSolve(struct ExtRelCache *cache, int setup)
132     {
133     struct ExternalFunc *efunc;
134     struct Slv_Interp slv_interp;
135     int ninputs,noutputs;
136     double *inputs, *outputs, *jacobian;
137     int (*init_func)(struct Slv_Interp *,
138     struct Instance *,struct gl_list_t *);
139     int nok = 0;
140    
141     if (!cache) return 1;
142     efunc = cache->efunc;
143     init_func = GetInitFunc(efunc);
144     Init_Slv_Interp(&slv_interp);
145     slv_interp.nodestamp = cache->nodestamp;
146     slv_interp.user_data = cache->user_data;
147     if (setup) {
148     slv_interp.first_call = (unsigned)1;
149     slv_interp.last_call = (unsigned)0;
150     slv_interp.check_args = (unsigned)1;
151     }
152     else{
153     slv_interp.first_call = (unsigned)0;
154     slv_interp.last_call = (unsigned)1;
155     slv_interp.check_args = (unsigned)0;
156     }
157     nok = (*init_func)(&slv_interp,cache->data,cache->arglist);
158     if (nok)
159     return 1;
160    
161     /*
162     * Save the user's data and update our status.
163     */
164     cache->user_data = slv_interp.user_data;
165     cache->newcalc_done = (unsigned)1; /* force at least one calculation */
166     cache->first_func_eval = (unsigned)0;
167     return 0;
168     }
169    
170    
171     static int ArgsDifferent(double new, double old)
172     {
173     if (fabs(new - old) >= g_external_tolerance)
174     return 1;
175     else
176     return 0;
177     }
178    
179     real64 ExtRel_Evaluate_RHS(struct rel_relation *rel)
180     {
181     struct Slv_Interp slv_interp;
182     struct ExtRelCache *cache;
183     struct ExternalFunc *efunc;
184     struct Instance *arg;
185     struct gl_list_t *inputlist;
186     double value;
187     int c,ninputs;
188     int nok;
189     unsigned long whichvar;
190     int newcalc_reqd=0;
191    
192     int (*eval_func)(struct Slv_Interp *,
193     int ninputs, int noutputs,
194     double *inputs, double *outputs, double *jacobian);
195    
196     assert(rel_extnodeinfo(rel));
197     cache = rel_extcache(rel);
198     efunc = cache->efunc;
199     eval_func = GetValueFunc(efunc);
200     inputlist = cache->inputlist;
201     ninputs = cache->ninputs;
202     whichvar = (unsigned long)rel_extwhichvar(rel);
203    
204     /*
205     * The determination of whether a new calculation is required needs
206     * some more thought. One thing we should insist upon is that all
207     * the relations for an external relation are forced into the same
208     * block.
209     */
210     for (c=0;c<ninputs;c++) {
211     arg = (struct Instance *)gl_fetch(inputlist,c+1);
212     value = RealAtomValue(arg);
213     if (ArgsDifferent(value,cache->inputs[c])) {
214     newcalc_reqd = 1;
215     cache->inputs[c] = value;
216     }
217     }
218     value = 0;
219     /*
220     * Do the calculations if necessary. Note that we have to *ensure*
221     * that we send the user the information that he provided to us.
222     * We have to update our user_data after each call of the user's function
223     * as he might move information around (not smart but possible), on us.
224     * If a function call is made, mark a new calculation as having been,
225     * done, otherwise dont.
226     */
227     if (newcalc_reqd) {
228     Init_Slv_Interp(&slv_interp);
229     slv_interp.nodestamp = cache->nodestamp;
230     slv_interp.user_data = cache->user_data;
231     slv_interp.func_eval = (unsigned)1;
232    
233     nok = (*eval_func)(&slv_interp, ninputs, cache->noutputs,
234     cache->inputs, cache->outputs, cache->jacobian);
235     value = cache->outputs[whichvar - ninputs - 1];
236     cache->newcalc_done = (unsigned)1; /* newcalc done */
237     cache->user_data = slv_interp.user_data; /* update user_data */
238     }
239     else{
240     value = cache->outputs[whichvar - ninputs - 1];
241     cache->newcalc_done = (unsigned)0; /* a result was simply returned */
242     }
243    
244     #ifdef KAA_DEBUG
245     FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_RHS result ->%g\n",
246     result);
247     #endif
248    
249     return value;
250     }
251    
252     /******* FIX FIX FIX **********/
253     real64 ExtRel_Evaluate_LHS(rel)
254     struct rel_relation *rel;
255     {
256     real64 res;
257     res = 1.0; /******* FIX FIX FIX **********/
258     FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_LHS result ->%g\n",
259     res);
260     return res;
261     }
262    
263    
264     /*
265     *********************************************************************
266     * ExtRel_Gradient evaluation routines.
267     *
268     * The following code implements gradient evaluation routines for
269     * external relations. The routines here will prepare the arguements
270     * and call a user supplied derivative routine, is same is non-NULL.
271     * If it is the user supplied function evaluation routine will be
272     * used to compute the gradients via finite differencing.
273     * The current solver will not necessarily call for the derivative
274     * all at once. This makes it necessary to do the gradient
275     * computations (user supplied or finite difference), and to cache
276     * away the results. Based on calculation flags, the appropriate
277     * *row* of this cached jacobian will be extracted and mapped to the
278     * main solve matrix.
279     *
280     * The cached jacobian is a contiguous vector ninputs*noutputs long
281     * and is loaded row wise. Indexing starts from 0. Each row corresponds
282     * to the partial derivatives of the output variable (associated with
283     * that row, wrt to all its incident input variables.
284     *
285     * Careful attention needs to be paid to the way this jacobian is
286     * loaded/unloaded, because of the multiple indexing schemes in use.
287     * i.e, arglist's and inputlists index 1..nvars, whereas all numeric
288     * vectors number from 0.
289     *
290     *********************************************************************
291     */
292    
293     struct deriv_data {
294     var_filter_t *filter;
295     mtx_matrix_t mtx;
296     mtx_coord_t nz;
297     };
298    
299    
300     /*
301     *********************************************************************
302     * ExtRel_MapDataToMtx
303     *
304     * This function attempts to map the information from the contiguous
305     * jacobian back into the main matrix, based on the current row <=>
306     * whichvar. The jacobian assumed to have been calculated.
307     * Since we are operating a relation at a time, we have to find
308     * out where to index into our jacobian. This index is computed as
309     * follows:
310     *
311     * index = (whichvar - ninputs - 1) * ninputs
312     *
313     * Example: a problem with 4 inputs, 3 outputs and whichvar = 6.
314     * with the counting for vars 1..nvars, but the jacobian indexing
315     * starting from 0 (c-wise).
316     *
317     * V-------- first output variable
318     * 1 2 3 4 5 6 7
319     * ^---------------- whichvar
320     * ------------------- grads for whichvar = 6
321     * | | | |
322     * v v v v
323     * index = 0 1 2 3 4 5 6 7 8 9 10 11
324     * jacobian = 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
325     *
326     * Hence jacobian index = (6 - 4 - 1) * 4 = 4
327     *********************************************************************
328     */
329    
330     static void ExtRel_MapDataToMtx(struct gl_list_t *inputlist,
331     unsigned long whichvar,
332     int ninputs,
333     double *jacobian,
334     struct deriv_data *d)
335     {
336     struct Instance *inst;
337     struct var_variable *var;
338     double value, *ptr;
339     boolean used;
340     unsigned long c;
341     int index;
342    
343     index = ((int)whichvar - ninputs - 1) * ninputs;
344     ptr = &jacobian[index];
345    
346     /* this is totally broken, thanks to kirk making the var=instance assumption */
347     Asc_Panic(2, "ExtRel_MapDataToMtx",
348     "ExtRel_MapDataToMtx is totally broken\n"
349     "Makes a var=instance assumption");
350     for (c=0;c<ninputs;c++) {
351     inst = (struct Instance *)gl_fetch(inputlist,c+1);
352     var = var_instance(inst);
353     used = var_apply_filter(var,d->filter);
354     if (used) {
355     d->nz.col = mtx_org_to_col(d->mtx,var_sindex(var));
356     value = ptr[c] + mtx_value(d->mtx,&(d->nz));
357     mtx_set_value(d->mtx,&(d->nz), value);
358     }
359     }
360     }
361    
362    
363     /*
364     *********************************************************************
365     * ExtRel Finite Differencing.
366     *
367     * This routine actually does the finite differencing.
368     * The jacobian is a single contiguous vector. We load information
369     * in it *row* wise. If we have noutputs x ninputs = 3 x 4, variables,
370     * then jacobian entry 4,
371     * would correspond to jacobian[1][0], i.e., = 0.5 for this eg.
372     *
373     * 2.0 9.0 4.0 6.0
374     * 0.5 1.3 0.0 9.7
375     * 80 7.0 1.0 2.5
376     *
377     * 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
378     * [0][0] [1][0] [2][1]
379     *
380     * When we are finite differencing variable c, we will be loading
381     * jacobian positions c, c+ninputs, c+2*ninputs ....
382     *********************************************************************
383     */
384    
385     static double CalculateInterval(double var_value)
386     {
387     return (1.0e-05);
388     }
389    
390     static int ExtRel__FDiff(struct Slv_Interp *slv_interp,
391     int (*eval_func) (/* ARGS */),
392     int ninputs, int noutputs,
393     double *inputs, double *outputs,
394     double *jacobian)
395     {
396     int c1,c2, nok = 0;
397     double *tmp_vector;
398     double *ptr;
399     double old_x,interval,value;
400    
401     tmp_vector = (double *)asccalloc(noutputs,sizeof(double));
402     for (c1=0;c1<ninputs;c1++) {
403     old_x = inputs[c1]; /* perturb x */
404     interval = CalculateInterval(old_x);
405     inputs[c1] = old_x + interval;
406     nok = (*eval_func)(slv_interp, ninputs, noutputs, /* call routine */
407     inputs, tmp_vector, jacobian);
408     if (nok) break;
409     ptr = &jacobian[c1];
410     for (c2=0;c2<noutputs;c2++) { /* load jacobian */
411     value = (tmp_vector[c2] - outputs[c2])/interval;
412     *ptr = value;
413     ptr += ninputs;
414     }
415     inputs[c1] = old_x;
416     }
417     ascfree((char *)tmp_vector); /* cleanup */
418     return nok;
419     }
420    
421    
422     int ExtRel_CalcDeriv(struct rel_relation *rel, struct deriv_data *d)
423     {
424     int nok = 0;
425     struct Slv_Interp slv_interp;
426     struct ExtRelCache *cache;
427     struct ExternalFunc *efunc;
428     unsigned long whichvar;
429     double *deriv_vector;
430     int (*eval_func)();
431     int (*deriv_func)();
432    
433     assert(rel_extnodeinfo(rel));
434     cache = rel_extcache(rel);
435     whichvar = rel_extwhichvar(rel);
436     efunc = cache->efunc;
437    
438     /*
439     * Check and deal with the special case of the first
440     * computation.
441     */
442     if (cache->first_deriv_eval) {
443     cache->newcalc_done = (unsigned)1;
444     cache->first_deriv_eval = (unsigned)0;
445     }
446    
447     /*
448     * If a function evaluation was not recently done, then we
449     * can return the results from the cached jacobian.
450     */
451     if (!cache->newcalc_done) {
452     ExtRel_MapDataToMtx(cache->inputlist, whichvar,
453     cache->ninputs, cache->jacobian, d);
454     return 0;
455     }
456    
457     /*
458     * If we are here, we have to do the derivative calculation.
459     * The only thing to determine is whether we do analytical
460     * derivatives (deriv_func != NULL) or finite differencing.
461     * In any case init the interpreter.
462     */
463     Init_Slv_Interp(&slv_interp);
464     slv_interp.deriv_eval = (unsigned)1;
465     slv_interp.user_data = cache->user_data;
466     deriv_func = GetDerivFunc(efunc);
467     if (deriv_func) {
468     nok = (*deriv_func)(&slv_interp, cache->ninputs, cache->noutputs,
469     cache->inputs, cache->outputs, cache->jacobian);
470     if (nok) return nok;
471     }
472     else{
473     eval_func = GetValueFunc(efunc);
474     nok = ExtRel__FDiff(&slv_interp, eval_func,
475     cache->ninputs, cache->noutputs,
476     cache->inputs, cache->outputs, cache->jacobian);
477     if (nok) return nok;
478     }
479    
480     /*
481     * Cleanup. Ensure that we update the users data, and load
482     * the main matrix with the derivative information.
483     */
484     cache->user_data = slv_interp.user_data; /* save user info */
485     ExtRel_MapDataToMtx(cache->inputlist, whichvar,
486     cache->ninputs, cache->jacobian, d);
487     return 0;
488     }
489    
490    
491     /*
492     *********************************************************************
493     * ExtRel Deriv routines.
494     *
495     * This is the entry point for most cases. ExtRel_CalcDeriv depends
496     * on ExtRel_Evaluate being called immediately before it.
497     *
498     *********************************************************************
499     */
500    
501     real64 ExtRel_Diffs_RHS(struct rel_relation *rel,
502     var_filter_t *filter,
503     int32 row,
504     mtx_matrix_t mtx)
505     {
506     int nok;
507     real64 res;
508     struct deriv_data data;
509    
510     data.filter = filter;
511     data.mtx = mtx;
512     data.nz.row = row;
513    
514     res = ExtRel_Evaluate_RHS(rel);
515     nok = ExtRel_CalcDeriv(rel,&data);
516     if (nok)
517     return 0.0;
518     else
519     return res;
520     }
521    
522    
523    
524     /******* FIX FIX FIX **********/
525    
526     real64 ExtRel_Diffs_LHS(struct rel_relation *rel,
527     var_filter_t *filter,
528     int32 row,
529     mtx_matrix_t mtx)
530     {
531     real64 res=0.0;
532     return 1.0; /******* FIX FIX FIX **********/
533     }
534    
535    
536    
537    
538    
539    
540    
541    
542    
543    

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22