models/sensitivity/sensitivity.c

#include <math.h>

#include "sensitivity.h"
#include <packages/sensitivity.h>
#include <compiler/instquery.h>
#include <compiler/atomvalue.h>
#include <utilities/ascMalloc.h>


/**
    Build then presolve an instance
*/
slv_system_t asc_sens_presolve(struct Instance *inst){
  slv_system_t sys;
  slv_parameters_t parameters;
  struct var_variable **vp;
  struct rel_relation **rp;
  int ind,len;
  char *tmp=NULL;

  sys = system_build(inst);
  if (sys==NULL) {
    ERROR_REPORTER_HERE(ASC_PROG_ERR,
      "Something radically wrong in creating solver.");
    return NULL;
  }
  if (g_SlvNumberOfRegisteredClients == 0) {
    return NULL;
  }
  ind = 0;
  while (strcmp(slv_solver_name(ind),"QRSlv")) {
    if (ind >= g_SlvNumberOfRegisteredClients) {
      ERROR_REPORTER_HERE(ASC_PROG_ERR,
        "QRSlv must be registered client.");
      return NULL;
    }
    ++ind;
  }
  slv_select_solver(sys,ind);

  slv_get_parameters(sys,&parameters);
  parameters.partition = 0;
  slv_set_parameters(sys,&parameters);
  slv_presolve(sys);

#if DEBUG
  vp = slv_get_solvers_var_list(sys);
  len = slv_get_num_solvers_vars(sys);
  for (ind=0 ; ind<len; ind++) {
    tmp = var_make_name(sys,vp[ind]);
    FPRINTF(stderr,"%s  %d\n",tmp,var_sindex(vp[ind]));
    if (tmp!=NULL) ascfree(tmp);
  }
  rp = slv_get_solvers_rel_list(sys);
  len = slv_get_num_solvers_rels(sys);
  for (ind=0 ; ind<len ; ind++) {
    tmp = rel_make_name(sys,rp[ind]);
    FPRINTF(stderr,"%s  %d\n",tmp,rel_sindex(rp[ind]));
    if (tmp) ascfree(tmp);
  }
#endif
  return sys;
}


int sensitivity_anal(
             struct Instance *inst, /* not used but will be */
             struct gl_list_t *arglist){
  struct Instance *which_instance,*tmp_inst, *atominst;
  struct gl_list_t *branch;
  struct var_variable **vlist = NULL;
  int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
  real64 **dy_dx = NULL;
  slv_system_t sys = NULL;
  int c;
  int noutputs = 0;
  int ninputs;
  int i,j;
  int offset;
  dof_t *dof;
  int num_vars,ind,found;

  linsolqr_system_t lqr_sys;    /* stuff for the linear system & matrix */
  mtx_matrix_t mtx;
  int32 capacity;
  real64 *scratch_vector = NULL;
  int result=0;

  /* Ignore unused params */
  (void) inst;

  (void)NumberFreeVars(NULL);       /* used to re-init the system */
  (void)NumberRels(NULL);       /* used to re-init the system */
  which_instance = FetchElement(arglist,1,1);
  sys = asc_sens_presolve(which_instance);
  if (!sys) {
    FPRINTF(stderr,"Early termination due to failure in Presolve\n");
    result = 1;
    goto error;
  }

  dof = slv_get_dofdata(sys);
  if (!(dof->n_rows == dof->n_cols &&
    dof->n_rows == dof->structural_rank)) {
    FPRINTF(stderr,"Early termination: non square system\n");
    result = 1;
    goto error;
  }
  /*
   * prepare the inputs list
   */
  vlist = slv_get_solvers_var_list(sys);

  branch = gl_fetch(arglist,2); /* input args */
  ninputs = gl_length(branch);
  inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);

  num_vars = slv_get_num_solvers_vars(sys);
  for (c=0;c<ninputs;c++) {
    atominst = (struct Instance *)gl_fetch(branch,c+1);
    found = 0;
    ind = num_vars - 1;
    /* search backwards because fixed vars are at the end of the
       var list */
    while (!found && ind >= 0){
      if (var_instance(vlist[ind]) == atominst) {
    inputs_ndx_list[c] = var_sindex(vlist[ind]);
    found = 1;
      }
      --ind;
    }
    if (!found) {
      FPRINTF(stderr,"Sensitivity input variable not found\n");
      result = 1;
      goto error;
    }
  }

  /*
   * prepare the outputs list
   */
  branch = gl_fetch(arglist,3); /* output args */
  noutputs = gl_length(branch);
  outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
  for (c=0;c<noutputs;c++) {
    atominst = (struct Instance *)gl_fetch(branch,c+1);
    found = 0;
    ind = 0;
    while (!found && ind < num_vars){
      if (var_instance(vlist[ind]) == atominst) {
    outputs_ndx_list[c] = var_sindex(vlist[ind]);
    found = 1;
      }
      ++ind;
    }
    if (!found) {
      FPRINTF(stderr,"Sensitivity ouput variable not found\n");
      result = 1;
      goto error;
    }
  }

  /*
   * prepare the results dy_dx.
   */
  dy_dx = make_matrix(noutputs,ninputs);


  result = Compute_J(sys);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
    goto error;
  }

  /*
   * Note: the RHS *has* to added here. We will construct the vector
   * of size matrix capacity and add it. It will be removed after
   * we finish computing dy_dx.
   */
  lqr_sys = slv_get_linsolqr_sys(sys);  /* get the linear system */
  mtx = linsolqr_get_matrix(lqr_sys);   /* get the matrix */
  capacity = mtx_capacity(mtx);
  scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
  linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
  result = LUFactorJacobian1(sys,dof->structural_rank);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
    goto error;
  }
  result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
                   inputs_ndx_list, ninputs,
                   outputs_ndx_list, noutputs);

  linsolqr_remove_rhs(lqr_sys,scratch_vector);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
    goto error;
  }

  /*
   * Write the results back to the partials array in the
   * instance tree
   */
  offset = 0;
  for (i=0;i<noutputs;i++) {
    for (j=0;j<ninputs;j++) {
      tmp_inst = FetchElement(arglist,4,offset+j+1);
      SetRealAtomValue(tmp_inst,dy_dx[i][j],(unsigned)0);
#if DEBUG
      FPRINTF(stderr,"%12.8f   i%d j%d",dy_dx[i][j],i,j);
#endif
    }
#if DEBUG
    FPRINTF(stderr,"\n");
#endif
    offset += ninputs;
  }
