/[ascend]/trunk/tcltk98/generic/interface/Sensitivity.c
ViewVC logotype

Contents of /trunk/tcltk98/generic/interface/Sensitivity.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 389 - (show annotations) (download) (as text)
Thu Mar 30 06:24:10 2006 UTC (16 years, 2 months ago) by johnpye
File MIME type: text/x-csrc
File size: 16390 byte(s)
Cleaning up #includes in the Tcl/Tk interface. Doing this
all as a group so that it can be reversed out if necessary.
1 /*
2 * Sensitivity.c
3 * by Kirk Abbott and Ben Allan
4 * Created: 1/94
5 * Version: $Revision: 1.31 $
6 * Version control file: $RCSfile: Sensitivity.c,v $
7 * Date last modified: $Date: 2003/08/23 18:43:08 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the ASCEND Tcl/Tk interface
11 *
12 * Copyright 1997, Carnegie Mellon University
13 *
14 * The ASCEND Tcl/Tk interface is free software; you can redistribute
15 * it and/or modify it under the terms of the GNU General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * The ASCEND Tcl/Tk interface is distributed in hope that it will be
20 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with the program; if not, write to the Free Software Foundation,
26 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
27 * COPYING. COPYING is found in ../compiler.
28 */
29
30
31 /*
32 * Sensititvity analysis code.
33 */
34
35 #include <math.h>
36 #include <tcl.h>
37 #include <utilities/ascConfig.h>
38 #include <utilities/ascMalloc.h>
39 #include <general/tm_time.h>
40 #include <general/list.h>
41 #include <solver/slv_types.h>
42 #include <solver/mtx.h>
43 #include <solver/var.h>
44 #include <solver/rel.h>
45 #include <solver/discrete.h>
46 #include <solver/conditional.h>
47 #include <solver/logrel.h>
48 #include <solver/bnd.h>
49 #include <solver/slv_common.h>
50 #include <solver/linsol.h>
51 #include <solver/linsolqr.h>
52 #include <solver/linutils.h>
53 #include <solver/calc.h>
54 #include <solver/relman.h>
55 #include <solver/slv_client.h>
56 #include <solver/system.h>
57 #include "HelpProc.h"
58 #include "Sensitivity.h"
59 #include "HelpProc.h"
60 #include "SolverGlobals.h"
61 #include <general/mathmacros.h>
62
63 #ifndef lint
64 static CONST char SensitivityID[] = "$Id: Sensitivity.c,v 1.31 2003/08/23 18:43:08 ballan Exp $";
65 #endif
66
67
68 #define SAFE_FIX_ME 0
69
70 #define DOTIME FALSE
71 #define DEBUG 1
72 /*
73 * This file attempts to implement the extraction of dy_dx from
74 * a system of equations. If one considers a black-box where x are
75 * the input variables, u are inuut parameters, y are the output
76 * variables, f(x,y,u) is the system of equations that will be solved
77 * for given x and u, then one can extract the sensitivity information
78 * of y wrt x.
79 * One crude but simple way of doing this is to finite-difference the
80 * given black box, i.e, perturb x, n\subx times by delta x, resolve
81 * the system of equations at the new value of x, and compute
82 * dy/dx = (f(x\sub1) - f(x\sub2))/(x\sub1 - x\sub2).
83 * This is known to be very expensive.
84 *
85 * The solution that will be adopted here is to formulate the Jacobian J of
86 * the system, (or the system is already solved, to grab the jacobian at
87 * the solution. Develop the sensitivity matrix df/dx by exact numnerical
88 * differentiation; this will be a n x n\subx matrix. Compute the LU factors
89 * for J, and to solve column by column to : LU dz/dx = df/dx. Here
90 * z, represents all the internal variables to the system and the strictly
91 * output variables y. In other words J = df/dz.
92 *
93 * Given the solution of the above problem we then simply extract the rows
94 * of dz/dx, that pertain to the y variables that we are interested in;
95 * this will give us dy/dx.
96 */
97
98 #ifdef THIS_MAY_BE_UNUSED_CODE
99 static real64 *zero_vector(real64 *vec, int size)
100 {
101 int c;
102 for (c=0;c<size;c++) {
103 vec[c] = 0.0;
104 }
105 return vec;
106 }
107 #endif
108
109 static int NumberIncludedRels(slv_system_t sys)
110 {
111 static int nrels = -1;
112 rel_filter_t rfilter;
113 int tmp;
114
115 if (!sys) { /* a NULL system may be used to reinitialise the count */
116 nrels = -1;
117 return -1;
118 }
119 if (nrels < 0) {
120 /*get used (included) relation count -- equalities only !*/
121 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
122 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
123 tmp = slv_count_solvers_rels(sys,&rfilter);
124 nrels = tmp;
125 return tmp;
126 } else {
127 return nrels;
128 }
129 }
130
131 static int NumberFreeVars(slv_system_t sys)
132 {
133 static int nvars = -1;
134 var_filter_t vfilter;
135 int tmp;
136
137 if (!sys) {
138 nvars = -1;
139 return -1;
140 }
141 if (nvars < 0) {
142 /*get used (free and incident) variable count */
143
144 vfilter.matchbits = (VAR_FIXED |VAR_INCIDENT| VAR_SVAR | VAR_ACTIVE);
145 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
146 tmp = slv_count_solvers_vars(sys,&vfilter);
147 nvars = tmp;
148 return tmp;
149 } else {
150 return nvars;
151 }
152 }
153
154
155 /*
156 *********************************************************************
157 * The following bit of code does the computation of the matrix
158 * dz/dx. It accepts a slv_system_t and a list of 'input' variables.
159 * The matrix df/dx now exists and sits to the right of the main
160 * square region of the Jacobian. We will do the following in turn
161 * for each variable in this list:
162 *
163 * 1) Get the variable index, and use it to extract the required column
164 * from the main gradient matrix, this will be stored in a temporary
165 * vector.
166 * 2) We will then clear the column in the original matrix, as we want to
167 * store the caluclated results back in place.
168 * 3) Add the vector extracted in 1) as rhs to the main matrix.
169 * 4) Call linsol solve on this rhs to yield a solution which represents
170 * a column of dx/dx.
171 * 6) Store this vector back in the location cleared out in 2).
172 * 7) Goto 1.
173 *
174 * At the end of this procedure we would have calculated all the columns of
175 * dz/dx, and left them back in the main matrix.
176 *********************************************************************
177 */
178
179
180 /*
181 * At this point we should have an empty jacobian. We now
182 * need to call relman_diff over the *entire* matrix.
183 * Calculates the entire jacobian. It is initially unscaled.
184 *
185 * Note: this assumes the sys given is using one of the ascend solvers
186 * with either linsol or linsolqr. Both are now allowed. baa 7/95
187 */
188 static int Compute_J(slv_system_t sys)
189 {
190 int32 row;
191 var_filter_t vfilter;
192 linsol_system_t lin_sys = NULL;
193 linsolqr_system_t lqr_sys = NULL;
194 mtx_matrix_t mtx;
195 struct rel_relation **rlist;
196 int nrows;
197 int qr=0;
198 #if DOTIME
199 double time1;
200 #endif
201
202 #if DOTIME
203 time1 = tm_cpu_time();
204 #endif
205 /*
206 * Get the linear system from the solve system.
207 * Get the matrix from the linear system.
208 * Get the relations list from the solve system.
209 */
210 lin_sys = slv_get_linsol_sys(sys);
211 if (lin_sys==NULL) {
212 qr=1;
213 lqr_sys=slv_get_linsolqr_sys(sys);
214 }
215 mtx = slv_get_sys_mtx(sys);
216 if (mtx==NULL || (lin_sys==NULL && lqr_sys==NULL)) {
217 FPRINTF(stderr,"Compute_dy_dx: error, found NULL.\n");
218 return 1;
219 }
220 rlist = slv_get_solvers_rel_list(sys);
221 nrows = NumberIncludedRels(sys);
222
223 calc_ok = TRUE;
224 vfilter.matchbits = VAR_SVAR;
225 vfilter.matchvalue = VAR_SVAR;
226 /*
227 * Clear the entire matrix and then compute
228 * the gradients at the current point.
229 */
230 mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
231 for(row=0; row<nrows; row++) {
232 struct rel_relation *rel;
233 /* added */
234 double resid;
235
236 rel = rlist[mtx_row_to_org(mtx,row)];
237 (void)relman_diffs(rel,&vfilter,mtx,&resid,SAFE_FIX_ME);
238
239 /* added */
240 rel_set_residual(rel,resid);
241
242 }
243 if (qr) {
244 linsolqr_matrix_was_changed(lqr_sys);
245 } else {
246 linsol_matrix_was_changed(lin_sys);
247 }
248 #if DOTIME
249 time1 = tm_cpu_time() - time1;
250 FPRINTF(stderr,"Time taken for Compute_J = %g\n",time1);
251 #endif
252 return(!calc_ok);
253 }
254
255
256 /*
257 * Note a rhs would have been previously added
258 * to keep the system happy.
259 * Now handles both linsol and linsolqr as needed.
260 * Assumes region to be factored is in upper left corner.
261 * Region is reordered using the last reorder method used in
262 * the case of linsolqr, so all linsolqr methods are available
263 * to this routine.
264 */
265 static int LUFactorJacobian(slv_system_t sys)
266 {
267 linsol_system_t lin_sys=NULL;
268 linsolqr_system_t lqr_sys=NULL;
269 mtx_region_t region;
270 int size,nvars,nrels;
271 #if DOTIME
272 double time1;
273 #endif
274
275 #if DOTIME
276 time1 = tm_cpu_time();
277 #endif
278
279 nvars = NumberFreeVars(sys);
280 nrels = NumberIncludedRels(sys);
281 size = MAX(nvars,nrels); /* get size of matrix */
282 mtx_region(&region,0,size-1,0,size-1); /* set the region */
283 lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
284 if (lin_sys!=NULL) {
285 linsol_matrix_was_changed(lin_sys);
286 linsol_reorder(lin_sys,&region); /* This was entire_MATRIX */
287 linsol_invert(lin_sys,&region); /* invert the given matrix over
288 * the given region */
289 } else {
290 enum factor_method fm;
291 int oldtiming;
292 lqr_sys = slv_get_linsolqr_sys(sys);
293 /* WE are ASSUMING that the system has been qr_preped before now! */
294 linsolqr_matrix_was_changed(lqr_sys);
295 linsolqr_reorder(lqr_sys,&region,natural); /* should retain original */
296 fm = linsolqr_fmethod(lqr_sys);
297 if (fm == unknown_f) {
298 fm = ranki_kw2; /* make sure somebody set it */
299 }
300 oldtiming = g_linsolqr_timing;
301 g_linsolqr_timing = 0;
302 linsolqr_factor(lqr_sys,fm);
303 g_linsolqr_timing = oldtiming;
304 }
305
306 #if DOTIME
307 time1 = tm_cpu_time() - time1;
308 FPRINTF(stderr,"Time taken for LUFactorJacobian = %g\n",time1);
309 #endif
310 return 0;
311 }
312
313
314 /*
315 * At this point the rhs would have already been added.
316 *
317 * Extended to handle either factorization code:
318 * linsol or linsolqr.
319 */
320 static int Compute_dy_dx_smart(slv_system_t sys,
321 real64 *rhs,
322 real64 **dy_dx,
323 int *inputs, int ninputs,
324 int *outputs, int noutputs)
325 {
326 linsol_system_t lin_sys=NULL;
327 linsolqr_system_t lqr_sys=NULL;
328 mtx_matrix_t mtx;
329 int col,current_col;
330 int row;
331 int capacity;
332 real64 *solution = NULL;
333 int i,j;
334 #if DOTIME
335 double time1;
336 #endif
337
338 #if DOTIME
339 time1 = tm_cpu_time();
340 #endif
341
342 lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
343 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
344 mtx = slv_get_sys_mtx(sys); /* get the matrix */
345
346 capacity = mtx_capacity(mtx);
347 solution = (real64 *)asccalloc(capacity,sizeof(real64));
348
349 /*
350 * The array inputs is a list of original indexes, of the variables
351 * that we are trying to obtain the sensitivity to. We have to
352 * get the *current* column from the matrix based on those indices.
353 * Hence we use mtx_org_to_col. This little manipulation is not
354 * necessary for the computed solution as the solve routine returns
355 * the results in the *original* order rather than the *current* order.
356 */
357 if (lin_sys!=NULL) {
358 for (j=0;j<ninputs;j++) {
359 col = inputs[j];
360 current_col = mtx_org_to_col(mtx,col);
361 mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS); /* get the column */
362
363 linsol_rhs_was_changed(lin_sys,rhs);
364 linsol_solve(lin_sys,rhs); /* solve */
365 linsol_copy_solution(lin_sys,rhs,solution); /* get the solution */
366
367 for (i=0;i<noutputs;i++) { /* copy the solution to dy_dx */
368 row = outputs[i];
369 dy_dx[i][j] = -1.0*solution[row]; /* the -1 comes from the lin alg */
370 }
371 /*
372 * zero the vector using the incidence pattern of our column.
373 * This is faster than the naiive approach of zeroing
374 * everything especially if the vector is large.
375 */
376 mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
377 }
378 } else {
379 for (j=0;j<ninputs;j++) {
380 col = inputs[j];
381 current_col = mtx_org_to_col(mtx,col);
382 mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS);
383
384 linsolqr_rhs_was_changed(lqr_sys,rhs);
385 linsolqr_solve(lqr_sys,rhs);
386 linsolqr_copy_solution(lqr_sys,rhs,solution);
387
388 for (i=0;i<noutputs;i++) {
389 row = outputs[i];
390 dy_dx[i][j] = -1.0*solution[row];
391 }
392 mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
393 }
394 }
395 if (solution) {
396 ascfree((char *)solution);
397 }
398
399 #if DOTIME
400 time1 = tm_cpu_time() - time1;
401 FPRINTF(stderr,"Time for Compute_dy_dx_smart = %g\n",time1);
402 #endif
403 return 0;
404 }
405
406
407 /*
408 *********************************************************************
409 * This code is provided for the benefit of a temporary
410 * fix for the derivative problem in Lsode.
411 * The proper permanent fix for lsode is to dump it in favor of
412 * cvode or dassl.
413 * Extended 7/95 baa to deal with linsolqr and blsode.
414 * It is assumed the system has been solved at the current point.
415 *********************************************************************
416 */
417 int Asc_BLsodeDerivatives(slv_system_t sys, double **dy_dx,
418 int *inputs_ndx_list, int ninputs,
419 int *outputs_ndx_list, int noutputs)
420 {
421 static int n_calls = 0;
422 linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
423 mtx_matrix_t mtx;
424 int32 capacity;
425 real64 *scratch_vector = NULL;
426 int result=0;
427
428 (void)NumberFreeVars(NULL); /* used to re-init the system */
429 (void)NumberIncludedRels(NULL); /* used to re-init the system */
430 if (!sys) {
431 FPRINTF(stderr,"The solve system does not exist !\n");
432 return 1;
433 }
434
435 result = Compute_J(sys);
436 if (result) {
437 FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
438 return 1;
439 }
440
441 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
442 if (lqr_sys==NULL) {
443 FPRINTF(stderr,"Early termination due to missing linsolqr system.\n");
444 return 1;
445 }
446 mtx = slv_get_sys_mtx(sys); /* get the matrix */
447 if (mtx==NULL) {
448 FPRINTF(stderr,"Early termination due to missing mtx in linsolqr.\n");
449 return 1;
450 }
451 capacity = mtx_capacity(mtx);
452 scratch_vector = (real64 *)asccalloc(capacity,sizeof(real64));
453 linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
454
455 result = LUFactorJacobian(sys);
456 if (result) {
457 FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
458 goto error;
459 }
460 result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
461 inputs_ndx_list, ninputs,
462 outputs_ndx_list, noutputs);
463
464 linsolqr_remove_rhs(lqr_sys,scratch_vector);
465 if (result) {
466 FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
467 goto error;
468 }
469
470 error:
471 n_calls++;
472 if (scratch_vector) {
473 ascfree((char *)scratch_vector);
474 }
475 return result;
476 }
477
478
479 int Asc_MtxNormsCmd(ClientData cdata, Tcl_Interp *interp,
480 int argc, CONST84 char *argv[])
481 {
482 slv_system_t sys;
483 linsol_system_t lin_sys = NULL;
484 linsolqr_system_t linqr_sys = NULL;
485 mtx_matrix_t mtx;
486 mtx_region_t reg;
487 double norm;
488 int solver;
489
490 (void)cdata; /* stop gcc whine about unused parameter */
491 (void)argv; /* stop gcc whine about unused parameter */
492
493 if (argc!=1) {
494 Tcl_SetResult(interp, "wrong # args: Usage __mtx_norms", TCL_STATIC);
495 return TCL_ERROR;
496 }
497 sys = g_solvsys_cur;
498 if (sys==NULL) {
499 Tcl_SetResult(interp, "__mtx_norms called with slv_system", TCL_STATIC);
500 return TCL_ERROR;
501 }
502 solver = slv_get_selected_solver(sys);
503 switch(solver) {
504 default:
505 case 0:
506 lin_sys = slv_get_linsol_sys(sys);
507 mtx = linsol_get_matrix(lin_sys);
508 reg.col.low = reg.row.low = 0;
509 reg.col.high = reg.row.high = mtx_symbolic_rank(mtx);
510 norm = linutils_A_1_norm(mtx,&reg);
511 FPRINTF(stderr,"A_1_norm = %g\n",norm);
512 norm = linutils_A_infinity_norm(mtx,&reg);
513 FPRINTF(stderr,"A_infinity_norm = %g\n",norm);
514 norm = linutils_A_Frobenius_norm(mtx,&reg);
515 FPRINTF(stderr,"A_Frobenius_norm = %g\n",norm);
516 norm = linutils_A_cond_kaa(lin_sys,mtx,NULL);
517 FPRINTF(stderr,"A_condition # = %g\n",norm);
518 break;
519 case 3:
520 case 5:
521 linqr_sys = slv_get_linsolqr_sys(sys);
522 mtx = linsolqr_get_matrix(linqr_sys);
523 reg.col.low = reg.row.low = 0;
524 reg.col.high = reg.row.high = mtx_symbolic_rank(mtx);
525 norm = linutils_A_1_norm(mtx,&reg);
526 FPRINTF(stderr,"A_1_norm = %g\n",norm);
527 norm = linutils_A_infinity_norm(mtx,&reg);
528 FPRINTF(stderr,"A_infinity_norm = %g\n",norm);
529 norm = linutils_A_Frobenius_norm(mtx,&reg);
530 FPRINTF(stderr,"A_Frobenius_norm = %g\n",norm);
531 norm = linutils_A_condqr_kaa(linqr_sys,mtx,NULL);
532 FPRINTF(stderr,"A_condition # = %g\n",norm);
533 break;
534 }
535 return TCL_OK;
536 }

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