dummy/solver/extrel.c

/*
 *  External Relations Cache for solvers.
 *  by Kirk Andre Abbott
 *  Created: 08/10/94
 *  Version: $Revision: 1.10 $
 *  Version control file: $RCSfile: extrel.c,v $
 *  Date last modified: $Date: 1997/07/18 12:14:14 $
 *  Last modified by: $Author: mthomas $
 *
 *  This file is part of the SLV solver.
 *
 *  Copyright (C) 1994 Kirk Abbott
 *
 *  The SLV solver is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public License as
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version.
 *
 *  The SLV solver is distributed in hope that it will be
 *  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along with
 *  the program; if not, write to the Free Software Foundation, Inc., 675
 *  Mass Ave, Cambridge, MA 02139 USA.  Check the file named COPYING.
 *  COPYING is found in ../compiler.
 */
#include <utilities/ascConfig.h>
#include <compiler/packages.h>
/* #include "solver/exprman.h" */
#include "rel.h"
#include "extrel.h"
#include <compiler/extcall.h>

double g_external_tolerance = 1.0e-12;

struct ExtRelCache *CreateExtRelCache(struct ExtCallNode *ext)
{
  struct ExtRelCache *result;
  unsigned long n_input_args, n_output_args;
  int ninputs, noutputs;

  assert(ext!=NULL);
  result = (struct ExtRelCache *)ascmalloc(sizeof(struct ExtRelCache));
  result->nodestamp = ExternalCallNodeStamp(ext);
  result->efunc = ExternalCallExtFunc(ext);
  result->data = ExternalCallDataInstance(ext);
  result->arglist = ExternalCallArgList(ext);
  result->user_data = NULL;     /* absolutely important to be
                     * initialized to NULL ! */

  n_input_args = NumberInputArgs(result->efunc);
  n_output_args = NumberOutputArgs(result->efunc);

  /*
   * We own the result of the LinearizeArgList call.
   */
  result->inputlist = LinearizeArgList(result->arglist,1,n_input_args);
  ninputs = (int)gl_length(result->inputlist);
  noutputs = (int)CountNumberOfArgs(result->arglist,n_input_args+1,
                    n_input_args+n_output_args);
  result->ninputs = ninputs;
  result->noutputs = noutputs;

  result->inputs = (double *)asccalloc(ninputs,sizeof(double));
  result->outputs = (double *)asccalloc(noutputs,sizeof(double));
  result->jacobian = (double *)asccalloc(ninputs*noutputs,sizeof(double));
  /*
   * Setup default flags for controlling calculations.
   */
  result->newcalc_done = (unsigned)1;
  return result;
}


struct ExtRelCache *CreateCacheFromInstance(struct Instance *relinst)
{
  struct ExtCallNode *ext;
  struct ExtRelCache *cache;
  CONST struct relation *reln;
  enum Expr_enum type;

  reln = GetInstanceRelation(relinst,&type);
  if (type!=e_blackbox) {
    FPRINTF(stderr,"Invalid relation type in (CreateCacheFromInstance)\n");
    return NULL;
  }
  ext = BlackBoxExtCall(reln);
  cache = CreateExtRelCache(ext);
  return cache;
}

void ExtRel_DestroyCache(struct ExtRelCache *cache)
{
  if (cache) {
    cache->nodestamp=-1;
    cache->efunc = NULL;
    cache->data = NULL;
    gl_destroy(cache->inputlist);   /* we own the list */
    ascfree(cache->inputs);
    ascfree(cache->outputs);
    ascfree(cache->jacobian);
    cache->inputlist = NULL;
    cache->inputs = NULL;
    cache->outputs = NULL;
    cache->jacobian = NULL;
    ascfree(cache);
  }
}


/*
 *********************************************************************
 * ExtRel_PreSolve:
 *
 * To deal with the first time we also want to do arguement
 * checking, and then turn off the first_func_eval flag.
 * Turn on the newcalc_done flag. The rationale behind this is
 * as follows:
 * The solver at the moment does not treat an external relation
 * specially, i.e., as a block. It also calls for its functions
 * a relation at a time. However the external relations compute
 * all their outputs at once. So as not to do unnecessary
 * recalculations we use these flag bits. We set newcalc_done
 * initially to true, so as to force *at least* one calculation
 * of the external relations. By similar reasoning first_func_eval (done)
 * is set to false.
 *********************************************************************
 */
