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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 122 - (hide annotations) (download) (as text)
Mon Dec 19 06:12:40 2005 UTC (16 years, 6 months ago) by johnpye
File MIME type: text/x-csrc
File size: 16430 byte(s)
Refactoring all MAX, MIN, ABS calls to general/mathmacros.
Adding a GCC optimisation for these macros.
1 aw0a 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 "interface/HelpProc.h"
58     #include "interface/Sensitivity.h"
59     #include "interface/HelpProc.h"
60     #include "interface/SolverGlobals.h"
61 johnpye 122 #include "general/mathmacros.h"
62 aw0a 1
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