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

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