/[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 1034 - (show annotations) (download) (as text)
Thu Jan 4 05:37:55 2007 UTC (13 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 16452 byte(s)
Switch 'void slv_*' function to output a 0-on-success error flag.
Added some debug output to ensure that such output gets through.
Fixed TestBlackBox.testfail1 to work, using above changes.
Simulation::solve now throws an exception on failure (will need to modify PyGTK GUI accordingly)
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/pairlist.h>
10 #include <general/dstring.h>
11 #include "compiler.h"
12 #include "symtab.h"
13 #include "functype.h"
14 #include "safe.h"
15 #include "fractions.h"
16 #include "dimen.h"
17 #include "expr_types.h"
18 #include "vlist.h"
19 #include "dimen_io.h"
20 #include "instance_enum.h"
21 #include "find.h"
22 #include "atomvalue.h"
23 #include "instance_name.h"
24 #include "extfunc.h"
25 #include "relation_type.h"
26 #include "rel_blackbox.h"
27 #include "relation.h"
28 #include "relation_util.h"
29 #include "instance_io.h"
30 #include "relation_io.h"
31 #include "instquery.h"
32 #include "visitinst.h"
33 #include "mathinst.h"
34 #include "extcall.h" /* for copy/destroy speciallist */
35 #include "packages.h" /* for init slv interp */
36 #include "name.h" /* for copy/destroy name */
37
38 #define BBDEBUG 0 /* set 0 if not wanting spew */
39
40 static int32 ArgsDifferent(double new, double old, double tol)
41 {
42 if (fabs(new - old) > fabs(tol)) {
43 return 1;
44 } else {
45 return 0;
46 }
47 }
48
49 /*------------------------------------------------------------------------------
50 RESIDUAL AND GRADIENT EVALUATION ROUTINES
51 */
52
53 /*
54 Note:
55 Assume a function in c/fortran that computes yhat(x),
56 aka func(inputs_x, outputs_yhat);.
57 This is fit into the ascend world by associating
58 inputs_x directly to ascend variables and defining
59 the ascend variables y in the outer model.
60 The residual of the blackbox equations in the ascend model
61 is y-yhat. The gradient is I - dyhat/dx, where dyhat/dx
62 is the reduced jacobian of the blackbox.
63 */
64 int BlackBoxCalcResidGrad(struct Instance *i, double *res, double *gradient, struct relation *r)
65 {
66 int residErr = 0;
67 int gradErr = 0;
68
69 residErr = BlackBoxCalcResidual(i, res, r);
70 gradErr = BlackBoxCalcGradient(i, gradient, r);
71 return (int)(fabs(residErr) + fabs(gradErr));
72 }
73
74 /*
75 Note:
76 Assume a function in c/fortran that computes yhat(x),
77 aka func(inputs_x, outputs_yhat);.
78 This is fit into the ascend world by associating
79 inputs_x directly to ascend variables and defining
80 the ascend variables y in the outer model.
81 The residual of the blackbox equations in the ascend model
82 is y-yhat. The gradient is I - dyhat/dx, where dyhat/dx
83 is the reduced jacobian of the blackbox.
84 */
85 int BlackBoxCalcResidual(struct Instance *i, double *res, struct relation *r)
86 {
87 /* decls */
88 unsigned long *argToVar;
89 unsigned long c;
90 unsigned long argToVarLen;
91 struct BlackBoxCache *common;
92 struct ExternalFunc *efunc;
93 ExtBBoxFunc *evalFunc;
94 int lhsVar;
95 int outputIndex;
96 struct Instance *arg;
97 double value;
98 double inputTolerance;
99 int updateNeeded;
100 int nok = 0;
101 static int warnexpt;
102
103 if(!warnexpt){
104 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox evaluation is experimental (%s)",__FUNCTION__);
105 warnexpt = 1;
106 }
107
108 (void)i;
109 efunc = RelationBlackBoxExtFunc(r);
110 common = RelationBlackBoxCache(r);
111 assert(efunc && common);
112 argToVarLen = common->inputsLen;
113 argToVar = RBBOX(r).inputArgs;
114 lhsVar = RBBOX(r).lhsvar;
115 outputIndex = RBBOX(r).lhsindex;
116 evalFunc = GetValueFunc(efunc);
117 inputTolerance = GetValueFuncTolerance(efunc);
118 updateNeeded = 0;
119
120 /* impl:
121 check input values changed.
122 if changed, recompute output values.
123 compute y - yhat for lhsvar.
124 */
125 /* check: */
126 if (common->residCount < 1) {
127 updateNeeded = 1;
128 } else {
129 for (c=0; c < argToVarLen ; c++) {
130 arg = RelationVariable(r,argToVar[c]);
131 value = RealAtomValue(arg);
132 if (ArgsDifferent(value, common->inputs[c], inputTolerance)) {
133 updateNeeded = 1;
134 break;
135 }
136 }
137 }
138 /* recompute */
139 if (updateNeeded) {
140 for (c=0; c < argToVarLen ;c++) {
141 arg = RelationVariable(r,argToVar[c]);
142 value = RealAtomValue(arg);
143 common->inputs[c] = value;
144 }
145 common->interp.task = bb_func_eval;
146
147 nok = (*evalFunc)(&(common->interp),
148 common->inputsLen,
149 common->outputsLen,
150 common->inputs,
151 common->outputs,
152 common->jacobian);
153 if(nok)CONSOLE_DEBUG("blackbox residual function returned error %d",nok);
154 common->residCount++;
155 }
156 value = common->outputs[outputIndex];
157
158 /* return lhsval - value */
159 arg = RelationVariable(r,lhsVar);
160 *res = RealAtomValue(arg) - value;
161 return nok;
162
163 }
164
165 /**
166 Calculate the gradient (slice of the overall jacobian) for the blackbox.
167
168 The tricky bit in this is that if input or output args are merged,
169 their partial derivatives get summed and the solver at least gets
170 what it deserves, which may or may not be sane in a modeling sense.
171
172 Implementation:
173 # check input values changed.
174 # if changed, recompute gradient in bbox.
175 # compute gradient per varlist from bbox row.
176 */
177
178 int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
179 , int ninputs, int noutputs
180 , double *inputs, double *outputs, double *jac
181 );
182
183 int BlackBoxCalcGradient(struct Instance *i, double *gradient
184 , struct relation *r
185 ){
186 /* decls */
187 unsigned long *argToVar;
188 unsigned long c, varlistLen;
189 unsigned long argToVarLen;
190 struct BlackBoxCache *common;
191 struct ExternalFunc *efunc;
192 ExtBBoxFunc *derivFunc;
193 int lhsVar;
194 int outputIndex;
195 struct Instance *arg;
196 double value;
197 double inputTolerance;
198 double *bboxRow;
199 int updateNeeded, offset;
200 unsigned int k;
201 int nok = 0;
202 static int warnexpt, warnfdiff;
203
204 if(!warnexpt){
205 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox gradient is experimental (%s)",__FUNCTION__);
206 warnexpt = 1;
207 }
208
209 /* prepare */
210 (void)i;
211 efunc = RelationBlackBoxExtFunc(r);
212 common = RelationBlackBoxCache(r);
213 argToVarLen = common->inputsLen;
214 argToVar = RBBOX(r).inputArgs;
215 lhsVar = RBBOX(r).lhsvar;
216 outputIndex = RBBOX(r).lhsindex;
217 derivFunc = GetDerivFunc(efunc);
218 inputTolerance = GetValueFuncTolerance(efunc);
219 updateNeeded = 0;
220 varlistLen = NumberVariables(r);
221
222 /* do we need to calculate? */
223 if(common->gradCount < 1){
224 updateNeeded = 1;
225 }else{
226 for(c=0; c<argToVarLen; c++){
227 arg = RelationVariable(r,argToVar[c]);
228 value = RealAtomValue(arg);
229 if (ArgsDifferent(value, common->inputsJac[c], inputTolerance)) {
230 updateNeeded = 1;
231 break;
232 }
233 }
234 }
235
236 /* if inputs have changed more than allowed, or it's first time: evaluate */
237 if (updateNeeded) {
238 for (c=0; c < argToVarLen ;c++) {
239 arg = RelationVariable(r,argToVar[c]);
240 value = RealAtomValue(arg);
241 common->inputsJac[c] = value;
242 }
243 common->interp.task = bb_deriv_eval;
244
245 if(derivFunc){
246 nok = (*derivFunc)(
247 &(common->interp)
248 , common->inputsLen, common->outputsLen
249 , common->inputsJac, common->outputs, common->jacobian
250 );
251 }else{
252 if(!warnfdiff){
253 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Using finite-difference "
254 " to compute derivatives for one or more black boxes in"
255 " this model."
256 );
257 warnfdiff = 1;
258 }
259
260 nok = blackbox_fdiff(GetValueFunc(efunc), &(common->interp)
261 , common->inputsLen, common->outputsLen
262 , common->inputsJac, common->outputs, common->jacobian
263 );
264 }
265 common->gradCount++;
266 }
267
268 for (k = 0; k < varlistLen; k++) {
269 gradient[k] = 0.0;
270 }
271
272 /* now compute d(y-yhat)/dx for this row as ( I - dyhat/dx ) */
273 k = lhsVar-1;
274 gradient[k] = 1.0; /* I */
275 offset = common->inputsLen * outputIndex; /* offset to row needed. */
276 bboxRow = common->jacobian + offset; /* pointer to row by ptr arithmetic */
277 for (c=0; c < argToVarLen; c++) {
278 k = ((int)(argToVar[c])) - 1; /* varnum is 1..n, not 0 .*/
279 gradient[k] -= bboxRow[c];
280 }
281
282 return nok;
283 }
284
285 struct Instance *BlackBoxGetOutputVar(CONST struct relation *r)
286 {
287 assert(r != NULL);
288 unsigned long lhsVarNumber;
289 struct Instance *i;
290
291 lhsVarNumber = RBBOX(r).lhsvar;
292 i = RelationVariable(r, lhsVarNumber);
293 return i;
294 }
295
296 static int g_cbbdcount=0;
297 struct BlackBoxData *CreateBlackBoxData(struct BlackBoxCache *common)
298 {
299 struct BlackBoxData *b = (struct BlackBoxData *)malloc(sizeof(struct BlackBoxData));
300 g_cbbdcount++;
301 b->count = g_cbbdcount;
302 assert(common!=NULL);
303 b->common = common;
304 AddRefBlackBoxCache(common);
305 #if BBDEBUG
306 FPRINTF(ASCERR,"CreateBlackBoxData(%p) made BBD#%d (%p)\n",common,b->count, b);
307 #endif
308 return b;
309 }
310
311 static
312 struct BlackBoxCache *CloneBlackBoxCache(struct BlackBoxCache *b)
313 {
314 return CreateBlackBoxCache(b->inputsLen , b->outputsLen, b->argListNames, b->dataName, b->efunc);
315 }
316
317 void CopyBlackBoxDataByReference(struct relation *src, struct relation *dest, void * bboxtable_p)
318 {
319 unsigned long entry;
320 struct pairlist_t *bboxtable = (struct pairlist_t *)bboxtable_p;
321 struct BlackBoxCache * newCache = NULL, *oldCache;
322 struct BlackBoxData *b;
323 struct BlackBoxData *a;
324 assert(src);
325 assert(dest);
326 a = RelationBlackBoxData(src);
327 oldCache = RelationBlackBoxCache(src);
328 assert(a != NULL);
329
330 /* at moment, key is src cache and val is first dest bbd */
331 entry = pairlist_contains(bboxtable, oldCache);
332 if (entry) {
333 /* everyone gets a unique bbd and shared cache */
334 b = (struct BlackBoxData *)pairlist_valueAt(bboxtable, entry);
335 #if BBDEBUG
336 FPRINTF(ASCERR,"CopyBlackBoxDataByReference(%p, %p): found already cloned cache C#%d (%p)in BBD#%d (%p)\n",src,dest,b->common->count, b->common, b->count, b);
337 #endif
338 newCache = b->common;
339 b = CreateBlackBoxData(newCache);
340 } else {
341 /* dup common */
342 /* This is common across a single box output array of rels,
343 not across all box instances sharing relation data.
344 All the array pointers, etc are different.
345 */
346 newCache = CloneBlackBoxCache(oldCache); /* starts refc 1 for our scope. */
347 b = CreateBlackBoxData(newCache);
348 DeleteRefBlackBoxCache(dest, &newCache); /* left our scope. now null */
349 pairlist_append(bboxtable, oldCache, b); /* next relation in same array will find its data in the table now, because it has the same a. */
350 #if BBDEBUG
351 FPRINTF(ASCERR,"CopyBlackBoxDataByReference(%p, %p): cloned cache C#%d (%p)into BBD#%d (%p) with new cache C#%d (%p)\n",src,dest,oldCache->count, oldCache, b->count, b,b->common->count, b->common);
352 #endif
353 }
354
355 dest->externalData = b;
356 }
357
358 void DestroyBlackBoxData(struct relation *rel, struct BlackBoxData *b)
359 {
360 #if BBDEBUG
361 FPRINTF(ASCERR,"DestroyBlackBoxData(%p): destroying bbd %p BBD#%d\n",rel,b,b->count);
362 #endif
363 DeleteRefBlackBoxCache(rel, &(b->common));
364 b->count *= -1;
365 ascfree(b);
366 }
367
368 /*------------------------------------------------------------------------------
369 BLACKBOX GRADIENT BY FINITE-DIFFERENCE
370 */
371
372 /**
373 This function works out what the peturbed value for a variable should be.
374
375 @TODO add some awareness of boundaries in this, hopefully.
376 */
377 static inline double blackbox_peturbation(double varvalue){
378 return (1.0e-05);
379 }
380
381 /**
382 Blackbox derivatives estimated by finite difference (by evaluation at
383 peturbed value of each input in turn)
384
385 Call signature as for ExtBBoxFunc.
386 */
387 int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
388 , int ninputs, int noutputs
389 , double *inputs, double *outputs, double *jac
390 ){
391 long c,r;
392 int nok = 0;
393 double *tmp_outputs;
394 double *ptr;
395 double old_x,deltax,value;
396
397 /* CONSOLE_DEBUG("NUMERICAL DERIVATIVE..."); */
398
399 tmp_outputs = ASC_NEW_ARRAY_CLEAR(double,noutputs);
400
401 for (c=0;c<ninputs;c++){
402 /* perturb each input in turn */
403 old_x = inputs[c];
404 deltax = blackbox_peturbation(old_x);
405 inputs[c] = old_x + deltax;
406 /* CONSOLE_DEBUG("PETURBATED VALUE of input[%ld] = %f",c,inputs[c]); */
407
408 /* call routine. note that the 'jac' parameter is just along for the ride */
409 nok = (*resfn)(interp, ninputs, noutputs, inputs, tmp_outputs, jac);
410 if(nok){
411 CONSOLE_DEBUG("External evaluation error (%d)",nok);
412 break;
413 }
414
415 /* fill load jacobian */
416 ptr = &jac[c];
417 for(r=0;r<noutputs;r++){
418 value = (tmp_outputs[r] - outputs[r]) / deltax;
419 /* CONSOLE_DEBUG("output[%ld]: value = %f, gradient = %f",r,tmp_outputs[r],value); */
420 *ptr = value;
421 ptr += ninputs;
422 }
423
424 /* now return this input to its old value */
425 inputs[c] = old_x;
426 }
427 ASC_FREE(tmp_outputs);
428 if(nok){
429 CONSOLE_DEBUG("External evaluation error");
430 }
431 return nok;
432
433 }
434
435 /*------------------------------------------------------------------------------
436 BLACK BOX CACHE
437 */
438
439 static int g_cbbccount = 0;
440
441 #define JACMAGIC -3.141592071828
442 struct BlackBoxCache *CreateBlackBoxCache(
443 int32 inputsLen,
444 int32 outputsLen,
445 struct gl_list_t *argListNames,
446 struct Name *dataName,
447 struct ExternalFunc *efunc
448 )
449 {
450 struct BlackBoxCache *b = (struct BlackBoxCache *)malloc(sizeof(struct BlackBoxCache));
451 g_cbbccount++;
452 b->count = g_cbbccount;
453 b->interp.task = bb_none;
454 b->interp.status = calc_all_ok;
455 b->interp.user_data = NULL;
456 b->argListNames = DeepCopySpecialList(argListNames,(CopyFunc)CopyName);
457 b->dataName = CopyName(dataName);
458 b->inputsLen = inputsLen;
459 b->outputsLen = outputsLen;
460 b->jacobianLen = inputsLen*outputsLen;
461 b->hessianLen = 0; /* FIXME. */
462 b->inputs = (double *)ascmalloc(inputsLen*sizeof(double));
463 b->inputsJac = (double *)ascmalloc(inputsLen*sizeof(double));
464 b->outputs = (double *)ascmalloc(outputsLen*sizeof(double));
465 b->jacobian = (double*)ascmalloc(outputsLen*inputsLen*sizeof(double)+sizeof(double));
466 b->jacobian[outputsLen*inputsLen] = JACMAGIC;
467 b->hessian = NULL;
468 b->residCount = 0;
469 b->gradCount = 0;
470 b->refCount = 1;
471 b->efunc = efunc;
472 return b;
473 }
474
475 void InitBBox(struct Instance *context, struct BlackBoxCache *b)
476 {
477 ExtBBoxInitFunc * init;
478 enum find_errors ferr = correct_instance;
479 unsigned long nbr, br, errpos;
480 struct gl_list_t *tmp;
481
482 struct gl_list_t *arglist;
483 struct Instance *data;
484
485 /* fish up data from name. */
486 if (b->dataName != NULL) {
487 tmp = FindInstances(context,b->dataName,&ferr);
488 assert(tmp != NULL);
489 assert(gl_length(tmp) == 1);
490 data = (struct Instance *)gl_fetch(tmp,1);
491 gl_destroy(tmp);
492 assert(data != NULL);
493 } else {
494 data = NULL;
495 }
496
497 /* convert arglistnames to arglist. */
498 WriteNamesInList2D(ASCERR, b->argListNames, ", ", "\n");
499 nbr = gl_length(b->argListNames);
500 arglist = gl_create(nbr);
501 for (br = 1; br <= nbr; br++) {
502 tmp = (struct gl_list_t *)gl_fetch(b->argListNames,br);
503 tmp = FindInstancesFromNames(context, tmp, &ferr, &errpos);
504 assert(tmp != NULL);
505 gl_append_ptr(arglist,tmp);
506 }
507
508 /* now do the init */
509 init = GetInitFunc(b->efunc);
510 b->interp.task = bb_first_call;
511 (*init)( &(b->interp), data, arglist);
512 b->interp.task = bb_none;
513 }
514
515 int32 BlackBoxCacheInputsLen(struct BlackBoxCache *b)
516 {
517 assert(b != NULL);
518 return b->inputsLen;
519 }
520
521 void AddRefBlackBoxCache(struct BlackBoxCache *b)
522 {
523 assert(b != NULL);
524 (b->refCount)++;
525 }
526
527 static void DestroyBlackBoxCache(struct relation *rel, struct BlackBoxCache *b)
528 {
529 struct ExternalFunc *efunc;
530 ExtBBoxFinalFunc *final;
531 assert(b != NULL);
532 if (b->jacobian[b->outputsLen*b->inputsLen] != JACMAGIC)
533 {
534 ERROR_REPORTER_HERE(ASC_PROG_WARNING, "Jacobian overrun detected while destroying BlackBoxCache. Debug the user blackbox implementation.\n");
535 }
536 #if BBDEBUG
537 FPRINTF(ASCERR,"DestroyBlackBoxCache(%p,%p): destroying %d\n",rel,b,b->count);
538 #endif
539 ascfree(b->inputs);
540 ascfree(b->inputsJac);
541 ascfree(b->outputs);
542 ascfree(b->jacobian);
543 DeepDestroySpecialList(b->argListNames,(DestroyFunc)DestroyName);
544 b->argListNames = NULL;
545 DestroyName(b->dataName);
546 b->dataName = NULL;
547 b->inputsLen = -(b->inputsLen);
548 b->outputsLen = -(b->outputsLen);
549 b->jacobianLen = -(b->jacobianLen);
550 b->hessianLen = -(b->hessianLen);
551 b->inputs = NULL;
552 b->inputsJac = NULL;
553 b->outputs = NULL;
554 b->jacobian = NULL;
555 b->hessian = NULL;
556 b->residCount = -(b->residCount);
557 b->gradCount = -(b->gradCount);
558 b->refCount = -3;
559 if (rel != NULL) {
560 /* rel is null with refcount 0 only when there's a bbox array
561 instantiation failure or we're in a weird copy process.
562 Either way, we don't need final. */
563 efunc = b->efunc;
564 final = GetFinalFunc(efunc);
565 b->interp.task = bb_last_call;
566 (*final)(&(b->interp));
567 b->efunc = NULL;
568 }
569 b->count *= -1;
570 ascfree(b);
571 }
572
573 void DeleteRefBlackBoxCache(struct relation *rel, struct BlackBoxCache **b)
574 {
575 struct BlackBoxCache * d = *b;
576 assert(b != NULL && *b != NULL);
577 if (d->refCount > 0) {
578 (d->refCount)--;
579 #if BBDEBUG
580 FPRINTF(ASCERR,"DeleteRefBlackBoxCache(%p,%p): C#%d refcount reduced to %d in bbox %p\n",rel,d,d->count,d->refCount,b);
581 #endif
582 *b = NULL;
583 }
584 if (d->refCount == 0) {
585 #if BBDEBUG
586 FPRINTF(ASCERR,"DeleteRefBlackBoxCache(%p,%p): bbc refcount reduced to %d in bbox %p\n",rel,d,d->count,b);
587 #endif
588 DestroyBlackBoxCache(rel,d);
589 *b = NULL;
590 return;
591 }
592 if (d->refCount < 0) {
593 #if BBDEBUG
594 FPRINTF(ASCERR, "Attempt to delete BlackBoxCache (%p) too often.\n", *b);
595 #endif
596 *b = NULL;
597 }
598 }

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