/[ascend]/trunk/models/sensitivity/sensitivity.c
ViewVC logotype

Contents of /trunk/models/sensitivity/sensitivity.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1162 - (show annotations) (download) (as text)
Wed Jan 17 04:34:59 2007 UTC (15 years, 6 months ago) by johnpye
File MIME type: text/x-csrc
File size: 17473 byte(s)
Fixed sensitivity module. It now runs as mostly external code, with just a few relics required by LSODE still in packages/sensitivity.c.
1 #include <math.h>
2
3 #include "sensitivity.h"
4 #include <packages/sensitivity.h>
5 #include <compiler/instquery.h>
6 #include <compiler/atomvalue.h>
7 #include <utilities/ascMalloc.h>
8
9
10 /**
11 Build then presolve an instance
12 */
13 slv_system_t asc_sens_presolve(struct Instance *inst){
14 slv_system_t sys;
15 slv_parameters_t parameters;
16 struct var_variable **vp;
17 struct rel_relation **rp;
18 int ind,len;
19 char *tmp=NULL;
20
21 sys = system_build(inst);
22 if (sys==NULL) {
23 ERROR_REPORTER_HERE(ASC_PROG_ERR,
24 "Something radically wrong in creating solver.");
25 return NULL;
26 }
27 if (g_SlvNumberOfRegisteredClients == 0) {
28 return NULL;
29 }
30 ind = 0;
31 while (strcmp(slv_solver_name(ind),"QRSlv")) {
32 if (ind >= g_SlvNumberOfRegisteredClients) {
33 ERROR_REPORTER_HERE(ASC_PROG_ERR,
34 "QRSlv must be registered client.");
35 return NULL;
36 }
37 ++ind;
38 }
39 slv_select_solver(sys,ind);
40
41 slv_get_parameters(sys,&parameters);
42 parameters.partition = 0;
43 slv_set_parameters(sys,&parameters);
44 slv_presolve(sys);
45
46 #if DEBUG
47 vp = slv_get_solvers_var_list(sys);
48 len = slv_get_num_solvers_vars(sys);
49 for (ind=0 ; ind<len; ind++) {
50 tmp = var_make_name(sys,vp[ind]);
51 FPRINTF(stderr,"%s %d\n",tmp,var_sindex(vp[ind]));
52 if (tmp!=NULL) ascfree(tmp);
53 }
54 rp = slv_get_solvers_rel_list(sys);
55 len = slv_get_num_solvers_rels(sys);
56 for (ind=0 ; ind<len ; ind++) {
57 tmp = rel_make_name(sys,rp[ind]);
58 FPRINTF(stderr,"%s %d\n",tmp,rel_sindex(rp[ind]));
59 if (tmp) ascfree(tmp);
60 }
61 #endif
62 return sys;
63 }
64
65
66 int sensitivity_anal(
67 struct Instance *inst, /* not used but will be */
68 struct gl_list_t *arglist){
69 struct Instance *which_instance,*tmp_inst, *atominst;
70 struct gl_list_t *branch;
71 struct var_variable **vlist = NULL;
72 int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
73 real64 **dy_dx = NULL;
74 slv_system_t sys = NULL;
75 int c;
76 int noutputs = 0;
77 int ninputs;
78 int i,j;
79 int offset;
80 dof_t *dof;
81 int num_vars,ind,found;
82
83 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
84 mtx_matrix_t mtx;
85 int32 capacity;
86 real64 *scratch_vector = NULL;
87 int result=0;
88
89 /* Ignore unused params */
90 (void) inst;
91
92 (void)NumberFreeVars(NULL); /* used to re-init the system */
93 (void)NumberRels(NULL); /* used to re-init the system */
94 which_instance = FetchElement(arglist,1,1);
95 sys = asc_sens_presolve(which_instance);
96 if (!sys) {
97 FPRINTF(stderr,"Early termination due to failure in Presolve\n");
98 result = 1;
99 goto error;
100 }
101
102 dof = slv_get_dofdata(sys);
103 if (!(dof->n_rows == dof->n_cols &&
104 dof->n_rows == dof->structural_rank)) {
105 FPRINTF(stderr,"Early termination: non square system\n");
106 result = 1;
107 goto error;
108 }
109 /*
110 * prepare the inputs list
111 */
112 vlist = slv_get_solvers_var_list(sys);
113
114 branch = gl_fetch(arglist,2); /* input args */
115 ninputs = gl_length(branch);
116 inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
117
118 num_vars = slv_get_num_solvers_vars(sys);
119 for (c=0;c<ninputs;c++) {
120 atominst = (struct Instance *)gl_fetch(branch,c+1);
121 found = 0;
122 ind = num_vars - 1;
123 /* search backwards because fixed vars are at the end of the
124 var list */
125 while (!found && ind >= 0){
126 if (var_instance(vlist[ind]) == atominst) {
127 inputs_ndx_list[c] = var_sindex(vlist[ind]);
128 found = 1;
129 }
130 --ind;
131 }
132 if (!found) {
133 FPRINTF(stderr,"Sensitivity input variable not found\n");
134 result = 1;
135 goto error;
136 }
137 }
138
139 /*
140 * prepare the outputs list
141 */
142 branch = gl_fetch(arglist,3); /* output args */
143 noutputs = gl_length(branch);
144 outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
145 for (c=0;c<noutputs;c++) {
146 atominst = (struct Instance *)gl_fetch(branch,c+1);
147 found = 0;
148 ind = 0;
149 while (!found && ind < num_vars){
150 if (var_instance(vlist[ind]) == atominst) {
151 outputs_ndx_list[c] = var_sindex(vlist[ind]);
152 found = 1;
153 }
154 ++ind;
155 }
156 if (!found) {
157 FPRINTF(stderr,"Sensitivity ouput variable not found\n");
158 result = 1;
159 goto error;
160 }
161 }
162
163 /*
164 * prepare the results dy_dx.
165 */
166 dy_dx = make_matrix(noutputs,ninputs);
167
168
169 result = Compute_J(sys);
170 if (result) {
171 FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
172 goto error;
173 }
174
175 /*
176 * Note: the RHS *has* to added here. We will construct the vector
177 * of size matrix capacity and add it. It will be removed after
178 * we finish computing dy_dx.
179 */
180 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
181 mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
182 capacity = mtx_capacity(mtx);
183 scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
184 linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
185 result = LUFactorJacobian1(sys,dof->structural_rank);
186 if (result) {
187 FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
188 goto error;
189 }
190 result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
191 inputs_ndx_list, ninputs,
192 outputs_ndx_list, noutputs);
193
194 linsolqr_remove_rhs(lqr_sys,scratch_vector);
195 if (result) {
196 FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
197 goto error;
198 }
199
200 /*
201 * Write the results back to the partials array in the
202 * instance tree
203 */
204 offset = 0;
205 for (i=0;i<noutputs;i++) {
206 for (j=0;j<ninputs;j++) {
207 tmp_inst = FetchElement(arglist,4,offset+j+1);
208 SetRealAtomValue(tmp_inst,dy_dx[i][j],(unsigned)0);
209 #if DEBUG
210 FPRINTF(stderr,"%12.8f i%d j%d",dy_dx[i][j],i,j);
211 #endif
212 }
213 #if DEBUG
214 FPRINTF(stderr,"\n");
215 #endif
216 offset += ninputs;
217 }
218 #if DEBUG
219 FPRINTF(stderr,"\n");
220 #endif
221
222 error:
223 if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
224 if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
225 if (dy_dx) free_matrix(dy_dx,noutputs);
226 if (scratch_vector) ascfree((char *)scratch_vector);
227 if (sys) system_destroy(sys);
228 return result;
229 }
230
231 /**
232 Finite-difference Check Arguments...?
233
234 @param arglist list of arguments
235
236 Argument list contains
237 . arg1 - model inst to be solved
238 . arg2 - array of input instances
239 . arg3 - array of output instances
240 . arg4 - matrix of partials to be written to
241 */
242 static int FiniteDiffCheckArgs(struct gl_list_t *arglist)
243 {
244 struct Instance *inst;
245 unsigned long len;
246 unsigned long ninputs, noutputs;
247
248 len = gl_length(arglist);
249 if (len != 4) {
250 FPRINTF(stderr,"wrong number of args -- 4 expected\n");
251 return 1;
252 }
253 inst = FetchElement(arglist,1,1);
254 if (InstanceKind(inst)!=MODEL_INST) {
255 FPRINTF(stderr,"Arg #1 is not a model instance\n");
256 return 1;
257 }
258 ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
259 /* input args */
260 noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
261 /* output args */
262 len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
263 /* partials matrix */
264 if (len != (ninputs*noutputs)) {
265 FPRINTF(stderr,
266 "The array of partials is inconsistent with the args given\n");
267 return 1;
268 }
269 return 0;
270 }
271
272
273
274 /*-------------------------------------------------
275 SENSITIVITY ANALYSIS CODE
276 */
277
278 /**
279 Sensitivity Analysis
280
281 @param arglist List of arguments
282
283 Argument list contains the following:
284 . arg1 - model inst for which the sensitivity analysis is to be done.
285 . arg2 - array of input instances.
286 . arg3 - array of output instances.
287 . arg4 matrix of partials to be written to.
288 */
289 int SensitivityCheckArgs(struct gl_list_t *arglist)
290 {
291 struct Instance *inst;
292 unsigned long len;
293 unsigned long ninputs, noutputs;
294
295 len = gl_length(arglist);
296 if (len != 4) {
297 FPRINTF(stderr,"wrong number of args -- 4 expected\n");
298 return 1;
299 }
300 inst = FetchElement(arglist,1,1);
301 if (InstanceKind(inst)!=MODEL_INST) {
302 FPRINTF(stderr,"Arg #1 is not a model instance\n");
303 return 1;
304 }
305 ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
306 /* input args */
307 noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
308 /* output args */
309 len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
310 /* partials matrix */
311 if (len != (ninputs*noutputs)) {
312 FPRINTF(stderr,
313 "The array of partials is inconsistent with the args given\n");
314 return 1;
315 }
316 return 0;
317 }
318
319
320 /**
321 @param arglist List of arguments
322 @param step_length ...?
323
324 The list of arguments should chontain
325
326 . arg1 - model inst for which the sensitivity analysis is to be done.
327 . arg2 - array of input instances.
328 . arg3 - new_input instances, for variable projection.
329 . arg4 - instance representing the step_length for projection.
330
331 The result will be written to standard out.
332 */
333 int SensitivityAllCheckArgs(struct gl_list_t *arglist, double *step_length)
334 {
335 struct Instance *inst;
336 unsigned long len;
337
338 /*
339
340 */
341
342 len = gl_length(arglist);
343 if (len != 4) {
344 FPRINTF(stderr,"wrong number of args -- 4 expected\n");
345 return 1;
346 }
347 inst = FetchElement(arglist,1,1);
348 if (InstanceKind(inst)!=MODEL_INST) {
349 FPRINTF(stderr,"Arg #1 is not a model instance\n");
350 return 1;
351 }
352 /*
353 * we should be checking that arg2 list contains solver vars
354 * and that they are fixed. The same applies to arglist 3... later.
355 * We will check and return the steplength though. 0 means dont do
356 * the variable projection.
357 */
358 inst = FetchElement(arglist,4,1);
359 *step_length = RealAtomValue(inst);
360 if (fabs(*step_length) < 1e-08)
361 *step_length = 0.0;
362 return 0;
363 }
364
365
366 /**
367 This function is very similar to sensitivity_anal, execept,
368 that it perform the analysis on the entire system, based on
369 the given inputs. Nothing is returned. As such the call from
370 ASCEND is:
371
372 <pre>
373 EXTERN sensitivity_anal_all(
374 this_instance,
375 u_old[1..n_inputs],
376 u_new[1..n_inputs],
377 step_length
378 );</pre>
379
380 The results will be witten to standard out.
381 This function is more expensive from a a memory point of view,
382 as we keep aroung a dense matrix n_outputs * n_inputs, but here
383 n_outputs may be *much* larger depending on problem size.
384 */
385 int sensitivity_anal_all( struct Instance *inst, /* not used but will be */
386 struct gl_list_t *arglist,
387 real64 step_length)
388 {
389 struct Instance *which_instance;
390 struct gl_list_t *branch2, *branch3;
391 dof_t *dof;
392 struct var_variable **inputs = NULL, **outputs = NULL;
393 int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
394 real64 **dy_dx = NULL;
395 struct var_variable **vp,**ptr;
396 slv_system_t sys = NULL;
397 long c;
398 int noutputs=0, ninputs;
399 var_filter_t vfilter;
400
401 struct var_variable **new_inputs = NULL; /* optional stuff for variable
402 * projection */
403
404 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
405 mtx_matrix_t mtx;
406 int32 capacity;
407 real64 *scratch_vector = NULL;
408 int result=0;
409
410 /* Ignore unused params */
411 (void)inst; (void)step_length;
412
413 /*
414 * Call the presolve for the system. This should number variables
415 * and relations as well create and order the main Jacobian. The
416 * only potential problem that I see here is that presolve using
417 * the slv0 solver *only* recognizes solver vars. So that if one
418 * wanted to see the sensitivity of a *parameter*, it would not
419 * be possible. We will have to trap this in CheckArgs.
420 *
421 * Also the current version of ascend is fucked in how the var module
422 * handles variables and their numbering through the interface ptr
423 * crap.
424 */
425
426 (void)NumberFreeVars(NULL); /* used to re-init the system */
427 (void)NumberRels(NULL); /* used to re-init the system */
428 which_instance = FetchElement(arglist,1,1);
429 sys = asc_sens_presolve(which_instance);
430 if (!sys) {
431 FPRINTF(stderr,"Early termination due to failure in asc_sens_presolve\n");
432 result = 1;
433 goto error;
434 }
435 dof = slv_get_dofdata(sys);
436
437 /*
438 * prepare the inputs list. We dont need the index of the new_inputs
439 * list. We will grab them later if necessary.
440 */
441 branch2 = gl_fetch(arglist,2); /* input args -- old u_values */
442 branch3 = gl_fetch(arglist,3); /* input args -- new u_values */
443 ninputs = gl_length(branch2);
444 inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
445 new_inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
446
447 inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
448 for (c=0;c<ninputs;c++) {
449 inputs[c] = (struct var_variable *)gl_fetch(branch2,c+1);
450 inputs_ndx_list[c] = var_mindex(inputs[c]);
451 new_inputs[c] = (struct var_variable *)gl_fetch(branch3,c+1);
452 }
453 inputs[ninputs] = NULL; /* null terminate the list */
454 new_inputs[ninputs] = NULL; /* null terminate the list */
455
456 /*
457 * prepare the outputs list. This is where we differ from
458 * the other function. The noutputs, and the indexes of these
459 * outputs is obtained from the entire solve system.
460 */
461 vfilter.matchbits = 0;
462 noutputs = slv_count_solvers_vars(sys,&vfilter);
463 outputs = ASC_NEW_ARRAY(struct var_variable *,noutputs+1);
464 outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
465 ptr = vp = slv_get_solvers_var_list(sys);
466 for (c=0;c<noutputs;c++) {
467 outputs[c] = *ptr;
468 outputs_ndx_list[c] = var_sindex(*ptr);
469 ptr++;
470 }
471 outputs[noutputs] = NULL; /* null terminate the list */
472
473 /*
474 * prepare the results dy_dx. This is the expensive part from a
475 * memory point of view. However I would like to have the entire
476 * noutputs * ninputs matrix even for a short while so that I
477 * can compute a number of different types of norms.
478 */
479 dy_dx = make_matrix(noutputs,ninputs);
480
481 result = Compute_J(sys);
482 if (result) {
483 FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
484 goto error;
485 }
486
487 /*
488 * Note: the RHS *has* to added here. We will construct the vector
489 * of size matrix capacity and add it. It will be removed after
490 * we finish computing dy_dx.
491 */
492 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
493 mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
494 capacity = mtx_capacity(mtx);
495 scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
496 linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
497 result = LUFactorJacobian1(sys,dof->structural_rank);
498 if (result) {
499 FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
500 goto error;
501 }
502 result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
503 inputs_ndx_list, ninputs,
504 outputs_ndx_list, noutputs);
505
506 linsolqr_remove_rhs(lqr_sys,scratch_vector);
507 if (result) {
508 FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
509 goto error;
510 }
511
512 /*
513 * Do some analysis on the results, inclusive of
514 * writing them to someplace useful.
515 */
516 if (DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx))
517 result = 1;
518
519 /*
520 * Some experimental projection stuff -- not used now.
521 * if (DoProject_X(inputs, new_inputs, step_length,
522 * outputs, ninputs, noutputs, dy_dx))
523 * result = 1;
524 */
525
526 error:
527 if (inputs) ascfree((char *)inputs);
528 if (new_inputs) ascfree((char *)new_inputs);
529 if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
530 if (outputs) ascfree((char *)outputs);
531 if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
532 if (dy_dx) free_matrix(dy_dx,noutputs);
533 if (scratch_vector) ascfree((char *)scratch_vector);
534 if (sys) system_destroy(sys);
535 return result;
536 }
537
538
539
540 /**
541 Calls 'DoSolve'
542
543 @see DoSolve
544 */
545 int do_solve_eval( struct Instance *i,
546 struct gl_list_t *arglist, void *user_data)
547 {
548 unsigned long len;
549 int result;
550 struct Instance *inst;
551 len = gl_length(arglist);
552
553 /* Ignore unused params */
554 (void)i;
555
556 if (len!=2) {
557 ERROR_REPORTER_HERE(ASC_USER_ERROR,
558 "Wrong number of args to do_solve_eval.");
559 return 1;
560 }
561 inst = FetchElement(arglist,1,1);
562 if (!inst)
563 return 1;
564 result = DoSolve(inst);
565 return result;
566 }
567
568
569 /**
570 Finite different evaluate...
571 */
572 int do_finite_diff_eval( struct Instance *i,
573 struct gl_list_t *arglist, void *user_data)
574 {
575 int result;
576
577 /* Ignore unused params */
578 (void)i;
579
580 if (FiniteDiffCheckArgs(arglist))
581 return 1;
582 result = finite_difference(arglist);
583 return result;
584 }
585
586
587 int do_sensitivity_eval( struct Instance *i,
588 struct gl_list_t *arglist, void *user_data)
589 {
590 int result;
591 if (SensitivityCheckArgs(arglist)) {
592 return 1;
593 }
594 result = sensitivity_anal(i,arglist);
595 return result;
596 }
597
598 int do_sensitivity_eval_all( struct Instance *i,
599 struct gl_list_t *arglist, void *user_data)
600 {
601 int result;
602 double step_length = 0.0;
603 if (SensitivityAllCheckArgs(arglist,&step_length)) {
604 return 1;
605 }
606 result = sensitivity_anal_all(i,arglist,step_length);
607 return result;
608 }
609
610
611 char sensitivity_help[] =
612 "This function does sensitivity analysis dy/dx. It requires 4 args:\n"
613 " 1. name: name of a reference instance or SELF.\n"
614 " 2. x: x, where x is an array of > solver_var.\n"
615 " 3. y: where y is an array of > solver_var.\n"
616 " 4. dy/dx: which dy_dx[1..n_y][1..n_x].";
617
618 int sensitivity_register(void){
619
620 int result=0;
621
622
623 result = CreateUserFunctionMethod("do_solve",
624 do_solve_eval,
625 2,NULL,NULL,NULL); /* was 2,0,null */
626 result += CreateUserFunctionMethod("do_finite_difference",
627 do_finite_diff_eval,
628 4,NULL,NULL,NULL); /* 4,0,null */
629 result += CreateUserFunctionMethod("do_sensitivity",
630 do_sensitivity_eval,
631 4,sensitivity_help,NULL,NULL);
632 result += CreateUserFunctionMethod("do_sensitivity_all",
633 do_sensitivity_eval_all,
634 4,"See do_sensitivity_eval for details",NULL,NULL);
635
636 return result;
637 }
638

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