/[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 91 - (show annotations) (download) (as text)
Wed Dec 7 16:17:03 2005 UTC (14 years, 7 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 /*********************************************************************\
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 #if 0
67 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 #endif
76
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
167 #if 0
168 static int ReSolve(slv_system_t sys)
169 {
170 if (!sys)
171 return 1;
172 slv_solve(sys);
173 return 0;
174 }
175 #endif
176
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
203 /* Ignore unused params */
204 (void)slv_interp; (void)i; (void)whichvar;
205
206 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 /* Ignore unused params */
324 (void)slv_interp; (void)i; (void)whichvar;
325
326 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
620 #if 0
621 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 #endif
654
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 real64 **dy_dx = NULL;
664 slv_system_t sys = NULL;
665 int c;
666 int noutputs = 0;
667 int ninputs;
668 int i,j;
669 int offset;
670 dof_t *dof;
671 int num_vars,ind,found;
672
673 /* Ignore unused params */
674 (void)slv_interp; (void) inst;
675
676 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
683 (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
887 #if 0
888 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 # if DEBUG
921 FPRINTF(stderr,"Old_y = %12.8g; Nex_y = %12.8g\n",old_y,new_y);
922 # endif
923 }
924 ascfree((char *)delta_x);
925 return 0;
926 }
927 #endif
928
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 real64 **dy_dx = NULL;
960 struct var_variable **vp,**ptr;
961 slv_system_t sys = NULL;
962 long c;
963 int noutputs=0, ninputs;
964 var_filter_t vfilter;
965
966 /* Ignore unused params */
967 (void)slv_interp; (void)inst; (void)step_length;
968
969 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