#if DEBUG
  FPRINTF(stderr,"\n");
#endif

error:
  if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
  if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
  if (dy_dx) free_matrix(dy_dx,noutputs);
  if (scratch_vector) ascfree((char *)scratch_vector);
  if (sys) system_destroy(sys);
  return result;
}

/**
    Finite-difference Check Arguments...?

    @param arglist list of arguments

    Argument list contains
    . arg1 - model inst to be solved
    . arg2 - array of input instances
    . arg3 - array of output instances
    . arg4 - matrix of partials to be written to
*/
static int FiniteDiffCheckArgs(struct gl_list_t *arglist)
{
  struct Instance *inst;
  unsigned long len;
  unsigned long ninputs, noutputs;

  len = gl_length(arglist);
  if (len != 4) {
    FPRINTF(stderr,"wrong number of args -- 4 expected\n");
    return 1;
  }
  inst = FetchElement(arglist,1,1);
  if (InstanceKind(inst)!=MODEL_INST) {
    FPRINTF(stderr,"Arg #1 is not a model instance\n");
    return 1;
  }
  ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
    /* input args */
  noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
    /* output args */
  len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
    /* partials matrix */
  if (len != (ninputs*noutputs)) {
    FPRINTF(stderr,
        "The array of partials is inconsistent with the args given\n");
    return 1;
  }
  return 0;
}


/*-------------------------------------------------
  SENSITIVITY ANALYSIS CODE
*/

/**
    Sensitivity Analysis

    @param arglist List of arguments

    Argument list contains the following:
      . arg1 - model inst for which the sensitivity analysis is to be done.
      . arg2 - array of input instances.
      . arg3 - array of output instances.
      . arg4 matrix of partials to be written to.
*/
int SensitivityCheckArgs(struct gl_list_t *arglist)
{
  struct Instance *inst;
  unsigned long len;
  unsigned long ninputs, noutputs;

  len = gl_length(arglist);
  if (len != 4) {
    FPRINTF(stderr,"wrong number of args -- 4 expected\n");
    return 1;
  }
  inst = FetchElement(arglist,1,1);
  if (InstanceKind(inst)!=MODEL_INST) {
    FPRINTF(stderr,"Arg #1 is not a model instance\n");
    return 1;
  }
  ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
    /* input args */
  noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
   /* output args */
  len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
        /* partials matrix */
  if (len != (ninputs*noutputs)) {
    FPRINTF(stderr,
        "The array of partials is inconsistent with the args given\n");
    return 1;
  }
  return 0;
}


/**
    @param arglist List of arguments
    @param step_length ...?

    The list of arguments should chontain

      . arg1 - model inst for which the sensitivity analysis is to be done.
      . arg2 - array of input instances.
      . arg3 - new_input instances, for variable projection.
      . arg4 - instance representing the step_length for projection.

    The result will be written to standard out.
*/
int SensitivityAllCheckArgs(struct gl_list_t *arglist, double *step_length)
{
  struct Instance *inst;
  unsigned long len;

  /*

   */

  len = gl_length(arglist);
  if (len != 4) {
    FPRINTF(stderr,"wrong number of args -- 4 expected\n");
    return 1;
  }
  inst = FetchElement(arglist,1,1);
  if (InstanceKind(inst)!=MODEL_INST) {
    FPRINTF(stderr,"Arg #1 is not a model instance\n");
    return 1;
  }
  /*
   * we should be checking that arg2 list contains solver vars
   * and that they are fixed. The same applies to arglist 3... later.
   * We will check and return the steplength though. 0 means dont do
   * the variable projection.
   */
  inst = FetchElement(arglist,4,1);
  *step_length = RealAtomValue(inst);
  if (fabs(*step_length) < 1e-08)
    *step_length = 0.0;
  return 0;
}


