/[ascend]/trunk/base/generic/packages/sensitivity.c
ViewVC logotype

Contents of /trunk/base/generic/packages/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: 11589 byte(s)
Refactored sensitivity module a little further, added self_test.
1 /* ASCEND modelling environment
2 Copyright (C) 1990-2006 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 Sensititvity analysis code
20
21 @see documentation in sensitivity.h.
22 *//*
23 by Kirk Abbott
24 */
25
26 #include "sensitivity.h"
27
28 #include <stdio.h>
29 #include <string.h>
30 #include <math.h>
31 #include <utilities/error.h>
32 #include <utilities/ascMalloc.h>
33 #include <utilities/ascPanic.h>
34 #include <compiler/packages.h>
35 #include <compiler/fractions.h>
36 #include <compiler/dimen.h>
37 #include <compiler/atomvalue.h>
38 #include <compiler/instquery.h>
39 #include <general/mathmacros.h>
40
41 #define DEBUG 1
42
43 #if 0
44 static real64 *zero_vector(real64 *vec, int size)
45 {
46 int c;
47 for (c=0;c<size;c++) {
48 vec[c] = 0.0;
49 }
50 return vec;
51 }
52 #endif
53
54 /*----------------------------------------------------
55 UTILITY ROUTINES
56 */
57
58 /**
59 Allocate memory for a matrix
60 @param nrows Number of rows
61 @param ncols Number of colums
62 @return Pointer to the allocated matrix memory location
63 */
64 real64 **make_matrix(int nrows, int ncols){
65 real64 **result;
66 int i;
67 result = (real64 **)calloc(nrows,sizeof(real64*));
68 for (i=0;i<nrows;i++) {
69 result[i] = (real64 *)calloc(ncols,sizeof(real64));
70 }
71 return result;
72 }
73
74 /**
75 Free a matrix from memory
76 @param matrix Memory location for the matrix
77 @param nrows Number of rows in the matrix
78 */
79 void free_matrix(real64 **matrix, int nrows){
80 int i;
81 if (!matrix)
82 return;
83 for (i=0;i<nrows;i++) {
84 if (matrix[i]) {
85 free(matrix[i]);
86 matrix[i] = NULL;
87 }
88 }
89 free(matrix);
90 }
91
92 /**
93 Fetch an element from a branch of an arglist...?
94 */
95 struct Instance *FetchElement(struct gl_list_t *arglist,
96 unsigned long whichbranch,
97 unsigned long whichelement){
98 struct gl_list_t *branch;
99 struct Instance *element;
100
101 if (!arglist) return NULL;
102 branch = (struct gl_list_t *)gl_fetch(arglist,whichbranch);
103 element = (struct Instance *)gl_fetch(branch,whichelement);
104 return element;
105 }
106
107 /**
108 Looks like it returns the number of relations in a solver's system
109
110 @param sys The system in question
111 @return int Number of relations in the system
112 @see slv_count_solvers_rels
113 */
114 int NumberRels(slv_system_t sys){
115 static int nrels = -1;
116 rel_filter_t rfilter;
117 int tmp;
118
119 if (!sys) { /* a NULL system may be used to reinitialise the count */
120 nrels = -1;
121 return -1;
122 }
123 if (nrels < 0) {
124 /*get used (included) relation count -- equalities only !*/
125 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
126 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
127 tmp = slv_count_solvers_rels(sys,&rfilter);
128 nrels = tmp;
129 return tmp;
130 }
131 else return nrels;
132 }
133
134 int NumberIncludedRels(slv_system_t sys){
135 static int nrels = -1;
136 rel_filter_t rfilter;
137 int tmp;
138
139 if (!sys) { /* a NULL system may be used to reinitialise the count */
140 nrels = -1;
141 return -1;
142 }
143 if (nrels < 0) {
144 /*get used (included) relation count -- equalities only !*/
145 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
146 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
147 tmp = slv_count_solvers_rels(sys,&rfilter);
148 nrels = tmp;
149 return tmp;
150 } else {
151 return nrels;
152 }
153 }
154
155 /**
156 Looks like it returns the number of free variables in a solver's system
157
158 @param sys The system in question
159 @return the number of free variables
160
161 @see slv_count_solvers_vars
162 */
163 int NumberFreeVars(slv_system_t sys){
164 static int nvars = -1;
165 var_filter_t vfilter;
166 int tmp;
167
168 if (!sys) {
169 /* no system: return error */
170 nvars = -1;
171 return -1;
172 }
173
174 if (nvars >=0){
175 /* return stored value */
176 return nvars;
177 }
178
179 /* get used (free and incident) variable count */
180 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
181 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
182 tmp = slv_count_solvers_vars(sys,&vfilter);
183 nvars = tmp;
184 return tmp;
185 }
186
187 /*
188 * The following bit of code does the computation of the matrix
189 * dz/dx. It accepts a slv_system_t and a list of 'input' variables.
190 * The matrix df/dx now exists and sits to the right of the main
191 * square region of the Jacobian. We will do the following in turn
192 * for each variable in this list:
193 *
194 * 1) Get the variable index, and use it to extract the required column
195 * from the main gradient matrix, this will be stored in a temporary
196 * vector.
197 * 2) We will then clear the column in the original matrix, as we want to
198 * store the caluclated results back in place.
199 * 3) Add the vector extracted in 1) as rhs to the main matrix.
200 * 4) Call linsol solve on this rhs to yield a solution which represents
201 * a column of dx/dx.
202 * 6) Store this vector back in the location cleared out in 2).
203 * 7) Goto 1.
204 *
205 * At the end of this procedure we would have calculated all the columns of
206 * dz/dx, and left them back in the main matrix.
207 */
208
209
210 /*
211 * At this point we should have an empty jacobian. We now
212 * need to call relman_diff over the *entire* matrix.
213 * fixed and free.
214 *
215 * Calculates the entire jacobian.
216 * It is initially unscaled.
217 */
218 int Compute_J(slv_system_t sys)
219 {
220 int32 row;
221 var_filter_t vfilter;
222 linsolqr_system_t lqr_sys;
223 mtx_matrix_t mtx;
224 struct rel_relation **rlist;
225 int nrows;
226 real64 resid;
227
228 /*
229 * Get the linear system from the solve system.
230 * Get the matrix from the linear system.
231 * Get the relations list from the solve system.
232 */
233 lqr_sys = slv_get_linsolqr_sys(sys);
234 mtx = linsolqr_get_matrix(lqr_sys);
235 rlist = slv_get_solvers_rel_list(sys);
236 nrows = slv_get_num_solvers_rels(sys);
237
238 calc_ok = TRUE;
239 vfilter.matchbits = (VAR_SVAR | VAR_ACTIVE) ;
240 vfilter.matchvalue = vfilter.matchbits ;
241
242 /*
243 * Clear the entire matrix and then compute
244 * the gradients at the current point.
245 */
246 mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
247 for(row=0; row<nrows; row++) {
248 struct rel_relation *rel;
249 rel = rlist[mtx_row_to_org(mtx,row)];
250 asc_assert(rel!=NULL);
251 (void)relman_diffs(rel,&vfilter,mtx,&resid,1);
252 }
253 linsolqr_matrix_was_changed(lqr_sys);
254
255 return(!calc_ok);
256 }
257
258 /**
259 @NOTE there is another (probably better?) version of this in models/sensitivity
260 that uses sparse matrices. This one got pulled out of some files in tcltk dir
261 and is used by the LSODE integrator.
262
263 * Note a rhs would have been previously added
264 * to keep the system happy.
265 * Now handles both linsol and linsolqr as needed.
266 * Assumes region to be factored is in upper left corner.
267 * Region is reordered using the last reorder method used in
268 * the case of linsolqr, so all linsolqr methods are available
269 * to this routine.
270 */
271 int LUFactorJacobian(slv_system_t sys){
272 linsol_system_t lin_sys=NULL;
273 linsolqr_system_t lqr_sys=NULL;
274 mtx_region_t region;
275 int size,nvars,nrels;
276 #if DOTIME
277 double time1;
278 #endif
279
280 #if DOTIME
281 time1 = tm_cpu_time();
282 #endif
283
284 nvars = NumberFreeVars(sys);
285 nrels = NumberIncludedRels(sys);
286 size = MAX(nvars,nrels); /* get size of matrix */
287 mtx_region(&region,0,size-1,0,size-1); /* set the region */
288 lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
289 if (lin_sys!=NULL) {
290 linsol_matrix_was_changed(lin_sys);
291 linsol_reorder(lin_sys,&region); /* This was entire_MATRIX */
292 linsol_invert(lin_sys,&region); /* invert the given matrix over
293 * the given region */
294 } else {
295 enum factor_method fm;
296 int oldtiming;
297 lqr_sys = slv_get_linsolqr_sys(sys);
298 /* WE are ASSUMING that the system has been qr_preped before now! */
299 linsolqr_matrix_was_changed(lqr_sys);
300 linsolqr_reorder(lqr_sys,&region,natural); /* should retain original */
301 fm = linsolqr_fmethod(lqr_sys);
302 if (fm == unknown_f) {
303 fm = ranki_kw2; /* make sure somebody set it */
304 }
305 oldtiming = g_linsolqr_timing;
306 g_linsolqr_timing = 0;
307 linsolqr_factor(lqr_sys,fm);
308 g_linsolqr_timing = oldtiming;
309 }
310
311 #if DOTIME
312 time1 = tm_cpu_time() - time1;
313 CONSOLE_DEBUG("Time taken for LUFactorJacobian = %g\n",time1);
314 #endif
315 return 0;
316 }
317
318
319 /**
320 At this point the rhs would have already been added.
321
322 Extended to handle either factorization code:
323 linsol or linsolqr.
324
325 This routine is part of the 'temporary solution' for derivatives in lsode.c
326 */
327 int Compute_dy_dx_smart(slv_system_t sys,
328 real64 *rhs,
329 real64 **dy_dx,
330 int *inputs, int ninputs,
331 int *outputs, int noutputs)
332 {
333 linsol_system_t lin_sys=NULL;
334 linsolqr_system_t lqr_sys=NULL;
335 mtx_matrix_t mtx;
336 int col,current_col;
337 int row;
338 int capacity;
339 real64 *solution = NULL;
340 int i,j;
341 #if DOTIME
342 double time1;
343 #endif
344
345 #if DOTIME
346 time1 = tm_cpu_time();
347 #endif
348
349 lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
350 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
351 mtx = slv_get_sys_mtx(sys); /* get the matrix */
352
353 capacity = mtx_capacity(mtx);
354 solution = ASC_NEW_ARRAY_CLEAR(real64,capacity);
355
356 /*
357 * The array inputs is a list of original indexes, of the variables
358 * that we are trying to obtain the sensitivity to. We have to
359 * get the *current* column from the matrix based on those indices.
360 * Hence we use mtx_org_to_col. This little manipulation is not
361 * necessary for the computed solution as the solve routine returns
362 * the results in the *original* order rather than the *current* order.
363 */
364 if (lin_sys!=NULL) {
365 for (j=0;j<ninputs;j++) {
366 col = inputs[j];
367 current_col = mtx_org_to_col(mtx,col);
368 mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS); /* get the column */
369
370 linsol_rhs_was_changed(lin_sys,rhs);
371 linsol_solve(lin_sys,rhs); /* solve */
372 linsol_copy_solution(lin_sys,rhs,solution); /* get the solution */
373
374 for (i=0;i<noutputs;i++) { /* copy the solution to dy_dx */
375 row = outputs[i];
376 dy_dx[i][j] = -1.0*solution[row]; /* the -1 comes from the lin alg */
377 }
378 /*
379 * zero the vector using the incidence pattern of our column.
380 * This is faster than the naiive approach of zeroing
381 * everything especially if the vector is large.
382 */
383 mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
384 }
385 } else {
386 for (j=0;j<ninputs;j++) {
387 col = inputs[j];
388 current_col = mtx_org_to_col(mtx,col);
389 mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS);
390
391 linsolqr_rhs_was_changed(lqr_sys,rhs);
392 linsolqr_solve(lqr_sys,rhs);
393 linsolqr_copy_solution(lqr_sys,rhs,solution);
394
395 for (i=0;i<noutputs;i++) {
396 row = outputs[i];
397 dy_dx[i][j] = -1.0*solution[row];
398 }
399 mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
400 }
401 }
402 if (solution) {
403 ascfree((char *)solution);
404 }
405
406 #if DOTIME
407 time1 = tm_cpu_time() - time1;
408 CONSOLE_DEBUG("Time for Compute_dy_dx_smart = %g\n",time1);
409 #endif
410 return 0;
411 }
412
413 #undef DEBUG

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