/[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 1162 - (show annotations) (download) (as text)
Wed Jan 17 04:34:59 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 20158 byte(s)
Fixed sensitivity module. It now runs as mostly external code, with just a few relics required by LSODE still in packages/sensitivity.c.
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 <string.h>
29 #include <math.h>
30 #include <utilities/error.h>
31 #include <utilities/ascMalloc.h>
32 #include <utilities/ascPanic.h>
33 #include <compiler/packages.h>
34 #include <compiler/fractions.h>
35 #include <compiler/dimen.h>
36 #include <compiler/atomvalue.h>
37 #include <compiler/instquery.h>
38 #include <general/mathmacros.h>
39
40 #define DEBUG 1
41
42 #if 0
43 static real64 *zero_vector(real64 *vec, int size)
44 {
45 int c;
46 for (c=0;c<size;c++) {
47 vec[c] = 0.0;
48 }
49 return vec;
50 }
51 #endif
52
53 /*----------------------------------------------------
54 UTILITY ROUTINES
55 */
56
57 /**
58 Allocate memory for a matrix
59 @param nrows Number of rows
60 @param ncols Number of colums
61 @return Pointer to the allocated matrix memory location
62 */
63 real64 **make_matrix(int nrows, int ncols){
64 real64 **result;
65 int i;
66 result = (real64 **)calloc(nrows,sizeof(real64*));
67 for (i=0;i<nrows;i++) {
68 result[i] = (real64 *)calloc(ncols,sizeof(real64));
69 }
70 return result;
71 }
72
73 /**
74 Free a matrix from memory
75 @param matrix Memory location for the matrix
76 @param nrows Number of rows in the matrix
77 */
78 void free_matrix(real64 **matrix, int nrows){
79 int i;
80 if (!matrix)
81 return;
82 for (i=0;i<nrows;i++) {
83 if (matrix[i]) {
84 free(matrix[i]);
85 matrix[i] = NULL;
86 }
87 }
88 free(matrix);
89 }
90
91 /**
92 Fetch an element from a branch of an arglist...?
93 */
94 struct Instance *FetchElement(struct gl_list_t *arglist,
95 unsigned long whichbranch,
96 unsigned long whichelement){
97 struct gl_list_t *branch;
98 struct Instance *element;
99
100 if (!arglist) return NULL;
101 branch = (struct gl_list_t *)gl_fetch(arglist,whichbranch);
102 element = (struct Instance *)gl_fetch(branch,whichelement);
103 return element;
104 }
105
106 #if 0
107 static int ReSolve(slv_system_t sys)
108 {
109 if (!sys)
110 return 1;
111 slv_solve(sys);
112 return 0;
113 }
114 #endif
115
116 /**
117 Build then presolve the solve an instance...
118 */
119 int DoSolve(struct Instance *inst){
120 slv_system_t sys;
121
122 sys = system_build(inst);
123 if (!sys) {
124 ERROR_REPORTER_HERE(ASC_PROG_ERR,
125 "Something radically wrong in creating solver.");
126 return 1;
127 }
128 (void)slv_select_solver(sys,0);
129 slv_presolve(sys);
130 slv_solve(sys);
131 system_destroy(sys);
132 return 0;
133 }
134
135 /**
136 Finite difference...
137 */
138 int finite_difference(struct gl_list_t *arglist){
139 struct Instance *model_inst, *xinst, *inst;
140 slv_system_t sys;
141 int ninputs,noutputs;
142 int i,j,offset;
143 real64 **partials;
144 real64 *y_old, *y_new;
145 real64 x;
146 real64 interval = 1e-6;
147 int result=0;
148
149 model_inst = FetchElement(arglist,1,1);
150 sys = system_build(model_inst);
151 if (!sys) {
152 FPRINTF(stdout,"Something radically wrong in creating solver\n");
153 return 1;
154 }
155 (void)slv_select_solver(sys,0);
156 slv_presolve(sys);
157 slv_solve(sys);
158
159 /* Make the necessary vectors */
160
161 ninputs = (int)gl_length((struct gl_list_t *)gl_fetch(arglist,2));
162 noutputs = (int)gl_length((struct gl_list_t *)gl_fetch(arglist,3));
163 y_old = (real64 *)calloc(noutputs,sizeof(real64));
164 y_new = (real64 *)calloc(noutputs,sizeof(real64));
165 partials = make_matrix(noutputs,ninputs);
166 for (i=0;i<noutputs;i++) { /* get old yvalues */
167 inst = FetchElement(arglist,3,i+1);
168 y_old[i] = RealAtomValue(inst);
169 }
170 for (j=0;j<ninputs;j++) {
171 xinst = FetchElement(arglist,2,j+1);
172 x = RealAtomValue(xinst);
173 SetRealAtomValue(xinst,x+interval,(unsigned)0); /* perturb system */
174 slv_presolve(sys);
175 slv_solve(sys);
176
177 for (i=0;i<noutputs;i++) { /* get new yvalues */
178 inst = FetchElement(arglist,3,i+1);
179 y_new[i] = RealAtomValue(inst);
180 partials[i][j] = -1.0 * (y_old[i] - y_new[i])/interval;
181 PRINTF("y_old = %20.12g y_new = %20.12g\n", y_old[i],y_new[i]);
182 }
183 SetRealAtomValue(xinst,x,(unsigned)0); /* unperturb system */
184 }
185 offset = 0;
186 for (i=0;i<noutputs;i++) {
187 for (j=0;j<ninputs;j++) {
188 inst = FetchElement(arglist,4,offset+j+1);
189 SetRealAtomValue(inst,partials[i][j],(unsigned)0);
190 PRINTF("%12.6f %s",partials[i][j], (j==(ninputs-1)) ? "\n" : " ");
191 }
192 offset += ninputs;
193 }
194 /* error: */
195 free(y_old);
196 free(y_new);
197 free_matrix(partials,noutputs);
198 system_destroy(sys);
199 return result;
200 }
201
202
203
204 /**
205 Looks like it returns the number of relations in a solver's system
206
207 @param sys The system in question
208 @return int Number of relations in the system
209 @see slv_count_solvers_rels
210 */
211 int NumberRels(slv_system_t sys){
212 static int nrels = -1;
213 rel_filter_t rfilter;
214 int tmp;
215
216 if (!sys) { /* a NULL system may be used to reinitialise the count */
217 nrels = -1;
218 return -1;
219 }
220 if (nrels < 0) {
221 /*get used (included) relation count -- equalities only !*/
222 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
223 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
224 tmp = slv_count_solvers_rels(sys,&rfilter);
225 nrels = tmp;
226 return tmp;
227 }
228 else return nrels;
229 }
230
231 int NumberIncludedRels(slv_system_t sys){
232 static int nrels = -1;
233 rel_filter_t rfilter;
234 int tmp;
235
236 if (!sys) { /* a NULL system may be used to reinitialise the count */
237 nrels = -1;
238 return -1;
239 }
240 if (nrels < 0) {
241 /*get used (included) relation count -- equalities only !*/
242 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
243 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
244 tmp = slv_count_solvers_rels(sys,&rfilter);
245 nrels = tmp;
246 return tmp;
247 } else {
248 return nrels;
249 }
250 }
251
252 /**
253 Looks like it returns the number of free variables in a solver's system
254
255 @param sys The system in question
256 @return the number of free variables
257
258 @see slv_count_solvers_vars
259 */
260 int NumberFreeVars(slv_system_t sys){
261 static int nvars = -1;
262 var_filter_t vfilter;
263 int tmp;
264
265 if (!sys) {
266 /* no system: return error */
267 nvars = -1;
268 return -1;
269 }
270
271 if (nvars >=0){
272 /* return stored value */
273 return nvars;
274 }
275
276 /* get used (free and incident) variable count */
277 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
278 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
279 tmp = slv_count_solvers_vars(sys,&vfilter);
280 nvars = tmp;
281 return tmp;
282 }
283
284 /*
285 * The following bit of code does the computation of the matrix
286 * dz/dx. It accepts a slv_system_t and a list of 'input' variables.
287 * The matrix df/dx now exists and sits to the right of the main
288 * square region of the Jacobian. We will do the following in turn
289 * for each variable in this list:
290 *
291 * 1) Get the variable index, and use it to extract the required column
292 * from the main gradient matrix, this will be stored in a temporary
293 * vector.
294 * 2) We will then clear the column in the original matrix, as we want to
295 * store the caluclated results back in place.
296 * 3) Add the vector extracted in 1) as rhs to the main matrix.
297 * 4) Call linsol solve on this rhs to yield a solution which represents
298 * a column of dx/dx.
299 * 6) Store this vector back in the location cleared out in 2).
300 * 7) Goto 1.
301 *
302 * At the end of this procedure we would have calculated all the columns of
303 * dz/dx, and left them back in the main matrix.
304 */
305
306
307 /*
308 * At this point we should have an empty jacobian. We now
309 * need to call relman_diff over the *entire* matrix.
310 * fixed and free.
311 *
312 * Calculates the entire jacobian.
313 * It is initially unscaled.
314 */
315 int Compute_J(slv_system_t sys)
316 {
317 int32 row;
318 var_filter_t vfilter;
319 linsolqr_system_t lqr_sys;
320 mtx_matrix_t mtx;
321 struct rel_relation **rlist;
322 int nrows;
323 real64 resid;
324
325 /*
326 * Get the linear system from the solve system.
327 * Get the matrix from the linear system.
328 * Get the relations list from the solve system.
329 */
330 lqr_sys = slv_get_linsolqr_sys(sys);
331 mtx = linsolqr_get_matrix(lqr_sys);
332 rlist = slv_get_solvers_rel_list(sys);
333 nrows = slv_get_num_solvers_rels(sys);
334
335 calc_ok = TRUE;
336 vfilter.matchbits = (VAR_SVAR | VAR_ACTIVE) ;
337 vfilter.matchvalue = vfilter.matchbits ;
338
339 /*
340 * Clear the entire matrix and then compute
341 * the gradients at the current point.
342 */
343 mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
344 for(row=0; row<nrows; row++) {
345 struct rel_relation *rel;
346 rel = rlist[mtx_row_to_org(mtx,row)];
347 asc_assert(rel!=NULL);
348 (void)relman_diffs(rel,&vfilter,mtx,&resid,1);
349 }
350 linsolqr_matrix_was_changed(lqr_sys);
351
352 return(!calc_ok);
353 }
354
355
356 #if 0
357 /**
358 * At this point we should have an empty jacobian. We now
359 * need to call relman_diff over the *entire* matrix.
360 * Calculates the entire jacobian. It is initially unscaled.
361 *
362 * Note: this assumes the sys given is using one of the ascend solvers
363 * with either linsol or linsolqr. Both are now allowed. baa 7/95
364 */
365 #define SAFE_FIX_ME 0
366 static int Compute_J_OLD(slv_system_t sys)
367 {
368 int32 row;
369 var_filter_t vfilter;
370 linsol_system_t lin_sys = NULL;
371 linsolqr_system_t lqr_sys = NULL;
372 mtx_matrix_t mtx;
373 struct rel_relation **rlist;
374 int nrows;
375 int qr=0;
376 #\if DOTIME
377 double time1;
378 #\endif
379
380 #\if DOTIME
381 time1 = tm_cpu_time();
382 #\endif
383 /*
384 * Get the linear system from the solve system.
385 * Get the matrix from the linear system.
386 * Get the relations list from the solve system.
387 */
388 lin_sys = slv_get_linsol_sys(sys);
389 if (lin_sys==NULL) {
390 qr=1;
391 lqr_sys=slv_get_linsolqr_sys(sys);
392 }
393 mtx = slv_get_sys_mtx(sys);
394 if (mtx==NULL || (lin_sys==NULL && lqr_sys==NULL)) {
395 FPRINTF(stderr,"Compute_dy_dx: error, found NULL.\n");
396 return 1;
397 }
398 rlist = slv_get_solvers_rel_list(sys);
399 nrows = NumberIncludedRels(sys);
400
401 calc_ok = TRUE;
402 vfilter.matchbits = VAR_SVAR;
403 vfilter.matchvalue = VAR_SVAR;
404 /*
405 * Clear the entire matrix and then compute
406 * the gradients at the current point.
407 */
408 mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
409 for(row=0; row<nrows; row++) {
410 struct rel_relation *rel;
411 /* added */
412 double resid;
413
414 rel = rlist[mtx_row_to_org(mtx,row)];
415 (void)relman_diffs(rel,&vfilter,mtx,&resid,SAFE_FIX_ME);
416
417 /* added */
418 rel_set_residual(rel,resid);
419
420 }
421 if (qr) {
422 linsolqr_matrix_was_changed(lqr_sys);
423 } else {
424 linsol_matrix_was_changed(lin_sys);
425 }
426 #\if DOTIME
427 time1 = tm_cpu_time() - time1;
428 FPRINTF(stderr,"Time taken for Compute_J = %g\n",time1);
429 #\endif
430 return(!calc_ok);
431 }
432 #endif
433
434 /*
435 LU Factor Jacobian
436
437 @note The RHS will have been will already have been added before
438 calling this function.
439
440 THERE SEEM TO BE TWO VERSIONS OF THIS... WHAT'S THE STORY?
441 */
442 int LUFactorJacobian1(slv_system_t sys,int rank){
443 linsolqr_system_t lqr_sys;
444 mtx_region_t region;
445 enum factor_method fm;
446
447 mtx_region(&region,0,rank-1,0,rank-1); /* set the region */
448 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
449
450 linsolqr_matrix_was_changed(lqr_sys);
451 linsolqr_reorder(lqr_sys,&region,natural);
452
453 fm = linsolqr_fmethod(lqr_sys);
454 if (fm == unknown_f) fm = ranki_kw2; /* make sure somebody set it */
455 linsolqr_factor(lqr_sys,fm);
456
457 return 0;
458 }
459
460 /**
461 This is the other version of this, pulled in from Tcl/Tk files.
462
463 * Note a rhs would have been previously added
464 * to keep the system happy.
465 * Now handles both linsol and linsolqr as needed.
466 * Assumes region to be factored is in upper left corner.
467 * Region is reordered using the last reorder method used in
468 * the case of linsolqr, so all linsolqr methods are available
469 * to this routine.
470 */
471 int LUFactorJacobian(slv_system_t sys){
472 linsol_system_t lin_sys=NULL;
473 linsolqr_system_t lqr_sys=NULL;
474 mtx_region_t region;
475 int size,nvars,nrels;
476 #if DOTIME
477 double time1;
478 #endif
479
480 #if DOTIME
481 time1 = tm_cpu_time();
482 #endif
483
484 nvars = NumberFreeVars(sys);
485 nrels = NumberIncludedRels(sys);
486 size = MAX(nvars,nrels); /* get size of matrix */
487 mtx_region(&region,0,size-1,0,size-1); /* set the region */
488 lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
489 if (lin_sys!=NULL) {
490 linsol_matrix_was_changed(lin_sys);
491 linsol_reorder(lin_sys,&region); /* This was entire_MATRIX */
492 linsol_invert(lin_sys,&region); /* invert the given matrix over
493 * the given region */
494 } else {
495 enum factor_method fm;
496 int oldtiming;
497 lqr_sys = slv_get_linsolqr_sys(sys);
498 /* WE are ASSUMING that the system has been qr_preped before now! */
499 linsolqr_matrix_was_changed(lqr_sys);
500 linsolqr_reorder(lqr_sys,&region,natural); /* should retain original */
501 fm = linsolqr_fmethod(lqr_sys);
502 if (fm == unknown_f) {
503 fm = ranki_kw2; /* make sure somebody set it */
504 }
505 oldtiming = g_linsolqr_timing;
506 g_linsolqr_timing = 0;
507 linsolqr_factor(lqr_sys,fm);
508 g_linsolqr_timing = oldtiming;
509 }
510
511 #if DOTIME
512 time1 = tm_cpu_time() - time1;
513 FPRINTF(stderr,"Time taken for LUFactorJacobian = %g\n",time1);
514 #endif
515 return 0;
516 }
517
518
519 /**
520 At this point the rhs would have already been added.
521
522 Extended to handle either factorization code:
523 linsol or linsolqr.
524
525 This routine is part of the 'temporary solution' for derivatives in lsode.c
526 */
527 int Compute_dy_dx_smart(slv_system_t sys,
528 real64 *rhs,
529 real64 **dy_dx,
530 int *inputs, int ninputs,
531 int *outputs, int noutputs)
532 {
533 linsol_system_t lin_sys=NULL;
534 linsolqr_system_t lqr_sys=NULL;
535 mtx_matrix_t mtx;
536 int col,current_col;
537 int row;
538 int capacity;
539 real64 *solution = NULL;
540 int i,j;
541 #if DOTIME
542 double time1;
543 #endif
544
545 #if DOTIME
546 time1 = tm_cpu_time();
547 #endif
548
549 lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
550 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
551 mtx = slv_get_sys_mtx(sys); /* get the matrix */
552
553 capacity = mtx_capacity(mtx);
554 solution = ASC_NEW_ARRAY_CLEAR(real64,capacity);
555
556 /*
557 * The array inputs is a list of original indexes, of the variables
558 * that we are trying to obtain the sensitivity to. We have to
559 * get the *current* column from the matrix based on those indices.
560 * Hence we use mtx_org_to_col. This little manipulation is not
561 * necessary for the computed solution as the solve routine returns
562 * the results in the *original* order rather than the *current* order.
563 */
564 if (lin_sys!=NULL) {
565 for (j=0;j<ninputs;j++) {
566 col = inputs[j];
567 current_col = mtx_org_to_col(mtx,col);
568 mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS); /* get the column */
569
570 linsol_rhs_was_changed(lin_sys,rhs);
571 linsol_solve(lin_sys,rhs); /* solve */
572 linsol_copy_solution(lin_sys,rhs,solution); /* get the solution */
573
574 for (i=0;i<noutputs;i++) { /* copy the solution to dy_dx */
575 row = outputs[i];
576 dy_dx[i][j] = -1.0*solution[row]; /* the -1 comes from the lin alg */
577 }
578 /*
579 * zero the vector using the incidence pattern of our column.
580 * This is faster than the naiive approach of zeroing
581 * everything especially if the vector is large.
582 */
583 mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
584 }
585 } else {
586 for (j=0;j<ninputs;j++) {
587 col = inputs[j];
588 current_col = mtx_org_to_col(mtx,col);
589 mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS);
590
591 linsolqr_rhs_was_changed(lqr_sys,rhs);
592 linsolqr_solve(lqr_sys,rhs);
593 linsolqr_copy_solution(lqr_sys,rhs,solution);
594
595 for (i=0;i<noutputs;i++) {
596 row = outputs[i];
597 dy_dx[i][j] = -1.0*solution[row];
598 }
599 mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
600 }
601 }
602 if (solution) {
603 ascfree((char *)solution);
604 }
605
606 #if DOTIME
607 time1 = tm_cpu_time() - time1;
608 FPRINTF(stderr,"Time for Compute_dy_dx_smart = %g\n",time1);
609 #endif
610 return 0;
611 }
612
613 #if 0
614 static int ComputeInverse(slv_system_t sys,
615 real64 *rhs)
616 {
617 linsolqr_system_t lqr_sys;
618 mtx_matrix_t mtx;
619 int capacity,order;
620 real64 *solution = NULL;
621 int j,k;
622
623 lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
624 mtx = linsolqr_get_matrix(lqr_sys); /* get the matrix */
625
626 capacity = mtx_capacity(mtx);
627 zero_vector(rhs,capacity); /* zero the rhs */
628 solution = ASC_NEW_ARRAY_CLEAR(real64,capacity);
629
630 order = mtx_order(mtx);
631 for (j=0;j<order;j++) {
632 rhs[j] = 1.0;
633 linsolqr_rhs_was_changed(lqr_sys,rhs);
634 linsolqr_solve(lqr_sys,rhs); /* solve */
635 linsolqr_copy_solution(lqr_sys,rhs,solution); /* get the solution */
636
637 FPRINTF(stderr,"This is the rhs and solution for rhs #%d\n",j);
638 for (k=0;k<order;k++) {
639 FPRINTF(stderr,"%12.8g %12.8g\n",rhs[k],solution[k]);
640 }
641 rhs[j] = 0.0;
642 }
643 if (solution) ascfree((char *)solution);
644 return 0;
645 }
646 #endif
647
648 /**
649 Do Data Analysis
650 */
651 int DoDataAnalysis(struct var_variable **inputs,
652 struct var_variable **outputs,
653 int ninputs, int noutputs,
654 real64 **dy_dx)
655 {
656 FILE *fp;
657 double *norm_2, *norm_1;
658 double input_nominal,maxvalue,sum;
659 int i,j;
660
661 norm_1 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
662 norm_2 = ASC_NEW_ARRAY_CLEAR(double,ninputs);
663
664 fp = fopen("sensitivity.out","w+");
665 if (!fp) return 0;
666
667 /*
668 * calculate the 1 and 2 norms; cache them away so that we
669 * can pretty print them. Style is everything !.
670 *
671 */
672 for (j=0;j<ninputs;j++) {
673 input_nominal = var_nominal(inputs[j]);
674 maxvalue = sum = 0;
675 for (i=0;i<noutputs;i++) {
676 dy_dx[i][j] *= input_nominal/var_nominal(outputs[i]);
677 maxvalue = MAX(fabs(dy_dx[i][j]),maxvalue);
678 sum += dy_dx[i][j]*dy_dx[i][j];
679 }
680 norm_1[j] = maxvalue;
681 norm_2[j] = sum;
682 }
683
684 for (j=0;j<ninputs;j++) { /* print the var_index */
685 FPRINTF(fp,"%8d ",var_mindex(inputs[j]));
686 }
687 FPRINTF(fp,"\n");
688
689 for (j=0;j<ninputs;j++) { /* print the 1 norms */
690 FPRINTF(fp,"%-#18.8f ",norm_1[j]);
691 }
692 FPRINTF(fp,"\n");
693
694 for (j=0;j<ninputs;j++) { /* print the 2 norms */
695 FPRINTF(fp,"%-#18.8f ",norm_2[j]);
696 }
697 FPRINTF(fp,"\n\n");
698 ascfree((char *)norm_1);
699 ascfree((char *)norm_2);
700
701 for (i=0;i<noutputs;i++) { /* print the scaled data */
702 for (j=0;j<ninputs;j++) {
703 FPRINTF(fp,"%-#18.8f %-4d",dy_dx[i][j],i);
704 }
705 if (var_fixed(outputs[i]))
706 FPRINTF(fp," **fixed*** \n");
707 else
708 PUTC('\n',fp);
709 }
710 FPRINTF(fp,"\n");
711 fclose(fp);
712 return 0;
713 }
714
715 #if 0
716 static int DoProject_X(struct var_variable **old_inputs,
717 struct var_variable **new_inputs, /* new values of u */
718 double step_length,
719 struct var_variable **outputs,
720 int ninputs, int noutputs,
721 real64 **dy_dx)
722 {
723 struct var_variable *var;
724 real64 old_y, new_y, tmp;
725 real64 *delta_x;
726 int i,j;
727
728 delta_x = ASC_NEW_ARRAY_CLEAR(real64,ninputs);
729
730 for (j=0;j<ninputs;j++) {
731 delta_x[j] = var_value(new_inputs[j]) - var_value(old_inputs[j]);
732 /* delta_x[j] = RealAtomValue(new_inputs[j]) - RealAtomValue(old_inputs[j]); */
733 }
734
735 for (i=0;i<noutputs;i++) {
736 var = outputs[i];
737 if (var_fixed(var) || !var_active(var)) /* project only the free vars */
738 continue;
739 tmp = 0.0;
740 for (j=0;j<ninputs;j++) {
741 tmp += (dy_dx[i][j] * delta_x[j]);
742 }
743 /* old_y = RealAtomValue(var); */
744 old_y = var_value(var);
745 new_y = old_y + step_length*tmp;
746 /* SetRealAtomValue(var,new_y,(unsigned)0); */
747 var_set_value(var,new_y);
748 # if DEBUG
749 FPRINTF(stderr,"Old_y = %12.8g; Nex_y = %12.8g\n",old_y,new_y);
750 # endif
751 }
752 ascfree((char *)delta_x);
753 return 0;
754 }
755 #endif
756
757
758
759
760
761
762 #undef DEBUG

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