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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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