generic/compiler/rel_blackbox.c


#include <math.h>
#include <errno.h>
#include <stdarg.h>
#include <utilities/ascConfig.h>
#include <utilities/ascMalloc.h>
#include <utilities/ascPanic.h>
#include <general/list.h>
#include <general/dstring.h>
#include "compiler.h"
#include "symtab.h"
#include "functype.h"
#include "safe.h"
#include "fractions.h"
#include "dimen.h"
#include "expr_types.h"
#include "vlist.h"
#include "dimen_io.h"
#include "instance_enum.h"
#include "find.h"
#include "atomvalue.h"
#include "instance_name.h"
#include "extfunc.h"
#include "relation_type.h"
#include "rel_blackbox.h"
#include "relation.h"
#include "relation_util.h"
#include "instance_io.h"
#include "instquery.h"
#include "visitinst.h"
#include "mathinst.h"
#include "extcall.h" /* for copy/destroy speciallist */
#include "packages.h" /* for init slv interp */

static int32 ArgsDifferent(double new, double old, double tol)
{
    if (fabs(new - old) > fabs(tol)) {
        return 1;
    } else {
        return 0;
    }
}

/*------------------------------------------------------------------------------
  RESIDUAL AND GRADIENT EVALUATION ROUTINES
*/

/* 
    Note:
    Assume a function in c/fortran that computes yhat(x),
    aka func(inputs_x, outputs_yhat);.
    This is fit into the ascend world by associating
    inputs_x directly to ascend variables and defining
    the ascend variables y in the outer model.
    The residual of the blackbox equations in the ascend model
    is y-yhat. The gradient is I - dyhat/dx, where dyhat/dx
    is the reduced jacobian of the blackbox.
*/
int BlackBoxCalcResidGrad(struct Instance *i, double *res, double *gradient, struct relation *r)
{
    int residErr = 0;
    int gradErr = 0;

    residErr = BlackBoxCalcResidual(i, res, r);
    gradErr = BlackBoxCalcGradient(i, gradient, r);
    return (int)(fabs(residErr) + fabs(gradErr));
}

/* 
    Note:
    Assume a function in c/fortran that computes yhat(x),
    aka func(inputs_x, outputs_yhat);.
    This is fit into the ascend world by associating
    inputs_x directly to ascend variables and defining
    the ascend variables y in the outer model.
    The residual of the blackbox equations in the ascend model
    is y-yhat. The gradient is I - dyhat/dx, where dyhat/dx
    is the reduced jacobian of the blackbox.
*/
int BlackBoxCalcResidual(struct Instance *i, double *res, struct relation *r)
{
/* decls */
    unsigned long *argToVar;
    unsigned long c;
    unsigned long argToVarLen;
    struct BlackBoxData *unique;
    struct BlackBoxCache *common;
    struct ExternalFunc *efunc;
    ExtBBoxFunc *evalFunc;
    int lhsVar;
    int outputIndex;
    struct Instance *arg;
    double value;
    double inputTolerance;
    int updateNeeded;
    int nok = 0;
    static int warnexpt;

    if(!warnexpt){
        ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox evaluation is experimental (%s)",__FUNCTION__);
        warnexpt = 1;
    }

    (void)i;
    efunc = RelationBlackBoxExtFunc(r);
    unique = RelationBlackBoxData(r);
    common = RelationBlackBoxCache(r);
    argToVar = unique->inputArgs;
    argToVarLen = common->inputsLen;
    lhsVar = unique->lhsvar;
    outputIndex = unique->lhsindex;
    evalFunc = GetValueFunc(efunc);
    inputTolerance = GetValueFuncTolerance(efunc);
    updateNeeded = 0;

/* impl:
    check input values changed.
    if changed, recompute output values.
    compute y - yhat for lhsvar.
*/
    /* check: */
    if (common->residCount < 1) {
        updateNeeded = 1;
    } else {
        for (c=0; c < argToVarLen ; c++) {
            arg = RelationVariable(r,argToVar[c]); 
            value = RealAtomValue(arg);
            if (ArgsDifferent(value, common->inputs[c], inputTolerance)) {
                updateNeeded = 1;
                break;
            }
        }
    }
    /* recompute */
    if (updateNeeded) {
        for (c=0; c < argToVarLen ;c++) {
            arg = RelationVariable(r,argToVar[c]);
            value = RealAtomValue(arg);
            common->inputs[c] = value;
        }
        common->interp.task = bb_func_eval;

        nok = (*evalFunc)(&(common->interp),
                common->inputsLen,
                common->outputsLen,
                common->inputs,
                common->outputs,
                common->jacobian);
        common->residCount++;
    }
    value = common->outputs[outputIndex];

    /* return lhsval - value */
    arg = RelationVariable(r,lhsVar);
    *res = RealAtomValue(arg) - value;
    return nok;

}

