/[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 910 - (show annotations) (download) (as text)
Thu Oct 26 13:35:25 2006 UTC (13 years, 8 months ago) by johnpye
File MIME type: text/x-csrc
File size: 12377 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
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 static int warnexpt;
98
99 if(!warnexpt){
100 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox evaluation is experimental (%s)",__FUNCTION__);
101 warnexpt = 1;
102 }
103
104 (void)i;
105 efunc = RelationBlackBoxExtFunc(r);
106 unique = RelationBlackBoxData(r);
107 common = RelationBlackBoxCache(r);
108 argToVar = unique->inputArgs;
109 argToVarLen = common->inputsLen;
110 lhsVar = unique->lhsvar;
111 outputIndex = unique->lhsindex;
112 evalFunc = GetValueFunc(efunc);
113 inputTolerance = GetValueFuncTolerance(efunc);
114 updateNeeded = 0;
115
116 /* impl:
117 check input values changed.
118 if changed, recompute output values.
119 compute y - yhat for lhsvar.
120 */
121 /* check: */
122 if (common->residCount < 1) {
123 updateNeeded = 1;
124 } else {
125 for (c=0; c < argToVarLen ; c++) {
126 arg = RelationVariable(r,argToVar[c]);
127 value = RealAtomValue(arg);
128 if (ArgsDifferent(value, common->inputs[c], inputTolerance)) {
129 updateNeeded = 1;
130 break;
131 }
132 }
133 }
134 /* recompute */
135 if (updateNeeded) {
136 for (c=0; c < argToVarLen ;c++) {
137 arg = RelationVariable(r,argToVar[c]);
138 value = RealAtomValue(arg);
139 common->inputs[c] = value;
140 }
141 common->interp.task = bb_func_eval;
142
143 nok = (*evalFunc)(&(common->interp),
144 common->inputsLen,
145 common->outputsLen,
146 common->inputs,
147 common->outputs,
148 common->jacobian);
149 common->residCount++;
150 }
151 value = common->outputs[outputIndex];
152
153 /* return lhsval - value */
154 arg = RelationVariable(r,lhsVar);
155 *res = RealAtomValue(arg) - value;
156 return nok;
157
158 }
159
160 /**
161 Calculate the gradient (slice of the overall jacobian) for the blackbox.
162
163 The tricky bit in this is that if input or output args are merged,
164 their partial derivatives get summed and the solver at least gets
165 what it deserves, which may or may not be sane in a modeling sense.
166
167 Implementation:
168 # check input values changed.
169 # if changed, recompute gradient in bbox.
170 # compute gradient per varlist from bbox row.
171 */
172
173 int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
174 , int ninputs, int noutputs
175 , double *inputs, double *outputs, double *jac
176 );
177
178 int BlackBoxCalcGradient(struct Instance *i, double *gradient
179 , struct relation *r
180 ){
181 /* decls */
182 unsigned long *argToVar;
183 unsigned long c, varlistLen;
184 unsigned long argToVarLen;
185 struct BlackBoxData *unique;
186 struct BlackBoxCache *common;
187 struct ExternalFunc *efunc;
188 ExtBBoxFunc *derivFunc;
189 int lhsVar;
190 int outputIndex;
191 struct Instance *arg;
192 double value;
193 double inputTolerance;
194 double *bboxRow;
195 int updateNeeded, offset;
196 unsigned int k;
197 int nok = 0;
198 static int warnexpt;
199
200 if(!warnexpt){
201 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox gradient is experimental (%s)",__FUNCTION__);
202 warnexpt = 1;
203 }
204
205 /* prepare */
206 (void)i;
207 efunc = RelationBlackBoxExtFunc(r);
208 unique = RelationBlackBoxData(r);
209 common = RelationBlackBoxCache(r);
210 argToVar = unique->inputArgs;
211 argToVarLen = common->inputsLen;
212 lhsVar = unique->lhsvar;
213 outputIndex = unique->lhsindex;
214 derivFunc = GetDerivFunc(efunc);
215 inputTolerance = GetValueFuncTolerance(efunc);
216 updateNeeded = 0;
217 varlistLen = NumberVariables(r);
218
219 /* do we need to calculate? */
220 if(common->gradCount < 1){
221 updateNeeded = 1;
222 }else{
223 for(c=0; c<argToVarLen; c++){
224 arg = RelationVariable(r,argToVar[c]);
225 value = RealAtomValue(arg);
226 if (ArgsDifferent(value, common->inputsJac[c], inputTolerance)) {
227 updateNeeded = 1;
228 break;
229 }
230 }
231 }
232
233 /* if inputs have changed more than allowed, or it's first time: evaluate */
234 if (updateNeeded) {
235 for (c=0; c < argToVarLen ;c++) {
236 arg = RelationVariable(r,argToVar[c]);
237 value = RealAtomValue(arg);
238 common->inputsJac[c] = value;
239 }
240 common->interp.task = bb_deriv_eval;
241
242 if(derivFunc){
243 nok = (*derivFunc)(
244 &(common->interp)
245 , common->inputsLen, common->outputsLen
246 , common->inputsJac, common->outputs, common->jacobian
247 );
248 }else{
249 CONSOLE_DEBUG("COMPUTING FINITE-DIFFERENCE BLACK-BOX DERIVATIVE");
250
251 nok = blackbox_fdiff(GetValueFunc(efunc), &(common->interp)
252 , common->inputsLen, common->outputsLen
253 , common->inputsJac, common->outputs, common->jacobian
254 );
255 }
256 common->gradCount++;
257 }
258
259 for (k = 0; k < varlistLen; k++) {
260 gradient[k] = 0.0;
261 }
262
263 /* now compute d(y-yhat)/dx for this row as ( I - dyhat/dx ) */
264 k = lhsVar-1;
265 gradient[k] = 1.0; /* I */
266 offset = common->inputsLen * outputIndex; /* offset to row needed. */
267 bboxRow = common->jacobian + offset; /* pointer to row by ptr arithmetic */
268 for (c=0; c < argToVarLen; c++) {
269 k = ((int)(argToVar[c])) - 1; /* varnum is 1..n, not 0 .*/
270 gradient[k] -= bboxRow[c];
271 }
272
273 return nok;
274 }
275
276 struct Instance *BlackBoxGetOutputVar(struct relation *r)
277 {
278 assert(r != NULL);
279 /* unsigned long lhsVarNumber; */
280 /** @TODO FIXME BlackBoxGetOutputVar */
281 return NULL;
282 }
283
284 struct BlackBoxData *CreateBlackBoxData(struct BlackBoxCache *common,
285 unsigned long lhsindex,
286 unsigned long lhsVarNumber)
287 {
288 struct BlackBoxData *b = (struct BlackBoxData *)malloc(sizeof(struct BlackBoxData));
289 assert(common!=NULL);
290 b->common = common;
291 AddRefBlackBoxCache(common);
292 b->inputArgs = (unsigned long *)ascmalloc(sizeof(long) * (common->inputsLen));
293 b->lhsindex = lhsindex;
294 b->lhsvar = lhsVarNumber;
295 return b;
296 }
297
298 void DestroyBlackBoxData(struct relation *rel, struct BlackBoxData *b)
299 {
300 b->lhsindex = -(b->lhsindex);
301 b->lhsvar = 0;
302 ascfree(b->inputArgs);
303 b->inputArgs = NULL;
304 DeleteRefBlackBoxCache(rel, &(b->common));
305 ascfree(b);
306 }
307
308 /*------------------------------------------------------------------------------
309 BLACKBOX GRADIENT BY FINITE-DIFFERENCE
310 */
311
312 /**
313 This function works out what the peturbed value for a variable should be.
314
315 @TODO add some awareness of boundaries in this, hopefully.
316 */
317 static inline double blackbox_peturbation(double varvalue){
318 return (1.0e-05);
319 }
320
321 /**
322 Blackbox derivatives estimated by finite difference (by evaluation at
323 peturbed value of each input in turn)
324
325 Call signature as for ExtBBoxFunc.
326 */
327 int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
328 , int ninputs, int noutputs
329 , double *inputs, double *outputs, double *jac
330 ){
331 long c,r;
332 int nok = 0;
333 double *tmp_outputs;
334 double *ptr;
335 double old_x,deltax,value;
336
337 /* CONSOLE_DEBUG("NUMERICAL DERIVATIVE..."); */
338
339 tmp_outputs = ASC_NEW_ARRAY_CLEAR(double,noutputs);
340
341 for (c=0;c<ninputs;c++){
342 /* perturb each input in turn */
343 old_x = inputs[c];
344 deltax = blackbox_peturbation(old_x);
345 inputs[c] = old_x + deltax;
346 /* CONSOLE_DEBUG("PETURBATED VALUE of input[%ld] = %f",c,inputs[c]); */
347
348 /* call routine. note that the 'jac' parameter is just along for the ride */
349 nok = (*resfn)(interp, ninputs, noutputs, inputs, tmp_outputs, jac);
350 if(nok){
351 CONSOLE_DEBUG("External evaluation error (%d)",nok);
352 break;
353 }
354
355 /* fill load jacobian */
356 ptr = &jac[c];
357 for(r=0;r<noutputs;r++){
358 value = (tmp_outputs[r] - outputs[r]) / deltax;
359 /* CONSOLE_DEBUG("output[%ld]: value = %f, gradient = %f",r,tmp_outputs[r],value); */
360 *ptr = value;
361 ptr += ninputs;
362 }
363
364 /* now return this input to its old value */
365 inputs[c] = old_x;
366 }
367 ASC_FREE(tmp_outputs);
368 if(nok){
369 CONSOLE_DEBUG("External evaluation error");
370 }
371 return nok;
372
373 }
374
375 /*------------------------------------------------------------------------------
376 BLACK BOX CACHE
377 */
378
379 #define JACMAGIC -3.141592071828
380 struct BlackBoxCache *CreateBlackBoxCache(
381 int32 inputsLen,
382 int32 outputsLen,
383 struct gl_list_t *formalArgs
384 )
385 {
386 struct BlackBoxCache *b = (struct BlackBoxCache *)malloc(sizeof(struct BlackBoxCache));
387 b->interp.task = bb_none;
388 b->interp.status = calc_all_ok;
389 b->interp.user_data = NULL;
390 b->formalArgList = CopySpecialList(formalArgs);
391 b->inputsLen = inputsLen;
392 b->outputsLen = outputsLen;
393 b->jacobianLen = inputsLen*outputsLen;
394 b->hessianLen = 0; /* FIXME. */
395 b->inputs = (double *)ascmalloc(inputsLen*sizeof(double));
396 b->inputsJac = (double *)ascmalloc(inputsLen*sizeof(double));
397 b->outputs = (double *)ascmalloc(outputsLen*sizeof(double));
398 b->jacobian = (double*)ascmalloc(outputsLen*inputsLen*sizeof(double)+sizeof(double));
399 b->jacobian[outputsLen*inputsLen] = JACMAGIC;
400 b->hessian = NULL;
401 b->residCount = 0;
402 b->gradCount = 0;
403 b->refCount = 1;
404 return b;
405 }
406
407 void AddRefBlackBoxCache(struct BlackBoxCache *b)
408 {
409 assert(b != NULL);
410 (b->refCount)++;
411 }
412
413 static void DestroyBlackBoxCache(struct relation *rel, struct BlackBoxCache *b)
414 {
415 struct ExternalFunc *efunc;
416 ExtBBoxFinalFunc *final;
417 assert(b != NULL);
418 if (b->jacobian[b->outputsLen*b->inputsLen] != JACMAGIC)
419 {
420 FPRINTF(ASCERR, "Jacobian overrun detected while destroying BlackBoxCache. Debug the user blackbox implementation.\n");
421 }
422 ascfree(b->inputs);
423 ascfree(b->inputsJac);
424 ascfree(b->outputs);
425 ascfree(b->jacobian);
426 DestroySpecialList(b->formalArgList);
427 b->formalArgList = NULL;
428 b->inputsLen = -(b->inputsLen);
429 b->outputsLen = -(b->outputsLen);
430 b->jacobianLen = -(b->jacobianLen);
431 b->hessianLen = -(b->hessianLen);
432 b->inputs = NULL;
433 b->inputsJac = NULL;
434 b->outputs = NULL;
435 b->jacobian = NULL;
436 b->hessian = NULL;
437 b->residCount = -(b->residCount);
438 b->gradCount = -(b->gradCount);
439 b->refCount = -3;
440 if (rel != NULL) {
441 /* rel is null with refcount 0 only when there's a bbox array
442 instantiation failure. If failure, we don't need final. */
443 efunc = RelationBlackBoxExtFunc(rel);
444 final = GetFinalFunc(efunc);
445 b->interp.task = bb_last_call;
446 (*final)(&(b->interp));
447 }
448 ascfree(b);
449 }
450
451 void DeleteRefBlackBoxCache(struct relation *rel, struct BlackBoxCache **b)
452 {
453 struct BlackBoxCache * d = *b;
454 assert(b != NULL && *b != NULL);
455 *b = NULL;
456 if (d->refCount > 0) {
457 (d->refCount)--;
458 }
459 if (d->refCount == 0) {
460 DestroyBlackBoxCache(rel,d);
461 return;
462 }
463 if (d->refCount < 0) {
464 FPRINTF(ASCERR, "Attempt to delete BlackBoxCache (%p) too often.\n", *b);
465 }
466 }

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