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

Annotation of /trunk/models/sensitivity/sensitivity.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2064 - (hide annotations) (download) (as text)
Mon Jul 6 12:03:36 2009 UTC (11 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 25378 byte(s)
Trying to solve some bugs with sensitivity external method.
1 johnpye 1163 /* 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 johnpye 1174
14 johnpye 1163 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 johnpye 1164 The module permits sensitivity analysis in an ASCEND model. In your model
21     file, you need to create some additional arrays like so:
22    
23     X[1..m] IS_A solver_var; (* outputs *)
24     U[1..n] IS_A solver_var; (* inputs *)
25     dx_du[1..m][1..n] IS_A solver_var; (* data space *)
26    
27     You must then 'ARE_THE_SAME' your problem variables into the spaces in the
28     first two of the above arrays.
29    
30     Then, run EXTERNAL do_sensitivity(SELF,U[1..n],X[1..m],dx_du[1..m][1..n]);
31    
32     In the array dx_du[1..m][1..n], values of the corresponding sensitivity
33     derivatives will be calculated.
34    
35     A test/example model is available in the models directory, called
36     sensitivity_test.a4c.
37    
38     --- implementation notes ---
39    
40 johnpye 1163 This file attempts to implement the extraction of dy_dx from
41     a system of equations. If one considers a black-box where x are
42     the input variables, u are input parameters, y are the output
43     variables, f(x,y,u) is the system of equations that will be solved
44     for given x and u, then one can extract the sensitivity information
45     of y wrt x.
46    
47     One crude but simple way of doing this is to finite-difference the
48     given black box, i.e, perturb x, n_x times by delta x, resolve
49     the system of equations at the new value of x, and compute
50    
51     dy/dx = (f(x_1) - f(x_2))/(x_1 - x_2).
52    
53     This is known to be very expensive.
54    
55     The solution that will be adopted here is to formulate the Jacobian J of
56     the system, (or the system is already solved, to grab the jacobian at
57     the solution. Develop the sensitivity matrix df/dx by exact numnerical
58     differentiation; this will be a n x n_x matrix. Compute the LU factors
59     for J, and to solve column by column to : LU dz/dx = df/dx. Here
60     z, represents all the internal variables to the system and the strictly
61     output variables y. In other words J = df/dz.
62    
63     Given the solution of the above problem we then simply extract the rows
64     of dz/dx, that pertain to the y variables that we are interested in;
65     this will give us dy/dx.
66    
67 johnpye 1164 @todo There are files in tcltk called Sensitivity.[ch]. Do we need them?
68 johnpye 1163 */
69    
70 johnpye 1162 #include <math.h>
71    
72 jpye 2018 #include <ascend/utilities/ascMalloc.h>
73     #include <ascend/general/mathmacros.h>
74 johnpye 1183
75 jpye 2018 #include <ascend/compiler/instquery.h>
76     #include <ascend/compiler/atomvalue.h>
77     #include <ascend/compiler/extfunc.h>
78 johnpye 1162
79 jpye 2018 #include <ascend/linear/densemtx.h>
80 johnpye 1183
81 jpye 2018 #include <ascend/packages/sensitivity.h>
82     #include <ascend/system/system.h>
83     #include <ascend/solver/solver.h>
84 johnpye 1183
85 jpye 2064 #define SENSITIVITY_DEBUG
86 johnpye 1162
87 johnpye 1163 ASC_EXPORT int sensitivity_register(void);
88    
89     ExtMethodRun do_sensitivity_eval;
90     ExtMethodRun do_sensitivity_eval_all;
91    
92 johnpye 1162 /**
93     Build then presolve an instance
94     */
95 johnpye 1164 slv_system_t sens_presolve(struct Instance *inst){
96 johnpye 1162 slv_system_t sys;
97     slv_parameters_t parameters;
98 jpye 1487 const SlvFunctionsT *S;
99 johnpye 1163 #ifdef SENSITIVITY_DEBUG
100 johnpye 1162 struct var_variable **vp;
101     struct rel_relation **rp;
102 jpye 2064 int len, ind;
103 johnpye 1162 char *tmp=NULL;
104 johnpye 1163 #endif
105 johnpye 1162 sys = system_build(inst);
106     if (sys==NULL) {
107 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to build system.");
108 johnpye 1162 return NULL;
109     }
110 jpye 1487 S = solver_engine_named("QRSlv");
111     if(!S){
112     ERROR_REPORTER_HERE(ASC_PROG_ERR,"QRSlv solver not found (required for sensitivity)");
113     return NULL;
114 johnpye 1162 }
115 jpye 1487 slv_select_solver(sys,S->number);
116 johnpye 1162
117     slv_get_parameters(sys,&parameters);
118     parameters.partition = 0;
119     slv_set_parameters(sys,&parameters);
120     slv_presolve(sys);
121    
122 johnpye 1163 #ifdef SENSITIVITY_DEBUG
123 johnpye 1162 vp = slv_get_solvers_var_list(sys);
124     len = slv_get_num_solvers_vars(sys);
125     for (ind=0 ; ind<len; ind++) {
126     tmp = var_make_name(sys,vp[ind]);
127 johnpye 1163 CONSOLE_DEBUG("%s %d\n",tmp,var_sindex(vp[ind]));
128 johnpye 1162 if (tmp!=NULL) ascfree(tmp);
129     }
130     rp = slv_get_solvers_rel_list(sys);
131     len = slv_get_num_solvers_rels(sys);
132     for (ind=0 ; ind<len ; ind++) {
133     tmp = rel_make_name(sys,rp[ind]);
134 johnpye 1163 CONSOLE_DEBUG("%s %d\n",tmp,rel_sindex(rp[ind]));
135 johnpye 1162 if (tmp) ascfree(tmp);
136     }
137     #endif
138     return sys;
139     }
140    
141 johnpye 1163 /*
142     LU Factor Jacobian
143 johnpye 1162
144 johnpye 1163 @note The RHS will have been will already have been added before
145     calling this function.
146    
147     @NOTE there is another version of this floating around in packages/senstivity.c. The
148     other one uses dense matrices so probably this one's better?
149     */
150     int LUFactorJacobian1(slv_system_t sys,int rank){
151     linsolqr_system_t lqr_sys;
152     mtx_region_t region;
153     enum factor_method fm;
154    
155     mtx_region(&region,0,rank-1,0,rank-1); /* set the region */
156     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
157    
158     linsolqr_matrix_was_changed(lqr_sys);
159     linsolqr_reorder(lqr_sys,&region,natural);
160    
161     fm = linsolqr_fmethod(lqr_sys);
162     if (fm == unknown_f) fm = ranki_kw2; /* make sure somebody set it */
163     linsolqr_factor(lqr_sys,fm);
164    
165     return 0;
166     }
167    
168 johnpye 1162 int sensitivity_anal(
169 johnpye 1163 struct Instance *inst, /* not used but will be */
170     struct gl_list_t *arglist
171     ){
172     struct Instance *which_instance,*tmp_inst, *atominst;
173     struct gl_list_t *branch;
174     struct var_variable **vlist = NULL;
175     int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
176 jpye 2064 DenseMatrix dy_dx = densematrix_create_empty();
177 johnpye 1163 slv_system_t sys = NULL;
178     int c;
179     int noutputs = 0;
180     int ninputs;
181     int i,j;
182     int offset;
183     dof_t *dof;
184     int num_vars,ind,found;
185 johnpye 1162
186 johnpye 1163 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
187     mtx_matrix_t mtx;
188     int32 capacity;
189     real64 *scratch_vector = NULL;
190     int result=0;
191 johnpye 1162
192 johnpye 1163 /* Ignore unused params */
193     (void) inst;
194 johnpye 1162
195 johnpye 1163 (void)NumberFreeVars(NULL); /* used to re-init the system */
196     (void)NumberRels(NULL); /* used to re-init the system */
197     which_instance = FetchElement(arglist,1,1);
198 johnpye 1164 sys = sens_presolve(which_instance);
199 johnpye 1163 if (!sys) {
200 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to presolve");
201 johnpye 1163 result = 1;
202 johnpye 1164 goto finish;
203 johnpye 1163 }
204 johnpye 1162
205 johnpye 1163 dof = slv_get_dofdata(sys);
206     if (!(dof->n_rows == dof->n_cols &&
207     dof->n_rows == dof->structural_rank)) {
208 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"System is not square");
209 johnpye 1163 result = 1;
210 johnpye 1164 goto finish;
211 johnpye 1163 }
212 johnpye 1162
213 johnpye 1163 CONSOLE_DEBUG("Presolved, square");
214 johnpye 1162
215 johnpye 1163 /*
216     prepare the inputs list
217     */
218     vlist = slv_get_solvers_var_list(sys);
219 johnpye 1162
220 johnpye 1163 branch = gl_fetch(arglist,2); /* input args */
221     ninputs = gl_length(branch);
222     inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
223 johnpye 1162
224 johnpye 1163 num_vars = slv_get_num_solvers_vars(sys);
225     for (c=0;c<ninputs;c++) {
226     atominst = (struct Instance *)gl_fetch(branch,c+1);
227     found = 0;
228     ind = num_vars - 1;
229     /* search backwards because fixed vars are at the end of the
230     var list */
231     while (!found && ind >= 0){
232     if (var_instance(vlist[ind]) == atominst) {
233     inputs_ndx_list[c] = var_sindex(vlist[ind]);
234     found = 1;
235     }
236     --ind;
237     }
238 johnpye 1164 if(!found){
239     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to find sensitivity input variable");
240 johnpye 1163 result = 1;
241 johnpye 1164 goto finish;
242 johnpye 1163 }
243     }
244 johnpye 1162
245 johnpye 1163 CONSOLE_DEBUG("%d inputs",ninputs);
246 johnpye 1162
247 johnpye 1163 /*
248     prepare the outputs list
249     */
250     branch = gl_fetch(arglist,3); /* output args */
251     noutputs = gl_length(branch);
252     outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
253     for (c=0;c<noutputs;c++) {
254     atominst = (struct Instance *)gl_fetch(branch,c+1);
255     found = 0;
256     ind = 0;
257     while (!found && ind < num_vars){
258     if (var_instance(vlist[ind]) == atominst) {
259     outputs_ndx_list[c] = var_sindex(vlist[ind]);
260     found = 1;
261     }
262     ++ind;
263     }
264     if (!found) {
265 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to find sensitivity ouput variable");
266 johnpye 1163 result = 1;
267 johnpye 1164 goto finish;
268 johnpye 1163 }
269     }
270 johnpye 1174
271 johnpye 1163 CONSOLE_DEBUG("%d outputs",noutputs);
272 johnpye 1162
273 johnpye 1163 /*
274     prepare the results dy_dx.
275     */
276 johnpye 1174 dy_dx = densematrix_create(noutputs,ninputs);
277 johnpye 1162
278 johnpye 1163 result = Compute_J(sys);
279     if (result) {
280 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to calculate Jacobian failure in calc Jacobian\n");
281     goto finish;
282 johnpye 1163 }
283 johnpye 1162
284 johnpye 1163 CONSOLE_DEBUG("Computed Jacobian");
285    
286     /*
287     Note: the RHS *has* to added here. We will construct the vector
288     of size matrix capacity and add it. It will be removed after
289     we finish computing dy_dx.
290     */
291     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
292     mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
293     capacity = mtx_capacity(mtx);
294     scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
295     linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
296     result = LUFactorJacobian1(sys,dof->structural_rank);
297     if(result){
298 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to calculate LUFactorJacobian");
299     goto finish;
300 johnpye 1163 }
301     result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
302     inputs_ndx_list, ninputs,
303     outputs_ndx_list, noutputs
304     );
305    
306     linsolqr_remove_rhs(lqr_sys,scratch_vector);
307     if (result) {
308 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed Compute_dy_dx");
309     goto finish;
310 johnpye 1163 }
311 johnpye 1174
312 johnpye 1163 CONSOLE_DEBUG("Computed dy/dx");
313    
314     /*
315     Write the results back to the partials array in the
316     instance tree
317     */
318     offset = 0;
319     for (i=0;i<noutputs;i++) {
320     for (j=0;j<ninputs;j++) {
321     tmp_inst = FetchElement(arglist,4,offset+j+1);
322 johnpye 1174 SetRealAtomValue(tmp_inst,DENSEMATRIX_ELEM(dy_dx,i,j),(unsigned)0);
323 johnpye 1163 #ifdef SENSITIVITY_DEBUG
324 johnpye 1174 CONSOLE_DEBUG("%12.8f i%d j%d",DENSEMATRIX_ELEM(dy_dx,i,j),i,j);
325 johnpye 1162 #endif
326 johnpye 1163 }
327     #ifdef SENSITIVITY_DEBUG
328     CONSOLE_DEBUG("\n");
329 johnpye 1162 #endif
330 johnpye 1163 offset += ninputs;
331     }
332     #ifdef SENSITIVITY_DEBUG
333     CONSOLE_DEBUG("\n");
334 johnpye 1162 #endif
335    
336 johnpye 1163 ERROR_REPORTER_HERE(ASC_USER_SUCCESS
337     ,"Sensitivity results for %d vars were written back to the model"
338     ,noutputs
339     );
340    
341 johnpye 1164 finish:
342 johnpye 1163 if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
343     if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
344 johnpye 1174 densematrix_destroy(dy_dx);
345 johnpye 1163 if (scratch_vector) ascfree((char *)scratch_vector);
346     if (sys) system_destroy(sys);
347     return result;
348 johnpye 1162 }
349    
350     /**
351 johnpye 1164 Do Data Analysis. Used by sensitivity_anal_all.
352 johnpye 1162 */
353 johnpye 1163 int DoDataAnalysis(struct var_variable **inputs,
354     struct var_variable **outputs,
355     int ninputs, int noutputs,
356 johnpye 1174 DenseMatrix dy_dx)
357 johnpye 1163 {
358     FILE *fp;
359     double *norm_2, *norm_1;
360     double input_nominal,maxvalue,sum;
361     int i,j;
362 johnpye 1162
363 johnpye 1163 norm_1 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
364     norm_2 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
365    
366     fp = fopen("sensitivity.out","w+");
367     if (!fp) return 0;
368    
369     /*
370     * calculate the 1 and 2 norms; cache them away so that we
371     * can pretty print them. Style is everything !.
372     *
373     */
374     for (j=0;j<ninputs;j++) {
375     input_nominal = var_nominal(inputs[j]);
376     maxvalue = sum = 0;
377     for (i=0;i<noutputs;i++) {
378 johnpye 1174 DENSEMATRIX_ELEM(dy_dx,i,j) *= input_nominal/var_nominal(outputs[i]);
379     maxvalue = MAX(fabs(DENSEMATRIX_ELEM(dy_dx,i,j)),maxvalue);
380     sum += DENSEMATRIX_ELEM(dy_dx,i,j)*DENSEMATRIX_ELEM(dy_dx,i,j);
381 johnpye 1163 }
382     norm_1[j] = maxvalue;
383     norm_2[j] = sum;
384     }
385    
386     for (j=0;j<ninputs;j++) { /* print the var_index */
387     fprintf(fp,"%8d ",var_mindex(inputs[j]));
388     }
389     fprintf(fp,"\n");
390    
391     for (j=0;j<ninputs;j++) { /* print the 1 norms */
392     fprintf(fp,"%-#18.8f ",norm_1[j]);
393     }
394     fprintf(fp,"\n");
395    
396     for (j=0;j<ninputs;j++) { /* print the 2 norms */
397     fprintf(fp,"%-#18.8f ",norm_2[j]);
398     }
399     fprintf(fp,"\n\n");
400     ascfree((char *)norm_1);
401     ascfree((char *)norm_2);
402    
403     for (i=0;i<noutputs;i++) { /* print the scaled data */
404     for (j=0;j<ninputs;j++) {
405 johnpye 1174 fprintf(fp,"%-#18.8f %-4d",DENSEMATRIX_ELEM(dy_dx,i,j),i);
406 johnpye 1163 }
407     if (var_fixed(outputs[i]))
408     fprintf(fp," **fixed*** \n");
409     else
410     PUTC('\n',fp);
411     }
412     fprintf(fp,"\n");
413     fclose(fp);
414     return 0;
415     }
416    
417    
418 johnpye 1162 /**
419     This function is very similar to sensitivity_anal, execept,
420     that it perform the analysis on the entire system, based on
421     the given inputs. Nothing is returned. As such the call from
422     ASCEND is:
423    
424     <pre>
425     EXTERN sensitivity_anal_all(
426     this_instance,
427     u_old[1..n_inputs],
428     u_new[1..n_inputs],
429     step_length
430     );</pre>
431    
432     The results will be witten to standard out.
433     This function is more expensive from a a memory point of view,
434 johnpye 1163 as we keep around a dense matrix n_outputs * n_inputs, but here
435 johnpye 1162 n_outputs may be *much* larger depending on problem size.
436     */
437     int sensitivity_anal_all( struct Instance *inst, /* not used but will be */
438 johnpye 1163 struct gl_list_t *arglist,
439     real64 step_length
440     ){
441     struct Instance *which_instance;
442     struct gl_list_t *branch2, *branch3;
443     dof_t *dof;
444     struct var_variable **inputs = NULL, **outputs = NULL;
445     int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
446 johnpye 1174 DenseMatrix dy_dx;
447 johnpye 1163 struct var_variable **vp,**ptr;
448     slv_system_t sys = NULL;
449     long c;
450     int noutputs=0, ninputs;
451     var_filter_t vfilter;
452 johnpye 1162
453 johnpye 1163 struct var_variable **new_inputs = NULL; /* optional stuff for variable
454 johnpye 1162 * projection */
455    
456 johnpye 1163 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
457     mtx_matrix_t mtx;
458     int32 capacity;
459     real64 *scratch_vector = NULL;
460     int result=0;
461    
462     /* Ignore unused params */
463     (void)inst; (void)step_length;
464    
465     /*
466     * Call the presolve for the system. This should number variables
467     * and relations as well create and order the main Jacobian. The
468     * only potential problem that I see here is that presolve using
469     * the slv0 solver *only* recognizes solver vars. So that if one
470     * wanted to see the sensitivity of a *parameter*, it would not
471     * be possible. We will have to trap this in CheckArgs.
472     *
473     * Also the current version of ascend is fucked in how the var module
474     * handles variables and their numbering through the interface ptr
475     * crap.
476     */
477    
478     (void)NumberFreeVars(NULL); /* used to re-init the system */
479     (void)NumberRels(NULL); /* used to re-init the system */
480     which_instance = FetchElement(arglist,1,1);
481 johnpye 1164 sys = sens_presolve(which_instance);
482 johnpye 1163 if (!sys) {
483 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to presolve");
484 johnpye 1163 result = 1;
485 johnpye 1164 goto finish;
486 johnpye 1163 }
487     dof = slv_get_dofdata(sys);
488    
489     /*
490     * prepare the inputs list. We dont need the index of the new_inputs
491     * list. We will grab them later if necessary.
492     */
493     branch2 = gl_fetch(arglist,2); /* input args -- old u_values */
494     branch3 = gl_fetch(arglist,3); /* input args -- new u_values */
495     ninputs = gl_length(branch2);
496     inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
497     new_inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
498    
499     inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
500     for (c=0;c<ninputs;c++) {
501     inputs[c] = (struct var_variable *)gl_fetch(branch2,c+1);
502     inputs_ndx_list[c] = var_mindex(inputs[c]);
503     new_inputs[c] = (struct var_variable *)gl_fetch(branch3,c+1);
504     }
505     inputs[ninputs] = NULL; /* null terminate the list */
506     new_inputs[ninputs] = NULL; /* null terminate the list */
507    
508     /*
509     * prepare the outputs list. This is where we differ from
510     * the other function. The noutputs, and the indexes of these
511     * outputs is obtained from the entire solve system.
512     */
513     vfilter.matchbits = 0;
514     noutputs = slv_count_solvers_vars(sys,&vfilter);
515     outputs = ASC_NEW_ARRAY(struct var_variable *,noutputs+1);
516     outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
517     ptr = vp = slv_get_solvers_var_list(sys);
518     for (c=0;c<noutputs;c++) {
519     outputs[c] = *ptr;
520     outputs_ndx_list[c] = var_sindex(*ptr);
521     ptr++;
522     }
523     outputs[noutputs] = NULL; /* null terminate the list */
524    
525     /*
526     * prepare the results dy_dx. This is the expensive part from a
527     * memory point of view. However I would like to have the entire
528     * noutputs * ninputs matrix even for a short while so that I
529     * can compute a number of different types of norms.
530     */
531 johnpye 1174 dy_dx = densematrix_create(noutputs,ninputs);
532 johnpye 1163
533     result = Compute_J(sys);
534     if (result) {
535 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to compute Jacobian");
536     goto finish;
537 johnpye 1163 }
538    
539     /*
540     Note: the RHS *has* to added here. We will construct the vector
541     of size matrix capacity and add it. It will be removed after
542     we finish computing dy_dx.
543     */
544     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
545     mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
546     capacity = mtx_capacity(mtx);
547     scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
548     linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
549     result = LUFactorJacobian1(sys,dof->structural_rank);
550    
551     if (result) {
552 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failure in LUFactorJacobian");
553     goto finish;
554 johnpye 1163 }
555    
556     result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
557     inputs_ndx_list, ninputs,
558     outputs_ndx_list, noutputs
559     );
560    
561     linsolqr_remove_rhs(lqr_sys,scratch_vector);
562     if (result) {
563 johnpye 1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failure in Compute_dy_dx");
564     goto finish;
565 johnpye 1163 }
566    
567     /*
568     * Do some analysis on the results, inclusive of
569     * writing them to someplace useful.
570     */
571     if(DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx)){
572     result = 1;
573     }
574     /*
575     * Some experimental projection stuff -- not used now.
576     * if (DoProject_X(inputs, new_inputs, step_length,
577     * outputs, ninputs, noutputs, dy_dx))
578     * result = 1;
579     */
580    
581 johnpye 1164 finish:
582 johnpye 1163 if (inputs) ascfree((char *)inputs);
583     if (new_inputs) ascfree((char *)new_inputs);
584     if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
585     if (outputs) ascfree((char *)outputs);
586     if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
587 johnpye 1174 densematrix_destroy(dy_dx);;
588 johnpye 1163 if (scratch_vector) ascfree((char *)scratch_vector);
589     if (sys) system_destroy(sys);
590     return result;
591     }
592    
593    
594    
595    
596     #if 0
597     static int ComputeInverse(slv_system_t sys,
598     real64 *rhs)
599     {
600     linsolqr_system_t lqr_sys;
601 johnpye 1162 mtx_matrix_t mtx;
602 johnpye 1163 int capacity,order;
603     real64 *solution = NULL;
604     int j,k;
605 johnpye 1162
606 johnpye 1163 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
607     mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
608 johnpye 1162
609 johnpye 1163 capacity = mtx_capacity(mtx);
610     zero_vector(rhs,capacity); /* zero the rhs */
611     solution = ASC_NEW_ARRAY_CLEAR(real64,capacity);
612 johnpye 1162
613 johnpye 1163 order = mtx_order(mtx);
614     for (j=0;j<order;j++) {
615     rhs[j] = 1.0;
616     linsolqr_rhs_was_changed(lqr_sys,rhs);
617     linsolqr_solve(lqr_sys,rhs); /* solve */
618     linsolqr_copy_solution(lqr_sys,rhs,solution); /* get the solution */
619    
620     CONSOLE_DEBUG("This is the rhs and solution for rhs #%d\n",j);
621     for (k=0;k<order;k++) {
622     CONSOLE_DEBUG("%12.8g %12.8g\n",rhs[k],solution[k]);
623     }
624     rhs[j] = 0.0;
625 johnpye 1162 }
626 johnpye 1163 if (solution) ascfree((char *)solution);
627     return 0;
628     }
629     #endif
630 johnpye 1162
631 johnpye 1163 #if 0
632     static int DoProject_X(struct var_variable **old_inputs,
633     struct var_variable **new_inputs, /* new values of u */
634     double step_length,
635     struct var_variable **outputs,
636     int ninputs, int noutputs,
637 johnpye 1174 DenseMatrix dy_dx)
638 johnpye 1163 {
639     struct var_variable *var;
640     real64 old_y, new_y, tmp;
641     real64 *delta_x;
642     int i,j;
643 johnpye 1162
644 johnpye 1163 delta_x = ASC_NEW_ARRAY_CLEAR(real64,ninputs);
645    
646     for (j=0;j<ninputs;j++) {
647     delta_x[j] = var_value(new_inputs[j]) - var_value(old_inputs[j]);
648     /* delta_x[j] = RealAtomValue(new_inputs[j]) - RealAtomValue(old_inputs[j]); */
649 johnpye 1162 }
650    
651 johnpye 1163 for (i=0;i<noutputs;i++) {
652     var = outputs[i];
653     if (var_fixed(var) || !var_active(var)) /* project only the free vars */
654     continue;
655     tmp = 0.0;
656     for (j=0;j<ninputs;j++) {
657 johnpye 1174 tmp += (DENSEMATRIX_ELEM(dy_dx,i,j)* delta_x[j]);
658 johnpye 1163 }
659     /* old_y = RealAtomValue(var); */
660     old_y = var_value(var);
661     new_y = old_y + step_length*tmp;
662     /* SetRealAtomValue(var,new_y,(unsigned)0); */
663     var_set_value(var,new_y);
664     # \if DEBUG
665     CONSOLE_DEBUG("Old_y = %12.8g; Nex_y = %12.8g\n",old_y,new_y);
666     # \endif
667 johnpye 1162 }
668 johnpye 1163 ascfree((char *)delta_x);
669     return 0;
670     }
671     #endif
672 johnpye 1162
673    
674 johnpye 1163 #if 0
675     /**
676     * At this point we should have an empty jacobian. We now
677     * need to call relman_diff over the *entire* matrix.
678     * Calculates the entire jacobian. It is initially unscaled.
679     *
680     * Note: this assumes the sys given is using one of the ascend solvers
681     * with either linsol or linsolqr. Both are now allowed. baa 7/95
682     */
683     #define SAFE_FIX_ME 0
684     static int Compute_J_OLD(slv_system_t sys)
685     {
686     int32 row;
687     var_filter_t vfilter;
688     linsol_system_t lin_sys = NULL;
689     linsolqr_system_t lqr_sys = NULL;
690     mtx_matrix_t mtx;
691     struct rel_relation **rlist;
692     int nrows;
693     int qr=0;
694     #\if DOTIME
695     double time1;
696     #\endif
697 johnpye 1162
698 johnpye 1163 #\if DOTIME
699     time1 = tm_cpu_time();
700     #\endif
701 johnpye 1162 /*
702 johnpye 1163 * Get the linear system from the solve system.
703     * Get the matrix from the linear system.
704     * Get the relations list from the solve system.
705 johnpye 1162 */
706 johnpye 1163 lin_sys = slv_get_linsol_sys(sys);
707     if (lin_sys==NULL) {
708     qr=1;
709     lqr_sys=slv_get_linsolqr_sys(sys);
710 johnpye 1162 }
711 johnpye 1163 mtx = slv_get_sys_mtx(sys);
712     if (mtx==NULL || (lin_sys==NULL && lqr_sys==NULL)) {
713     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Compute_dy_dx: error, found NULL.\n");
714     return 1;
715 johnpye 1162 }
716 johnpye 1163 rlist = slv_get_solvers_rel_list(sys);
717     nrows = NumberIncludedRels(sys);
718 johnpye 1162
719 johnpye 1163 calc_ok = TRUE;
720     vfilter.matchbits = VAR_SVAR;
721     vfilter.matchvalue = VAR_SVAR;
722 johnpye 1162 /*
723 johnpye 1163 * Clear the entire matrix and then compute
724     * the gradients at the current point.
725 johnpye 1162 */
726 johnpye 1163 mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
727     for(row=0; row<nrows; row++) {
728     struct rel_relation *rel;
729     /* added */
730     double resid;
731 johnpye 1162
732 johnpye 1163 rel = rlist[mtx_row_to_org(mtx,row)];
733     (void)relman_diffs(rel,&vfilter,mtx,&resid,SAFE_FIX_ME);
734 johnpye 1162
735 johnpye 1163 /* added */
736     rel_set_residual(rel,resid);
737    
738     }
739     if (qr) {
740     linsolqr_matrix_was_changed(lqr_sys);
741     } else {
742     linsol_matrix_was_changed(lin_sys);
743     }
744     #\if DOTIME
745     time1 = tm_cpu_time() - time1;
746     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Time taken for Compute_J = %g\n",time1);
747     #\endif
748     return(!calc_ok);
749 johnpye 1162 }
750 johnpye 1163 #endif
751 johnpye 1162
752    
753 johnpye 1164 /**
754     Sensitivity Analysis
755 johnpye 1162
756 johnpye 1164 @param arglist List of arguments
757    
758     Argument list contains the following:
759     . arg1 - model inst for which the sensitivity analysis is to be done.
760     . arg2 - array of input instances.
761     . arg3 - array of output instances.
762     . arg4 matrix of partials to be written to.
763 johnpye 1163 */
764 johnpye 1164 int SensitivityCheckArgs(struct gl_list_t *arglist){
765     struct Instance *inst;
766     unsigned long len;
767     unsigned long ninputs, noutputs;
768 johnpye 1163
769 johnpye 1164 len = gl_length(arglist);
770     if (len != 4) {
771     ERROR_REPORTER_HERE(ASC_PROG_ERR,"wrong number of args -- 4 expected\n");
772 johnpye 1163 return 1;
773     }
774 johnpye 1164 inst = FetchElement(arglist,1,1);
775     if (InstanceKind(inst)!=MODEL_INST) {
776     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Arg #1 is not a model instance\n");
777     return 1;
778     }
779     ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
780     /* input args */
781     noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
782     /* output args */
783     len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
784     /* partials matrix */
785     if (len != (ninputs*noutputs)) {
786     ERROR_REPORTER_HERE(ASC_PROG_ERR,
787     "The array of partials is inconsistent with the args given\n");
788     return 1;
789     }
790 johnpye 1163 return 0;
791     }
792    
793 johnpye 1162
794 johnpye 1164 int do_sensitivity_eval( struct Instance *i,
795 johnpye 1163 struct gl_list_t *arglist, void *user_data
796     ){
797 johnpye 1164 CONSOLE_DEBUG("Starting sensitivity analysis...");
798     if(SensitivityCheckArgs(arglist))return 1;
799 johnpye 1162
800 johnpye 1164 return sensitivity_anal(i,arglist);
801 johnpye 1162 }
802    
803    
804     /**
805 johnpye 1164 @param arglist List of arguments
806     @param step_length ...?
807 johnpye 1162
808 johnpye 1164 The list of arguments should chontain
809 johnpye 1162
810 johnpye 1164 . arg1 - model inst for which the sensitivity analysis is to be done.
811     . arg2 - array of input instances.
812     . arg3 - new_input instances, for variable projection.
813     . arg4 - instance representing the step_length for projection.
814 johnpye 1163
815 johnpye 1164 The result will be written to standard out.
816     */
817     int SensitivityAllCheckArgs(struct gl_list_t *arglist, double *step_length)
818     {
819     struct Instance *inst;
820     unsigned long len;
821 johnpye 1163
822 johnpye 1164 len = gl_length(arglist);
823     if (len != 4) {
824     ERROR_REPORTER_HERE(ASC_PROG_ERR,"wrong number of args -- 4 expected\n");
825     return 1;
826 johnpye 1163 }
827 johnpye 1164 inst = FetchElement(arglist,1,1);
828     if (InstanceKind(inst)!=MODEL_INST) {
829     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Arg #1 is not a model instance\n");
830     return 1;
831 johnpye 1163 }
832 johnpye 1164 /*
833     * we should be checking that arg2 list contains solver vars
834     * and that they are fixed. The same applies to arglist 3... later.
835     * We will check and return the steplength though. 0 means dont do
836     * the variable projection.
837     */
838     inst = FetchElement(arglist,4,1);
839     *step_length = RealAtomValue(inst);
840     if (fabs(*step_length) < 1e-08)
841     *step_length = 0.0;
842     return 0;
843 johnpye 1162 }
844    
845    
846     int do_sensitivity_eval_all( struct Instance *i,
847 johnpye 1163 struct gl_list_t *arglist, void *user_data
848     ){
849     int result;
850     double step_length = 0.0;
851     if (SensitivityAllCheckArgs(arglist,&step_length)) {
852     return 1;
853     }
854     result = sensitivity_anal_all(i,arglist,step_length);
855     return result;
856 johnpye 1162 }
857    
858    
859 johnpye 1163 const char sensitivity_help[] =
860 johnpye 1167 "This function does sensitivity analysis dy/dx. It requires 4 args:\n\n"
861 johnpye 1162 " 1. name: name of a reference instance or SELF.\n"
862     " 2. x: x, where x is an array of > solver_var.\n"
863     " 3. y: where y is an array of > solver_var.\n"
864 johnpye 1167 " 4. dy/dx: which dy_dx[1..n_y][1..n_x].\n\n"
865     "See also sensitivity_anal_all.";
866 johnpye 1162
867 johnpye 1167 /** @TODO document what 'u_new' is all about...? */
868    
869     const char sensitivity_all_help[] =
870     "Analyse the sensitivity of *all* variables in the system with respect\n"
871     "to the specific set of inputs 'u'. Instead of returning values to a\n"
872     "a special array inside the model, the results are written to the\n"
873     "console. Usage example:\n\n"
874     "EXTERN sensitivity_anal_all(\n"
875     " this_instance,\n"
876     " u_old[1..n_inputs],\n"
877     " u_new[1..n_inputs],\n"
878     " step_length\n"
879     ");\n\n"
880     "See also sensitivity_anal.";
881    
882 johnpye 1162 int sensitivity_register(void){
883 johnpye 1163 int result=0;
884 johnpye 1162
885 johnpye 1163 result += CreateUserFunctionMethod("do_sensitivity",
886     do_sensitivity_eval,
887     4,sensitivity_help,NULL,NULL
888     );
889     result += CreateUserFunctionMethod("do_sensitivity_all",
890     do_sensitivity_eval_all,
891 johnpye 1167 4,sensitivity_all_help,NULL,NULL
892 johnpye 1163 );
893 johnpye 1162
894 johnpye 1163 return result;
895 johnpye 1162 }
896    
897 johnpye 1163 /* :ex: set ts=4: */

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