/**
    This function is very similar to sensitivity_anal, execept,
    that it perform the analysis on the entire system, based on
    the given inputs. Nothing is returned. As such the call from
    ASCEND is:

<pre>
EXTERN sensitivity_anal_all(
    this_instance,
    u_old[1..n_inputs],
    u_new[1..n_inputs],
    step_length
);</pre>

    The results will be witten to standard out.
    This function is more expensive from a a memory point of view,
    as we keep aroung a dense matrix n_outputs * n_inputs, but here
    n_outputs may be *much* larger depending on problem size.
*/
int sensitivity_anal_all( struct Instance *inst,  /* not used but will be */
             struct gl_list_t *arglist,
             real64 step_length)
{
  struct Instance *which_instance;
  struct gl_list_t *branch2, *branch3;
  dof_t *dof;
  struct var_variable **inputs = NULL, **outputs = NULL;
  int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
  real64 **dy_dx = NULL;
  struct var_variable **vp,**ptr;
  slv_system_t sys = NULL;
  long c;
  int noutputs=0, ninputs;
  var_filter_t vfilter;

  struct var_variable **new_inputs = NULL; /* optional stuff for variable
                      * projection */

  linsolqr_system_t lqr_sys;    /* stuff for the linear system & matrix */
  mtx_matrix_t mtx;
  int32 capacity;
  real64 *scratch_vector = NULL;
  int result=0;

  /* Ignore unused params */
  (void)inst; (void)step_length;

  /*
   * Call the presolve for the system. This should number variables
   * and relations as well create and order the main Jacobian. The
   * only potential problem that I see here is that presolve using
   * the slv0 solver *only* recognizes solver vars. So that if one
   * wanted to see the sensitivity of a *parameter*, it would not
   * be possible. We will have to trap this in CheckArgs.
   *
   * Also the current version of ascend is fucked in how the var module
   * handles variables and their numbering through the interface ptr
   * crap.
   */

  (void)NumberFreeVars(NULL);       /* used to re-init the system */
  (void)NumberRels(NULL);       /* used to re-init the system */
  which_instance = FetchElement(arglist,1,1);
  sys = asc_sens_presolve(which_instance);
  if (!sys) {
    FPRINTF(stderr,"Early termination due to failure in asc_sens_presolve\n");
    result = 1;
    goto error;
  }
  dof = slv_get_dofdata(sys);

  /*
   * prepare the inputs list. We dont need the index of the new_inputs
   * list. We will grab them later if necessary.
   */
  branch2 = gl_fetch(arglist,2); /* input args -- old u_values */
  branch3 = gl_fetch(arglist,3); /* input args -- new u_values */
  ninputs = gl_length(branch2);
  inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
  new_inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);

  inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
  for (c=0;c<ninputs;c++) {
    inputs[c] = (struct var_variable *)gl_fetch(branch2,c+1);
    inputs_ndx_list[c] = var_mindex(inputs[c]);
    new_inputs[c] = (struct var_variable *)gl_fetch(branch3,c+1);
  }
  inputs[ninputs] = NULL;   /* null terminate the list */
  new_inputs[ninputs] = NULL;   /* null terminate the list */

  /*
   * prepare the outputs list. This is where we differ from
   * the other function. The noutputs, and the indexes of these
   * outputs is obtained from the entire solve system.
   */
  vfilter.matchbits = 0;
  noutputs = slv_count_solvers_vars(sys,&vfilter);
  outputs = ASC_NEW_ARRAY(struct var_variable *,noutputs+1);
  outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
  ptr = vp = slv_get_solvers_var_list(sys);
  for (c=0;c<noutputs;c++) {
    outputs[c] = *ptr;
    outputs_ndx_list[c] = var_sindex(*ptr);
    ptr++;
  }
  outputs[noutputs] = NULL; /* null terminate the list */

  /*
   * prepare the results dy_dx. This is the expensive part from a
   * memory point of view. However I would like to have the entire
   * noutputs * ninputs matrix even for a short while so that I
   * can compute a number of  different types of norms.
   */
  dy_dx = make_matrix(noutputs,ninputs);

  result = Compute_J(sys);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
    goto error;
  }

  /*
   * Note: the RHS *has* to added here. We will construct the vector
   * of size matrix capacity and add it. It will be removed after
   * we finish computing dy_dx.
   */
  lqr_sys = slv_get_linsolqr_sys(sys);  /* get the linear system */
  mtx = linsolqr_get_matrix(lqr_sys);   /* get the matrix */
  capacity = mtx_capacity(mtx);
  scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
  linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
  result = LUFactorJacobian1(sys,dof->structural_rank);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
    goto error;
  }
  result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
                   inputs_ndx_list, ninputs,
                   outputs_ndx_list, noutputs);

  linsolqr_remove_rhs(lqr_sys,scratch_vector);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
    goto error;
  }

  /*
   * Do some analysis on the results, inclusive of
   * writing them to someplace useful.
   */
  if (DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx))
    result = 1;

  /*
   * Some experimental projection stuff -- not used now.
   * if (DoProject_X(inputs, new_inputs, step_length,
   *     outputs, ninputs, noutputs, dy_dx))
   * result = 1;
   */

 error:
  if (inputs) ascfree((char *)inputs);
  if (new_inputs) ascfree((char *)new_inputs);
  if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
  if (outputs) ascfree((char *)outputs);
  if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
  if (dy_dx) free_matrix(dy_dx,noutputs);
  if (scratch_vector) ascfree((char *)scratch_vector);
  if (sys) system_destroy(sys);
  return result;
}


