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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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