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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 16242 byte(s)
Setting up web subdirectory in repository
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