/**
    Calls 'DoSolve'

    @see DoSolve
*/
int do_solve_eval( struct Instance *i,
          struct gl_list_t *arglist, void *user_data)
{
  unsigned long len;
  int result;
  struct Instance *inst;
  len = gl_length(arglist);

  /* Ignore unused params */
  (void)i;

  if (len!=2) {
    ERROR_REPORTER_HERE(ASC_USER_ERROR,
      "Wrong number of args to do_solve_eval.");
    return 1;
  }
  inst = FetchElement(arglist,1,1);
  if (!inst)
    return 1;
  result = DoSolve(inst);
  return result;
}


/**
    Finite different evaluate...
*/
int do_finite_diff_eval( struct Instance *i,
             struct gl_list_t *arglist, void *user_data)
{
  int result;

  /* Ignore unused params */
  (void)i;

  if (FiniteDiffCheckArgs(arglist))
    return 1;
  result = finite_difference(arglist);
  return result;
}


int do_sensitivity_eval( struct Instance *i,
             struct gl_list_t *arglist, void *user_data)
{
  int result;
  if (SensitivityCheckArgs(arglist)) {
    return 1;
  }
  result = sensitivity_anal(i,arglist);
  return result;
}

int do_sensitivity_eval_all( struct Instance *i,
                struct gl_list_t *arglist, void *user_data)
{
  int result;
  double step_length = 0.0;
  if (SensitivityAllCheckArgs(arglist,&step_length)) {
    return 1;
  }
  result = sensitivity_anal_all(i,arglist,step_length);
  return result;
}


char sensitivity_help[] =
    "This function does sensitivity analysis dy/dx. It requires 4 args:\n"
    "  1. name: name of a reference instance or SELF.\n"
    "  2. x: x, where x is an array of > solver_var.\n"
    "  3. y: where y is an array of > solver_var.\n"
    "  4. dy/dx: which dy_dx[1..n_y][1..n_x].";

int sensitivity_register(void){

  int result=0;


  result = CreateUserFunctionMethod("do_solve",
                  do_solve_eval,
                  2,NULL,NULL,NULL); /* was 2,0,null */
  result += CreateUserFunctionMethod("do_finite_difference",
                   do_finite_diff_eval,
                   4,NULL,NULL,NULL); /* 4,0,null */
  result += CreateUserFunctionMethod("do_sensitivity",
                   do_sensitivity_eval,
                   4,sensitivity_help,NULL,NULL);
  result += CreateUserFunctionMethod("do_sensitivity_all",
                   do_sensitivity_eval_all,
                   4,"See do_sensitivity_eval for details",NULL,NULL);

  return result;
}