int ExtRel_PreSolve(struct ExtRelCache *cache, int setup)
{
  struct ExternalFunc *efunc;
  struct Slv_Interp slv_interp;
  int ninputs,noutputs;
  double *inputs, *outputs, *jacobian;
  int (*init_func)(struct Slv_Interp *,
           struct Instance *,struct gl_list_t *);
  int nok = 0;

  if (!cache) return 1;
  efunc = cache->efunc;
  init_func = GetInitFunc(efunc);
  Init_Slv_Interp(&slv_interp);
  slv_interp.nodestamp = cache->nodestamp;
  slv_interp.user_data = cache->user_data;
  if (setup) {
    slv_interp.first_call = (unsigned)1;
    slv_interp.last_call = (unsigned)0;
    slv_interp.check_args = (unsigned)1;
  }
  else{
    slv_interp.first_call = (unsigned)0;
    slv_interp.last_call = (unsigned)1;
    slv_interp.check_args = (unsigned)0;
  }
  nok = (*init_func)(&slv_interp,cache->data,cache->arglist);
  if (nok)
    return 1;

  /*
   * Save the user's data and update our status.
   */
  cache->user_data = slv_interp.user_data;
  cache->newcalc_done = (unsigned)1;    /* force at least one calculation */
  cache->first_func_eval = (unsigned)0;
  return 0;
}


static int ArgsDifferent(double new, double old)
{
  if (fabs(new - old) >= g_external_tolerance)
    return 1;
  else
    return 0;
}

real64 ExtRel_Evaluate_RHS(struct rel_relation *rel)
{
  struct Slv_Interp slv_interp;
  struct ExtRelCache *cache;
  struct ExternalFunc *efunc;
  struct Instance *arg;
  struct gl_list_t *inputlist;
  double value;
  int c,ninputs;
  int nok;
  unsigned long whichvar;
  int newcalc_reqd=0;

  int (*eval_func)(struct Slv_Interp *,
           int ninputs, int noutputs,
           double *inputs, double *outputs, double *jacobian);

  assert(rel_extnodeinfo(rel));
  cache = rel_extcache(rel);
  efunc = cache->efunc;
  eval_func = GetValueFunc(efunc);
  inputlist = cache->inputlist;
  ninputs = cache->ninputs;
  whichvar = (unsigned long)rel_extwhichvar(rel);

  /*
   * The determination of whether a new calculation is required needs
   * some more thought. One thing we should insist upon is that all
   * the relations for an external relation are forced into the same
   * block.
   */
  for (c=0;c<ninputs;c++) {
    arg = (struct Instance *)gl_fetch(inputlist,c+1);
    value = RealAtomValue(arg);
    if (ArgsDifferent(value,cache->inputs[c])) {
      newcalc_reqd = 1;
      cache->inputs[c] = value;
    }
  }
  value = 0;
  /*
   * Do the calculations if necessary. Note that we have to *ensure*
   * that we send the user the information that he provided to us.
   * We have to update our user_data after each call of the user's function
   * as he might move information around (not smart but possible), on us.
   * If a function call is made, mark a new calculation as having been,
   * done, otherwise dont.
   */
  if (newcalc_reqd) {
    Init_Slv_Interp(&slv_interp);
    slv_interp.nodestamp = cache->nodestamp;
    slv_interp.user_data = cache->user_data;
    slv_interp.func_eval = (unsigned)1;

    nok = (*eval_func)(&slv_interp, ninputs, cache->noutputs,
               cache->inputs, cache->outputs, cache->jacobian);
    value = cache->outputs[whichvar - ninputs - 1];
    cache->newcalc_done = (unsigned)1;          /* newcalc done */
    cache->user_data = slv_interp.user_data;        /* update user_data */
  }
  else{
    value = cache->outputs[whichvar - ninputs - 1];
    cache->newcalc_done = (unsigned)0; /* a result was simply returned */
  }

#ifdef KAA_DEBUG
  FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_RHS result ->%g\n",
      result);
#endif

  return value;
}

    /******* FIX FIX FIX **********/
real64 ExtRel_Evaluate_LHS(rel)
struct rel_relation *rel;
{
  real64 res;
  res = 1.0;        /******* FIX FIX FIX **********/
  FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_LHS result ->%g\n",
      res);
  return res;
}


