/[ascend]/trunk/base/generic/packages/sensitivity.c
ViewVC logotype

Annotation of /trunk/base/generic/packages/sensitivity.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 91 - (hide annotations) (download) (as text)
Wed Dec 7 16:17:03 2005 UTC (14 years, 8 months ago) by johnpye
File MIME type: text/x-csrc
File size: 32439 byte(s)
Removed rootfind.c from solver lib in jamfile
Tidied up error messages in safe.c
Fixed a missing #endif in sensitivity.c
1 aw0a 27 /*********************************************************************\
2     Sensititvity analysis code. Kirk Abbott.
3     \*********************************************************************/
4    
5     #include <string.h>
6     #include <math.h>
7     #include "utilities/ascConfig.h"
8     #include "compiler/instance_enum.h"
9     #include "compiler/compiler.h"
10     #include "general/list.h"
11     #include "utilities/ascMalloc.h"
12     #include "compiler/extfunc.h"
13     #include "compiler/packages.h"
14     #include "compiler/fractions.h"
15     #include "compiler/dimen.h"
16     #include "compiler/atomvalue.h"
17     #include "compiler/instquery.h"
18     #include "solver/mtx.h"
19     #include "solver/mtx_basic.h"
20     #include "solver/mtx_perms.h"
21     #include "solver/mtx_query.h"
22     #include "solver/linsol.h"
23     #include "solver/linsolqr.h"
24     #include "solver/slv_types.h"
25     #include "solver/var.h"
26     #include "solver/rel.h"
27     #include "solver/discrete.h"
28     #include "solver/conditional.h"
29     #include "solver/logrel.h"
30     #include "solver/bnd.h"
31     #include "solver/calc.h"
32     #include "solver/relman.h"
33     #include "solver/slv_common.h"
34     #include "solver/slv_stdcalls.h"
35     #include "solver/system.h"
36     #include "solver/slv_client.h"
37    
38     #define DEBUG 1
39    
40     /*
41     * This file attempts to implement the extraction of dy_dx from
42     * a system of equations. If one considers a black-box where x are
43     * the input variables, u are inuut parameters, y are the output
44     * variables, f(x,y,u) is the system of equations that will be solved
45     * for given x and u, then one can extract the sensitivity information
46     * of y wrt x.
47     * One crude but simple way of doing this is to finite-difference the
48     * given black box, i.e, perturb x, n\subx times by delta x, resolve
49     * the system of equations at the new value of x, and compute
50     * dy/dx = (f(x\sub1) - f(x\sub2))/(x\sub1 - x\sub2).
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\subx 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    
66 johnpye 89 #if 0
67 aw0a 27 static real64 *zero_vector(real64 *vec, int size)
68     {
69     int c;
70     for (c=0;c<size;c++) {
71     vec[c] = 0.0;
72     }
73     return vec;
74     }
75 johnpye 89 #endif
76 aw0a 27
77     static real64 **make_matrix(int nrows, int ncols)
78     {
79     real64 **result;
80     int i;
81     result = (real64 **)calloc(nrows,sizeof(real64*));
82     for (i=0;i<nrows;i++) {
83     result[i] = (real64 *)calloc(ncols,sizeof(real64));
84     }
85     return result;
86     }
87    
88     static void free_matrix(real64 **matrix, int nrows)
89     {
90     int i;
91     if (!matrix)
92     return;
93     for (i=0;i<nrows;i++) {
94     if (matrix[i]) {
95     free(matrix[i]);
96     matrix[i] = NULL;
97     }
98     }
99     free(matrix);
100     }
101    
102     static struct Instance *FetchElement(struct gl_list_t *arglist,
103     unsigned long whichbranch,
104     unsigned long whichelement)
105     {
106     struct gl_list_t *branch;
107     struct Instance *element;
108    
109     if (!arglist) return NULL;
110     branch = (struct gl_list_t *)gl_fetch(arglist,whichbranch);
111     element = (struct Instance *)gl_fetch(branch,whichelement);
112     return element;
113     }
114    
115     static slv_system_t PreSolve(struct Instance *inst)
116     {
117     slv_system_t sys;
118     slv_parameters_t parameters;
119     struct var_variable **vp;
120     struct rel_relation **rp;
121     int ind,len;
122     char *tmp=NULL;
123    
124     sys = system_build(inst);
125     if (sys==NULL) {
126     FPRINTF(stdout,
127     "(sensitivity.c): Something radically wrong in creating solver\n");
128     return NULL;
129     }
130     if (g_SlvNumberOfRegisteredClients == 0) {
131     return NULL;
132     }
133     ind = 0;
134     while (strcmp(slv_solver_name(ind),"QRSlv")) {
135     if (ind >= g_SlvNumberOfRegisteredClients) {
136     FPRINTF(stderr,"(sensitivity.c): QRSlv must be registered client\n");
137     return NULL;
138     }
139     ++ind;
140     }
141     slv_select_solver(sys,ind);
142    
143     slv_get_parameters(sys,&parameters);
144     parameters.partition = 0;
145     slv_set_parameters(sys,&parameters);
146     slv_presolve(sys);
147    
148     #if DEBUG
149     vp = slv_get_solvers_var_list(sys);
150     len = slv_get_num_solvers_vars(sys);
151     for (ind=0 ; ind<len; ind++) {
152     tmp = var_make_name(sys,vp[ind]);
153     FPRINTF(stderr,"%s %d\n",tmp,var_sindex(vp[ind]));
154     if (tmp!=NULL) ascfree(tmp);
155     }
156     rp = slv_get_solvers_rel_list(sys);
157     len = slv_get_num_solvers_rels(sys);
158     for (ind=0 ; ind<len ; ind++) {
159     tmp = rel_make_name(sys,rp[ind]);
160     FPRINTF(stderr,"%s %d\n",tmp,rel_sindex(rp[ind]));
161     if (tmp) ascfree(tmp);
162     }
163     #endif
164     return sys;
165     }
166 johnpye 89
167     #if 0
168 aw0a 27 static int ReSolve(slv_system_t sys)
169     {
170     if (!sys)
171     return 1;
172     slv_solve(sys);
173     return 0;
174     }
175 johnpye 89 #endif
176 aw0a 27
177     static int DoSolve(struct Instance *inst)
178     {
179     slv_system_t sys;
180    
181     sys = system_build(inst);
182     if (!sys) {
183     FPRINTF(stdout,"Something radically wrong in creating solver\n");
184     return 1;
185     }
186     (void)slv_select_solver(sys,0);
187     slv_presolve(sys);
188     slv_solve(sys);
189     system_destroy(sys);
190     return 0;
191     }
192    
193     int do_solve_eval(struct Slv_Interp *slv_interp,
194     struct Instance *i,
195     struct gl_list_t *arglist,
196     unsigned long whichvar)
197     {
198     unsigned long len;
199     int result;
200     struct Instance *inst;
201     len = gl_length(arglist);
202 johnpye 89
203     /* Ignore unused params */
204     (void)slv_interp; (void)i; (void)whichvar;
205    
206 aw0a 27 if (len!=2) {
207     FPRINTF(stdout,"Wrong number of args to do_solve_eval\n");
208     return 1;
209     }
210     inst = FetchElement(arglist,1,1);
211     if (!inst)
212     return 1;
213     result = DoSolve(inst);
214     return result;
215     }
216    
217     static int FiniteDiffCheckArgs(struct gl_list_t *arglist)
218     {
219     struct Instance *inst;
220     unsigned long len;
221     unsigned long ninputs, noutputs;
222    
223     /*
224     * arg1 - model inst to be solved.
225     * arg2 - array of input instances.
226     * arg3 - array of output instances.
227     * matrix of partials to be written to.
228     */
229     len = gl_length(arglist);
230     if (len != 4) {
231     FPRINTF(stderr,"wrong number of args -- 4 expected\n");
232     return 1;
233     }
234     inst = FetchElement(arglist,1,1);
235     if (InstanceKind(inst)!=MODEL_INST) {
236     FPRINTF(stderr,"Arg #1 is not a model instance\n");
237     return 1;
238     }
239     ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
240     /* input args */
241     noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
242     /* output args */
243     len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
244     /* partials matrix */
245     if (len != (ninputs*noutputs)) {
246     FPRINTF(stderr,
247     "The array of partials is inconsistent with the args given\n");
248     return 1;
249     }
250     return 0;
251     }
252    
253     static int finite_difference(struct gl_list_t *arglist)
254     {
255     struct Instance *model_inst, *xinst, *inst;
256     slv_system_t sys;
257     int ninputs,noutputs;
258     int i,j,offset;
259     real64 **partials;
260     real64 *y_old, *y_new;
261     real64 x;
262     real64 interval = 1e-6;
263     int result=0;
264    
265     model_inst = FetchElement(arglist,1,1);
266     sys = system_build(model_inst);
267     if (!sys) {
268     FPRINTF(stdout,"Something radically wrong in creating solver\n");
269     return 1;
270     }
271     (void)slv_select_solver(sys,0);
272     slv_presolve(sys);
273     slv_solve(sys);
274     /*
275     * Make the necessary vectors.
276     */
277     ninputs = (int)gl_length((struct gl_list_t *)gl_fetch(arglist,2));
278     noutputs = (int)gl_length((struct gl_list_t *)gl_fetch(arglist,3));
279     y_old = (real64 *)calloc(noutputs,sizeof(real64));
280     y_new = (real64 *)calloc(noutputs,sizeof(real64));
281     partials = make_matrix(noutputs,ninputs);
282     for (i=0;i<noutputs;i++) { /* get old yvalues */
283     inst = FetchElement(arglist,3,i+1);
284     y_old[i] = RealAtomValue(inst);
285     }
286     for (j=0;j<ninputs;j++) {
287     xinst = FetchElement(arglist,2,j+1);
288     x = RealAtomValue(xinst);
289     SetRealAtomValue(xinst,x+interval,(unsigned)0); /* perturb system */
290     slv_presolve(sys);
291     slv_solve(sys);
292    
293     for (i=0;i<noutputs;i++) { /* get new yvalues */
294     inst = FetchElement(arglist,3,i+1);
295     y_new[i] = RealAtomValue(inst);
296     partials[i][j] = -1.0 * (y_old[i] - y_new[i])/interval;
297     PRINTF("y_old = %20.12g y_new = %20.12g\n", y_old[i],y_new[i]);
298     }
299     SetRealAtomValue(xinst,x,(unsigned)0); /* unperturb system */
300     }
301     offset = 0;
302     for (i=0;i<noutputs;i++) {
303     for (j=0;j<ninputs;j++) {
304     inst = FetchElement(arglist,4,offset+j+1);
305     SetRealAtomValue(inst,partials[i][j],(unsigned)0);
306     PRINTF("%12.6f %s",partials[i][j], (j==(ninputs-1)) ? "\n" : " ");
307     }
308     offset += ninputs;
309     }
310     /* error: */
311     free(y_old);
312     free(y_new);
313     free_matrix(partials,noutputs);
314     system_destroy(sys);
315     return result;
316     }
317    
318     int do_finite_diff_eval(struct Slv_Interp *slv_interp,
319     struct Instance *i,
320     struct gl_list_t *arglist,
321     unsigned long whichvar)
322     {
323 johnpye 89 /* Ignore unused params */
324     (void)slv_interp; (void)i; (void)whichvar;
325    
326 aw0a 27 int result;
327     if (FiniteDiffCheckArgs(arglist))
328     return 1;
329     result = finite_difference(arglist);
330     return result;
331     }
332    
333     /*********************************************************************\
334     Sensititvity analysis code.
335     \*********************************************************************/
336    
337     int SensitivityCheckArgs(struct gl_list_t *arglist)
338     {
339     struct Instance *inst;
340     unsigned long len;
341     unsigned long ninputs, noutputs;
342    
343     /*
344     * arg1 - model inst for which the sensitivity analysis is to be done.
345     * arg2 - array of input instances.
346     * arg3 - array of output instances.
347     * arg4 matrix of partials to be written to.
348     */
349    
350     len = gl_length(arglist);
351     if (len != 4) {
352     FPRINTF(stderr,"wrong number of args -- 4 expected\n");
353     return 1;
354     }
355     inst = FetchElement(arglist,1,1);
356     if (InstanceKind(inst)!=MODEL_INST) {
357     FPRINTF(stderr,"Arg #1 is not a model instance\n");
358     return 1;
359     }
360     ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
361     /* input args */
362     noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
363     /* output args */
364     len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
365     /* partials matrix */
366     if (len != (ninputs*noutputs)) {
367     FPRINTF(stderr,
368     "The array of partials is inconsistent with the args given\n");
369     return 1;
370     }
371     return 0;
372     }
373    
374     int SensitivityAllCheckArgs(struct gl_list_t *arglist, double *step_length)
375     {
376     struct Instance *inst;
377     unsigned long len;
378    
379     /*
380     * arg1 - model inst for which the sensitivity analysis is to be done.
381     * arg2 - array of input instances.
382     * arg3 - new_input instances, for variable projection.
383     * arg4 - instance representing the step_length for projection.
384     * The result will be written to standard out.
385     */
386    
387     len = gl_length(arglist);
388     if (len != 4) {
389     FPRINTF(stderr,"wrong number of args -- 4 expected\n");
390     return 1;
391     }
392     inst = FetchElement(arglist,1,1);
393     if (InstanceKind(inst)!=MODEL_INST) {
394     FPRINTF(stderr,"Arg #1 is not a model instance\n");
395     return 1;
396     }
397     /*
398     * we should be checking that arg2 list contains solver vars
399     * and that they are fixed. The same applies to arglist 3... later.
400     * We will check and return the steplength though. 0 means dont do
401     * the variable projection.
402     */
403     inst = FetchElement(arglist,4,1);
404     *step_length = RealAtomValue(inst);
405     if (fabs(*step_length) < 1e-08)
406     *step_length = 0.0;
407     return 0;
408     }
409    
410    
411     static int NumberRels(slv_system_t sys)
412     {
413     static int nrels = -1;
414     rel_filter_t rfilter;
415     int tmp;
416    
417     if (!sys) { /* a NULL system may be used to reinitialise the count */
418     nrels = -1;
419     return -1;
420     }
421     if (nrels < 0) {
422     /*get used (included) relation count -- equalities only !*/
423     rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
424     rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
425     tmp = slv_count_solvers_rels(sys,&rfilter);
426     nrels = tmp;
427     return tmp;
428     }
429     else return nrels;
430     }
431    
432     static int NumberFreeVars(slv_system_t sys)
433     {
434     static int nvars = -1;
435     var_filter_t vfilter;
436     int tmp;
437    
438     if (!sys) {
439     nvars = -1;
440     return -1;
441     }
442     if (nvars < 0) {
443     /*get used (free and incident) variable count */
444     vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
445     vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
446     tmp = slv_count_solvers_vars(sys,&vfilter);
447     nvars = tmp;
448     return tmp;
449     }
450     else return nvars;
451     }
452    
453    
454     /***************************************************************/
455     /***************************************************************/
456     /*
457     * The following bit of code does the computation of the matrix
458     * dz/dx. It accepts a slv_system_t and a list of 'input' variables.
459     * The matrix df/dx now exists and sits to the right of the main
460     * square region of the Jacobian. We will do the following in turn
461     * for each variable in this list:
462     *
463     * 1) Get the variable index, and use it to extract the required column
464     * from the main gradient matrix, this will be stored in a temporary
465     * vector.
466     * 2) We will then clear the column in the original matrix, as we want to
467     * store the caluclated results back in place.
468     * 3) Add the vector extracted in 1) as rhs to the main matrix.
469     * 4) Call linsol solve on this rhs to yield a solution which represents
470     * a column of dx/dx.
471     * 6) Store this vector back in the location cleared out in 2).
472     * 7) Goto 1.
473     *
474     * At the end of this procedure we would have calculated all the columns of
475     * dz/dx, and left them back in the main matrix.
476     */
477    
478     /*******************************************************************/
479     /*******************************************************************/
480    
481     /*
482     * At this point we should have an empty jacobian. We now
483     * need to call relman_diff over the *entire* matrix.
484     * fixed and free.
485     *
486     * Calculates the entire jacobian.
487     * It is initially unscaled.
488     */
489     static int Compute_J(slv_system_t sys)
490     {
491     int32 row;
492     var_filter_t vfilter;
493     linsolqr_system_t lqr_sys;
494     mtx_matrix_t mtx;
495     struct rel_relation **rlist;
496     int nrows;
497     real64 resid;
498    
499     /*
500     * Get the linear system from the solve system.
501     * Get the matrix from the linear system.
502     * Get the relations list from the solve system.
503     */
504     lqr_sys = slv_get_linsolqr_sys(sys);
505     mtx = linsolqr_get_matrix(lqr_sys);
506     rlist = slv_get_solvers_rel_list(sys);
507     nrows = NumberRels(sys);
508    
509     calc_ok = TRUE;
510     vfilter.matchbits = (VAR_SVAR | VAR_ACTIVE) ;
511     vfilter.matchvalue = vfilter.matchbits ;
512    
513     /*
514     * Clear the entire matrix and then compute
515     * the gradients at the current point.
516     */
517     mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
518     for(row=0; row<nrows; row++) {
519     struct rel_relation *rel;
520     rel = rlist[mtx_row_to_org(mtx,row)];
521     (void)relman_diffs(rel,&vfilter,mtx,&resid,1);
522     }
523     linsolqr_matrix_was_changed(lqr_sys);
524    
525     return(!calc_ok);
526     }
527    
528     /*
529     * Note a rhs would have been previously added
530     * to keep the system happy.
531     */
532    
533     static int LUFactorJacobian(slv_system_t sys,int rank)
534     {
535     linsolqr_system_t lqr_sys;
536     mtx_region_t region;
537     enum factor_method fm;
538    
539     mtx_region(&region,0,rank-1,0,rank-1); /* set the region */
540     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
541    
542     linsolqr_matrix_was_changed(lqr_sys);
543     linsolqr_reorder(lqr_sys,&region,natural);
544    
545     fm = linsolqr_fmethod(lqr_sys);
546     if (fm == unknown_f) fm = ranki_kw2; /* make sure somebody set it */
547     linsolqr_factor(lqr_sys,fm);
548    
549     return 0;
550     }
551    
552    
553     /*
554     * At this point the rhs would have already been added.
555     *
556     */
557    
558    
559     static int Compute_dy_dx_smart(slv_system_t sys,
560     real64 *rhs,
561     real64 **dy_dx,
562     int *inputs, int ninputs,
563     int *outputs, int noutputs)
564     {
565     linsolqr_system_t lqr_sys;
566     mtx_matrix_t mtx;
567     int col,current_col;
568     int row;
569     int capacity;
570     real64 *solution = NULL;
571     int i,j,k;
572    
573     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
574     mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
575    
576     capacity = mtx_capacity(mtx);
577     solution = (real64 *)asccalloc(capacity,sizeof(real64));
578    
579     /*
580     * The array inputs is a list of original indexes, of the variables
581     * that we are trying to obtain the sensitivity to. We have to
582     * get the *current* column from the matrix based on those indices.
583     * Hence we use mtx_org_to_col. This little manipulation is not
584     * necessary for the computed solution as the solve routine returns
585     * the results in the *original* order rather than the *current* order.
586     */
587     for (j=0;j<ninputs;j++) {
588     col = inputs[j];
589     current_col = mtx_org_to_col(mtx,col);
590     mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS);
591     /* get the column in org row indexed order */
592    
593     linsolqr_rhs_was_changed(lqr_sys,rhs);
594     linsolqr_solve(lqr_sys,rhs); /* solve */
595     linsolqr_copy_solution(lqr_sys,rhs,solution); /* get the solution */
596    
597     #if DEBUG
598     FPRINTF(stderr,"This is the rhs and solution for rhs #%d orgcol%d\n",j,col);
599     for (k=0;k<capacity;k++) {
600     FPRINTF(stderr,"%12.8g %12.8g\n",rhs[k],solution[k]);
601     }
602     #endif
603    
604     for (i=0;i<noutputs;i++) { /* copy the solution to dy_dx */
605     row = outputs[i];
606     dy_dx[i][j] = -1.0*solution[row]; /* the -1 comes from the lin algebra */
607     }
608     /*
609     * zero the vector using the incidence pattern of our column.
610     * This is faster than the naiive approach of zeroing
611     * everything especially if the vector is large.
612     */
613     mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
614     }
615    
616     if (solution) ascfree((char *)solution);
617     return 0;
618     }
619 johnpye 89
620     #if 0
621 aw0a 27 static int ComputeInverse(slv_system_t sys,
622     real64 *rhs)
623     {
624     linsolqr_system_t lqr_sys;
625     mtx_matrix_t mtx;
626     int capacity,order;
627     real64 *solution = NULL;
628     int j,k;
629    
630     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
631     mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
632    
633     capacity = mtx_capacity(mtx);
634     zero_vector(rhs,capacity); /* zero the rhs */
635     solution = (real64 *)asccalloc(capacity,sizeof(real64));
636    
637     order = mtx_order(mtx);
638     for (j=0;j<order;j++) {
639     rhs[j] = 1.0;
640     linsolqr_rhs_was_changed(lqr_sys,rhs);
641     linsolqr_solve(lqr_sys,rhs); /* solve */
642     linsolqr_copy_solution(lqr_sys,rhs,solution); /* get the solution */
643    
644     FPRINTF(stderr,"This is the rhs and solution for rhs #%d\n",j);
645     for (k=0;k<order;k++) {
646     FPRINTF(stderr,"%12.8g %12.8g\n",rhs[k],solution[k]);
647     }
648     rhs[j] = 0.0;
649     }
650     if (solution) ascfree((char *)solution);
651     return 0;
652     }
653 johnpye 89 #endif
654 aw0a 27
655     int sensitivity_anal(struct Slv_Interp *slv_interp,
656     struct Instance *inst, /* not used but will be */
657     struct gl_list_t *arglist)
658     {
659     struct Instance *which_instance,*tmp_inst, *atominst;
660     struct gl_list_t *branch;
661     struct var_variable **vlist = NULL;
662     int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
663 johnpye 89 real64 **dy_dx = NULL;
664 aw0a 27 slv_system_t sys = NULL;
665 johnpye 89 int c;
666     int noutputs = 0;
667     int ninputs;
668 aw0a 27 int i,j;
669     int offset;
670     dof_t *dof;
671     int num_vars,ind,found;
672    
673 johnpye 89 /* Ignore unused params */
674     (void)slv_interp; (void) inst;
675    
676 aw0a 27 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
677     mtx_matrix_t mtx;
678     int32 capacity;
679     real64 *scratch_vector = NULL;
680     int result=0;
681    
682 johnpye 89
683 aw0a 27 (void)NumberFreeVars(NULL); /* used to re-init the system */
684     (void)NumberRels(NULL); /* used to re-init the system */
685     which_instance = FetchElement(arglist,1,1);
686     sys = PreSolve(which_instance);
687     if (!sys) {
688     FPRINTF(stderr,"Early termination due to failure in Presolve\n");
689     result = 1;
690     goto error;
691     }
692    
693     dof = slv_get_dofdata(sys);
694     if (!(dof->n_rows == dof->n_cols &&
695     dof->n_rows == dof->structural_rank)) {
696     FPRINTF(stderr,"Early termination: non square system\n");
697     result = 1;
698     goto error;
699     }
700     /*
701     * prepare the inputs list
702     */
703     vlist = slv_get_solvers_var_list(sys);
704    
705     branch = gl_fetch(arglist,2); /* input args */
706     ninputs = gl_length(branch);
707     inputs_ndx_list = (int *)ascmalloc(ninputs*sizeof(int));
708    
709     num_vars = slv_get_num_solvers_vars(sys);
710     for (c=0;c<ninputs;c++) {
711     atominst = (struct Instance *)gl_fetch(branch,c+1);
712     found = 0;
713     ind = num_vars - 1;
714     /* search backwards because fixed vars are at the end of the
715     var list */
716     while (!found && ind >= 0){
717     if (var_instance(vlist[ind]) == atominst) {
718     inputs_ndx_list[c] = var_sindex(vlist[ind]);
719     found = 1;
720     }
721     --ind;
722     }
723     if (!found) {
724     FPRINTF(stderr,"Sensitivity input variable not found\n");
725     result = 1;
726     goto error;
727     }
728     }
729    
730     /*
731     * prepare the outputs list
732     */
733     branch = gl_fetch(arglist,3); /* output args */
734     noutputs = gl_length(branch);
735     outputs_ndx_list = (int *)ascmalloc(noutputs*sizeof(int));
736     for (c=0;c<noutputs;c++) {
737     atominst = (struct Instance *)gl_fetch(branch,c+1);
738     found = 0;
739     ind = 0;
740     while (!found && ind < num_vars){
741     if (var_instance(vlist[ind]) == atominst) {
742     outputs_ndx_list[c] = var_sindex(vlist[ind]);
743     found = 1;
744     }
745     ++ind;
746     }
747     if (!found) {
748     FPRINTF(stderr,"Sensitivity ouput variable not found\n");
749     result = 1;
750     goto error;
751     }
752     }
753    
754     /*
755     * prepare the results dy_dx.
756     */
757     dy_dx = make_matrix(noutputs,ninputs);
758    
759    
760     result = Compute_J(sys);
761     if (result) {
762     FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
763     goto error;
764     }
765    
766     /*
767     * Note: the RHS *has* to added here. We will construct the vector
768     * of size matrix capacity and add it. It will be removed after
769     * we finish computing dy_dx.
770     */
771     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
772     mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
773     capacity = mtx_capacity(mtx);
774     scratch_vector = (real64 *)asccalloc(capacity,sizeof(real64));
775     linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
776     result = LUFactorJacobian(sys,dof->structural_rank);
777     if (result) {
778     FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
779     goto error;
780     }
781     result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
782     inputs_ndx_list, ninputs,
783     outputs_ndx_list, noutputs);
784    
785     linsolqr_remove_rhs(lqr_sys,scratch_vector);
786     if (result) {
787     FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
788     goto error;
789     }
790    
791     /*
792     * Write the results back to the partials array in the
793     * instance tree
794     */
795     offset = 0;
796     for (i=0;i<noutputs;i++) {
797     for (j=0;j<ninputs;j++) {
798     tmp_inst = FetchElement(arglist,4,offset+j+1);
799     SetRealAtomValue(tmp_inst,dy_dx[i][j],(unsigned)0);
800     #if DEBUG
801     FPRINTF(stderr,"%12.8f i%d j%d",dy_dx[i][j],i,j);
802     #endif
803     }
804     #if DEBUG
805     FPRINTF(stderr,"\n");
806     #endif
807     offset += ninputs;
808     }
809     #if DEBUG
810     FPRINTF(stderr,"\n");
811     #endif
812    
813     error:
814     if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
815     if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
816     if (dy_dx) free_matrix(dy_dx,noutputs);
817     if (scratch_vector) ascfree((char *)scratch_vector);
818     if (sys) system_destroy(sys);
819     return result;
820     }
821    
822    
823     static int DoDataAnalysis(struct var_variable **inputs,
824     struct var_variable **outputs,
825     int ninputs, int noutputs,
826     real64 **dy_dx)
827     {
828     FILE *fp;
829     double *norm_2, *norm_1;
830     double input_nominal,maxvalue,sum;
831     int i,j;
832    
833     norm_1 = (double *)asccalloc(ninputs,sizeof(double));
834     norm_2 = (double *)asccalloc(ninputs,sizeof(double));
835    
836     fp = fopen("sensitivity.out","w+");
837     if (!fp) return 0;
838    
839     /*
840     * calculate the 1 and 2 norms; cache them away so that we
841     * can pretty print them. Style is everything !.
842     *
843     */
844     for (j=0;j<ninputs;j++) {
845     input_nominal = var_nominal(inputs[j]);
846     maxvalue = sum = 0;
847     for (i=0;i<noutputs;i++) {
848     dy_dx[i][j] *= input_nominal/var_nominal(outputs[i]);
849     maxvalue = MAX(fabs(dy_dx[i][j]),maxvalue);
850     sum += dy_dx[i][j]*dy_dx[i][j];
851     }
852     norm_1[j] = maxvalue;
853     norm_2[j] = sum;
854     }
855    
856     for (j=0;j<ninputs;j++) { /* print the var_index */
857     FPRINTF(fp,"%8d ",var_mindex(inputs[j]));
858     }
859     FPRINTF(fp,"\n");
860    
861     for (j=0;j<ninputs;j++) { /* print the 1 norms */
862     FPRINTF(fp,"%-#18.8f ",norm_1[j]);
863     }
864     FPRINTF(fp,"\n");
865    
866     for (j=0;j<ninputs;j++) { /* print the 2 norms */
867     FPRINTF(fp,"%-#18.8f ",norm_2[j]);
868     }
869     FPRINTF(fp,"\n\n");
870     ascfree((char *)norm_1);
871     ascfree((char *)norm_2);
872    
873     for (i=0;i<noutputs;i++) { /* print the scaled data */
874     for (j=0;j<ninputs;j++) {
875     FPRINTF(fp,"%-#18.8f %-4d",dy_dx[i][j],i);
876     }
877     if (var_fixed(outputs[i]))
878     FPRINTF(fp," **fixed*** \n");
879     else
880     PUTC('\n',fp);
881     }
882     FPRINTF(fp,"\n");
883     fclose(fp);
884     return 0;
885     }
886 johnpye 89
887     #if 0
888 aw0a 27 static int DoProject_X(struct var_variable **old_inputs,
889     struct var_variable **new_inputs, /* new values of u */
890     double step_length,
891     struct var_variable **outputs,
892     int ninputs, int noutputs,
893     real64 **dy_dx)
894     {
895     struct var_variable *var;
896     real64 old_y, new_y, tmp;
897     real64 *delta_x;
898     int i,j;
899    
900     delta_x = (real64 *)asccalloc(ninputs,sizeof(real64));
901    
902     for (j=0;j<ninputs;j++) {
903     delta_x[j] = var_value(new_inputs[j]) - var_value(old_inputs[j]);
904     /* delta_x[j] = RealAtomValue(new_inputs[j]) - RealAtomValue(old_inputs[j]); */
905     }
906    
907     for (i=0;i<noutputs;i++) {
908     var = outputs[i];
909     if (var_fixed(var) || !var_active(var)) /* project only the free vars */
910     continue;
911     tmp = 0.0;
912     for (j=0;j<ninputs;j++) {
913     tmp += (dy_dx[i][j] * delta_x[j]);
914     }
915     /* old_y = RealAtomValue(var); */
916     old_y = var_value(var);
917     new_y = old_y + step_length*tmp;
918     /* SetRealAtomValue(var,new_y,(unsigned)0); */
919     var_set_value(var,new_y);
920 johnpye 89 # if DEBUG
921 aw0a 27 FPRINTF(stderr,"Old_y = %12.8g; Nex_y = %12.8g\n",old_y,new_y);
922 johnpye 91 # endif
923 aw0a 27 }
924     ascfree((char *)delta_x);
925     return 0;
926     }
927 johnpye 89 #endif
928 aw0a 27
929    
930     /*
931     * This function is very similar to sensitivity_anal, execept,
932     * that it perform the analysis on the entire system, based on
933     * the given inputs. Nothing is returned. As such the call from
934     * ASCEND is :
935     *
936     * EXTERN sensitivity_anal_all(
937     * this_instance,
938     * u_old[1..n_inputs],
939     * u_new[1..n_inputs],
940     * step_length
941     * );
942     *
943     * The results will be witten to standard out.
944     * This function is more expensive from a a memory point of view,
945     * as we keep aroung a dense matrix n_outputs * n_inputs, but here
946     * n_outputs may be *much* larger depending on problem size.
947     */
948    
949     int sensitivity_anal_all(struct Slv_Interp *slv_interp,
950     struct Instance *inst, /* not used but will be */
951     struct gl_list_t *arglist,
952     real64 step_length)
953     {
954     struct Instance *which_instance;
955     struct gl_list_t *branch2, *branch3;
956     dof_t *dof;
957     struct var_variable **inputs = NULL, **outputs = NULL;
958     int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
959 johnpye 89 real64 **dy_dx = NULL;
960 aw0a 27 struct var_variable **vp,**ptr;
961     slv_system_t sys = NULL;
962 johnpye 89 long c;
963     int noutputs=0, ninputs;
964 aw0a 27 var_filter_t vfilter;
965    
966 johnpye 89 /* Ignore unused params */
967     (void)slv_interp; (void)inst; (void)step_length;
968    
969 aw0a 27 struct var_variable **new_inputs = NULL; /* optional stuff for variable
970     * projection */
971    
972     linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
973     mtx_matrix_t mtx;
974     int32 capacity;
975     real64 *scratch_vector = NULL;
976     int result=0;
977    
978     /*
979     * Call the presolve for the system. This should number variables
980     * and relations as well create and order the main Jacobian. The
981     * only potential problem that I see here is that presolve using
982     * the slv0 solver *only* recognizes solver vars. So that if one
983     * wanted to see the sensitivity of a *parameter*, it would not
984     * be possible. We will have to trap this in CheckArgs.
985     *
986     * Also the current version of ascend is fucked in how the var module
987     * handles variables and their numbering through the interface ptr
988     * crap.
989     */
990    
991     (void)NumberFreeVars(NULL); /* used to re-init the system */
992     (void)NumberRels(NULL); /* used to re-init the system */
993     which_instance = FetchElement(arglist,1,1);
994     sys = PreSolve(which_instance);
995     if (!sys) {
996     FPRINTF(stderr,"Early termination due to failure in Presolve\n");
997     result = 1;
998     goto error;
999     }
1000     dof = slv_get_dofdata(sys);
1001    
1002     /*
1003     * prepare the inputs list. We dont need the index of the new_inputs
1004     * list. We will grab them later if necessary.
1005     */
1006     branch2 = gl_fetch(arglist,2); /* input args -- old u_values */
1007     branch3 = gl_fetch(arglist,3); /* input args -- new u_values */
1008     ninputs = gl_length(branch2);
1009     inputs = (struct var_variable **)ascmalloc((ninputs+1)*sizeof(struct var_variable *));
1010     new_inputs = (struct var_variable **)ascmalloc((ninputs+1)*sizeof(struct var_variable *));
1011    
1012     inputs_ndx_list = (int *)ascmalloc(ninputs*sizeof(int));
1013     for (c=0;c<ninputs;c++) {
1014     inputs[c] = (struct var_variable *)gl_fetch(branch2,c+1);
1015     inputs_ndx_list[c] = var_mindex(inputs[c]);
1016     new_inputs[c] = (struct var_variable *)gl_fetch(branch3,c+1);
1017     }
1018     inputs[ninputs] = NULL; /* null terminate the list */
1019     new_inputs[ninputs] = NULL; /* null terminate the list */
1020    
1021     /*
1022     * prepare the outputs list. This is where we differ from
1023     * the other function. The noutputs, and the indexes of these
1024     * outputs is obtained from the entire solve system.
1025     */
1026     vfilter.matchbits = 0;
1027     noutputs = slv_count_solvers_vars(sys,&vfilter);
1028     outputs = (struct var_variable **)ascmalloc((noutputs+1)*sizeof(struct var_variable *));
1029     outputs_ndx_list = (int *)ascmalloc(noutputs*sizeof(int));
1030     ptr = vp = slv_get_solvers_var_list(sys);
1031     for (c=0;c<noutputs;c++) {
1032     outputs[c] = *ptr;
1033     outputs_ndx_list[c] = var_sindex(*ptr);
1034     ptr++;
1035     }
1036     outputs[noutputs] = NULL; /* null terminate the list */
1037    
1038     /*
1039     * prepare the results dy_dx. This is the expensive part from a
1040     * memory point of view. However I would like to have the entire
1041     * noutputs * ninputs matrix even for a short while so that I
1042     * can compute a number of different types of norms.
1043     */
1044     dy_dx = make_matrix(noutputs,ninputs);
1045    
1046     result = Compute_J(sys);
1047     if (result) {
1048     FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
1049     goto error;
1050     }
1051    
1052     /*
1053     * Note: the RHS *has* to added here. We will construct the vector
1054     * of size matrix capacity and add it. It will be removed after
1055     * we finish computing dy_dx.
1056     */
1057     lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
1058     mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
1059     capacity = mtx_capacity(mtx);
1060     scratch_vector = (real64 *)asccalloc(capacity,sizeof(real64));
1061     linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
1062     result = LUFactorJacobian(sys,dof->structural_rank);
1063     if (result) {
1064     FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
1065     goto error;
1066     }
1067     result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
1068     inputs_ndx_list, ninputs,
1069     outputs_ndx_list, noutputs);
1070    
1071     linsolqr_remove_rhs(lqr_sys,scratch_vector);
1072     if (result) {
1073     FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
1074     goto error;
1075     }
1076    
1077     /*
1078     * Do some analysis on the results, inclusive of
1079     * writing them to someplace useful.
1080     */
1081     if (DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx))
1082     result = 1;
1083    
1084     /*
1085     * Some experimental projection stuff -- not used now.
1086     * if (DoProject_X(inputs, new_inputs, step_length,
1087     * outputs, ninputs, noutputs, dy_dx))
1088     * result = 1;
1089     */
1090    
1091     error:
1092     if (inputs) ascfree((char *)inputs);
1093     if (new_inputs) ascfree((char *)new_inputs);
1094     if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
1095     if (outputs) ascfree((char *)outputs);
1096     if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
1097     if (dy_dx) free_matrix(dy_dx,noutputs);
1098     if (scratch_vector) ascfree((char *)scratch_vector);
1099     if (sys) system_destroy(sys);
1100     return result;
1101     }
1102    
1103    
1104     int do_sensitivity_eval(struct Slv_Interp *slv_interp,
1105     struct Instance *i,
1106     struct gl_list_t *arglist)
1107     {
1108     int result;
1109     if (SensitivityCheckArgs(arglist)) {
1110     return 1;
1111     }
1112     result = sensitivity_anal(slv_interp,i,arglist);
1113     return result;
1114     }
1115    
1116     int do_sensitivity_eval_all(struct Slv_Interp *slv_interp,
1117     struct Instance *i,
1118     struct gl_list_t *arglist)
1119     {
1120     int result;
1121     double step_length = 0.0;
1122     if (SensitivityAllCheckArgs(arglist,&step_length)) {
1123     return 1;
1124     }
1125     result = sensitivity_anal_all(slv_interp,i,arglist,step_length);
1126     return result;
1127     }
1128    
1129     char sensitivity_help[] =
1130     "This function does sensitivity analysis dy/dx. It requires 4 args.\n\
1131     The first arg is the name of a reference instance or SELF.\n\
1132     The second arg is x, where x is an array of > solver_var\n.\
1133     The third arg y, where y is an array of > solver_var\n. \
1134     The fourth arg is dy/dx which dy_dx[1..n_y][1..n_x].\n";
1135    
1136     #undef DEBUG

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