/**
    Calculate the gradient (slice of the overall jacobian) for the blackbox.

    The tricky bit in this is that if input or output args are merged,
    their partial derivatives get summed and the solver at least gets
    what it deserves, which may or may not be sane in a modeling sense.

    Implementation:
        # check input values changed.
        # if changed, recompute gradient in bbox.
        # compute gradient per varlist from bbox row.
*/

int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
    , int ninputs, int noutputs
    , double *inputs, double *outputs, double *jac
);

int BlackBoxCalcGradient(struct Instance *i, double *gradient
        , struct relation *r
){
/* decls */
    unsigned long *argToVar;
    unsigned long c, varlistLen;
    unsigned long argToVarLen;
    struct BlackBoxData *unique;
    struct BlackBoxCache *common;
    struct ExternalFunc *efunc;
    ExtBBoxFunc *derivFunc;
    int lhsVar;
    int outputIndex;
    struct Instance *arg;
    double value;
    double inputTolerance;
    double *bboxRow;
    int updateNeeded, offset;
    unsigned int k;
    int nok = 0;
    static int warnexpt, warnfdiff;

    if(!warnexpt){ 
        ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Blackbox gradient is experimental (%s)",__FUNCTION__);
        warnexpt = 1;
    }

    /* prepare */
    (void)i;
    efunc = RelationBlackBoxExtFunc(r);
    unique = RelationBlackBoxData(r);
    common = RelationBlackBoxCache(r);
    argToVar = unique->inputArgs;
    argToVarLen = common->inputsLen;
    lhsVar = unique->lhsvar;
    outputIndex = unique->lhsindex;
    derivFunc = GetDerivFunc(efunc);
    inputTolerance = GetValueFuncTolerance(efunc);
    updateNeeded = 0;
    varlistLen = NumberVariables(r);

    /* do we need to calculate? */
    if(common->gradCount < 1){
        updateNeeded = 1;
    }else{
        for(c=0; c<argToVarLen; c++){
            arg = RelationVariable(r,argToVar[c]);
            value = RealAtomValue(arg);
            if (ArgsDifferent(value, common->inputsJac[c], inputTolerance)) {
                updateNeeded = 1;
                break;
            }
        }
    }
    
    /* if inputs have changed more than allowed, or it's first time: evaluate */
    if (updateNeeded) {
        for (c=0; c < argToVarLen ;c++) {
            arg = RelationVariable(r,argToVar[c]);
            value = RealAtomValue(arg);
            common->inputsJac[c] = value;
        }
        common->interp.task = bb_deriv_eval;

        if(derivFunc){          
            nok = (*derivFunc)(
                &(common->interp)
                , common->inputsLen, common->outputsLen
                , common->inputsJac, common->outputs, common->jacobian
            );
        }else{
            if(!warnfdiff){
                ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Using finite-difference "
                    " to compute derivatives for one or more black boxes in"
                    " this model."
                );
                warnfdiff = 1;
            }
            
            nok = blackbox_fdiff(GetValueFunc(efunc), &(common->interp)
                , common->inputsLen, common->outputsLen
                , common->inputsJac, common->outputs, common->jacobian
            );
        }
        common->gradCount++;
    }

    for (k = 0; k < varlistLen; k++) {
        gradient[k] = 0.0;
    }

    /* now compute d(y-yhat)/dx for this row as ( I - dyhat/dx ) */
    k = lhsVar-1;
    gradient[k] = 1.0; /* I */
    offset = common->inputsLen * outputIndex; /* offset to row needed. */
    bboxRow = common->jacobian + offset; /* pointer to row by ptr arithmetic */
    for (c=0; c < argToVarLen; c++) {
        k = ((int)(argToVar[c])) - 1; /* varnum is 1..n, not 0 .*/
        gradient[k] -= bboxRow[c];
    }

    return nok;
}