/*
 *********************************************************************
 * ExtRel_Gradient evaluation routines.
 *
 * The following code implements gradient evaluation routines for
 * external relations. The routines here will prepare the arguements
 * and call a user supplied derivative routine, is same is non-NULL.
 * If it is the user supplied function evaluation routine will be
 * used to compute the gradients via finite differencing.
 * The current solver will not necessarily call for the derivative
 * all at once. This makes it necessary to do the gradient
 * computations (user supplied or finite difference), and to cache
 * away the results. Based on calculation flags, the appropriate
 * *row* of this cached jacobian will be extracted and mapped to the
 * main solve matrix.
 *
 * The cached jacobian is a contiguous vector ninputs*noutputs long
 * and is loaded row wise. Indexing starts from 0. Each row corresponds
 * to the partial derivatives of the output variable (associated with
 * that row, wrt to all its incident input variables.
 *
 * Careful attention needs to be paid to the way this jacobian is
 * loaded/unloaded, because of the multiple indexing schemes in use.
 * i.e, arglist's and inputlists index 1..nvars, whereas all numeric
 * vectors number from 0.
 *
 *********************************************************************
 */

struct deriv_data {
  var_filter_t *filter;
  mtx_matrix_t mtx;
  mtx_coord_t nz;
};


/*
 *********************************************************************
 * ExtRel_MapDataToMtx
 *
 * This function attempts to map the information from the contiguous
 * jacobian back into the main matrix, based on the current row <=>
 * whichvar. The jacobian assumed to have been calculated.
 * Since we are operating a relation at a time, we have to find
 * out where to index into our jacobian. This index is computed as
 * follows:
 *
 * index = (whichvar - ninputs - 1) * ninputs
 *
 * Example: a problem with 4 inputs, 3 outputs and whichvar = 6.
 * with the counting for vars 1..nvars, but the jacobian indexing
 * starting from 0 (c-wise).
 *
 *          V-------- first output variable
 *  1 2 3 4 5 6 7
 *            ^---------------- whichvar
 *                  ------------------- grads for whichvar = 6
 *                  |    |    |    |
 *                  v    v    v    v
 *  index    =  0    1    2    3    4    5    6    7   8    9   10    11
 *  jacobian = 2.0  9.0  4.0  6.0  0.5  1.3  0.0  9.7  80  7.0  1.0  2.5
 *
 * Hence jacobian index = (6 - 4 - 1) * 4 = 4
 *********************************************************************
 */

static void ExtRel_MapDataToMtx(struct gl_list_t *inputlist,
                unsigned long whichvar,
                int ninputs,
                double *jacobian,
                struct deriv_data *d)
{
  struct Instance *inst;
  struct var_variable *var;
  double value, *ptr;
  boolean used;
  unsigned long c;
  int index;

  index = ((int)whichvar - ninputs - 1) * ninputs;
  ptr = &jacobian[index];

/* this is totally broken, thanks to kirk making the var=instance assumption */
  Asc_Panic(2, "ExtRel_MapDataToMtx",
            "ExtRel_MapDataToMtx is totally broken\n"
            "Makes a var=instance assumption");
  for (c=0;c<ninputs;c++) {
    inst = (struct Instance *)gl_fetch(inputlist,c+1);
    var = var_instance(inst);
    used = var_apply_filter(var,d->filter);
    if (used) {
      d->nz.col = mtx_org_to_col(d->mtx,var_sindex(var));
      value = ptr[c] + mtx_value(d->mtx,&(d->nz));
      mtx_set_value(d->mtx,&(d->nz), value);
    }
  }
}


/*
 *********************************************************************
 * ExtRel Finite Differencing.
 *
 * This routine actually does the finite differencing.
 * The jacobian is a single contiguous vector. We load information
 * in it *row* wise. If we have noutputs x ninputs = 3 x 4, variables,
 * then jacobian entry 4,
 * would correspond to jacobian[1][0], i.e., = 0.5 for this eg.
 *
 *    2.0  9.0   4.0  6.0
 *    0.5  1.3   0.0  9.7
 *    80   7.0   1.0  2.5
 *
 *  2.0  9.0  4.0  6.0  0.5  1.3  0.0  9.7  80  7.0  1.0  2.5
 * [0][0]              [1][0]             [2][1]
 *
 * When we are finite differencing variable c, we will be loading
 * jacobian positions c, c+ninputs, c+2*ninputs ....
 *********************************************************************
 */

static double CalculateInterval(double var_value)
{
  return (1.0e-05);
}

static int ExtRel__FDiff(struct Slv_Interp *slv_interp,
             int (*eval_func) (/* ARGS */),
             int ninputs, int noutputs,
             double *inputs, double *outputs,
             double *jacobian)
{
  int c1,c2, nok = 0;
  double *tmp_vector;
  double *ptr;
  double old_x,interval,value;

  tmp_vector = (double *)asccalloc(noutputs,sizeof(double));
  for (c1=0;c1<ninputs;c1++) {
    old_x = inputs[c1];                 /* perturb x */
    interval = CalculateInterval(old_x);
    inputs[c1] = old_x + interval;
    nok = (*eval_func)(slv_interp, ninputs, noutputs,   /* call routine */
               inputs, tmp_vector, jacobian);
    if (nok) break;
    ptr = &jacobian[c1];
    for (c2=0;c2<noutputs;c2++) {           /* load jacobian */
      value = (tmp_vector[c2] - outputs[c2])/interval;
      *ptr = value;
      ptr += ninputs;
    }
    inputs[c1] = old_x;
  }
  ascfree((char *)tmp_vector);              /* cleanup */
  return nok;
}