1	johnpye	1162	#include <math.h>
2
3			#include "sensitivity.h"
4			#include <packages/sensitivity.h>
5			#include <compiler/instquery.h>
6			#include <compiler/atomvalue.h>
7			#include <utilities/ascMalloc.h>
8
9
10			/**
11			Build then presolve an instance
12			*/
13			slv_system_t asc_sens_presolve(struct Instance *inst){
14			slv_system_t sys;
15			slv_parameters_t parameters;
16			struct var_variable **vp;
17			struct rel_relation **rp;
18			int ind,len;
19			char *tmp=NULL;
20
21			sys = system_build(inst);
22			if (sys==NULL) {
23			ERROR_REPORTER_HERE(ASC_PROG_ERR,
24			"Something radically wrong in creating solver.");
25			return NULL;
26			}
27			if (g_SlvNumberOfRegisteredClients == 0) {
28			return NULL;
29			}
30			ind = 0;
31			while (strcmp(slv_solver_name(ind),"QRSlv")) {
32			if (ind >= g_SlvNumberOfRegisteredClients) {
33			ERROR_REPORTER_HERE(ASC_PROG_ERR,
34			"QRSlv must be registered client.");
35			return NULL;
36			}
37			++ind;
38			}
39			slv_select_solver(sys,ind);
40
41			slv_get_parameters(sys,&parameters);
42			parameters.partition = 0;
43			slv_set_parameters(sys,&parameters);
44			slv_presolve(sys);
45
46			#if DEBUG
47			vp = slv_get_solvers_var_list(sys);
48			len = slv_get_num_solvers_vars(sys);
49			for (ind=0 ; ind<len; ind++) {
50			tmp = var_make_name(sys,vp[ind]);
51			FPRINTF(stderr,"%s %d\n",tmp,var_sindex(vp[ind]));
52			if (tmp!=NULL) ascfree(tmp);
53			}
54			rp = slv_get_solvers_rel_list(sys);
55			len = slv_get_num_solvers_rels(sys);
56			for (ind=0 ; ind<len ; ind++) {
57			tmp = rel_make_name(sys,rp[ind]);
58			FPRINTF(stderr,"%s %d\n",tmp,rel_sindex(rp[ind]));
59			if (tmp) ascfree(tmp);
60			}
61			#endif
62			return sys;
63			}
64
65
66			int sensitivity_anal(
67			struct Instance inst, / not used but will be */
68			struct gl_list_t *arglist){
69			struct Instance which_instance,tmp_inst, *atominst;
70			struct gl_list_t *branch;
71			struct var_variable **vlist = NULL;
72			int inputs_ndx_list = NULL, outputs_ndx_list = NULL;
73			real64 **dy_dx = NULL;
74			slv_system_t sys = NULL;
75			int c;
76			int noutputs = 0;
77			int ninputs;
78			int i,j;
79			int offset;
80			dof_t *dof;
81			int num_vars,ind,found;
82
83			linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
84			mtx_matrix_t mtx;
85			int32 capacity;
86			real64 *scratch_vector = NULL;
87			int result=0;
88
89			/* Ignore unused params */
90			(void) inst;
91
92			(void)NumberFreeVars(NULL); /* used to re-init the system */
93			(void)NumberRels(NULL); /* used to re-init the system */
94			which_instance = FetchElement(arglist,1,1);
95			sys = asc_sens_presolve(which_instance);
96			if (!sys) {
97			FPRINTF(stderr,"Early termination due to failure in Presolve\n");
98			result = 1;
99			goto error;
100			}
101
102			dof = slv_get_dofdata(sys);
103			if (!(dof->n_rows == dof->n_cols &&
104			dof->n_rows == dof->structural_rank)) {
105			FPRINTF(stderr,"Early termination: non square system\n");
106			result = 1;
107			goto error;
108			}
109			/*
110			* prepare the inputs list
111			*/
112			vlist = slv_get_solvers_var_list(sys);
113
114			branch = gl_fetch(arglist,2); /* input args */
115			ninputs = gl_length(branch);
116			inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
117
118			num_vars = slv_get_num_solvers_vars(sys);
119			for (c=0;c<ninputs;c++) {
120			atominst = (struct Instance *)gl_fetch(branch,c+1);
121			found = 0;
122			ind = num_vars - 1;
123			/* search backwards because fixed vars are at the end of the
124			var list */
125			while (!found && ind >= 0){
126			if (var_instance(vlist[ind]) == atominst) {
127			inputs_ndx_list[c] = var_sindex(vlist[ind]);
128			found = 1;
129			}
130			--ind;
131			}
132			if (!found) {
133			FPRINTF(stderr,"Sensitivity input variable not found\n");
134			result = 1;
135			goto error;
136			}
137			}
138
139			/*
140			* prepare the outputs list
141			*/
142			branch = gl_fetch(arglist,3); /* output args */
143			noutputs = gl_length(branch);
144			outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
145			for (c=0;c<noutputs;c++) {
146			atominst = (struct Instance *)gl_fetch(branch,c+1);
147			found = 0;
148			ind = 0;
149			while (!found && ind < num_vars){
150			if (var_instance(vlist[ind]) == atominst) {
151			outputs_ndx_list[c] = var_sindex(vlist[ind]);
152			found = 1;
153			}
154			++ind;
155			}
156			if (!found) {
157			FPRINTF(stderr,"Sensitivity ouput variable not found\n");
158			result = 1;
159			goto error;
160			}
161			}
162
163			/*
164			* prepare the results dy_dx.
165			*/
166			dy_dx = make_matrix(noutputs,ninputs);
167
168
169			result = Compute_J(sys);
170			if (result) {
171			FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
172			goto error;
173			}
174
175			/*
176			* Note: the RHS has to added here. We will construct the vector
177			* of size matrix capacity and add it. It will be removed after
178			* we finish computing dy_dx.
179			*/
180			lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
181			mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
182			capacity = mtx_capacity(mtx);
183			scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
184			linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
185			result = LUFactorJacobian1(sys,dof->structural_rank);
186			if (result) {
187			FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
188			goto error;
189			}
190			result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
191			inputs_ndx_list, ninputs,
192			outputs_ndx_list, noutputs);
193
194			linsolqr_remove_rhs(lqr_sys,scratch_vector);
195			if (result) {
196			FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
197			goto error;
198			}
199
200			/*
201			* Write the results back to the partials array in the
202			* instance tree
203			*/
204			offset = 0;
205			for (i=0;i<noutputs;i++) {
206			for (j=0;j<ninputs;j++) {
207			tmp_inst = FetchElement(arglist,4,offset+j+1);
208			SetRealAtomValue(tmp_inst,dy_dx[i][j],(unsigned)0);
209			#if DEBUG
210			FPRINTF(stderr,"%12.8f i%d j%d",dy_dx[i][j],i,j);
211			#endif
212			}
213			#if DEBUG
214			FPRINTF(stderr,"\n");
215			#endif
216			offset += ninputs;
217			}
218			#if DEBUG
219			FPRINTF(stderr,"\n");
220			#endif
221
222			error:
223			if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
224			if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
225			if (dy_dx) free_matrix(dy_dx,noutputs);
226			if (scratch_vector) ascfree((char *)scratch_vector);
227			if (sys) system_destroy(sys);
228			return result;
229			}
230
231			/**
232			Finite-difference Check Arguments...?
233
234			@param arglist list of arguments
235
236			Argument list contains
237			. arg1 - model inst to be solved
238			. arg2 - array of input instances
239			. arg3 - array of output instances
240			. arg4 - matrix of partials to be written to
241			*/
242			static int FiniteDiffCheckArgs(struct gl_list_t *arglist)
243			{
244			struct Instance *inst;
245			unsigned long len;
246			unsigned long ninputs, noutputs;
247
248			len = gl_length(arglist);
249			if (len != 4) {
250			FPRINTF(stderr,"wrong number of args -- 4 expected\n");
251			return 1;
252			}
253			inst = FetchElement(arglist,1,1);
254			if (InstanceKind(inst)!=MODEL_INST) {
255			FPRINTF(stderr,"Arg #1 is not a model instance\n");
256			return 1;
257			}
258			ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
259			/* input args */
260			noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
261			/* output args */
262			len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
263			/* partials matrix */
264			if (len != (ninputs*noutputs)) {
265			FPRINTF(stderr,
266			"The array of partials is inconsistent with the args given\n");
267			return 1;
268			}
269			return 0;
270			}
271
272
273
274			/*-------------------------------------------------
275			SENSITIVITY ANALYSIS CODE
276			*/
277
278			/**
279			Sensitivity Analysis
280
281			@param arglist List of arguments
282
283			Argument list contains the following:
284			. arg1 - model inst for which the sensitivity analysis is to be done.
285			. arg2 - array of input instances.
286			. arg3 - array of output instances.
287			. arg4 matrix of partials to be written to.
288			*/
289			int SensitivityCheckArgs(struct gl_list_t *arglist)
290			{
291			struct Instance *inst;
292			unsigned long len;
293			unsigned long ninputs, noutputs;
294
295			len = gl_length(arglist);
296			if (len != 4) {
297			FPRINTF(stderr,"wrong number of args -- 4 expected\n");
298			return 1;
299			}
300			inst = FetchElement(arglist,1,1);
301			if (InstanceKind(inst)!=MODEL_INST) {
302			FPRINTF(stderr,"Arg #1 is not a model instance\n");
303			return 1;
304			}
305			ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
306			/* input args */
307			noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
308			/* output args */
309			len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
310			/* partials matrix */
311			if (len != (ninputs*noutputs)) {
312			FPRINTF(stderr,
313			"The array of partials is inconsistent with the args given\n");
314			return 1;
315			}
316			return 0;
317			}
318
319
320			/**
321			@param arglist List of arguments
322			@param step_length ...?
323
324			The list of arguments should chontain
325
326			. arg1 - model inst for which the sensitivity analysis is to be done.
327			. arg2 - array of input instances.
328			. arg3 - new_input instances, for variable projection.
329			. arg4 - instance representing the step_length for projection.
330
331			The result will be written to standard out.
332			*/
333			int SensitivityAllCheckArgs(struct gl_list_t arglist, double step_length)
334			{
335			struct Instance *inst;
336			unsigned long len;
337
338			/*
339
340			*/
341
342			len = gl_length(arglist);
343			if (len != 4) {
344			FPRINTF(stderr,"wrong number of args -- 4 expected\n");
345			return 1;
346			}
347			inst = FetchElement(arglist,1,1);
348			if (InstanceKind(inst)!=MODEL_INST) {
349			FPRINTF(stderr,"Arg #1 is not a model instance\n");
350			return 1;
351			}
352			/*
353			* we should be checking that arg2 list contains solver vars
354			* and that they are fixed. The same applies to arglist 3... later.
355			* We will check and return the steplength though. 0 means dont do
356			* the variable projection.
357			*/
358			inst = FetchElement(arglist,4,1);
359			*step_length = RealAtomValue(inst);
360			if (fabs(*step_length) < 1e-08)
361			*step_length = 0.0;
362			return 0;
363			}
364
365
366			/**
367			This function is very similar to sensitivity_anal, execept,
368			that it perform the analysis on the entire system, based on
369			the given inputs. Nothing is returned. As such the call from
370			ASCEND is:
371
372			<pre>
373			EXTERN sensitivity_anal_all(
374			this_instance,
375			u_old[1..n_inputs],
376			u_new[1..n_inputs],
377			step_length
378			);</pre>
379
380			The results will be witten to standard out.
381			This function is more expensive from a a memory point of view,
382			as we keep aroung a dense matrix n_outputs * n_inputs, but here
383			n_outputs may be much larger depending on problem size.
384			*/
385			int sensitivity_anal_all( struct Instance inst, / not used but will be */
386			struct gl_list_t *arglist,
387			real64 step_length)
388			{
389			struct Instance *which_instance;
390			struct gl_list_t branch2, branch3;
391			dof_t *dof;
392			struct var_variable inputs = NULL, outputs = NULL;
393			int inputs_ndx_list = NULL, outputs_ndx_list = NULL;
394			real64 **dy_dx = NULL;
395			struct var_variable vp,ptr;
396			slv_system_t sys = NULL;
397			long c;
398			int noutputs=0, ninputs;
399			var_filter_t vfilter;
400
401			struct var_variable *new_inputs = NULL; / optional stuff for variable
402			* projection */
403
404			linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
405			mtx_matrix_t mtx;
406			int32 capacity;
407			real64 *scratch_vector = NULL;
408			int result=0;
409
410			/* Ignore unused params */
411			(void)inst; (void)step_length;
412
413			/*
414			* Call the presolve for the system. This should number variables
415			* and relations as well create and order the main Jacobian. The
416			* only potential problem that I see here is that presolve using
417			* the slv0 solver only recognizes solver vars. So that if one
418			* wanted to see the sensitivity of a parameter, it would not
419			* be possible. We will have to trap this in CheckArgs.
420			*
421			* Also the current version of ascend is fucked in how the var module
422			* handles variables and their numbering through the interface ptr
423			* crap.
424			*/
425
426			(void)NumberFreeVars(NULL); /* used to re-init the system */
427			(void)NumberRels(NULL); /* used to re-init the system */
428			which_instance = FetchElement(arglist,1,1);
429			sys = asc_sens_presolve(which_instance);
430			if (!sys) {
431			FPRINTF(stderr,"Early termination due to failure in asc_sens_presolve\n");
432			result = 1;
433			goto error;
434			}
435			dof = slv_get_dofdata(sys);
436
437			/*
438			* prepare the inputs list. We dont need the index of the new_inputs
439			* list. We will grab them later if necessary.
440			*/
441			branch2 = gl_fetch(arglist,2); /* input args -- old u_values */
442			branch3 = gl_fetch(arglist,3); /* input args -- new u_values */
443			ninputs = gl_length(branch2);
444			inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
445			new_inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
446
447			inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
448			for (c=0;c<ninputs;c++) {
449			inputs[c] = (struct var_variable *)gl_fetch(branch2,c+1);
450			inputs_ndx_list[c] = var_mindex(inputs[c]);
451			new_inputs[c] = (struct var_variable *)gl_fetch(branch3,c+1);
452			}
453			inputs[ninputs] = NULL; /* null terminate the list */
454			new_inputs[ninputs] = NULL; /* null terminate the list */
455
456			/*
457			* prepare the outputs list. This is where we differ from
458			* the other function. The noutputs, and the indexes of these
459			* outputs is obtained from the entire solve system.
460			*/
461			vfilter.matchbits = 0;
462			noutputs = slv_count_solvers_vars(sys,&vfilter);
463			outputs = ASC_NEW_ARRAY(struct var_variable *,noutputs+1);
464			outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
465			ptr = vp = slv_get_solvers_var_list(sys);
466			for (c=0;c<noutputs;c++) {
467			outputs[c] = *ptr;
468			outputs_ndx_list[c] = var_sindex(*ptr);
469			ptr++;
470			}
471			outputs[noutputs] = NULL; /* null terminate the list */
472
473			/*
474			* prepare the results dy_dx. This is the expensive part from a
475			* memory point of view. However I would like to have the entire
476			* noutputs * ninputs matrix even for a short while so that I
477			* can compute a number of different types of norms.
478			*/
479			dy_dx = make_matrix(noutputs,ninputs);
480
481			result = Compute_J(sys);
482			if (result) {
483			FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
484			goto error;
485			}
486
487			/*
488			* Note: the RHS has to added here. We will construct the vector
489			* of size matrix capacity and add it. It will be removed after
490			* we finish computing dy_dx.
491			*/
492			lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
493			mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
494			capacity = mtx_capacity(mtx);
495			scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
496			linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
497			result = LUFactorJacobian1(sys,dof->structural_rank);
498			if (result) {
499			FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
500			goto error;
501			}
502			result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
503			inputs_ndx_list, ninputs,
504			outputs_ndx_list, noutputs);
505
506			linsolqr_remove_rhs(lqr_sys,scratch_vector);
507			if (result) {
508			FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
509			goto error;
510			}
511
512			/*
513			* Do some analysis on the results, inclusive of
514			* writing them to someplace useful.
515			*/
516			if (DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx))
517			result = 1;
518
519			/*
520			* Some experimental projection stuff -- not used now.
521			* if (DoProject_X(inputs, new_inputs, step_length,
522			* outputs, ninputs, noutputs, dy_dx))
523			* result = 1;
524			*/
525
526			error:
527			if (inputs) ascfree((char *)inputs);
528			if (new_inputs) ascfree((char *)new_inputs);
529			if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
530			if (outputs) ascfree((char *)outputs);
531			if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
532			if (dy_dx) free_matrix(dy_dx,noutputs);
533			if (scratch_vector) ascfree((char *)scratch_vector);
534			if (sys) system_destroy(sys);
535			return result;
536			}
537
538
539
540			/**
541			Calls 'DoSolve'
542
543			@see DoSolve
544			*/
545			int do_solve_eval( struct Instance *i,
546			struct gl_list_t arglist, void user_data)
547			{
548			unsigned long len;
549			int result;
550			struct Instance *inst;
551			len = gl_length(arglist);
552
553			/* Ignore unused params */
554			(void)i;
555
556			if (len!=2) {
557			ERROR_REPORTER_HERE(ASC_USER_ERROR,
558			"Wrong number of args to do_solve_eval.");
559			return 1;
560			}
561			inst = FetchElement(arglist,1,1);
562			if (!inst)
563			return 1;
564			result = DoSolve(inst);
565			return result;
566			}
567
568
569			/**
570			Finite different evaluate...
571			*/
572			int do_finite_diff_eval( struct Instance *i,
573			struct gl_list_t arglist, void user_data)
574			{
575			int result;
576
577			/* Ignore unused params */
578			(void)i;
579
580			if (FiniteDiffCheckArgs(arglist))
581			return 1;
582			result = finite_difference(arglist);
583			return result;
584			}
585
586
587			int do_sensitivity_eval( struct Instance *i,
588			struct gl_list_t arglist, void user_data)
589			{
590			int result;
591			if (SensitivityCheckArgs(arglist)) {
592			return 1;
593			}
594			result = sensitivity_anal(i,arglist);
595			return result;
596			}
597
598			int do_sensitivity_eval_all( struct Instance *i,
599			struct gl_list_t arglist, void user_data)
600			{
601			int result;
602			double step_length = 0.0;
603			if (SensitivityAllCheckArgs(arglist,&step_length)) {
604			return 1;
605			}
606			result = sensitivity_anal_all(i,arglist,step_length);
607			return result;
608			}
609
610
611			char sensitivity_help[] =
612			"This function does sensitivity analysis dy/dx. It requires 4 args:\n"
613			" 1. name: name of a reference instance or SELF.\n"
614			" 2. x: x, where x is an array of > solver_var.\n"
615			" 3. y: where y is an array of > solver_var.\n"
616			" 4. dy/dx: which dy_dx[1..n_y][1..n_x].";
617
618			int sensitivity_register(void){
619
620			int result=0;
621
622
623			result = CreateUserFunctionMethod("do_solve",
624			do_solve_eval,
625			2,NULL,NULL,NULL); /* was 2,0,null */
626			result += CreateUserFunctionMethod("do_finite_difference",
627			do_finite_diff_eval,
628			4,NULL,NULL,NULL); /* 4,0,null */
629			result += CreateUserFunctionMethod("do_sensitivity",
630			do_sensitivity_eval,
631			4,sensitivity_help,NULL,NULL);
632			result += CreateUserFunctionMethod("do_sensitivity_all",
633			do_sensitivity_eval_all,
634			4,"See do_sensitivity_eval for details",NULL,NULL);
635
636			return result;
637			}
638