/[ascend]/trunk/models/sensitivity/sensitivity.c
ViewVC logotype

Diff of /trunk/models/sensitivity/sensitivity.c

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

revision 1162 by johnpye, Wed Jan 17 04:34:59 2007 UTC revision 1163 by johnpye, Wed Jan 17 06:33:30 2007 UTC
# Line 1  Line 1 
1    /*  ASCEND modelling environment
2        Copyright (C) 1996-2007 Carnegie Mellon University
3    
4        This program is free software; you can redistribute it and/or modify
5        it under the terms of the GNU General Public License as published by
6        the Free Software Foundation; either version 2, or (at your option)
7        any later version.
8    
9        This program is distributed in the hope that it will be useful,
10        but WITHOUT ANY WARRANTY; without even the implied warranty of
11        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12        GNU General Public License for more details.
13        
14        You should have received a copy of the GNU General Public License
15        along with this program; if not, write to the Free Software
16        Foundation, Inc., 59 Temple Place - Suite 330,
17        Boston, MA 02111-1307, USA.
18    *//** @file
19    
20        This file attempts to implement the extraction of dy_dx from
21        a system of equations. If one considers a black-box where x are
22        the input variables, u are input parameters, y are the output
23        variables, f(x,y,u) is the system of equations that will be solved
24        for given x and u, then one can extract the sensitivity information
25        of y wrt x.
26    
27        One crude but simple way of doing this is to finite-difference the
28        given black box, i.e, perturb x, n_x times by delta x, resolve
29        the system of equations at the new value of x, and compute
30    
31            dy/dx = (f(x_1) - f(x_2))/(x_1 - x_2).
32    
33        This is known to be very expensive.
34    
35        The solution that will be adopted here is to formulate the Jacobian J of
36        the system, (or the system is already solved, to grab the jacobian at
37        the solution. Develop the sensitivity matrix df/dx by exact numnerical
38        differentiation; this will be a n x n_x matrix. Compute the LU factors
39        for J, and to solve column by column to : LU dz/dx = df/dx. Here
40        z, represents all the internal variables to the system and the strictly
41        output variables y. In other words J = df/dz.
42    
43        Given the solution of the above problem we then simply extract the rows
44        of dz/dx, that pertain to the y variables that we are interested in;
45        this will give us dy/dx.
46    
47        @todo There is are files in tcltk called Sensitivity.[ch]. Do we need them?
48    */
49    
50  #include <math.h>  #include <math.h>
51    
 #include "sensitivity.h"  
52  #include <packages/sensitivity.h>  #include <packages/sensitivity.h>
53  #include <compiler/instquery.h>  #include <compiler/instquery.h>
54  #include <compiler/atomvalue.h>  #include <compiler/atomvalue.h>
55  #include <utilities/ascMalloc.h>  #include <utilities/ascMalloc.h>
56    #include <compiler/extfunc.h>
57    #include <general/mathmacros.h>
58    
59    /* #define SENSITIVITY_DEBUG */
60    
61    ASC_EXPORT int sensitivity_register(void);
62    
63    ExtMethodRun do_solve_eval;
64    ExtMethodRun do_finite_diff_eval;
65    ExtMethodRun do_sensitivity_eval;
66    ExtMethodRun do_sensitivity_eval_all;
67    
68  /**  /**
69      Build then presolve an instance      Build then presolve an instance
# Line 13  Line 71 
71  slv_system_t asc_sens_presolve(struct Instance *inst){  slv_system_t asc_sens_presolve(struct Instance *inst){
72    slv_system_t sys;    slv_system_t sys;
73    slv_parameters_t parameters;    slv_parameters_t parameters;
74      int ind;
75    #ifdef SENSITIVITY_DEBUG
76    struct var_variable **vp;    struct var_variable **vp;
77    struct rel_relation **rp;    struct rel_relation **rp;
78    int ind,len;    int len;
79    char *tmp=NULL;    char *tmp=NULL;
80    #endif
81    sys = system_build(inst);    sys = system_build(inst);
82    if (sys==NULL) {    if (sys==NULL) {
83      ERROR_REPORTER_HERE(ASC_PROG_ERR,      ERROR_REPORTER_HERE(ASC_PROG_ERR,
# Line 43  slv_system_t asc_sens_presolve(struct In Line 103  slv_system_t asc_sens_presolve(struct In
103    slv_set_parameters(sys,&parameters);    slv_set_parameters(sys,&parameters);
104    slv_presolve(sys);    slv_presolve(sys);
105    
106  #if DEBUG  #ifdef SENSITIVITY_DEBUG
107    vp = slv_get_solvers_var_list(sys);    vp = slv_get_solvers_var_list(sys);
108    len = slv_get_num_solvers_vars(sys);    len = slv_get_num_solvers_vars(sys);
109    for (ind=0 ; ind<len; ind++) {    for (ind=0 ; ind<len; ind++) {
110      tmp = var_make_name(sys,vp[ind]);      tmp = var_make_name(sys,vp[ind]);
111      FPRINTF(stderr,"%s  %d\n",tmp,var_sindex(vp[ind]));      CONSOLE_DEBUG("%s  %d\n",tmp,var_sindex(vp[ind]));
112      if (tmp!=NULL) ascfree(tmp);      if (tmp!=NULL) ascfree(tmp);
113    }    }
114    rp = slv_get_solvers_rel_list(sys);    rp = slv_get_solvers_rel_list(sys);
115    len = slv_get_num_solvers_rels(sys);    len = slv_get_num_solvers_rels(sys);
116    for (ind=0 ; ind<len ; ind++) {    for (ind=0 ; ind<len ; ind++) {
117      tmp = rel_make_name(sys,rp[ind]);      tmp = rel_make_name(sys,rp[ind]);
118      FPRINTF(stderr,"%s  %d\n",tmp,rel_sindex(rp[ind]));      CONSOLE_DEBUG("%s  %d\n",tmp,rel_sindex(rp[ind]));
119      if (tmp) ascfree(tmp);      if (tmp) ascfree(tmp);
120    }    }
121  #endif  #endif
122    return sys;    return sys;
123  }  }
124    
125    /*
126        LU Factor Jacobian
127    
128  int sensitivity_anal(      @note The RHS will have been  will already have been added before
129               struct Instance *inst, /* not used but will be */          calling this function.
              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;  