int ExtRel_CalcDeriv(struct rel_relation *rel, struct deriv_data *d)
{
  int nok = 0;
  struct Slv_Interp slv_interp;
  struct ExtRelCache *cache;
  struct ExternalFunc *efunc;
  unsigned long whichvar;
  double *deriv_vector;
  int (*eval_func)();
  int (*deriv_func)();

  assert(rel_extnodeinfo(rel));
  cache = rel_extcache(rel);
  whichvar = rel_extwhichvar(rel);
  efunc = cache->efunc;

  /*
   * Check and deal with the special case of the first
   * computation.
   */
  if (cache->first_deriv_eval) {
    cache->newcalc_done = (unsigned)1;
    cache->first_deriv_eval = (unsigned)0;
  }

  /*
   * If a function evaluation was not recently done, then we
   * can return the results from the cached jacobian.
   */
  if (!cache->newcalc_done) {
    ExtRel_MapDataToMtx(cache->inputlist, whichvar,
            cache->ninputs, cache->jacobian, d);
    return 0;
  }

  /*
   * If we are here, we have to do the derivative calculation.
   * The only thing to determine is whether we do analytical
   * derivatives (deriv_func != NULL) or finite differencing.
   * In any case init the interpreter.
   */
  Init_Slv_Interp(&slv_interp);
  slv_interp.deriv_eval = (unsigned)1;
  slv_interp.user_data = cache->user_data;
  deriv_func = GetDerivFunc(efunc);
  if (deriv_func) {
    nok = (*deriv_func)(&slv_interp, cache->ninputs, cache->noutputs,
            cache->inputs, cache->outputs, cache->jacobian);
    if (nok) return nok;
  }
  else{
    eval_func = GetValueFunc(efunc);
    nok = ExtRel__FDiff(&slv_interp, eval_func,
            cache->ninputs, cache->noutputs,
            cache->inputs, cache->outputs, cache->jacobian);
    if (nok) return nok;
  }

  /*
   * Cleanup. Ensure that we update the users data, and load
   * the main matrix with the derivative information.
   */
  cache->user_data = slv_interp.user_data;  /* save user info */
  ExtRel_MapDataToMtx(cache->inputlist, whichvar,
              cache->ninputs, cache->jacobian, d);
  return 0;
}


/*
 *********************************************************************
 * ExtRel Deriv routines.
 *
 * This is the entry point for most cases. ExtRel_CalcDeriv depends
 * on ExtRel_Evaluate being called  immediately before it.
 *
 *********************************************************************
 */

real64 ExtRel_Diffs_RHS(struct rel_relation *rel,
                  var_filter_t *filter,
                  int32 row,
                  mtx_matrix_t mtx)
{
  int nok;
  real64 res;
  struct deriv_data data;

  data.filter = filter;
  data.mtx = mtx;
  data.nz.row = row;

  res = ExtRel_Evaluate_RHS(rel);
  nok = ExtRel_CalcDeriv(rel,&data);
  if (nok)
    return 0.0;
  else
    return res;
}


    /******* FIX FIX FIX **********/

real64 ExtRel_Diffs_LHS(struct rel_relation *rel,
                  var_filter_t *filter,
                  int32 row,
                  mtx_matrix_t mtx)
{
  real64 res=0.0;
  return 1.0;   /******* FIX FIX FIX **********/
}


