/[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 915 - (show annotations) (download) (as text)
Sat Oct 28 07:06:22 2006 UTC (13 years, 8 months ago) by johnpye
File MIME type: text/x-csrc
File size: 12518 byte(s)
Working on what the problem is with 'on_load' methods in the C++/Python code.
Seems that 'EXTERNAL solve(SELF);' puts the model into a state which 'works' but if that's
not done first, then that set values in the model don't work correctly...
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, warnfdiff;
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 if(!warnfdiff){
250 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Using finite-difference "
251 " to compute derivatives for one or more black boxes in"
252 " this model."
253 );
254 warnfdiff = 1;
255 }
256
257 nok = blackbox_fdiff(GetValueFunc(efunc), &(common->interp)
258 , common->inputsLen, common->outputsLen
259 , common->inputsJac, common->outputs, common->jacobian
260 );
261 }
262 common->gradCount++;
263 }
264
265 for (k = 0; k < varlistLen; k++) {
266 gradient[k] = 0.0;
267 }
268
269 /* now compute d(y-yhat)/dx for this row as ( I - dyhat/dx ) */
270 k = lhsVar-1;
271 gradient[k] = 1.0; /* I */
272 offset = common->inputsLen * outputIndex; /* offset to row needed. */
273 bboxRow = common->jacobian + offset; /* pointer to row by ptr arithmetic */
274 for (c=0; c < argToVarLen; c++) {
275 k = ((int)(argToVar[c])) - 1; /* varnum is 1..n, not 0 .*/
276 gradient[k] -= bboxRow[c];
277 }
278
279 return nok;
280 }
281
282 struct Instance *BlackBoxGetOutputVar(struct relation *r)
283 {
284 assert(r != NULL);
285 /* unsigned long lhsVarNumber; */
286 /** @TODO FIXME BlackBoxGetOutputVar */
287 return NULL;
288 }
289
290 struct BlackBoxData *CreateBlackBoxData(struct BlackBoxCache *common,
291 unsigned long lhsindex,
292 unsigned long lhsVarNumber)
293 {
294 struct BlackBoxData *b = (struct BlackBoxData *)malloc(sizeof(struct BlackBoxData));
295 assert(common!=NULL);
296 b->common = common;
297 AddRefBlackBoxCache(common);
298 b->inputArgs = (unsigned long *)ascmalloc(sizeof(long) * (common->inputsLen));
299 b->lhsindex = lhsindex;
300 b->lhsvar = lhsVarNumber;
301 return b;
302 }
303
304 void DestroyBlackBoxData(struct relation *rel, struct BlackBoxData *b)
305 {
306 b->lhsindex = -(b->lhsindex);
307 b->lhsvar = 0;
308 ascfree(b->inputArgs);
309 b->inputArgs = NULL;
310 DeleteRefBlackBoxCache(rel, &(b->common));
311 ascfree(b);
312 }
313
314 /*------------------------------------------------------------------------------
315 BLACKBOX GRADIENT BY FINITE-DIFFERENCE
316 */
317
318 /**
319 This function works out what the peturbed value for a variable should be.
320
321 @TODO add some awareness of boundaries in this, hopefully.
322 */
323 static inline double blackbox_peturbation(double varvalue){
324 return (1.0e-05);
325 }
326
327 /**
328 Blackbox derivatives estimated by finite difference (by evaluation at
329 peturbed value of each input in turn)
330
331 Call signature as for ExtBBoxFunc.
332 */
333 int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
334 , int ninputs, int noutputs
335 , double *inputs, double *outputs, double *jac
336 ){
337 long c,r;
338 int nok = 0;
339 double *tmp_outputs;
340 double *ptr;
341 double old_x,deltax,value;
342
343 /* CONSOLE_DEBUG("NUMERICAL DERIVATIVE..."); */
344
345 tmp_outputs = ASC_NEW_ARRAY_CLEAR(double,noutputs);
346
347 for (c=0;c<ninputs;c++){
348 /* perturb each input in turn */
349 old_x = inputs[c];
350 deltax = blackbox_peturbation(old_x);
351 inputs[c] = old_x + deltax;
352 /* CONSOLE_DEBUG("PETURBATED VALUE of input[%ld] = %f",c,inputs[c]); */
353
354 /* call routine. note that the 'jac' parameter is just along for the ride */
355 nok = (*resfn)(interp, ninputs, noutputs, inputs, tmp_outputs, jac);
356 if(nok){
357 CONSOLE_DEBUG("External evaluation error (%d)",nok);
358 break;
359 }
360
361 /* fill load jacobian */
362 ptr = &jac[c];
363 for(r=0;r<noutputs;r++){
364 value = (tmp_outputs[r] - outputs[r]) / deltax;
365 /* CONSOLE_DEBUG("output[%ld]: value = %f, gradient = %f",r,tmp_outputs[r],value); */
366 *ptr = value;
367 ptr += ninputs;
368 }
369
370 /* now return this input to its old value */
371 inputs[c] = old_x;
372 }
373 ASC_FREE(tmp_outputs);
374 if(nok){
375 CONSOLE_DEBUG("External evaluation error");
376 }
377 return nok;
378
379 }
380
381 /*------------------------------------------------------------------------------
382 BLACK BOX CACHE
383 */
384
385 #define JACMAGIC -3.141592071828
386 struct BlackBoxCache *CreateBlackBoxCache(
387 int32 inputsLen,
388 int32 outputsLen,
389 struct gl_list_t *formalArgs
390 )
391 {
392 struct BlackBoxCache *b = (struct BlackBoxCache *)malloc(sizeof(struct BlackBoxCache));
393 b->interp.task = bb_none;
394 b->interp.status = calc_all_ok;
395 b->interp.user_data = NULL;
396 b->formalArgList = CopySpecialList(formalArgs);
397 b->inputsLen = inputsLen;
398 b->outputsLen = outputsLen;
399 b->jacobianLen = inputsLen*outputsLen;
400 b->hessianLen = 0; /* FIXME. */
401 b->inputs = (double *)ascmalloc(inputsLen*sizeof(double));
402 b->inputsJac = (double *)ascmalloc(inputsLen*sizeof(double));
403 b->outputs = (double *)ascmalloc(outputsLen*sizeof(double));
404 b->jacobian = (double*)ascmalloc(outputsLen*inputsLen*sizeof(double)+sizeof(double));
405 b->jacobian[outputsLen*inputsLen] = JACMAGIC;
406 b->hessian = NULL;
407 b->residCount = 0;
408 b->gradCount = 0;
409 b->refCount = 1;
410 return b;
411 }
412
413 void AddRefBlackBoxCache(struct BlackBoxCache *b)
414 {
415 assert(b != NULL);
416 (b->refCount)++;
417 }
418
419 static void DestroyBlackBoxCache(struct relation *rel, struct BlackBoxCache *b)
420 {
421 struct ExternalFunc *efunc;
422 ExtBBoxFinalFunc *final;
423 assert(b != NULL);
424 if (b->jacobian[b->outputsLen*b->inputsLen] != JACMAGIC)
425 {
426 FPRINTF(ASCERR, "Jacobian overrun detected while destroying BlackBoxCache. Debug the user blackbox implementation.\n");
427 }
428 ascfree(b->inputs);
429 ascfree(b->inputsJac);
430 ascfree(b->outputs);
431 ascfree(b->jacobian);
432 DestroySpecialList(b->formalArgList);
433 b->formalArgList = NULL;
434 b->inputsLen = -(b->inputsLen);
435 b->outputsLen = -(b->outputsLen);
436 b->jacobianLen = -(b->jacobianLen);
437 b->hessianLen = -(b->hessianLen);
438 b->inputs = NULL;
439 b->inputsJac = NULL;
440 b->outputs = NULL;
441 b->jacobian = NULL;
442 b->hessian = NULL;
443 b->residCount = -(b->residCount);
444 b->gradCount = -(b->gradCount);
445 b->refCount = -3;
446 if (rel != NULL) {
447 /* rel is null with refcount 0 only when there's a bbox array
448 instantiation failure. If failure, we don't need final. */
449 efunc = RelationBlackBoxExtFunc(rel);
450 final = GetFinalFunc(efunc);
451 b->interp.task = bb_last_call;
452 (*final)(&(b->interp));
453 }
454 ascfree(b);
455 }
456
457 void DeleteRefBlackBoxCache(struct relation *rel, struct BlackBoxCache **b)
458 {
459 struct BlackBoxCache * d = *b;
460 assert(b != NULL && *b != NULL);
461 *b = NULL;
462 if (d->refCount > 0) {
463 (d->refCount)--;
464 }
465 if (d->refCount == 0) {
466 DestroyBlackBoxCache(rel,d);
467 return;
468 }
469 if (d->refCount < 0) {
470 FPRINTF(ASCERR, "Attempt to delete BlackBoxCache (%p) too often.\n", *b);
471 }
472 }

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