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

Contents of /trunk/base/generic/compiler/rel_blackbox.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 909 - (show annotations) (download) (as text)
Thu Oct 26 12:44:41 2006 UTC (15 years, 8 months ago) by johnpye
File MIME type: text/x-csrc
File size: 12255 byte(s)
Added finite-difference evaluation of gradients in blackboxes.
Some work on J*v evaluation with IDA.
1
2 #include <math.h>
3 #include <errno.h>
4 #include <stdarg.h>
5 #include <utilities/ascConfig.h>
6 #include <utilities/ascMalloc.h>
7 #include <utilities/ascPanic.h>
8 #include <general/list.h>
9 #include <general/dstring.h>
10 #include "compiler.h"
11 #include "symtab.h"
12 #include "functype.h"
13 #include "safe.h"
14 #include "fractions.h"
15 #include "dimen.h"
16 #include "expr_types.h"
17 #include "vlist.h"
18 #include "dimen_io.h"
19 #include "instance_enum.h"
20 #include "find.h"
21 #include "atomvalue.h"
22 #include "instance_name.h"
23 #include "extfunc.h"
24 #include "relation_type.h"
25 #include "rel_blackbox.h"
26 #include "relation.h"
27 #include "relation_util.h"
28 #include "instance_io.h"
29 #include "instquery.h"
30 #include "visitinst.h"
31 #include "mathinst.h"
32 #include "extcall.h" /* for copy/destroy speciallist */
33 #include "packages.h" /* for init slv interp */
34
35 static int32 ArgsDifferent(double new, double old, double tol)
36 {
37 if (fabs(new - old) > fabs(tol)) {
38 return 1;
39 } else {
40 return 0;
41 }
42 }
43
44 /*------------------------------------------------------------------------------
45 RESIDUAL AND GRADIENT EVALUATION ROUTINES
46 */
47
48 /*
49 Note:
50 Assume a function in c/fortran that computes yhat(x),
51 aka func(inputs_x, outputs_yhat);.
52 This is fit into the ascend world by associating
53 inputs_x directly to ascend variables and defining
54 the ascend variables y in the outer model.
55 The residual of the blackbox equations in the ascend model
56 is y-yhat. The gradient is I - dyhat/dx, where dyhat/dx
57 is the reduced jacobian of the blackbox.
58 */
59 int BlackBoxCalcResidGrad(struct Instance *i, double *res, double *gradient, struct relation *r)
60 {
61 int residErr = 0;
62 int gradErr = 0;
63
64 residErr = BlackBoxCalcResidual(i, res, r);
65 gradErr = BlackBoxCalcGradient(i, gradient, r);
66 return (int)(fabs(residErr) + fabs(gradErr));
67 }
68
69 /*
70 Note:
71 Assume a function in c/fortran that computes yhat(x),
72 aka func(inputs_x, outputs_yhat);.
73 This is fit into the ascend world by associating
74 inputs_x directly to ascend variables and defining
75 the ascend variables y in the outer model.
76 The residual of the blackbox equations in the ascend model
77 is y-yhat. The gradient is I - dyhat/dx, where dyhat/dx
78 is the reduced jacobian of the blackbox.
79 */
80 int BlackBoxCalcResidual(struct Instance *i, double *res, struct relation *r)
81 {
82 /* decls */
83 unsigned long *argToVar;
84 unsigned long c;
85 unsigned long argToVarLen;
86 struct BlackBoxData *unique;
87 struct BlackBoxCache *common;
88 struct ExternalFunc *efunc;
89 ExtBBoxFunc *evalFunc;
90 int lhsVar;
91 int outputIndex;
92 struct Instance *arg;
93 double value;
94 double inputTolerance;
95 int updateNeeded;
96 int nok = 0;
97
98 /* impl setup */
99 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox evaluation is experimental (%s)",__FUNCTION__);
100 (void)i;
101 efunc = RelationBlackBoxExtFunc(r);
102 unique = RelationBlackBoxData(r);
103 common = RelationBlackBoxCache(r);
104 argToVar = unique->inputArgs;
105 argToVarLen = common->inputsLen;
106 lhsVar = unique->lhsvar;
107 outputIndex = unique->lhsindex;
108 evalFunc = GetValueFunc(efunc);
109 inputTolerance = GetValueFuncTolerance(efunc);
110 updateNeeded = 0;
111
112 /* impl:
113 check input values changed.
114 if changed, recompute output values.
115 compute y - yhat for lhsvar.
116 */
117 /* check: */
118 if (common->residCount < 1) {
119 updateNeeded = 1;
120 } else {
121 for (c=0; c < argToVarLen ; c++) {
122 arg = RelationVariable(r,argToVar[c]);
123 value = RealAtomValue(arg);
124 if (ArgsDifferent(value, common->inputs[c], inputTolerance)) {
125 updateNeeded = 1;
126 break;
127 }
128 }
129 }
130 /* recompute */
131 if (updateNeeded) {
132 for (c=0; c < argToVarLen ;c++) {
133 arg = RelationVariable(r,argToVar[c]);
134 value = RealAtomValue(arg);
135 common->inputs[c] = value;
136 }
137 common->interp.task = bb_func_eval;
138
139 nok = (*evalFunc)(&(common->interp),
140 common->inputsLen,
141 common->outputsLen,
142 common->inputs,
143 common->outputs,
144 common->jacobian);
145 common->residCount++;
146 }
147 value = common->outputs[outputIndex];
148
149 /* return lhsval - value */
150 arg = RelationVariable(r,lhsVar);
151 *res = RealAtomValue(arg) - value;
152 return nok;
153
154 }
155
156 /**
157 Calculate the gradient (slice of the overall jacobian) for the blackbox.
158
159 The tricky bit in this is that if input or output args are merged,
160 their partial derivatives get summed and the solver at least gets
161 what it deserves, which may or may not be sane in a modeling sense.
162
163 Implementation:
164 # check input values changed.
165 # if changed, recompute gradient in bbox.
166 # compute gradient per varlist from bbox row.
167 */
168
169 int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
170 , int ninputs, int noutputs
171 , double *inputs, double *outputs, double *jac
172 );
173
174 int BlackBoxCalcGradient(struct Instance *i, double *gradient
175 , struct relation *r
176 ){
177 /* decls */
178 unsigned long *argToVar;
179 unsigned long c, varlistLen;
180 unsigned long argToVarLen;
181 struct BlackBoxData *unique;
182 struct BlackBoxCache *common;
183 struct ExternalFunc *efunc;
184 ExtBBoxFunc *derivFunc;
185 int lhsVar;
186 int outputIndex;
187 struct Instance *arg;
188 double value;
189 double inputTolerance;
190 double *bboxRow;
191 int updateNeeded, offset;
192 unsigned int k;
193 int nok = 0;
194
195 /* prepare */
196 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox gradient is experimental (%s)",__FUNCTION__);
197 (void)i;
198 efunc = RelationBlackBoxExtFunc(r);
199 unique = RelationBlackBoxData(r);
200 common = RelationBlackBoxCache(r);
201 argToVar = unique->inputArgs;
202 argToVarLen = common->inputsLen;
203 lhsVar = unique->lhsvar;
204 outputIndex = unique->lhsindex;
205 derivFunc = GetDerivFunc(efunc);
206 inputTolerance = GetValueFuncTolerance(efunc);
207 updateNeeded = 0;
208 varlistLen = NumberVariables(r);
209
210 /* do we need to calculate? */
211 if(common->gradCount < 1){
212 updateNeeded = 1;
213 }else{
214 for(c=0; c<argToVarLen; c++){
215 arg = RelationVariable(r,argToVar[c]);
216 value = RealAtomValue(arg);
217 if (ArgsDifferent(value, common->inputsJac[c], inputTolerance)) {
218 updateNeeded = 1;
219 break;
220 }
221 }
222 }
223
224 /* if inputs have changed more than allowed, or it's first time: evaluate */
225 if (updateNeeded) {
226 for (c=0; c < argToVarLen ;c++) {
227 arg = RelationVariable(r,argToVar[c]);
228 value = RealAtomValue(arg);
229 common->inputsJac[c] = value;
230 }
231 common->interp.task = bb_deriv_eval;
232
233 if(derivFunc){
234 nok = (*derivFunc)(
235 &(common->interp)
236 , common->inputsLen, common->outputsLen
237 , common->inputsJac, common->outputs, common->jacobian
238 );
239 }else{
240 CONSOLE_DEBUG("COMPUTING FINITE-DIFFERENCE BLACK-BOX DERIVATIVE");
241
242 nok = blackbox_fdiff(GetValueFunc(efunc), &(common->interp)
243 , common->inputsLen, common->outputsLen
244 , common->inputsJac, common->outputs, common->jacobian
245 );
246 }
247 common->gradCount++;
248 }
249
250 for (k = 0; k < varlistLen; k++) {
251 gradient[k] = 0.0;
252 }
253
254 /* now compute d(y-yhat)/dx for this row as ( I - dyhat/dx ) */
255 k = lhsVar-1;
256 gradient[k] = 1.0; /* I */
257 offset = common->inputsLen * outputIndex; /* offset to row needed. */
258 bboxRow = common->jacobian + offset; /* pointer to row by ptr arithmetic */
259 for (c=0; c < argToVarLen; c++) {
260 k = ((int)(argToVar[c])) - 1; /* varnum is 1..n, not 0 .*/
261 gradient[k] -= bboxRow[c];
262 }
263
264 return nok;
265 }
266
267 struct Instance *BlackBoxGetOutputVar(struct relation *r)
268 {
269 assert(r != NULL);
270 /* unsigned long lhsVarNumber; */
271 /** @TODO FIXME BlackBoxGetOutputVar */
272 return NULL;
273 }
274
275 struct BlackBoxData *CreateBlackBoxData(struct BlackBoxCache *common,
276 unsigned long lhsindex,
277 unsigned long lhsVarNumber)
278 {
279 struct BlackBoxData *b = (struct BlackBoxData *)malloc(sizeof(struct BlackBoxData));
280 assert(common!=NULL);
281 b->common = common;
282 AddRefBlackBoxCache(common);
283 b->inputArgs = (unsigned long *)ascmalloc(sizeof(long) * (common->inputsLen));
284 b->lhsindex = lhsindex;
285 b->lhsvar = lhsVarNumber;
286 return b;
287 }
288
289 void DestroyBlackBoxData(struct relation *rel, struct BlackBoxData *b)
290 {
291 b->lhsindex = -(b->lhsindex);
292 b->lhsvar = 0;
293 ascfree(b->inputArgs);
294 b->inputArgs = NULL;
295 DeleteRefBlackBoxCache(rel, &(b->common));
296 ascfree(b);
297 }
298
299 /*------------------------------------------------------------------------------
300 BLACKBOX GRADIENT BY FINITE-DIFFERENCE
301 */
302
303 /**
304 This function works out what the peturbed value for a variable should be.
305
306 @TODO add some awareness of boundaries in this, hopefully.
307 */
308 static inline double blackbox_peturbation(double varvalue){
309 return (1.0e-05);
310 }
311
312 /**
313 Blackbox derivatives estimated by finite difference (by evaluation at
314 peturbed value of each input in turn)
315
316 Call signature as for ExtBBoxFunc.
317 */
318 int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
319 , int ninputs, int noutputs
320 , double *inputs, double *outputs, double *jac
321 ){
322 long c,r;
323 int nok = 0;
324 double *tmp_outputs;
325 double *ptr;
326 double old_x,deltax,value;
327
328 CONSOLE_DEBUG("NUMERICAL DERIVATIVE...");
329
330 tmp_outputs = ASC_NEW_ARRAY_CLEAR(double,noutputs);
331
332 for (c=0;c<ninputs;c++){
333 /* perturb each input in turn */
334 old_x = inputs[c];
335 deltax = blackbox_peturbation(old_x);
336 inputs[c] = old_x + deltax;
337 CONSOLE_DEBUG("PETURBATED VALUE of input[%ld] = %f",c,inputs[c]);
338
339 /* call routine. note that the 'jac' parameter is just along for the ride */
340 nok = (*resfn)(interp, ninputs, noutputs, inputs, tmp_outputs, jac);
341 if(nok){
342 CONSOLE_DEBUG("External evaluation error (%d)",nok);
343 break;
344 }
345
346 /* fill load jacobian */
347 ptr = &jac[c];
348 for(r=0;r<noutputs;r++){
349 value = (tmp_outputs[r] - outputs[r]) / deltax;
350 CONSOLE_DEBUG("output[%ld]: value = %f, gradient = %f",r,tmp_outputs[r],value);
351 *ptr = value;
352 ptr += ninputs;
353 }
354
355 /* now return this input to its old value */
356 inputs[c] = old_x;
357 }
358 ASC_FREE(tmp_outputs);
359 if(nok){
360 CONSOLE_DEBUG("External evaluation error");
361 }
362 return nok;
363
364 }
365
366 /*------------------------------------------------------------------------------
367 BLACK BOX CACHE
368 */
369
370 #define JACMAGIC -3.141592071828
371 struct BlackBoxCache *CreateBlackBoxCache(
372 int32 inputsLen,
373 int32 outputsLen,
374 struct gl_list_t *formalArgs
375 )
376 {
377 struct BlackBoxCache *b = (struct BlackBoxCache *)malloc(sizeof(struct BlackBoxCache));
378 b->interp.task = bb_none;
379 b->interp.status = calc_all_ok;
380 b->interp.user_data = NULL;
381 b->formalArgList = CopySpecialList(formalArgs);
382 b->inputsLen = inputsLen;
383 b->outputsLen = outputsLen;
384 b->jacobianLen = inputsLen*outputsLen;
385 b->hessianLen = 0; /* FIXME. */
386 b->inputs = (double *)ascmalloc(inputsLen*sizeof(double));
387 b->inputsJac = (double *)ascmalloc(inputsLen*sizeof(double));
388 b->outputs = (double *)ascmalloc(outputsLen*sizeof(double));
389 b->jacobian = (double*)ascmalloc(outputsLen*inputsLen*sizeof(double)+sizeof(double));
390 b->jacobian[outputsLen*inputsLen] = JACMAGIC;
391 b->hessian = NULL;
392 b->residCount = 0;
393 b->gradCount = 0;
394 b->refCount = 1;
395 return b;
396 }
397
398 void AddRefBlackBoxCache(struct BlackBoxCache *b)
399 {
400 assert(b != NULL);
401 (b->refCount)++;
402 }
403
404 static void DestroyBlackBoxCache(struct relation *rel, struct BlackBoxCache *b)
405 {
406 struct ExternalFunc *efunc;
407 ExtBBoxFinalFunc *final;
408 assert(b != NULL);
409 if (b->jacobian[b->outputsLen*b->inputsLen] != JACMAGIC)
410 {
411 FPRINTF(ASCERR, "Jacobian overrun detected while destroying BlackBoxCache. Debug the user blackbox implementation.\n");
412 }
413 ascfree(b->inputs);
414 ascfree(b->inputsJac);
415 ascfree(b->outputs);
416 ascfree(b->jacobian);
417 DestroySpecialList(b->formalArgList);
418 b->formalArgList = NULL;
419 b->inputsLen = -(b->inputsLen);
420 b->outputsLen = -(b->outputsLen);
421 b->jacobianLen = -(b->jacobianLen);
422 b->hessianLen = -(b->hessianLen);
423 b->inputs = NULL;
424 b->inputsJac = NULL;
425 b->outputs = NULL;
426 b->jacobian = NULL;
427 b->hessian = NULL;
428 b->residCount = -(b->residCount);
429 b->gradCount = -(b->gradCount);
430 b->refCount = -3;
431 if (rel != NULL) {
432 /* rel is null with refcount 0 only when there's a bbox array
433 instantiation failure. If failure, we don't need final. */
434 efunc = RelationBlackBoxExtFunc(rel);
435 final = GetFinalFunc(efunc);
436 b->interp.task = bb_last_call;
437 (*final)(&(b->interp));
438 }
439 ascfree(b);
440 }
441
442 void DeleteRefBlackBoxCache(struct relation *rel, struct BlackBoxCache **b)
443 {
444 struct BlackBoxCache * d = *b;
445 assert(b != NULL && *b != NULL);
446 *b = NULL;
447 if (d->refCount > 0) {
448 (d->refCount)--;
449 }
450 if (d->refCount == 0) {
451 DestroyBlackBoxCache(rel,d);
452 return;
453 }
454 if (d->refCount < 0) {
455 FPRINTF(ASCERR, "Attempt to delete BlackBoxCache (%p) too often.\n", *b);
456 }
457 }

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