1	aw0a	1	/*
2			* External Relations Cache for solvers.
3			* by Kirk Andre Abbott
4			* Created: 08/10/94
5			* Version: $Revision: 1.10 $
6			* Version control file: $RCSfile: extrel.c,v $
7			* Date last modified: $Date: 1997/07/18 12:14:14 $
8			* Last modified by: $Author: mthomas $
9			*
10			* This file is part of the SLV solver.
11			*
12			* Copyright (C) 1994 Kirk Abbott
13			*
14			* The SLV solver is free software; you can redistribute
15			* it and/or modify it under the terms of the GNU General Public License as
16			* published by the Free Software Foundation; either version 2 of the
17			* License, or (at your option) any later version.
18			*
19			* The SLV solver is distributed in hope that it will be
20			* useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22			* General Public License for more details.
23			*
24			* You should have received a copy of the GNU General Public License along with
25			* the program; if not, write to the Free Software Foundation, Inc., 675
26			* Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
27			* COPYING is found in ../compiler.
28			*/
29	johnpye	399	#include <utilities/ascConfig.h>
30			#include <compiler/packages.h>
31	aw0a	1	/* #include "solver/exprman.h" */
32	johnpye	399	#include "rel.h"
33			#include "extrel.h"
34			#include <compiler/extcall.h>
35	aw0a	1
36			double g_external_tolerance = 1.0e-12;
37
38			struct ExtRelCache CreateExtRelCache(struct ExtCallNode ext)
39			{
40			struct ExtRelCache *result;
41			unsigned long n_input_args, n_output_args;
42			int ninputs, noutputs;
43
44			assert(ext!=NULL);
45			result = (struct ExtRelCache *)ascmalloc(sizeof(struct ExtRelCache));
46			result->nodestamp = ExternalCallNodeStamp(ext);
47			result->efunc = ExternalCallExtFunc(ext);
48			result->data = ExternalCallDataInstance(ext);
49			result->arglist = ExternalCallArgList(ext);
50			result->user_data = NULL; /* absolutely important to be
51			* initialized to NULL ! */
52
53			n_input_args = NumberInputArgs(result->efunc);
54			n_output_args = NumberOutputArgs(result->efunc);
55
56			/*
57			* We own the result of the LinearizeArgList call.
58			*/
59			result->inputlist = LinearizeArgList(result->arglist,1,n_input_args);
60			ninputs = (int)gl_length(result->inputlist);
61			noutputs = (int)CountNumberOfArgs(result->arglist,n_input_args+1,
62			n_input_args+n_output_args);
63			result->ninputs = ninputs;
64			result->noutputs = noutputs;
65
66			result->inputs = (double *)asccalloc(ninputs,sizeof(double));
67			result->outputs = (double *)asccalloc(noutputs,sizeof(double));
68			result->jacobian = (double )asccalloc(ninputsnoutputs,sizeof(double));
69			/*
70			* Setup default flags for controlling calculations.
71			*/
72			result->newcalc_done = (unsigned)1;
73			return result;
74			}
75
76
77			struct ExtRelCache CreateCacheFromInstance(struct Instance relinst)
78			{
79			struct ExtCallNode *ext;
80			struct ExtRelCache *cache;
81			CONST struct relation *reln;
82			enum Expr_enum type;
83
84			reln = GetInstanceRelation(relinst,&type);
85			if (type!=e_blackbox) {
86			FPRINTF(stderr,"Invalid relation type in (CreateCacheFromInstance)\n");
87			return NULL;
88			}
89			ext = BlackBoxExtCall(reln);
90			cache = CreateExtRelCache(ext);
91			return cache;
92			}
93
94			void ExtRel_DestroyCache(struct ExtRelCache *cache)
95			{
96			if (cache) {
97			cache->nodestamp=-1;
98			cache->efunc = NULL;
99			cache->data = NULL;
100			gl_destroy(cache->inputlist); /* we own the list */
101			ascfree(cache->inputs);
102			ascfree(cache->outputs);
103			ascfree(cache->jacobian);
104			cache->inputlist = NULL;
105			cache->inputs = NULL;
106			cache->outputs = NULL;
107			cache->jacobian = NULL;
108			ascfree(cache);
109			}
110			}
111
112
113			/*
114			*********************************************************************
115			* ExtRel_PreSolve:
116			*
117			* To deal with the first time we also want to do arguement
118			* checking, and then turn off the first_func_eval flag.
119			* Turn on the newcalc_done flag. The rationale behind this is
120			* as follows:
121			* The solver at the moment does not treat an external relation
122			* specially, i.e., as a block. It also calls for its functions
123			* a relation at a time. However the external relations compute
124			* all their outputs at once. So as not to do unnecessary
125			* recalculations we use these flag bits. We set newcalc_done
126			* initially to true, so as to force at least one calculation
127			* of the external relations. By similar reasoning first_func_eval (done)
128			* is set to false.
129			*********************************************************************
130			*/
131			int ExtRel_PreSolve(struct ExtRelCache *cache, int setup)
132			{
133			struct ExternalFunc *efunc;
134			struct Slv_Interp slv_interp;
135			int ninputs,noutputs;
136			double inputs, outputs, *jacobian;
137			int (init_func)(struct Slv_Interp ,
138			struct Instance ,struct gl_list_t );
139			int nok = 0;
140
141			if (!cache) return 1;
142			efunc = cache->efunc;
143			init_func = GetInitFunc(efunc);
144			Init_Slv_Interp(&slv_interp);
145			slv_interp.nodestamp = cache->nodestamp;
146			slv_interp.user_data = cache->user_data;
147			if (setup) {
148			slv_interp.first_call = (unsigned)1;
149			slv_interp.last_call = (unsigned)0;
150			slv_interp.check_args = (unsigned)1;
151			}
152			else{
153			slv_interp.first_call = (unsigned)0;
154			slv_interp.last_call = (unsigned)1;
155			slv_interp.check_args = (unsigned)0;
156			}
157			nok = (*init_func)(&slv_interp,cache->data,cache->arglist);
158			if (nok)
159			return 1;
160
161			/*
162			* Save the user's data and update our status.
163			*/
164			cache->user_data = slv_interp.user_data;
165			cache->newcalc_done = (unsigned)1; /* force at least one calculation */
166			cache->first_func_eval = (unsigned)0;
167			return 0;
168			}
169
170
171			static int ArgsDifferent(double new, double old)
172			{
173			if (fabs(new - old) >= g_external_tolerance)
174			return 1;
175			else
176			return 0;
177			}
178
179			real64 ExtRel_Evaluate_RHS(struct rel_relation *rel)
180			{
181			struct Slv_Interp slv_interp;
182			struct ExtRelCache *cache;
183			struct ExternalFunc *efunc;
184			struct Instance *arg;
185			struct gl_list_t *inputlist;
186			double value;
187			int c,ninputs;
188			int nok;
189			unsigned long whichvar;
190			int newcalc_reqd=0;
191
192			int (eval_func)(struct Slv_Interp ,
193			int ninputs, int noutputs,
194			double inputs, double outputs, double *jacobian);
195
196			assert(rel_extnodeinfo(rel));
197			cache = rel_extcache(rel);
198			efunc = cache->efunc;
199			eval_func = GetValueFunc(efunc);
200			inputlist = cache->inputlist;
201			ninputs = cache->ninputs;
202			whichvar = (unsigned long)rel_extwhichvar(rel);
203
204			/*
205			* The determination of whether a new calculation is required needs
206			* some more thought. One thing we should insist upon is that all
207			* the relations for an external relation are forced into the same
208			* block.
209			*/
210			for (c=0;c<ninputs;c++) {
211			arg = (struct Instance *)gl_fetch(inputlist,c+1);
212			value = RealAtomValue(arg);
213			if (ArgsDifferent(value,cache->inputs[c])) {
214			newcalc_reqd = 1;
215			cache->inputs[c] = value;
216			}
217			}
218			value = 0;
219			/*
220			* Do the calculations if necessary. Note that we have to ensure
221			* that we send the user the information that he provided to us.
222			* We have to update our user_data after each call of the user's function
223			* as he might move information around (not smart but possible), on us.
224			* If a function call is made, mark a new calculation as having been,
225			* done, otherwise dont.
226			*/
227			if (newcalc_reqd) {
228			Init_Slv_Interp(&slv_interp);
229			slv_interp.nodestamp = cache->nodestamp;
230			slv_interp.user_data = cache->user_data;
231			slv_interp.func_eval = (unsigned)1;
232
233			nok = (*eval_func)(&slv_interp, ninputs, cache->noutputs,
234			cache->inputs, cache->outputs, cache->jacobian);
235			value = cache->outputs[whichvar - ninputs - 1];
236			cache->newcalc_done = (unsigned)1; /* newcalc done */
237			cache->user_data = slv_interp.user_data; /* update user_data */
238			}
239			else{
240			value = cache->outputs[whichvar - ninputs - 1];
241			cache->newcalc_done = (unsigned)0; /* a result was simply returned */
242			}
243
244			#ifdef KAA_DEBUG
245			FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_RHS result ->%g\n",
246			result);
247			#endif
248
249			return value;
250			}
251
252			/***** FIX FIX FIX ********/
253			real64 ExtRel_Evaluate_LHS(rel)
254			struct rel_relation *rel;
255			{
256			real64 res;
257			res = 1.0; /***** FIX FIX FIX ********/
258			FPRINTF(stderr,"Finsished calling ExtRel_Evaluate_LHS result ->%g\n",
259			res);
260			return res;
261			}
262
263
264			/*
265			*********************************************************************
266			* ExtRel_Gradient evaluation routines.
267			*
268			* The following code implements gradient evaluation routines for
269			* external relations. The routines here will prepare the arguements
270			* and call a user supplied derivative routine, is same is non-NULL.
271			* If it is the user supplied function evaluation routine will be
272			* used to compute the gradients via finite differencing.
273			* The current solver will not necessarily call for the derivative
274			* all at once. This makes it necessary to do the gradient
275			* computations (user supplied or finite difference), and to cache
276			* away the results. Based on calculation flags, the appropriate
277			* row of this cached jacobian will be extracted and mapped to the
278			* main solve matrix.
279			*
280			* The cached jacobian is a contiguous vector ninputs*noutputs long
281			* and is loaded row wise. Indexing starts from 0. Each row corresponds
282			* to the partial derivatives of the output variable (associated with
283			* that row, wrt to all its incident input variables.
284			*
285			* Careful attention needs to be paid to the way this jacobian is
286			* loaded/unloaded, because of the multiple indexing schemes in use.
287			* i.e, arglist's and inputlists index 1..nvars, whereas all numeric
288			* vectors number from 0.
289			*
290			*********************************************************************
291			*/
292
293			struct deriv_data {
294			var_filter_t *filter;
295			mtx_matrix_t mtx;
296			mtx_coord_t nz;
297			};
298
299
300			/*
301			*********************************************************************
302			* ExtRel_MapDataToMtx
303			*
304			* This function attempts to map the information from the contiguous
305			* jacobian back into the main matrix, based on the current row <=>
306			* whichvar. The jacobian assumed to have been calculated.
307			* Since we are operating a relation at a time, we have to find
308			* out where to index into our jacobian. This index is computed as
309			* follows:
310			*
311			* index = (whichvar - ninputs - 1) * ninputs
312			*
313			* Example: a problem with 4 inputs, 3 outputs and whichvar = 6.
314			* with the counting for vars 1..nvars, but the jacobian indexing
315			* starting from 0 (c-wise).
316			*
317			* V-------- first output variable
318			* 1 2 3 4 5 6 7
319			* ^---------------- whichvar
320			* ------------------- grads for whichvar = 6
321			* \| \| \| \|
322			* v v v v
323			* index = 0 1 2 3 4 5 6 7 8 9 10 11
324			* jacobian = 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
325			*
326			* Hence jacobian index = (6 - 4 - 1) * 4 = 4
327			*********************************************************************
328			*/
329
330			static void ExtRel_MapDataToMtx(struct gl_list_t *inputlist,
331			unsigned long whichvar,
332			int ninputs,
333			double *jacobian,
334			struct deriv_data *d)
335			{
336			struct Instance *inst;
337			struct var_variable *var;
338			double value, *ptr;
339			boolean used;
340			unsigned long c;
341			int index;
342
343			index = ((int)whichvar - ninputs - 1) * ninputs;
344			ptr = &jacobian[index];
345
346			/* this is totally broken, thanks to kirk making the var=instance assumption */
347			Asc_Panic(2, "ExtRel_MapDataToMtx",
348			"ExtRel_MapDataToMtx is totally broken\n"
349			"Makes a var=instance assumption");
350			for (c=0;c<ninputs;c++) {
351			inst = (struct Instance *)gl_fetch(inputlist,c+1);
352			var = var_instance(inst);
353			used = var_apply_filter(var,d->filter);
354			if (used) {
355			d->nz.col = mtx_org_to_col(d->mtx,var_sindex(var));
356			value = ptr[c] + mtx_value(d->mtx,&(d->nz));
357			mtx_set_value(d->mtx,&(d->nz), value);
358			}
359			}
360			}
361
362
363			/*
364			*********************************************************************
365			* ExtRel Finite Differencing.
366			*
367			* This routine actually does the finite differencing.
368			* The jacobian is a single contiguous vector. We load information
369			* in it row wise. If we have noutputs x ninputs = 3 x 4, variables,
370			* then jacobian entry 4,
371			* would correspond to jacobian[1][0], i.e., = 0.5 for this eg.
372			*
373			* 2.0 9.0 4.0 6.0
374			* 0.5 1.3 0.0 9.7
375			* 80 7.0 1.0 2.5
376			*
377			* 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
378			* [0][0] [1][0] [2][1]
379			*
380			* When we are finite differencing variable c, we will be loading
381			* jacobian positions c, c+ninputs, c+2*ninputs ....
382			*********************************************************************
383			*/
384
385			static double CalculateInterval(double var_value)
386			{
387			return (1.0e-05);
388			}
389
390			static int ExtRel__FDiff(struct Slv_Interp *slv_interp,
391			int (eval_func) (/ ARGS */),
392			int ninputs, int noutputs,
393			double inputs, double outputs,
394			double *jacobian)
395			{
396			int c1,c2, nok = 0;
397			double *tmp_vector;
398			double *ptr;
399			double old_x,interval,value;
400
401			tmp_vector = (double *)asccalloc(noutputs,sizeof(double));
402			for (c1=0;c1<ninputs;c1++) {
403			old_x = inputs[c1]; /* perturb x */
404			interval = CalculateInterval(old_x);
405			inputs[c1] = old_x + interval;
406			nok = (eval_func)(slv_interp, ninputs, noutputs, / call routine */
407			inputs, tmp_vector, jacobian);
408			if (nok) break;
409			ptr = &jacobian[c1];
410			for (c2=0;c2<noutputs;c2++) { /* load jacobian */
411			value = (tmp_vector[c2] - outputs[c2])/interval;
412			*ptr = value;
413			ptr += ninputs;
414			}
415			inputs[c1] = old_x;
416			}
417			ascfree((char )tmp_vector); / cleanup */
418			return nok;
419			}
420
421
422			int ExtRel_CalcDeriv(struct rel_relation rel, struct deriv_data d)
423			{
424			int nok = 0;
425			struct Slv_Interp slv_interp;
426			struct ExtRelCache *cache;
427			struct ExternalFunc *efunc;
428			unsigned long whichvar;
429			double *deriv_vector;
430			int (*eval_func)();
431			int (*deriv_func)();
432
433			assert(rel_extnodeinfo(rel));
434			cache = rel_extcache(rel);
435			whichvar = rel_extwhichvar(rel);
436			efunc = cache->efunc;
437
438			/*
439			* Check and deal with the special case of the first
440			* computation.
441			*/
442			if (cache->first_deriv_eval) {
443			cache->newcalc_done = (unsigned)1;
444			cache->first_deriv_eval = (unsigned)0;
445			}
446
447			/*
448			* If a function evaluation was not recently done, then we
449			* can return the results from the cached jacobian.
450			*/
451			if (!cache->newcalc_done) {
452			ExtRel_MapDataToMtx(cache->inputlist, whichvar,
453			cache->ninputs, cache->jacobian, d);
454			return 0;
455			}
456
457			/*
458			* If we are here, we have to do the derivative calculation.
459			* The only thing to determine is whether we do analytical
460			* derivatives (deriv_func != NULL) or finite differencing.
461			* In any case init the interpreter.
462			*/
463			Init_Slv_Interp(&slv_interp);
464			slv_interp.deriv_eval = (unsigned)1;
465			slv_interp.user_data = cache->user_data;
466			deriv_func = GetDerivFunc(efunc);
467			if (deriv_func) {
468			nok = (*deriv_func)(&slv_interp, cache->ninputs, cache->noutputs,
469			cache->inputs, cache->outputs, cache->jacobian);
470			if (nok) return nok;
471			}
472			else{
473			eval_func = GetValueFunc(efunc);
474			nok = ExtRel__FDiff(&slv_interp, eval_func,
475			cache->ninputs, cache->noutputs,
476			cache->inputs, cache->outputs, cache->jacobian);
477			if (nok) return nok;
478			}
479
480			/*
481			* Cleanup. Ensure that we update the users data, and load
482			* the main matrix with the derivative information.
483			*/
484			cache->user_data = slv_interp.user_data; /* save user info */
485			ExtRel_MapDataToMtx(cache->inputlist, whichvar,
486			cache->ninputs, cache->jacobian, d);
487			return 0;
488			}
489
490
491			/*
492			*********************************************************************
493			* ExtRel Deriv routines.
494			*
495			* This is the entry point for most cases. ExtRel_CalcDeriv depends
496			* on ExtRel_Evaluate being called immediately before it.
497			*
498			*********************************************************************
499			*/
500
501			real64 ExtRel_Diffs_RHS(struct rel_relation *rel,
502			var_filter_t *filter,
503			int32 row,
504			mtx_matrix_t mtx)
505			{
506			int nok;
507			real64 res;
508			struct deriv_data data;
509
510			data.filter = filter;
511			data.mtx = mtx;
512			data.nz.row = row;
513
514			res = ExtRel_Evaluate_RHS(rel);
515			nok = ExtRel_CalcDeriv(rel,&data);
516			if (nok)
517			return 0.0;
518			else
519			return res;
520			}
521
522
523
524			/***** FIX FIX FIX ********/
525
526			real64 ExtRel_Diffs_LHS(struct rel_relation *rel,
527			var_filter_t *filter,
528			int32 row,
529			mtx_matrix_t mtx)
530			{
531			real64 res=0.0;
532			return 1.0; /***** FIX FIX FIX ********/
533			}
534
535
536
537
538
539
540
541
542
543