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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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