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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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