130    
131    /* Ignore unused params */      @NOTE there is another version of this floating around in packages/senstivity.c. The
132    (void) inst;      other one uses dense matrices so probably this one's better?
133    */
134    (void)NumberFreeVars(NULL);       /* used to re-init the system */  int LUFactorJacobian1(slv_system_t sys,int rank){
135    (void)NumberRels(NULL);       /* used to re-init the system */    linsolqr_system_t lqr_sys;
136    which_instance = FetchElement(arglist,1,1);    mtx_region_t region;
137    sys = asc_sens_presolve(which_instance);    enum factor_method fm;
   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;  
     }  
   }  
138    
139    /*    mtx_region(&region,0,rank-1,0,rank-1);    /* set the region */
140     * prepare the outputs list    lqr_sys = slv_get_linsolqr_sys(sys);  /* get the linear system */
    */  
   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;  
     }  
   }  
141    
142    /*    linsolqr_matrix_was_changed(lqr_sys);
143     * prepare the results dy_dx.    linsolqr_reorder(lqr_sys,&region,natural);
    */  
   dy_dx = make_matrix(noutputs,ninputs);  
144    
145      fm = linsolqr_fmethod(lqr_sys);
146      if (fm == unknown_f) fm = ranki_kw2; /* make sure somebody set it */
147      linsolqr_factor(lqr_sys,fm);
148    
149    result = Compute_J(sys);    return 0;
150    if (result) {  }
     FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");  
     goto error;  
   }  
