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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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  */  */
168  int BlackBoxCalcGradient(struct Instance *i, double *gradient, struct relation *r)  
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 */  /* decls */
178      unsigned long *argToVar;      unsigned long *argToVar;
179      unsigned long c, varlistLen;      unsigned long c, varlistLen;
# Line 173  int BlackBoxCalcGradient(struct Instance Line 192  int BlackBoxCalcGradient(struct Instance
192      unsigned int k;      unsigned int k;
193      int nok = 0;      int nok = 0;
194        
195  /* impl setup */      /* prepare */
196      ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox gradient is experimental (%s)",__FUNCTION__);      ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox gradient is experimental (%s)",__FUNCTION__);
197      (void)i;      (void)i;
198      efunc = RelationBlackBoxExtFunc(r);      efunc = RelationBlackBoxExtFunc(r);
# Line 188  int BlackBoxCalcGradient(struct Instance Line 207  int BlackBoxCalcGradient(struct Instance
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: */  
     if (common->gradCount < 1) {  
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;
# Line 206  int BlackBoxCalcGradient(struct Instance Line 220  int BlackBoxCalcGradient(struct Instance
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]);
# Line 215  int BlackBoxCalcGradient(struct Instance Line 230  int BlackBoxCalcGradient(struct Instance
230          }          }
231          common->interp.task = bb_deriv_eval;          common->interp.task = bb_deriv_eval;
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            }
247          common->gradCount++;          common->gradCount++;
248      }      }
249    
250      for (k = 0; k < varlistLen; k++) {      for (k = 0; k < varlistLen; k++) {
251          gradient[k] = 0.0;          gradient[k] = 0.0;
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. */
# Line 242  int BlackBoxCalcGradient(struct Instance Line 267  int BlackBoxCalcGradient(struct Instance
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    /*------------------------------------------------------------------------------
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  #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

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