struct Instance *BlackBoxGetOutputVar(struct relation *r)
{
    assert(r != NULL);
    /* unsigned long lhsVarNumber; */
    /** @TODO FIXME BlackBoxGetOutputVar */ 
    return NULL;
}

struct BlackBoxData *CreateBlackBoxData(struct BlackBoxCache *common,
                    unsigned long lhsindex,
                    unsigned long lhsVarNumber)
{
    struct BlackBoxData *b = (struct BlackBoxData *)malloc(sizeof(struct BlackBoxData));
    assert(common!=NULL);
    b->common = common;
    AddRefBlackBoxCache(common);
    b->inputArgs = (unsigned long *)ascmalloc(sizeof(long) * (common->inputsLen));
    b->lhsindex = lhsindex;
    b->lhsvar = lhsVarNumber;
    return b;
}

void DestroyBlackBoxData(struct relation *rel, struct BlackBoxData *b)
{
    b->lhsindex = -(b->lhsindex);
    b->lhsvar = 0;
    ascfree(b->inputArgs);
    b->inputArgs = NULL;
    DeleteRefBlackBoxCache(rel, &(b->common));
    ascfree(b);
}

/*------------------------------------------------------------------------------
  BLACKBOX GRADIENT BY FINITE-DIFFERENCE
*/

/**
    This function works out what the peturbed value for a variable should be.

    @TODO add some awareness of boundaries in this, hopefully.
*/
static inline double blackbox_peturbation(double varvalue){
  return (1.0e-05);
}

/**
    Blackbox derivatives estimated by finite difference (by evaluation at 
    peturbed value of each input in turn)

    Call signature as for ExtBBoxFunc.
*/
int blackbox_fdiff(ExtBBoxFunc *resfn, struct BBoxInterp *interp
    , int ninputs, int noutputs
    , double *inputs, double *outputs, double *jac
){
  long c,r;
  int nok = 0;
  double *tmp_outputs;
  double *ptr;
  double old_x,deltax,value;

  /* CONSOLE_DEBUG("NUMERICAL DERIVATIVE..."); */

  tmp_outputs = ASC_NEW_ARRAY_CLEAR(double,noutputs);

  for (c=0;c<ninputs;c++){
    /* perturb each input in turn */
    old_x = inputs[c];
    deltax = blackbox_peturbation(old_x);
    inputs[c] = old_x + deltax;
    /* CONSOLE_DEBUG("PETURBATED VALUE of input[%ld] = %f",c,inputs[c]); */

    /* call routine. note that the 'jac' parameter is just along for the ride */
    nok = (*resfn)(interp, ninputs, noutputs, inputs, tmp_outputs, jac);
    if(nok){
        CONSOLE_DEBUG("External evaluation error (%d)",nok);
        break;
    }

    /* fill load jacobian */
    ptr = &jac[c];
    for(r=0;r<noutputs;r++){
      value = (tmp_outputs[r] - outputs[r]) / deltax;
      /* CONSOLE_DEBUG("output[%ld]: value = %f, gradient = %f",r,tmp_outputs[r],value); */
      *ptr = value;
      ptr += ninputs;
    }

    /* now return this input to its old value */
    inputs[c] = old_x;
  }
  ASC_FREE(tmp_outputs);
  if(nok){
    CONSOLE_DEBUG("External evaluation error");
  }
  return nok;

}