151    
   /*  
    * 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;  
   }  
152    
153    /*  int sensitivity_anal(
154     * Write the results back to the partials array in the           struct Instance *inst, /* not used but will be */
155     * instance tree           struct gl_list_t *arglist
156     */  ){
157    offset = 0;      struct Instance *which_instance,*tmp_inst, *atominst;
158    for (i=0;i<noutputs;i++) {      struct gl_list_t *branch;
159      for (j=0;j<ninputs;j++) {      struct var_variable **vlist = NULL;
160        tmp_inst = FetchElement(arglist,4,offset+j+1);      int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
161        SetRealAtomValue(tmp_inst,dy_dx[i][j],(unsigned)0);      real64 **dy_dx = NULL;
162  #if DEBUG      slv_system_t sys = NULL;
163        FPRINTF(stderr,"%12.8f   i%d j%d",dy_dx[i][j],i,j);      int c;
164        int noutputs = 0;
165        int ninputs;
166        int i,j;
167        int offset;
168        dof_t *dof;
169        int num_vars,ind,found;
170    
171        linsolqr_system_t lqr_sys;  /* stuff for the linear system & matrix */
172        mtx_matrix_t mtx;
173        int32 capacity;
174        real64 *scratch_vector = NULL;
175        int result=0;
176    
177        /* Ignore unused params */
178        (void) inst;
179    
180        (void)NumberFreeVars(NULL);     /* used to re-init the system */
181        (void)NumberRels(NULL);     /* used to re-init the system */
182        which_instance = FetchElement(arglist,1,1);
183        sys = asc_sens_presolve(which_instance);
184        if (!sys) {
185            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in Presolve\n");
186            result = 1;
187            goto error;
188        }
189    
190        dof = slv_get_dofdata(sys);
191        if (!(dof->n_rows == dof->n_cols &&
192            dof->n_rows == dof->structural_rank)) {
193            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination: non square system\n");
194            result = 1;
195            goto error;
196        }
197    
198        CONSOLE_DEBUG("Presolved, square");
199    
200        /*
201            prepare the inputs list
202        */
203        vlist = slv_get_solvers_var_list(sys);
204    
205        branch = gl_fetch(arglist,2); /* input args */
206        ninputs = gl_length(branch);
207        inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
208    
209        num_vars = slv_get_num_solvers_vars(sys);
210        for (c=0;c<ninputs;c++) {
211            atominst = (struct Instance *)gl_fetch(branch,c+1);
212            found = 0;
213            ind = num_vars - 1;
214            /* search backwards because fixed vars are at the end of the
215               var list */
216            while (!found && ind >= 0){
217                if (var_instance(vlist[ind]) == atominst) {
218                    inputs_ndx_list[c] = var_sindex(vlist[ind]);
219                    found = 1;
220                }
221                --ind;
222            }
223            if (!found) {
224                ERROR_REPORTER_HERE(ASC_PROG_ERR,"Sensitivity input variable not found\n");
225                result = 1;
226                goto error;
227            }
228        }
229    
230        CONSOLE_DEBUG("%d inputs",ninputs);
231    
232        /*
233            prepare the outputs list
234        */
235        branch = gl_fetch(arglist,3); /* output args */
236        noutputs = gl_length(branch);
237        outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
238        for (c=0;c<noutputs;c++) {
239            atominst = (struct Instance *)gl_fetch(branch,c+1);
240            found = 0;
241            ind = 0;
242            while (!found && ind < num_vars){
243                if (var_instance(vlist[ind]) == atominst) {
244                    outputs_ndx_list[c] = var_sindex(vlist[ind]);
245                    found = 1;
246                }
247                ++ind;
248            }
249            if (!found) {
250                ERROR_REPORTER_HERE(ASC_PROG_ERR,"Sensitivity ouput variable not found\n");
251                result = 1;
252                goto error;
253            }
254        }
255        
256        CONSOLE_DEBUG("%d outputs",noutputs);
257    
258        /*
259            prepare the results dy_dx.
260        */
261        dy_dx = make_matrix(noutputs,ninputs);
262    
263        result = Compute_J(sys);
264        if (result) {
265            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in calc Jacobian\n");
266            goto error;
267        }
268    
269        CONSOLE_DEBUG("Computed Jacobian");
270    
271        /*
272            Note: the RHS *has* to added here. We will construct the vector
273            of size matrix capacity and add it. It will be removed after
274            we finish computing dy_dx.
275        */
276        lqr_sys = slv_get_linsolqr_sys(sys);    /* get the linear system */
277        mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
278        capacity = mtx_capacity(mtx);
279        scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
280        linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
281        result = LUFactorJacobian1(sys,dof->structural_rank);
282        if(result){
283            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in LUFactorJacobian\n");
284            goto error;
285        }
286        result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
287            inputs_ndx_list, ninputs,
288            outputs_ndx_list, noutputs
289        );
290    
291        linsolqr_remove_rhs(lqr_sys,scratch_vector);
292        if (result) {
293            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in Compute_dy_dx\n");
294            goto error;
295        }
296        
297        CONSOLE_DEBUG("Computed dy/dx");
298    
299        /*
300            Write the results back to the partials array in the
301            instance tree
302        */
303        offset = 0;
304        for (i=0;i<noutputs;i++) {
305            for (j=0;j<ninputs;j++) {
306                tmp_inst = FetchElement(arglist,4,offset+j+1);
307                SetRealAtomValue(tmp_inst,dy_dx[i][j],(unsigned)0);
308    #ifdef SENSITIVITY_DEBUG
309                CONSOLE_DEBUG("%12.8f   i%d j%d",dy_dx[i][j],i,j);
310  #endif  #endif
311      }          }
312  #if DEBUG  #ifdef SENSITIVITY_DEBUG
313      FPRINTF(stderr,"\n");          CONSOLE_DEBUG("\n");
314  #endif  #endif
315      offset += ninputs;          offset += ninputs;
316    }      }
317  #if DEBUG  #ifdef SENSITIVITY_DEBUG
318    FPRINTF(stderr,"\n");      CONSOLE_DEBUG("\n");
319  #endif  #endif
320    
321        ERROR_REPORTER_HERE(ASC_USER_SUCCESS
322            ,"Sensitivity results for %d vars were written back to the model"
323            ,noutputs
324        );
325    
326  error:  error:
327    if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);      if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
328    if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);      if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
329    if (dy_dx) free_matrix(dy_dx,noutputs);      if (dy_dx) free_matrix(dy_dx,noutputs);
330    if (scratch_vector) ascfree((char *)scratch_vector);      if (scratch_vector) ascfree((char *)scratch_vector);
331    if (sys) system_destroy(sys);      if (sys) system_destroy(sys);
332    return result;      return result;
333  }  }
334    
335  /**  /**
# Line 247  static int FiniteDiffCheckArgs(struct gl Line 351  static int FiniteDiffCheckArgs(struct gl
351    
352    len = gl_length(arglist);    len = gl_length(arglist);
353    if (len != 4) {    if (len != 4) {
354      FPRINTF(stderr,"wrong number of args -- 4 expected\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"wrong number of args -- 4 expected\n");
355      return 1;      return 1;
356    }    }
357    inst = FetchElement(arglist,1,1);    inst = FetchElement(arglist,1,1);
358    if (InstanceKind(inst)!=MODEL_INST) {    if (InstanceKind(inst)!=MODEL_INST) {
359      FPRINTF(stderr,"Arg #1 is not a model instance\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Arg #1 is not a model instance\n");
360      return 1;      return 1;
361    }    }
362    ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));    ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
# Line 262  static int FiniteDiffCheckArgs(struct gl Line 366  static int FiniteDiffCheckArgs(struct gl
366    len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));    len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
367      /* partials matrix */      /* partials matrix */
368    if (len != (ninputs*noutputs)) {    if (len != (ninputs*noutputs)) {
369      FPRINTF(stderr,      ERROR_REPORTER_HERE(ASC_PROG_ERR,
370          "The array of partials is inconsistent with the args given\n");          "The array of partials is inconsistent with the args given."
371        );
372      return 1;      return 1;
373    }    }
374    return 0;    return 0;
# Line 286  static int FiniteDiffCheckArgs(struct gl Line 391  static int FiniteDiffCheckArgs(struct gl
391        . arg3 - array of output instances.        . arg3 - array of output instances.
392        . arg4 matrix of partials to be written to.        . arg4 matrix of partials to be written to.
393  */  */
394  int SensitivityCheckArgs(struct gl_list_t *arglist)  int SensitivityCheckArgs(struct gl_list_t *arglist){
 {  
395    struct Instance *inst;    struct Instance *inst;
396    unsigned long len;    unsigned long len;
397    unsigned long ninputs, noutputs;    unsigned long ninputs, noutputs;
398    
399    len = gl_length(arglist);    len = gl_length(arglist);
400    if (len != 4) {    if (len != 4) {
401      FPRINTF(stderr,"wrong number of args -- 4 expected\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"wrong number of args -- 4 expected\n");
402      return 1;      return 1;
403    }    }
404    inst = FetchElement(arglist,1,1);    inst = FetchElement(arglist,1,1);
405    if (InstanceKind(inst)!=MODEL_INST) {    if (InstanceKind(inst)!=MODEL_INST) {
406      FPRINTF(stderr,"Arg #1 is not a model instance\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Arg #1 is not a model instance\n");
407      return 1;      return 1;
408    }    }
409    ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));    ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
# Line 309  int SensitivityCheckArgs(struct gl_list_ Line 413  int SensitivityCheckArgs(struct gl_list_
413    len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));    len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
414          /* partials matrix */          /* partials matrix */
415    if (len != (ninputs*noutputs)) {    if (len != (ninputs*noutputs)) {
416      FPRINTF(stderr,      ERROR_REPORTER_HERE(ASC_PROG_ERR,
417          "The array of partials is inconsistent with the args given\n");          "The array of partials is inconsistent with the args given\n");
418      return 1;      return 1;
419    }    }
# Line 335  int SensitivityAllCheckArgs(struct gl_li Line 439  int SensitivityAllCheckArgs(struct gl_li
439    struct Instance *inst;    struct Instance *inst;
440    unsigned long len;    unsigned long len;
441    
   /*  
   
    */  
   
