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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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