# Diff of /trunk/base/generic/compiler/rel_blackbox.c

revision 908 by johnpye, Thu Oct 26 10:18:53 2006 UTC revision 909 by johnpye, Thu Oct 26 12:44:41 2006 UTC
# Line 41  static int32 ArgsDifferent(double new, d Line 41  static int32 ArgsDifferent(double new, d
41      }      }
42  }  }
43
44    /*------------------------------------------------------------------------------
45      RESIDUAL AND GRADIENT EVALUATION ROUTINES
46    */
47
48  /*  /*
49      Note:      Note:
50      Assume a function in c/fortran that computes yhat(x),      Assume a function in c/fortran that computes yhat(x),
# Line 149  int BlackBoxCalcResidual(struct Instance Line 153  int BlackBoxCalcResidual(struct Instance
153
154  }  }
155
156  /* The tricky bit in this is that if input or output args are merged,  /**
157  their partial derivatives get summed and the solver at least gets      Calculate the gradient (slice of the overall jacobian) for the blackbox.
158  what it deserves, which may or may not be sane in a modeling sense.
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  */  */
169  {  int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
170        , int ninputs, int noutputs
171        , double *inputs, double *outputs, double *jac
172    );
173
175            , struct relation *r
176    ){
177  /* decls */  /* decls */
178      unsigned long *argToVar;      unsigned long *argToVar;
179      unsigned long c, varlistLen;      unsigned long c, varlistLen;
192      unsigned int k;      unsigned int k;
193      int nok = 0;      int nok = 0;
194
195  /* impl setup */      /* prepare */
197      (void)i;      (void)i;
198      efunc = RelationBlackBoxExtFunc(r);      efunc = RelationBlackBoxExtFunc(r);
207      updateNeeded = 0;      updateNeeded = 0;
208      varlistLen = NumberVariables(r);      varlistLen = NumberVariables(r);
209
210  /* impl:      /* do we need to calculate? */
211      check input values changed.      if(common->gradCount < 1){
if changed, recompute gradient in bbox.
compute gradient per varlist from bbox row.
*/
/* check: */
212          updateNeeded = 1;          updateNeeded = 1;
213      } else {      }else{
214          for (c=0; c < argToVarLen ; c++) {          for(c=0; c<argToVarLen; c++){
215              arg = RelationVariable(r,argToVar[c]);              arg = RelationVariable(r,argToVar[c]);
216              value = RealAtomValue(arg);              value = RealAtomValue(arg);
217              if (ArgsDifferent(value, common->inputsJac[c], inputTolerance)) {              if (ArgsDifferent(value, common->inputsJac[c], inputTolerance)) {
218                  updateNeeded = 1;                  updateNeeded = 1;
220              }              }
221          }          }
222      }      }
223      /* recompute */
224        /* if inputs have changed more than allowed, or it's first time: evaluate */
225      if (updateNeeded) {      if (updateNeeded) {
226          for (c=0; c < argToVarLen ;c++) {          for (c=0; c < argToVarLen ;c++) {
227              arg = RelationVariable(r,argToVar[c]);              arg = RelationVariable(r,argToVar[c]);
230          }          }
232
233          nok = (*derivFunc)(&(common->interp),          if(derivFunc){
234                  common->inputsLen,              nok = (*derivFunc)(
235                  common->outputsLen,                  &(common->interp)
236                  common->inputsJac,                  , common->inputsLen, common->outputsLen
237                  common->outputs,                  , common->inputsJac, common->outputs, common->jacobian
238                  common->jacobian);              );
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            }
248      }      }
249
250      for (k = 0; k < varlistLen; k++) {      for (k = 0; k < varlistLen; k++) {
252      }      }
253      /* now compute d(y-yhat)/dx for this row as I - dyhat/dx*/
254        /* now compute d(y-yhat)/dx for this row as ( I - dyhat/dx ) */
255      k = lhsVar-1;      k = lhsVar-1;
256      gradient[k] = 1.0; /* I */      gradient[k] = 1.0; /* I */
257      offset = common->inputsLen * outputIndex; /* offset to row needed. */      offset = common->inputsLen * outputIndex; /* offset to row needed. */
267  struct Instance *BlackBoxGetOutputVar(struct relation *r)  struct Instance *BlackBoxGetOutputVar(struct relation *r)
268  {  {
269      assert(r != NULL);      assert(r != NULL);
270      unsigned long lhsVarNumber;      /* unsigned long lhsVarNumber; */
271      /* FIXME BlackBoxGetOutputVar */          /** @TODO FIXME BlackBoxGetOutputVar */
272      return NULL;      return NULL;
273  }  }
274
# Line 271  void DestroyBlackBoxData(struct relation Line 296  void DestroyBlackBoxData(struct relation
296      ascfree(b);      ascfree(b);
297  }  }
298
299    /*------------------------------------------------------------------------------
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  #define JACMAGIC -3.141592071828
371  struct BlackBoxCache *CreateBlackBoxCache(  struct BlackBoxCache *CreateBlackBoxCache(
372      int32 inputsLen,      int32 inputsLen,

Legend:
 Removed from v.908 changed lines Added in v.909