442    len = gl_length(arglist);    len = gl_length(arglist);
443    if (len != 4) {    if (len != 4) {
444      FPRINTF(stderr,"wrong number of args -- 4 expected\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"wrong number of args -- 4 expected\n");
445      return 1;      return 1;
446    }    }
447    inst = FetchElement(arglist,1,1);    inst = FetchElement(arglist,1,1);
448    if (InstanceKind(inst)!=MODEL_INST) {    if (InstanceKind(inst)!=MODEL_INST) {
449      FPRINTF(stderr,"Arg #1 is not a model instance\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Arg #1 is not a model instance\n");
450      return 1;      return 1;
451    }    }
452    /*    /*
# Line 362  int SensitivityAllCheckArgs(struct gl_li Line 462  int SensitivityAllCheckArgs(struct gl_li
462    return 0;    return 0;
463  }  }
464    
465    /**
466        Do Data Analysis
467    */
468    int DoDataAnalysis(struct var_variable **inputs,
469                  struct var_variable **outputs,
470                  int ninputs, int noutputs,
471                  real64 **dy_dx)
472    {
473      FILE *fp;
474      double *norm_2, *norm_1;
475      double input_nominal,maxvalue,sum;
476      int i,j;
477    
478      norm_1 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
479      norm_2 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
480    
481      fp = fopen("sensitivity.out","w+");
482      if (!fp) return 0;
483    
484      /*
485       * calculate the 1 and 2 norms; cache them away so that we
486       * can pretty print them. Style is everything !.
487       *
488       */
489      for (j=0;j<ninputs;j++) {
490        input_nominal = var_nominal(inputs[j]);
491        maxvalue = sum = 0;
492        for (i=0;i<noutputs;i++) {
493          dy_dx[i][j] *= input_nominal/var_nominal(outputs[i]);
494          maxvalue = MAX(fabs(dy_dx[i][j]),maxvalue);
495          sum += dy_dx[i][j]*dy_dx[i][j];
496        }
497        norm_1[j] = maxvalue;
498        norm_2[j] = sum;
499      }
500    
501      for (j=0;j<ninputs;j++) {     /* print the var_index */
502        fprintf(fp,"%8d    ",var_mindex(inputs[j]));
503      }
504      fprintf(fp,"\n");
505    
506      for (j=0;j<ninputs;j++) {     /* print the 1 norms */
507        fprintf(fp,"%-#18.8f    ",norm_1[j]);
508      }
509      fprintf(fp,"\n");
510    
511      for (j=0;j<ninputs;j++) {     /* print the 2 norms */
512        fprintf(fp,"%-#18.8f    ",norm_2[j]);
513      }
514      fprintf(fp,"\n\n");
515      ascfree((char *)norm_1);
516      ascfree((char *)norm_2);
517    
518      for (i=0;i<noutputs;i++) {        /* print the scaled data */
519        for (j=0;j<ninputs;j++) {
520          fprintf(fp,"%-#18.8f   %-4d",dy_dx[i][j],i);
521        }
522        if (var_fixed(outputs[i]))
523          fprintf(fp,"    **fixed*** \n");
524        else
525          PUTC('\n',fp);
526      }
527      fprintf(fp,"\n");
528      fclose(fp);
529      return 0;
530    }
531    
532    
533    
534  /**  /**
535      This function is very similar to sensitivity_anal, execept,      This function is very similar to sensitivity_anal, execept,
# Line 379  EXTERN sensitivity_anal_all( Line 547  EXTERN sensitivity_anal_all(
547    
548      The results will be witten to standard out.      The results will be witten to standard out.
549      This function is more expensive from a a memory point of view,      This function is more expensive from a a memory point of view,
550      as we keep aroung a dense matrix n_outputs * n_inputs, but here      as we keep around a dense matrix n_outputs * n_inputs, but here
551      n_outputs may be *much* larger depending on problem size.      n_outputs may be *much* larger depending on problem size.
552  */  */
553  int sensitivity_anal_all( struct Instance *inst,  /* not used but will be */  int sensitivity_anal_all( struct Instance *inst,  /* not used but will be */
554               struct gl_list_t *arglist,          struct gl_list_t *arglist,
555               real64 step_length)          real64 step_length
556  {  ){
557    struct Instance *which_instance;      struct Instance *which_instance;
558    struct gl_list_t *branch2, *branch3;      struct gl_list_t *branch2, *branch3;
559    dof_t *dof;      dof_t *dof;
560    struct var_variable **inputs = NULL, **outputs = NULL;      struct var_variable **inputs = NULL, **outputs = NULL;
561    int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;      int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
562    real64 **dy_dx = NULL;      real64 **dy_dx = NULL;
563    struct var_variable **vp,**ptr;      struct var_variable **vp,**ptr;
564    slv_system_t sys = NULL;      slv_system_t sys = NULL;
565    long c;      long c;
566    int noutputs=0, ninputs;      int noutputs=0, ninputs;
567    var_filter_t vfilter;      var_filter_t vfilter;
568    
569    struct var_variable **new_inputs = NULL; /* optional stuff for variable      struct var_variable **new_inputs = NULL; /* optional stuff for variable
570                        * projection */                        * projection */
571    
572    linsolqr_system_t lqr_sys;    /* stuff for the linear system & matrix */      linsolqr_system_t lqr_sys;  /* stuff for the linear system & matrix */
573        mtx_matrix_t mtx;
574        int32 capacity;
575        real64 *scratch_vector = NULL;
576        int result=0;
577    
578        /* Ignore unused params */
579        (void)inst; (void)step_length;
580    
581        /*
582        * Call the presolve for the system. This should number variables
583        * and relations as well create and order the main Jacobian. The
584        * only potential problem that I see here is that presolve using
585        * the slv0 solver *only* recognizes solver vars. So that if one
586        * wanted to see the sensitivity of a *parameter*, it would not
587        * be possible. We will have to trap this in CheckArgs.
588        *
589        * Also the current version of ascend is fucked in how the var module
590        * handles variables and their numbering through the interface ptr
591        * crap.
592        */
593    
594        (void)NumberFreeVars(NULL);     /* used to re-init the system */
595        (void)NumberRels(NULL);     /* used to re-init the system */
596        which_instance = FetchElement(arglist,1,1);
597        sys = asc_sens_presolve(which_instance);
598        if (!sys) {
599            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in asc_sens_presolve\n");
600            result = 1;
601            goto error;
602        }
603        dof = slv_get_dofdata(sys);
604    
605        /*
606        * prepare the inputs list. We dont need the index of the new_inputs
607        * list. We will grab them later if necessary.
608        */
609        branch2 = gl_fetch(arglist,2); /* input args -- old u_values */
610        branch3 = gl_fetch(arglist,3); /* input args -- new u_values */
611        ninputs = gl_length(branch2);
612        inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
613        new_inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
614    
615        inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
616        for (c=0;c<ninputs;c++) {
617            inputs[c] = (struct var_variable *)gl_fetch(branch2,c+1);
618            inputs_ndx_list[c] = var_mindex(inputs[c]);
619            new_inputs[c] = (struct var_variable *)gl_fetch(branch3,c+1);
620        }
621        inputs[ninputs] = NULL;     /* null terminate the list */
622        new_inputs[ninputs] = NULL;     /* null terminate the list */
623    
624        /*
625        * prepare the outputs list. This is where we differ from
626        * the other function. The noutputs, and the indexes of these
627        * outputs is obtained from the entire solve system.
628        */
629        vfilter.matchbits = 0;
630        noutputs = slv_count_solvers_vars(sys,&vfilter);
631        outputs = ASC_NEW_ARRAY(struct var_variable *,noutputs+1);
632        outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
633        ptr = vp = slv_get_solvers_var_list(sys);
634        for (c=0;c<noutputs;c++) {
635            outputs[c] = *ptr;
636            outputs_ndx_list[c] = var_sindex(*ptr);
637            ptr++;
638        }
639        outputs[noutputs] = NULL; /* null terminate the list */
640    
641        /*
642        * prepare the results dy_dx. This is the expensive part from a
643        * memory point of view. However I would like to have the entire
644        * noutputs * ninputs matrix even for a short while so that I
645        * can compute a number of  different types of norms.
646        */
647        dy_dx = make_matrix(noutputs,ninputs);
648    
649        result = Compute_J(sys);
650        if (result) {
651            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in calc Jacobian\n");
652            goto error;
653        }
654    
655        /*
656            Note: the RHS *has* to added here. We will construct the vector
657            of size matrix capacity and add it. It will be removed after
658            we finish computing dy_dx.
659        */
660        lqr_sys = slv_get_linsolqr_sys(sys);    /* get the linear system */
661        mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
662        capacity = mtx_capacity(mtx);
663        scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
664        linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
665        result = LUFactorJacobian1(sys,dof->structural_rank);
666    
667        if (result) {
668            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in LUFactorJacobian\n");
669            goto error;
670        }
671    
672        result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
673            inputs_ndx_list, ninputs,
674            outputs_ndx_list, noutputs
675        );
676    
677        linsolqr_remove_rhs(lqr_sys,scratch_vector);
678        if (result) {
679            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Early termination due to failure in Compute_dy_dx\n");
680            goto error;
681        }
682    
683        /*
684        * Do some analysis on the results, inclusive of
685        * writing them to someplace useful.
686        */
687        if(DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx)){
688            result = 1;
689        }
690        /*
691        * Some experimental projection stuff -- not used now.
692        * if (DoProject_X(inputs, new_inputs, step_length,
693        *     outputs, ninputs, noutputs, dy_dx))
694        * result = 1;
695        */
696    
697    error:
698        if (inputs) ascfree((char *)inputs);
699        if (new_inputs) ascfree((char *)new_inputs);
700        if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
701        if (outputs) ascfree((char *)outputs);
702        if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
703        if (dy_dx) free_matrix(dy_dx,noutputs);
704        if (scratch_vector) ascfree((char *)scratch_vector);
705        if (sys) system_destroy(sys);
706        return result;
707    }
708    
709    
710    
711    
712    #if 0
713    static int ComputeInverse(slv_system_t sys,
714                  real64 *rhs)
715    {
716      linsolqr_system_t lqr_sys;
717    mtx_matrix_t mtx;    mtx_matrix_t mtx;
718    int32 capacity;    int capacity,order;
719    real64 *scratch_vector = NULL;    real64 *solution = NULL;
720    int result=0;    int j,k;
721    
722    /* Ignore unused params */    lqr_sys = slv_get_linsolqr_sys(sys);  /* get the linear system */
723    (void)inst; (void)step_length;    mtx = linsolqr_get_matrix(lqr_sys);       /* get the matrix */
724    
725    /*    capacity = mtx_capacity(mtx);
726     * Call the presolve for the system. This should number variables    zero_vector(rhs,capacity);        /* zero the rhs */
727     * and relations as well create and order the main Jacobian. The    solution = ASC_NEW_ARRAY_CLEAR(real64,capacity);
    * 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.  
    */  
728    
729    (void)NumberFreeVars(NULL);       /* used to re-init the system */    order = mtx_order(mtx);
730    (void)NumberRels(NULL);       /* used to re-init the system */    for (j=0;j<order;j++) {
731    which_instance = FetchElement(arglist,1,1);      rhs[j] = 1.0;
732    sys = asc_sens_presolve(which_instance);      linsolqr_rhs_was_changed(lqr_sys,rhs);
733    if (!sys) {      linsolqr_solve(lqr_sys,rhs);            /* solve */
734      FPRINTF(stderr,"Early termination due to failure in asc_sens_presolve\n");      linsolqr_copy_solution(lqr_sys,rhs,solution);   /* get the solution */
735      result = 1;  
736      goto error;      CONSOLE_DEBUG("This is the rhs and solution for rhs #%d\n",j);
737        for (k=0;k<order;k++) {
738          CONSOLE_DEBUG("%12.8g  %12.8g\n",rhs[k],solution[k]);
739        }
740        rhs[j] = 0.0;
741    }    }
742    dof = slv_get_dofdata(sys);    if (solution) ascfree((char *)solution);
743      return 0;
744    }
745    #endif
746    
747    /*  #if 0
748     * prepare the inputs list. We dont need the index of the new_inputs  static int DoProject_X(struct var_variable **old_inputs,
749     * list. We will grab them later if necessary.                 struct var_variable **new_inputs, /* new values of u */
750     */                 double step_length,
751    branch2 = gl_fetch(arglist,2); /* input args -- old u_values */                 struct var_variable **outputs,
752    branch3 = gl_fetch(arglist,3); /* input args -- new u_values */                 int ninputs, int noutputs,
753    ninputs = gl_length(branch2);                 real64 **dy_dx)
754    inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);  {
755    new_inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);    struct var_variable *var;
756      real64 old_y, new_y, tmp;
757    inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);    real64 *delta_x;
758    for (c=0;c<ninputs;c++) {    int i,j;
     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 */  
759    
760    /*    delta_x = ASC_NEW_ARRAY_CLEAR(real64,ninputs);
    * 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 */  
761    
762    /*    for (j=0;j<ninputs;j++) {
763     * prepare the results dy_dx. This is the expensive part from a      delta_x[j] = var_value(new_inputs[j]) - var_value(old_inputs[j]);
764     * memory point of view. However I would like to have the entire      /*    delta_x[j] = RealAtomValue(new_inputs[j]) - RealAtomValue(old_inputs[j]); */
765     * 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);  
766    
767    result = Compute_J(sys);    for (i=0;i<noutputs;i++) {
768    if (result) {      var = outputs[i];
769      FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");      if (var_fixed(var) || !var_active(var))    /* project only the free vars */
770      goto error;        continue;
771        tmp = 0.0;
772        for (j=0;j<ninputs;j++) {
773          tmp += (dy_dx[i][j] * delta_x[j]);
774        }
775        /*    old_y = RealAtomValue(var); */
776        old_y = var_value(var);
777        new_y = old_y + step_length*tmp;
778        /*    SetRealAtomValue(var,new_y,(unsigned)0);  */
779        var_set_value(var,new_y);
780    # \if  DEBUG
781        CONSOLE_DEBUG("Old_y = %12.8g; Nex_y = %12.8g\n",old_y,new_y);
782    # \endif
783    }    }
784      ascfree((char *)delta_x);
785      return 0;
786    }
787    #endif
788    
789    
790    #if 0
791    /**
792     * At this point we should have an empty jacobian. We now
793     * need to call relman_diff over the *entire* matrix.
794     * Calculates the entire jacobian. It is initially unscaled.
795     *
796     * Note: this assumes the sys given is using one of the ascend solvers
797     * with either linsol or linsolqr. Both are now allowed. baa 7/95
798     */
799    #define SAFE_FIX_ME 0
800    static int Compute_J_OLD(slv_system_t sys)
801    {
802      int32 row;
803      var_filter_t vfilter;
804      linsol_system_t lin_sys = NULL;
805      linsolqr_system_t lqr_sys = NULL;
806      mtx_matrix_t mtx;
807      struct rel_relation **rlist;
808      int nrows;
809      int qr=0;
810    #\if DOTIME
811      double time1;
812    #\endif
813    
814    #\if DOTIME
815      time1 = tm_cpu_time();
816    #\endif
817    /*    /*
818     * Note: the RHS *has* to added here. We will construct the vector     * Get the linear system from the solve system.
819     * of size matrix capacity and add it. It will be removed after     * Get the matrix from the linear system.
820     * we finish computing dy_dx.     * Get the relations list from the solve system.
821     */     */
822    lqr_sys = slv_get_linsolqr_sys(sys);  /* get the linear system */    lin_sys = slv_get_linsol_sys(sys);
823    mtx = linsolqr_get_matrix(lqr_sys);   /* get the matrix */    if (lin_sys==NULL) {
824    capacity = mtx_capacity(mtx);      qr=1;
825    scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);      lqr_sys=slv_get_linsolqr_sys(sys);
826    linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);    }
827    result = LUFactorJacobian1(sys,dof->structural_rank);    mtx = slv_get_sys_mtx(sys);
828    if (result) {    if (mtx==NULL || (lin_sys==NULL && lqr_sys==NULL)) {
829      FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Compute_dy_dx: error, found NULL.\n");
830      goto error;      return 1;
   }  
   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;  
831    }    }
832      rlist = slv_get_solvers_rel_list(sys);
833      nrows = NumberIncludedRels(sys);
834    
835      calc_ok = TRUE;
836      vfilter.matchbits = VAR_SVAR;
837      vfilter.matchvalue = VAR_SVAR;
838    /*    /*
839     * Do some analysis on the results, inclusive of     * Clear the entire matrix and then compute
840     * writing them to someplace useful.     * the gradients at the current point.
841     */     */
842    if (DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx))    mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
843      result = 1;    for(row=0; row<nrows; row++) {
844        struct rel_relation *rel;
845        /* added  */
846        double resid;
847    
848        rel = rlist[mtx_row_to_org(mtx,row)];
849        (void)relman_diffs(rel,&vfilter,mtx,&resid,SAFE_FIX_ME);
850    
851        /* added  */
852        rel_set_residual(rel,resid);
853    
854      }
855      if (qr) {
856        linsolqr_matrix_was_changed(lqr_sys);
857      } else {
858        linsol_matrix_was_changed(lin_sys);
859      }
860    #\if DOTIME
861      time1 = tm_cpu_time() - time1;
862      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Time taken for Compute_J = %g\n",time1);
863    #\endif
864      return(!calc_ok);
865    }
866    #endif
867    
   /*  
    * Some experimental projection stuff -- not used now.  
    * if (DoProject_X(inputs, new_inputs, step_length,  
    *     outputs, ninputs, noutputs, dy_dx))  
    * result = 1;  
    */  