/*------------------------------------------------------------------------------
  BLACK BOX CACHE
*/

#define JACMAGIC -3.141592071828
struct BlackBoxCache *CreateBlackBoxCache(
    int32 inputsLen,
    int32 outputsLen,
    struct gl_list_t *formalArgs
)
{
    struct BlackBoxCache *b = (struct BlackBoxCache *)malloc(sizeof(struct BlackBoxCache));
    b->interp.task = bb_none;
    b->interp.status = calc_all_ok;
    b->interp.user_data = NULL;
    b->formalArgList = CopySpecialList(formalArgs);
    b->inputsLen = inputsLen;
    b->outputsLen = outputsLen;
    b->jacobianLen = inputsLen*outputsLen;
    b->hessianLen = 0; /* FIXME. */
    b->inputs = (double *)ascmalloc(inputsLen*sizeof(double));
    b->inputsJac = (double *)ascmalloc(inputsLen*sizeof(double));
    b->outputs = (double *)ascmalloc(outputsLen*sizeof(double));
    b->jacobian = (double*)ascmalloc(outputsLen*inputsLen*sizeof(double)+sizeof(double));
    b->jacobian[outputsLen*inputsLen] = JACMAGIC;
    b->hessian = NULL;
    b->residCount = 0;
    b->gradCount = 0;
    b->refCount = 1;
    return b;
}

void AddRefBlackBoxCache(struct BlackBoxCache *b)
{
    assert(b != NULL);
    (b->refCount)++;
}

static void DestroyBlackBoxCache(struct relation *rel, struct BlackBoxCache *b)
{
    struct ExternalFunc *efunc;
    ExtBBoxFinalFunc *final;
    assert(b != NULL);
    if (b->jacobian[b->outputsLen*b->inputsLen] != JACMAGIC)
    {
        FPRINTF(ASCERR, "Jacobian overrun detected while destroying BlackBoxCache. Debug the user blackbox implementation.\n");
    }
    ascfree(b->inputs);
    ascfree(b->inputsJac);
    ascfree(b->outputs);
    ascfree(b->jacobian);
    DestroySpecialList(b->formalArgList);
    b->formalArgList = NULL;
    b->inputsLen = -(b->inputsLen);
    b->outputsLen = -(b->outputsLen);
    b->jacobianLen = -(b->jacobianLen);
    b->hessianLen = -(b->hessianLen);
    b->inputs = NULL;
    b->inputsJac = NULL;
    b->outputs = NULL;
    b->jacobian = NULL;
    b->hessian = NULL;
    b->residCount = -(b->residCount);
    b->gradCount = -(b->gradCount);
    b->refCount = -3;
    if (rel != NULL) {
        /* rel is null with refcount 0 only when there's a bbox array
        instantiation failure. If failure, we don't need final. */ 
        efunc = RelationBlackBoxExtFunc(rel);
        final = GetFinalFunc(efunc);
        b->interp.task = bb_last_call;
        (*final)(&(b->interp));
    }
    ascfree(b);
}

void DeleteRefBlackBoxCache(struct relation *rel, struct BlackBoxCache **b)
{
    struct BlackBoxCache * d = *b;
    assert(b != NULL && *b != NULL);
    *b = NULL;
    if (d->refCount > 0) {
        (d->refCount)--;
    }
    if (d->refCount == 0) {
        DestroyBlackBoxCache(rel,d);
        return;
    }
    if (d->refCount < 0) {
        FPRINTF(ASCERR, "Attempt to delete BlackBoxCache (%p) too often.\n", *b);
    }
}
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(inputsLensizeof(double));
402	b->inputsJac = (double )ascmalloc(inputsLensizeof(double));
403	b->outputs = (double )ascmalloc(outputsLensizeof(double));
404	b->jacobian = (double)ascmalloc(outputsLeninputsLen*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	}