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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2649 - (show annotations) (download) (as text)
Wed Dec 12 12:39:25 2012 UTC (9 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 25309 byte(s)
Fixing GPL header, removing postal address (rpmlint incorrect-fsf-address)
1 /* ASCEND modelling environment
2 Copyright (C) 1996-2007 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 *//** @file
17
18 The module permits sensitivity analysis in an ASCEND model. In your model
19 file, you need to create some additional arrays like so:
20
21 X[1..m] IS_A solver_var; (* outputs *)
22 U[1..n] IS_A solver_var; (* inputs *)
23 dx_du[1..m][1..n] IS_A solver_var; (* data space *)
24
25 You must then 'ARE_THE_SAME' your problem variables into the spaces in the
26 first two of the above arrays.
27
28 Then, run EXTERNAL do_sensitivity(SELF,U[1..n],X[1..m],dx_du[1..m][1..n]);
29
30 In the array dx_du[1..m][1..n], values of the corresponding sensitivity
31 derivatives will be calculated.
32
33 A test/example model is available in the models directory, called
34 sensitivity_test.a4c.
35
36 --- implementation notes ---
37
38 This file attempts to implement the extraction of dy_dx from
39 a system of equations. If one considers a black-box where x are
40 the input variables, u are input parameters, y are the output
41 variables, f(x,y,u) is the system of equations that will be solved
42 for given x and u, then one can extract the sensitivity information
43 of y wrt x.
44
45 One crude but simple way of doing this is to finite-difference the
46 given black box, i.e, perturb x, n_x times by delta x, resolve
47 the system of equations at the new value of x, and compute
48
49 dy/dx = (f(x_1) - f(x_2))/(x_1 - x_2).
50
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_x 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 @todo There are files in tcltk called Sensitivity.[ch]. Do we need them?
66 */
67
68 #include <math.h>
69
70 #include <ascend/general/ascMalloc.h>
71 #include <ascend/general/mathmacros.h>
72
73 #include <ascend/compiler/instquery.h>
74 #include <ascend/compiler/atomvalue.h>
75 #include <ascend/compiler/extfunc.h>
76
77 #include <ascend/linear/densemtx.h>
78
79 #include <ascend/packages/sensitivity.h>
80 #include <ascend/system/system.h>
81 #include <ascend/solver/solver.h>
82
83 #define SENSITIVITY_DEBUG
84
85 ASC_EXPORT int sensitivity_register(void);
86
87 ExtMethodRun do_sensitivity_eval;
88 ExtMethodRun do_sensitivity_eval_all;
89
90 /**
91 Build then presolve an instance
92 */
93 slv_system_t sens_presolve(struct Instance *inst){
94 slv_system_t sys;
95 slv_parameters_t parameters;
96 const SlvFunctionsT *S;
97 #ifdef SENSITIVITY_DEBUG
98 struct var_variable **vp;
99 struct rel_relation **rp;
100 int len, ind;
101 char *tmp=NULL;
102 #endif
103 sys = system_build(inst);
104 if (sys==NULL) {
105 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to build system.");
106 return NULL;
107 }
108 S = solver_engine_named("QRSlv");
109 if(!S){
110 ERROR_REPORTER_HERE(ASC_PROG_ERR,"QRSlv solver not found (required for sensitivity)");
111 return NULL;
112 }
113 slv_select_solver(sys,S->number);
114
115 slv_get_parameters(sys,&parameters);
116 parameters.partition = 0;
117 slv_set_parameters(sys,&parameters);
118 slv_presolve(sys);
119
120 #ifdef SENSITIVITY_DEBUG
121 vp = slv_get_solvers_var_list(sys);
122 len = slv_get_num_solvers_vars(sys);
123 for (ind=0 ; ind<len; ind++) {
124 tmp = var_make_name(sys,vp[ind]);
125 CONSOLE_DEBUG("%s %d\n",tmp,var_sindex(vp[ind]));
126 if (tmp!=NULL) ascfree(tmp);
127 }
128 rp = slv_get_solvers_rel_list(sys);
129 len = slv_get_num_solvers_rels(sys);
130 for (ind=0 ; ind<len ; ind++) {
131 tmp = rel_make_name(sys,rp[ind]);
132 CONSOLE_DEBUG("%s %d\n",tmp,rel_sindex(rp[ind]));
133 if (tmp) ascfree(tmp);
134 }
135 #endif
136 return sys;
137 }
138
139 /*
140 LU Factor Jacobian
141
142 @note The RHS will have been will already have been added before
143 calling this function.
144
145 @NOTE there is another version of this floating around in packages/senstivity.c. The
146 other one uses dense matrices so probably this one's better?
147 */
148 int LUFactorJacobian1(slv_system_t sys,int rank){
149 linsolqr_system_t lqr_sys;
150 mtx_region_t region;
151 enum factor_method fm;
152
153 mtx_region(&region,0,rank-1,0,rank-1); /* set the region */
154 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
155
156 linsolqr_matrix_was_changed(lqr_sys);
157 linsolqr_reorder(lqr_sys,&region,natural);
158
159 fm = linsolqr_fmethod(lqr_sys);
160 if (fm == unknown_f) fm = ranki_kw2; /* make sure somebody set it */
161 linsolqr_factor(lqr_sys,fm);
162
163 return 0;
164 }
165
166 int sensitivity_anal(
167 struct Instance *inst, /* not used but will be */
168 struct gl_list_t *arglist
169 ){
170 struct Instance *which_instance,*tmp_inst, *atominst;
171 struct gl_list_t *branch;
172 struct var_variable **vlist = NULL;
173 int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
174 DenseMatrix dy_dx = densematrix_create_empty();
175 slv_system_t sys = NULL;
176 int c;
177 int noutputs = 0;
178 int ninputs;
179 int i,j;
180 int offset;
181 dof_t *dof;
182 int num_vars,ind,found;
183
184 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
185 mtx_matrix_t mtx;
186 int32 capacity;
187 real64 *scratch_vector = NULL;
188 int result=0;
189
190 /* Ignore unused params */
191 (void) inst;
192
193 (void)NumberFreeVars(NULL); /* used to re-init the system */
194 (void)NumberRels(NULL); /* used to re-init the system */
195 which_instance = FetchElement(arglist,1,1);
196 sys = sens_presolve(which_instance);
197 if (!sys) {
198 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to presolve");
199 result = 1;
200 goto finish;
201 }
202
203 dof = slv_get_dofdata(sys);
204 if (!(dof->n_rows == dof->n_cols &&
205 dof->n_rows == dof->structural_rank)) {
206 ERROR_REPORTER_HERE(ASC_PROG_ERR,"System is not square");
207 result = 1;
208 goto finish;
209 }
210
211 CONSOLE_DEBUG("Presolved, square");
212
213 /*
214 prepare the inputs list
215 */
216 vlist = slv_get_solvers_var_list(sys);
217
218 branch = gl_fetch(arglist,2); /* input args */
219 ninputs = gl_length(branch);
220 inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
221
222 num_vars = slv_get_num_solvers_vars(sys);
223 for (c=0;c<ninputs;c++) {
224 atominst = (struct Instance *)gl_fetch(branch,c+1);
225 found = 0;
226 ind = num_vars - 1;
227 /* search backwards because fixed vars are at the end of the
228 var list */
229 while (!found && ind >= 0){
230 if (var_instance(vlist[ind]) == atominst) {
231 inputs_ndx_list[c] = var_sindex(vlist[ind]);
232 found = 1;
233 }
234 --ind;
235 }
236 if(!found){
237 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to find sensitivity input variable");
238 result = 1;
239 goto finish;
240 }
241 }
242
243 CONSOLE_DEBUG("%d inputs",ninputs);
244
245 /*
246 prepare the outputs list
247 */
248 branch = gl_fetch(arglist,3); /* output args */
249 noutputs = gl_length(branch);
250 outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
251 for (c=0;c<noutputs;c++) {
252 atominst = (struct Instance *)gl_fetch(branch,c+1);
253 found = 0;
254 ind = 0;
255 while (!found && ind < num_vars){
256 if (var_instance(vlist[ind]) == atominst) {
257 outputs_ndx_list[c] = var_sindex(vlist[ind]);
258 found = 1;
259 }
260 ++ind;
261 }
262 if (!found) {
263 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to find sensitivity ouput variable");
264 result = 1;
265 goto finish;
266 }
267 }
268
269 CONSOLE_DEBUG("%d outputs",noutputs);
270
271 /*
272 prepare the results dy_dx.
273 */
274 dy_dx = densematrix_create(noutputs,ninputs);
275
276 result = Compute_J(sys);
277 if (result) {
278 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to calculate Jacobian failure in calc Jacobian\n");
279 goto finish;
280 }
281
282 CONSOLE_DEBUG("Computed Jacobian");
283
284 /*
285 Note: the RHS *has* to added here. We will construct the vector
286 of size matrix capacity and add it. It will be removed after
287 we finish computing dy_dx.
288 */
289 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
290 mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
291 capacity = mtx_capacity(mtx);
292 scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
293 linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
294 result = LUFactorJacobian1(sys,dof->structural_rank);
295 if(result){
296 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to calculate LUFactorJacobian");
297 goto finish;
298 }
299 result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
300 inputs_ndx_list, ninputs,
301 outputs_ndx_list, noutputs
302 );
303
304 linsolqr_remove_rhs(lqr_sys,scratch_vector);
305 if (result) {
306 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed Compute_dy_dx");
307 goto finish;
308 }
309
310 CONSOLE_DEBUG("Computed dy/dx");
311
312 /*
313 Write the results back to the partials array in the
314 instance tree
315 */
316 offset = 0;
317 for (i=0;i<noutputs;i++) {
318 for (j=0;j<ninputs;j++) {
319 tmp_inst = FetchElement(arglist,4,offset+j+1);
320 SetRealAtomValue(tmp_inst,DENSEMATRIX_ELEM(dy_dx,i,j),(unsigned)0);
321 #ifdef SENSITIVITY_DEBUG
322 CONSOLE_DEBUG("%12.8f i%d j%d",DENSEMATRIX_ELEM(dy_dx,i,j),i,j);
323 #endif
324 }
325 #ifdef SENSITIVITY_DEBUG
326 CONSOLE_DEBUG("\n");
327 #endif
328 offset += ninputs;
329 }
330 #ifdef SENSITIVITY_DEBUG
331 CONSOLE_DEBUG("\n");
332 #endif
333
334 ERROR_REPORTER_HERE(ASC_USER_SUCCESS
335 ,"Sensitivity results for %d vars were written back to the model"
336 ,noutputs
337 );
338
339 finish:
340 if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
341 if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
342 densematrix_destroy(dy_dx);
343 if (scratch_vector) ascfree((char *)scratch_vector);
344 if (sys) system_destroy(sys);
345 return result;
346 }
347
348 /**
349 Do Data Analysis. Used by sensitivity_anal_all.
350 */
351 int DoDataAnalysis(struct var_variable **inputs,
352 struct var_variable **outputs,
353 int ninputs, int noutputs,
354 DenseMatrix dy_dx)
355 {
356 FILE *fp;
357 double *norm_2, *norm_1;
358 double input_nominal,maxvalue,sum;
359 int i,j;
360
361 norm_1 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
362 norm_2 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
363
364 fp = fopen("sensitivity.out","w+");
365 if (!fp) return 0;
366
367 /*
368 * calculate the 1 and 2 norms; cache them away so that we
369 * can pretty print them. Style is everything !.
370 *
371 */
372 for (j=0;j<ninputs;j++) {
373 input_nominal = var_nominal(inputs[j]);
374 maxvalue = sum = 0;
375 for (i=0;i<noutputs;i++) {
376 DENSEMATRIX_ELEM(dy_dx,i,j) *= input_nominal/var_nominal(outputs[i]);
377 maxvalue = MAX(fabs(DENSEMATRIX_ELEM(dy_dx,i,j)),maxvalue);
378 sum += DENSEMATRIX_ELEM(dy_dx,i,j)*DENSEMATRIX_ELEM(dy_dx,i,j);
379 }
380 norm_1[j] = maxvalue;
381 norm_2[j] = sum;
382 }
383
384 for (j=0;j<ninputs;j++) { /* print the var_index */
385 fprintf(fp,"%8d ",var_mindex(inputs[j]));
386 }
387 fprintf(fp,"\n");
388
389 for (j=0;j<ninputs;j++) { /* print the 1 norms */
390 fprintf(fp,"%-#18.8f ",norm_1[j]);
391 }
392 fprintf(fp,"\n");
393
394 for (j=0;j<ninputs;j++) { /* print the 2 norms */
395 fprintf(fp,"%-#18.8f ",norm_2[j]);
396 }
397 fprintf(fp,"\n\n");
398 ascfree((char *)norm_1);
399 ascfree((char *)norm_2);
400
401 for (i=0;i<noutputs;i++) { /* print the scaled data */
402 for (j=0;j<ninputs;j++) {
403 fprintf(fp,"%-#18.8f %-4d",DENSEMATRIX_ELEM(dy_dx,i,j),i);
404 }
405 if (var_fixed(outputs[i]))
406 fprintf(fp," **fixed*** \n");
407 else
408 PUTC('\n',fp);
409 }
410 fprintf(fp,"\n");
411 fclose(fp);
412 return 0;
413 }
414
415
416 /**
417 This function is very similar to sensitivity_anal, execept,
418 that it perform the analysis on the entire system, based on
419 the given inputs. Nothing is returned. As such the call from
420 ASCEND is:
421
422 <pre>
423 EXTERN sensitivity_anal_all(
424 this_instance,
425 u_old[1..n_inputs],
426 u_new[1..n_inputs],
427 step_length
428 );</pre>
429
430 The results will be witten to standard out.
431 This function is more expensive from a a memory point of view,
432 as we keep around a dense matrix n_outputs * n_inputs, but here
433 n_outputs may be *much* larger depending on problem size.
434 */
435 int sensitivity_anal_all( struct Instance *inst, /* not used but will be */
436 struct gl_list_t *arglist,
437 real64 step_length
438 ){
439 struct Instance *which_instance;
440 struct gl_list_t *branch2, *branch3;
441 dof_t *dof;
442 struct var_variable **inputs = NULL, **outputs = NULL;
443 int *inputs_ndx_list = NULL, *outputs_ndx_list = NULL;
444 DenseMatrix dy_dx;
445 struct var_variable **vp,**ptr;
446 slv_system_t sys = NULL;
447 long c;
448 int noutputs=0, ninputs;
449 var_filter_t vfilter;
450
451 struct var_variable **new_inputs = NULL; /* optional stuff for variable
452 * projection */
453
454 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
455 mtx_matrix_t mtx;
456 int32 capacity;
457 real64 *scratch_vector = NULL;
458 int result=0;
459
460 /* Ignore unused params */
461 (void)inst; (void)step_length;
462
463 /*
464 * Call the presolve for the system. This should number variables
465 * and relations as well create and order the main Jacobian. The
466 * only potential problem that I see here is that presolve using
467 * the slv0 solver *only* recognizes solver vars. So that if one
468 * wanted to see the sensitivity of a *parameter*, it would not
469 * be possible. We will have to trap this in CheckArgs.
470 *
471 * Also the current version of ascend is fucked in how the var module
472 * handles variables and their numbering through the interface ptr
473 * crap.
474 */
475
476 (void)NumberFreeVars(NULL); /* used to re-init the system */
477 (void)NumberRels(NULL); /* used to re-init the system */
478 which_instance = FetchElement(arglist,1,1);
479 sys = sens_presolve(which_instance);
480 if (!sys) {
481 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to presolve");
482 result = 1;
483 goto finish;
484 }
485 dof = slv_get_dofdata(sys);
486
487 /*
488 * prepare the inputs list. We dont need the index of the new_inputs
489 * list. We will grab them later if necessary.
490 */
491 branch2 = gl_fetch(arglist,2); /* input args -- old u_values */
492 branch3 = gl_fetch(arglist,3); /* input args -- new u_values */
493 ninputs = gl_length(branch2);
494 inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
495 new_inputs = ASC_NEW_ARRAY(struct var_variable *,ninputs+1);
496
497 inputs_ndx_list = ASC_NEW_ARRAY(int,ninputs);
498 for (c=0;c<ninputs;c++) {
499 inputs[c] = (struct var_variable *)gl_fetch(branch2,c+1);
500 inputs_ndx_list[c] = var_mindex(inputs[c]);
501 new_inputs[c] = (struct var_variable *)gl_fetch(branch3,c+1);
502 }
503 inputs[ninputs] = NULL; /* null terminate the list */
504 new_inputs[ninputs] = NULL; /* null terminate the list */
505
506 /*
507 * prepare the outputs list. This is where we differ from
508 * the other function. The noutputs, and the indexes of these
509 * outputs is obtained from the entire solve system.
510 */
511 vfilter.matchbits = 0;
512 noutputs = slv_count_solvers_vars(sys,&vfilter);
513 outputs = ASC_NEW_ARRAY(struct var_variable *,noutputs+1);
514 outputs_ndx_list = ASC_NEW_ARRAY(int,noutputs);
515 ptr = vp = slv_get_solvers_var_list(sys);
516 for (c=0;c<noutputs;c++) {
517 outputs[c] = *ptr;
518 outputs_ndx_list[c] = var_sindex(*ptr);
519 ptr++;
520 }
521 outputs[noutputs] = NULL; /* null terminate the list */
522
523 /*
524 * prepare the results dy_dx. This is the expensive part from a
525 * memory point of view. However I would like to have the entire
526 * noutputs * ninputs matrix even for a short while so that I
527 * can compute a number of different types of norms.
528 */
529 dy_dx = densematrix_create(noutputs,ninputs);
530
531 result = Compute_J(sys);
532 if (result) {
533 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to compute Jacobian");
534 goto finish;
535 }
536
537 /*
538 Note: the RHS *has* to added here. We will construct the vector
539 of size matrix capacity and add it. It will be removed after
540 we finish computing dy_dx.
541 */
542 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
543 mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
544 capacity = mtx_capacity(mtx);
545 scratch_vector = ASC_NEW_ARRAY_CLEAR(real64,capacity);
546 linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
547 result = LUFactorJacobian1(sys,dof->structural_rank);
548
549 if (result) {
550 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failure in LUFactorJacobian");
551 goto finish;
552 }
553
554 result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
555 inputs_ndx_list, ninputs,
556 outputs_ndx_list, noutputs
557 );
558
559 linsolqr_remove_rhs(lqr_sys,scratch_vector);
560 if (result) {
561 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failure in Compute_dy_dx");
562 goto finish;
563 }
564
565 /*
566 * Do some analysis on the results, inclusive of
567 * writing them to someplace useful.
568 */
569 if(DoDataAnalysis(inputs, outputs, ninputs, noutputs, dy_dx)){
570 result = 1;
571 }
572 /*
573 * Some experimental projection stuff -- not used now.
574 * if (DoProject_X(inputs, new_inputs, step_length,
575 * outputs, ninputs, noutputs, dy_dx))
576 * result = 1;
577 */
578
579 finish:
580 if (inputs) ascfree((char *)inputs);
581 if (new_inputs) ascfree((char *)new_inputs);
582 if (inputs_ndx_list) ascfree((char *)inputs_ndx_list);
583 if (outputs) ascfree((char *)outputs);
584 if (outputs_ndx_list) ascfree((char *)outputs_ndx_list);
585 densematrix_destroy(dy_dx);;
586 if (scratch_vector) ascfree((char *)scratch_vector);
587 if (sys) system_destroy(sys);
588 return result;
589 }
590
591
592
593
594 #if 0
595 static int ComputeInverse(slv_system_t sys,
596 real64 *rhs)
597 {
598 linsolqr_system_t lqr_sys;
599 mtx_matrix_t mtx;
600 int capacity,order;
601 real64 *solution = NULL;
602 int j,k;
603
604 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
605 mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
606
607 capacity = mtx_capacity(mtx);
608 zero_vector(rhs,capacity); /* zero the rhs */
609 solution = ASC_NEW_ARRAY_CLEAR(real64,capacity);
610
611 order = mtx_order(mtx);
612 for (j=0;j<order;j++) {
613 rhs[j] = 1.0;
614 linsolqr_rhs_was_changed(lqr_sys,rhs);
615 linsolqr_solve(lqr_sys,rhs); /* solve */
616 linsolqr_copy_solution(lqr_sys,rhs,solution); /* get the solution */
617
618 CONSOLE_DEBUG("This is the rhs and solution for rhs #%d\n",j);
619 for (k=0;k<order;k++) {
620 CONSOLE_DEBUG("%12.8g %12.8g\n",rhs[k],solution[k]);
621 }
622 rhs[j] = 0.0;
623 }
624 if (solution) ascfree((char *)solution);
625 return 0;
626 }
627 #endif
628
629 #if 0
630 static int DoProject_X(struct var_variable **old_inputs,
631 struct var_variable **new_inputs, /* new values of u */
632 double step_length,
633 struct var_variable **outputs,
634 int ninputs, int noutputs,
635 DenseMatrix dy_dx)
636 {
637 struct var_variable *var;
638 real64 old_y, new_y, tmp;
639 real64 *delta_x;
640 int i,j;
641
642 delta_x = ASC_NEW_ARRAY_CLEAR(real64,ninputs);
643
644 for (j=0;j<ninputs;j++) {
645 delta_x[j] = var_value(new_inputs[j]) - var_value(old_inputs[j]);
646 /* delta_x[j] = RealAtomValue(new_inputs[j]) - RealAtomValue(old_inputs[j]); */
647 }
648
649 for (i=0;i<noutputs;i++) {
650 var = outputs[i];
651 if (var_fixed(var) || !var_active(var)) /* project only the free vars */
652 continue;
653 tmp = 0.0;
654 for (j=0;j<ninputs;j++) {
655 tmp += (DENSEMATRIX_ELEM(dy_dx,i,j)* delta_x[j]);
656 }
657 /* old_y = RealAtomValue(var); */
658 old_y = var_value(var);
659 new_y = old_y + step_length*tmp;
660 /* SetRealAtomValue(var,new_y,(unsigned)0); */
661 var_set_value(var,new_y);
662 # \if DEBUG
663 CONSOLE_DEBUG("Old_y = %12.8g; Nex_y = %12.8g\n",old_y,new_y);
664 # \endif
665 }
666 ascfree((char *)delta_x);
667 return 0;
668 }
669 #endif
670
671
672 #if 0
673 /**
674 * At this point we should have an empty jacobian. We now
675 * need to call relman_diff over the *entire* matrix.
676 * Calculates the entire jacobian. It is initially unscaled.
677 *
678 * Note: this assumes the sys given is using one of the ascend solvers
679 * with either linsol or linsolqr. Both are now allowed. baa 7/95
680 */
681 #define SAFE_FIX_ME 0
682 static int Compute_J_OLD(slv_system_t sys)
683 {
684 int32 row;
685 var_filter_t vfilter;
686 linsol_system_t lin_sys = NULL;
687 linsolqr_system_t lqr_sys = NULL;
688 mtx_matrix_t mtx;
689 struct rel_relation **rlist;
690 int nrows;
691 int qr=0;
692 #\if DOTIME
693 double time1;
694 #\endif
695
696 #\if DOTIME
697 time1 = tm_cpu_time();
698 #\endif
699 /*
700 * Get the linear system from the solve system.
701 * Get the matrix from the linear system.
702 * Get the relations list from the solve system.
703 */
704 lin_sys = slv_get_linsol_sys(sys);
705 if (lin_sys==NULL) {
706 qr=1;
707 lqr_sys=slv_get_linsolqr_sys(sys);
708 }
709 mtx = slv_get_sys_mtx(sys);
710 if (mtx==NULL || (lin_sys==NULL && lqr_sys==NULL)) {
711 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Compute_dy_dx: error, found NULL.\n");
712 return 1;
713 }
714 rlist = slv_get_solvers_rel_list(sys);
715 nrows = NumberIncludedRels(sys);
716
717 calc_ok = TRUE;
718 vfilter.matchbits = VAR_SVAR;
719 vfilter.matchvalue = VAR_SVAR;
720 /*
721 * Clear the entire matrix and then compute
722 * the gradients at the current point.
723 */
724 mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
725 for(row=0; row<nrows; row++) {
726 struct rel_relation *rel;
727 /* added */
728 double resid;
729
730 rel = rlist[mtx_row_to_org(mtx,row)];
731 (void)relman_diffs(rel,&vfilter,mtx,&resid,SAFE_FIX_ME);
732
733 /* added */
734 rel_set_residual(rel,resid);
735
736 }
737 if (qr) {
738 linsolqr_matrix_was_changed(lqr_sys);
739 } else {
740 linsol_matrix_was_changed(lin_sys);
741 }
742 #\if DOTIME
743 time1 = tm_cpu_time() - time1;
744 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Time taken for Compute_J = %g\n",time1);
745 #\endif
746 return(!calc_ok);
747 }
748 #endif
749
750
751 /**
752 Sensitivity Analysis
753
754 @param arglist List of arguments
755
756 Argument list contains the following:
757 . arg1 - model inst for which the sensitivity analysis is to be done.
758 . arg2 - array of input instances.
759 . arg3 - array of output instances.
760 . arg4 matrix of partials to be written to.
761 */
762 int SensitivityCheckArgs(struct gl_list_t *arglist){
763 struct Instance *inst;
764 unsigned long len;
765 unsigned long ninputs, noutputs;
766
767 len = gl_length(arglist);
768 if (len != 4) {
769 ERROR_REPORTER_HERE(ASC_PROG_ERR,"wrong number of args -- 4 expected\n");
770 return 1;
771 }
772 inst = FetchElement(arglist,1,1);
773 if (InstanceKind(inst)!=MODEL_INST) {
774 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Arg #1 is not a model instance\n");
775 return 1;
776 }
777 ninputs = gl_length((struct gl_list_t *)gl_fetch(arglist,2));
778 /* input args */
779 noutputs = gl_length((struct gl_list_t *)gl_fetch(arglist,3));
780 /* output args */
781 len = gl_length((struct gl_list_t *)gl_fetch(arglist,4));
782 /* partials matrix */
783 if (len != (ninputs*noutputs)) {
784 ERROR_REPORTER_HERE(ASC_PROG_ERR,
785 "The array of partials is inconsistent with the args given\n");
786 return 1;
787 }
788 return 0;
789 }
790
791
792 int do_sensitivity_eval( struct Instance *i,
793 struct gl_list_t *arglist, void *user_data
794 ){
795 CONSOLE_DEBUG("Starting sensitivity analysis...");
796 if(SensitivityCheckArgs(arglist))return 1;
797
798 return sensitivity_anal(i,arglist);
799 }
800
801
802 /**
803 @param arglist List of arguments
804 @param step_length ...?
805
806 The list of arguments should chontain
807
808 . arg1 - model inst for which the sensitivity analysis is to be done.
809 . arg2 - array of input instances.
810 . arg3 - new_input instances, for variable projection.
811 . arg4 - instance representing the step_length for projection.
812
813 The result will be written to standard out.
814 */
815 int SensitivityAllCheckArgs(struct gl_list_t *arglist, double *step_length)
816 {
817 struct Instance *inst;
818 unsigned long len;
819
820 len = gl_length(arglist);
821 if (len != 4) {
822 ERROR_REPORTER_HERE(ASC_PROG_ERR,"wrong number of args -- 4 expected\n");
823 return 1;
824 }
825 inst = FetchElement(arglist,1,1);
826 if (InstanceKind(inst)!=MODEL_INST) {
827 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Arg #1 is not a model instance\n");
828 return 1;
829 }
830 /*
831 * we should be checking that arg2 list contains solver vars
832 * and that they are fixed. The same applies to arglist 3... later.
833 * We will check and return the steplength though. 0 means dont do
834 * the variable projection.
835 */
836 inst = FetchElement(arglist,4,1);
837 *step_length = RealAtomValue(inst);
838 if (fabs(*step_length) < 1e-08)
839 *step_length = 0.0;
840 return 0;
841 }
842
843
844 int do_sensitivity_eval_all( struct Instance *i,
845 struct gl_list_t *arglist, void *user_data
846 ){
847 int result;
848 double step_length = 0.0;
849 if (SensitivityAllCheckArgs(arglist,&step_length)) {
850 return 1;
851 }
852 result = sensitivity_anal_all(i,arglist,step_length);
853 return result;
854 }
855
856
857 const char sensitivity_help[] =
858 "This function does sensitivity analysis dy/dx. It requires 4 args:\n\n"
859 " 1. name: name of a reference instance or SELF.\n"
860 " 2. x: x, where x is an array of > solver_var.\n"
861 " 3. y: where y is an array of > solver_var.\n"
862 " 4. dy/dx: which dy_dx[1..n_y][1..n_x].\n\n"
863 "See also sensitivity_anal_all.";
864
865 /** @TODO document what 'u_new' is all about...? */
866
867 const char sensitivity_all_help[] =
868 "Analyse the sensitivity of *all* variables in the system with respect\n"
869 "to the specific set of inputs 'u'. Instead of returning values to a\n"
870 "a special array inside the model, the results are written to the\n"
871 "console. Usage example:\n\n"
872 "EXTERN sensitivity_anal_all(\n"
873 " this_instance,\n"
874 " u_old[1..n_inputs],\n"
875 " u_new[1..n_inputs],\n"
876 " step_length\n"
877 ");\n\n"
878 "See also sensitivity_anal.";
879
880 int sensitivity_register(void){
881 int result=0;
882
883 result += CreateUserFunctionMethod("do_sensitivity",
884 do_sensitivity_eval,
885 4,sensitivity_help,NULL,NULL
886 );
887 result += CreateUserFunctionMethod("do_sensitivity_all",
888 do_sensitivity_eval_all,
889 4,sensitivity_all_help,NULL,NULL
890 );
891
892 return result;
893 }
894
895 /* :ex: set ts=4: */

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