868    
869   error:  #if 0
870    if (inputs) ascfree((char *)inputs);  static int ReSolve(slv_system_t sys)
871    if (new_inputs) ascfree((char *)new_inputs);  {
872    if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);    if (!sys)
873    if (outputs) ascfree((char *)outputs);      return 1;
874    if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);    slv_solve(sys);
875    if (dy_dx) free_matrix(dy_dx,noutputs);    return 0;
   if (scratch_vector) ascfree((char *)scratch_vector);  
   if (sys) system_destroy(sys);  
   return result;  
876  }  }
877    #endif
878    
879    /**
880        Build then presolve the solve an instance...
881    */
882    int DoSolve(struct Instance *inst){
883      slv_system_t sys;
884    
885      sys = system_build(inst);
886      if (!sys) {
887        ERROR_REPORTER_HERE(ASC_PROG_ERR,
888          "Something radically wrong in creating solver.");
889        return 1;
890      }
891      (void)slv_select_solver(sys,0);
892      slv_presolve(sys);
893      slv_solve(sys);
894      system_destroy(sys);
895      return 0;
896    }
897    
898  /**  /**
899      Calls 'DoSolve'      Calls 'DoSolve'
# Line 543  int sensitivity_anal_all( struct Instanc Line 901  int sensitivity_anal_all( struct Instanc
901      @see DoSolve      @see DoSolve
902  */  */
903  int do_solve_eval( struct Instance *i,  int do_solve_eval( struct Instance *i,
904            struct gl_list_t *arglist, void *user_data)          struct gl_list_t *arglist, void *user_data
905  {  ){
906    unsigned long len;    unsigned long len;
907    int result;    int result;
908    struct Instance *inst;    struct Instance *inst;
909    len = gl_length(arglist);    len = gl_length(arglist);
910    
911    /* Ignore unused params */    (void)i; /* not used */
   (void)i;  
912    
913    if (len!=2) {    if (len!=2) {
914      ERROR_REPORTER_HERE(ASC_USER_ERROR,      ERROR_REPORTER_HERE(ASC_USER_ERROR,
# Line 567  int do_solve_eval( struct Instance *i, Line 924  int do_solve_eval( struct Instance *i,
924    
925    
926  /**  /**
927      Finite different evaluate...      Finite difference...
928  */  */
929  int do_finite_diff_eval( struct Instance *i,  int finite_difference(struct gl_list_t *arglist){
930               struct gl_list_t *arglist, void *user_data)    struct Instance *model_inst, *xinst, *inst;
931  {    slv_system_t sys;
932    int result;    int ninputs,noutputs;
933      int i,j,offset;
934    /* Ignore unused params */    real64 **partials;
935    (void)i;    real64 *y_old, *y_new;
936      real64 x;
937      real64 interval = 1e-6;
938      int result=0;
939    
940    if (FiniteDiffCheckArgs(arglist))    model_inst = FetchElement(arglist,1,1);
941      sys = system_build(model_inst);
942      if (!sys) {
943        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Something radically wrong in creating solver\n");
944      return 1;      return 1;
945    result = finite_difference(arglist);    }
946      (void)slv_select_solver(sys,0);
947      slv_presolve(sys);
948      slv_solve(sys);
949    
950      /* Make the necessary vectors */
951    
952      ninputs = (int)gl_length((struct gl_list_t *)gl_fetch(arglist,2));
953      noutputs = (int)gl_length((struct gl_list_t *)gl_fetch(arglist,3));
954      y_old = (real64 *)calloc(noutputs,sizeof(real64));
955      y_new = (real64 *)calloc(noutputs,sizeof(real64));
956      partials = make_matrix(noutputs,ninputs);
957      for (i=0;i<noutputs;i++) {        /* get old yvalues */
958        inst = FetchElement(arglist,3,i+1);
959        y_old[i] = RealAtomValue(inst);
960      }
961      for (j=0;j<ninputs;j++) {
962        xinst = FetchElement(arglist,2,j+1);
963        x = RealAtomValue(xinst);
964        SetRealAtomValue(xinst,x+interval,(unsigned)0); /* perturb system */
965        slv_presolve(sys);
966        slv_solve(sys);
967    
968        for (i=0;i<noutputs;i++) {      /* get new yvalues */
969          inst = FetchElement(arglist,3,i+1);
970          y_new[i] = RealAtomValue(inst);
971          partials[i][j] = -1.0 * (y_old[i] - y_new[i])/interval;
972          PRINTF("y_old = %20.12g  y_new = %20.12g\n", y_old[i],y_new[i]);
973        }
974        SetRealAtomValue(xinst,x,(unsigned)0); /* unperturb system */
975      }
976      offset = 0;
977      for (i=0;i<noutputs;i++) {
978        for (j=0;j<ninputs;j++) {
979          inst = FetchElement(arglist,4,offset+j+1);
980          SetRealAtomValue(inst,partials[i][j],(unsigned)0);
981          PRINTF("%12.6f %s",partials[i][j], (j==(ninputs-1)) ? "\n" : "    ");
982        }
983        offset += ninputs;
984      }
985      /* error: */
986      free(y_old);
987      free(y_new);
988      free_matrix(partials,noutputs);
989      system_destroy(sys);
990    return result;    return result;
991  }  }
992    
993    
994    
995    /**
996        Finite different evaluate...
997    */
998    int do_finite_diff_eval( struct Instance *i,
999            struct gl_list_t *arglist, void *user_data
1000    ){
1001        int result;
1002        (void)i; /* not used */
1003    
1004        if(FiniteDiffCheckArgs(arglist))
1005            return 1;
1006        result = finite_difference(arglist);
1007        return result;
1008    }
1009    
1010    
1011  int do_sensitivity_eval( struct Instance *i,  int do_sensitivity_eval( struct Instance *i,
1012               struct gl_list_t *arglist, void *user_data)          struct gl_list_t *arglist, void *user_data
1013  {  ){
1014    int result;      CONSOLE_DEBUG("Starting sensitivity analysis...");
1015    if (SensitivityCheckArgs(arglist)) {      if(SensitivityCheckArgs(arglist))return 1;
1016      return 1;  
1017    }      return sensitivity_anal(i,arglist);
   result = sensitivity_anal(i,arglist);  
   return result;  
1018  }  }
1019    
1020  int do_sensitivity_eval_all( struct Instance *i,  int do_sensitivity_eval_all( struct Instance *i,
1021                  struct gl_list_t *arglist, void *user_data)          struct gl_list_t *arglist, void *user_data
1022  {  ){
1023    int result;      int result;
1024    double step_length = 0.0;      double step_length = 0.0;
1025    if (SensitivityAllCheckArgs(arglist,&step_length)) {      if (SensitivityAllCheckArgs(arglist,&step_length)) {
1026      return 1;          return 1;
1027    }      }
1028    result = sensitivity_anal_all(i,arglist,step_length);      result = sensitivity_anal_all(i,arglist,step_length);
1029    return result;      return result;
1030  }  }
1031    
1032    
1033  char sensitivity_help[] =  const char sensitivity_help[] =
1034      "This function does sensitivity analysis dy/dx. It requires 4 args:\n"      "This function does sensitivity analysis dy/dx. It requires 4 args:\n"
1035      "  1. name: name of a reference instance or SELF.\n"      "  1. name: name of a reference instance or SELF.\n"
1036      "  2. x: x, where x is an array of > solver_var.\n"      "  2. x: x, where x is an array of > solver_var.\n"
# Line 616  char sensitivity_help[] = Line 1038  char sensitivity_help[] =
1038      "  4. dy/dx: which dy_dx[1..n_y][1..n_x].";      "  4. dy/dx: which dy_dx[1..n_y][1..n_x].";
1039    
1040  int sensitivity_register(void){  int sensitivity_register(void){
1041        int result=0;
1042    
1043    int result=0;      result = CreateUserFunctionMethod("do_solve",
1044            do_solve_eval,
1045            2,NULL,NULL,NULL
1046    result = CreateUserFunctionMethod("do_solve",      );
1047                    do_solve_eval,      result += CreateUserFunctionMethod("do_finite_difference",
1048                    2,NULL,NULL,NULL); /* was 2,0,null */          do_finite_diff_eval,
1049    result += CreateUserFunctionMethod("do_finite_difference",          4,NULL,NULL,NULL
1050                     do_finite_diff_eval,      );
1051                     4,NULL,NULL,NULL); /* 4,0,null */      result += CreateUserFunctionMethod("do_sensitivity",
1052    result += CreateUserFunctionMethod("do_sensitivity",          do_sensitivity_eval,
1053                     do_sensitivity_eval,          4,sensitivity_help,NULL,NULL
1054                     4,sensitivity_help,NULL,NULL);      );
1055    result += CreateUserFunctionMethod("do_sensitivity_all",      result += CreateUserFunctionMethod("do_sensitivity_all",
1056                     do_sensitivity_eval_all,          do_sensitivity_eval_all,
1057                     4,"See do_sensitivity_eval for details",NULL,NULL);          4,"See do_sensitivity_eval for details",NULL,NULL
1058        );
1059    
1060    return result;      return result;
1061  }  }
1062    
1063    /* :ex: set ts=4: */
Legend:
Removed from v.1162  
changed lines
  Added in v.1163
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22