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

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