/[ascend]/trunk/base/generic/solver/integrator.c
ViewVC logotype

Contents of /trunk/base/generic/solver/integrator.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 946 - (show annotations) (download) (as text)
Sat Nov 25 15:28:56 2006 UTC (13 years, 9 months ago) by johnpye
File MIME type: text/x-csrc
File size: 42984 byte(s)
Found bug in integrator_analyse_dae... ongoing
1 /* ASCEND modelling environment
2 Copyright (C) 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 *//**
19 @file
20 Integrator API for ASCEND, for solving systems of ODEs and/or DAEs.
21 *//*
22 by John Pye, May 2006
23 based on parts of Integrators.c in the Tcl/Tk interface directory, heavily
24 modified to provide a non-GUI-specific API and modularised for multiple
25 integration engines.
26 */
27 #include <time.h>
28 #include <string.h>
29 #include "integrator.h"
30 #include "lsode.h"
31 #include "ida.h"
32 #include "slv_common.h"
33 #include "slv_client.h"
34 #include <utilities/ascPanic.h>
35
36 #include "samplelist.h"
37
38 /**
39 Define as TRUE to enable debug message printing.
40 @TODO this needs to go away.
41 */
42 #define INTEG_DEBUG TRUE
43
44 /**
45 Print a debug message with value if INTEG_DEBUG is TRUE.
46 */
47 #define print_debug(msg, value) \
48 if(INTEG_DEBUG){ CONSOLE_DEBUG(msg, value); }
49
50 /**
51 Print a debug message string if INTEG_DEBUG is TRUE.
52 */
53 #define print_debugstring(msg) \
54 if(INTEG_DEBUG){ CONSOLE_DEBUG(msg); }
55
56 /*------------------------------------------------------------------------------
57 The following names are of solver_var children or attributes
58 * we support (at least temporarily) to determine who is a state and
59 * who matching derivative.
60 * These should be supported directly in a future solveratominst.
61 */
62
63 static symchar *g_symbols[3];
64
65 #define STATEFLAG g_symbols[0]
66 /*
67 Integer child. 0= algebraic, 1 = state, 2 = derivative, 3 = 2nd deriv etc
68 independent variable is -1.
69 */
70 #define INTEG_OTHER_VAR -1L
71 #define INTEG_ALGEBRAIC_VAR 0L
72 #define INTEG_STATE_VAR 1L
73
74 #define STATEINDEX g_symbols[1]
75 /* Integer child. all variables with the same STATEINDEX value are taken to
76 * be derivatives of the same state variable. We really need a compiler
77 * that maintains this info by backpointers, but oh well until then.
78 */
79 #define OBSINDEX g_symbols[2]
80 /* Integer child. All variables with OBSINDEX !=0 will be recorded in
81 * the blsode output file. Tis someone else's job to grok this output.
82 */
83
84
85 /** Temporary catcher of dynamic variable and observation variable data */
86 struct Integ_var_t {
87 long index;
88 long type;
89 struct Integ_var_t *derivative;
90 struct Integ_var_t *derivative_of;
91 struct var_variable *i;
92 int varindx; /**< index into slv_get_master_vars_list, or -1 if not there */
93 int isstate;
94 };
95
96 /*------------------------------------------------------------------------------
97 forward declarations
98 */
99
100 /* abstractions of setup/teardown procedures for the specific solvers */
101 void integrator_create_engine(IntegratorSystem *sys);
102 void integrator_free_engine(IntegratorSystem *sys);
103
104 static int integrator_analyse_ode(IntegratorSystem *sys);
105 static int integrator_analyse_dae(IntegratorSystem *sys);
106
107 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
108 static void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
109 IntegratorVarVisitorFn integrator_ode_classify_var;
110 IntegratorVarVisitorFn integrator_dae_classify_var;
111 IntegratorVarVisitorFn integrator_classify_indep_var;
112
113 static int integrator_sort_obs_vars(IntegratorSystem *sys);
114 static void integrator_print_var_stats(IntegratorSystem *sys);
115 static int integrator_check_indep_var(IntegratorSystem *sys);
116
117 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2);
118 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2);
119 static void Integ_SetObsId(struct var_variable *v, long oindex);
120
121 static long DynamicVarInfo(struct var_variable *v,long *vindex);
122 static struct var_variable *ObservationVar(struct var_variable *v, long *oindex);
123 static void IntegInitSymbols(void);
124
125 /*------------------------------------------------------------------------------
126 INSTANTIATION AND DESTRUCTION
127 */
128
129 /**
130 Create a new IntegratorSystem and assign a slv_system_t to it.
131 */
132 IntegratorSystem *integrator_new(slv_system_t slvsys, struct Instance *inst){
133 IntegratorSystem *sys;
134
135 if (slvsys == NULL) {
136 ERROR_REPORTER_HERE(ASC_PROG_ERR,"slvsys is NULL!");
137 return NULL;
138 }
139
140 sys = ASC_NEW_CLEAR(IntegratorSystem);
141 sys->system = slvsys;
142 sys->instance = inst;
143
144 sys->states = NULL; sys->derivs = NULL;
145 sys->dynvars = NULL; sys->obslist = NULL; sys->indepvars = NULL;
146
147 sys->y_id = NULL;
148 sys->obs_id = NULL;
149 sys->y = NULL;
150 sys->ydot = NULL;
151 sys->obs = NULL;
152 return sys;
153 }
154
155 /**
156 Carefully trash any data in the IntegratorSystem that we own,
157 then destroy the IntegratorSystem struct.
158
159 Note that the integrator doesn't own the samplelist.
160
161 @param sys will be destroyed and set to NULL.
162 */
163 void integrator_free(IntegratorSystem *sys){
164 if(sys==NULL)return;
165
166 integrator_free_engine(sys);
167
168 if(sys->states != NULL)gl_destroy(sys->states);
169 if(sys->derivs != NULL)gl_destroy(sys->derivs);
170
171 if(sys->dynvars != NULL)gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
172 if(sys->obslist != NULL)gl_free_and_destroy(sys->obslist); /* and obslist */
173 if (sys->indepvars != NULL)gl_free_and_destroy(sys->indepvars); /* and indepvars */
174
175 /* if(sys->y_id != NULL)ASC_FREE(sys->y_id); */
176 if(sys->obs_id != NULL)ASC_FREE(sys->obs_id);
177
178 if(sys->y != NULL && !sys->ycount)ASC_FREE(sys->y);
179 if(sys->ydot != NULL && !sys->ydotcount)ASC_FREE(sys->ydot);
180 if(sys->obs != NULL && !sys->obscount)ASC_FREE(sys->obs);
181
182 slv_destroy_parms(&(sys->params));
183
184 ASC_FREE(sys);
185 CONSOLE_DEBUG("Destroyed IntegratorSystem");
186 sys=NULL;
187 }
188
189 /**
190 Utility function to retreive pointers to the symbols we'll be looking for
191 in the instance hierarchy.
192 */
193 static void IntegInitSymbols(void){
194 STATEFLAG = AddSymbol("ode_type");
195 STATEINDEX = AddSymbol("ode_id");
196 OBSINDEX = AddSymbol("obs_id");
197 }
198
199 /*------------------------------------------------------------------------------
200 INTEGRATOR ENGINE
201 */
202
203 /* return 0 on success */
204 int integrator_set_engine(IntegratorSystem *sys, IntegratorEngine engine){
205
206 /* verify integrator type ok. always passes for nonNULL inst. */
207 if(engine==INTEG_UNKNOWN){
208 ERROR_REPORTER_NOLINE(ASC_USER_ERROR
209 ,"Integrator has not been specified (or is unknown)."
210 );
211 return 1;
212 }else{
213 /** @TODO other engine-specific tests */
214 }
215
216 if(engine==sys->engine){
217 return 0;
218 }
219 if(sys->engine!=INTEG_UNKNOWN){
220 integrator_free_engine(sys);
221 }
222 sys->engine = engine;
223 integrator_create_engine(sys);
224
225 return 0;
226 }
227
228 IntegratorEngine integrator_get_engine(const IntegratorSystem *sys){
229 return sys->engine;
230 }
231
232 /**
233 Free any engine-specific data that was required for the solution of
234 this system. Note that this data is pointed to by sys->enginedata.
235 */
236 void integrator_free_engine(IntegratorSystem *sys){
237 switch(sys->engine){
238 case INTEG_LSODE: integrator_lsode_free(sys->enginedata); break;
239 #ifdef ASC_WITH_IDA
240 case INTEG_IDA: integrator_ida_free(sys->enginedata); break;
241 #endif
242 default: break;
243 }
244 sys->enginedata=NULL;
245 }
246
247 /**
248 Create enginedata memory if required for this solver. This doesn't include
249 allocating computation space, since we assume that this stage all we know
250 is that we want to use a specified integrator engine, not the full details
251 of the problem at hand. Allocating space inside enginedata should be done
252 during the solve stage (and freed inside integrator_free_engine)
253 */
254 void integrator_create_engine(IntegratorSystem *sys){
255 if(sys->enginedata!=NULL)return;
256 switch(sys->engine){
257 case INTEG_LSODE: integrator_lsode_create(sys); break;
258 #ifdef ASC_WITH_IDA
259 case INTEG_IDA: integrator_ida_create(sys); break;
260 #endif
261 default: break;
262 }
263 }
264
265 /*------------------------------------------------------------------------------
266 PARAMETERS
267 */
268
269 /**
270 Reset the parameters in this IntegratorSystem to the default ones for this
271 Integrator.
272
273 @return 0 on success, 1 on error
274 */
275 static int integrator_params_default(IntegratorSystem *sys){
276 switch(sys->engine){
277 case INTEG_LSODE: return integrator_lsode_params_default(sys);
278 #ifdef ASC_WITH_IDA
279 case INTEG_IDA: return integrator_ida_params_default(sys);
280 #endif
281 default: return 0;
282 }
283 }
284
285 int integrator_params_get(const IntegratorSystem *sys, slv_parameters_t *parameters){
286 asc_assert(sys!=NULL);
287 asc_assert(sys->params.num_parms > 0);
288 memcpy(parameters,&(sys->params),sizeof(slv_parameters_t));
289 return 0;
290 }
291
292 int integrator_params_set(IntegratorSystem *sys, const slv_parameters_t *parameters){
293 asc_assert(sys!=NULL);
294 asc_assert(parameters!=NULL);
295 memcpy(&(sys->params),parameters,sizeof(slv_parameters_t));
296 return 0;
297 }
298
299 /*------------------------------------------------------------------------------
300 ANALYSIS
301
302 Provide two modes in order to provide analysis suitable for solution of both
303 ODEs (the previous technique) and DAEs (new code). These share a common
304 "visit" method that needs to eventually be integrated with the code in
305 <solver/analyze.c>. For the moment, we're just hacking in to the compiler.
306 */
307
308 /**
309 Locate the independent variable. For the purpose of GUI design, this needs
310 to work independent of the integration engine being used.
311 */
312 int integrator_find_indep_var(IntegratorSystem *sys){
313 int result = 0;
314
315 if(sys->x != NULL){
316 CONSOLE_DEBUG("sys->x already set");
317 return 1;
318 }
319 assert(sys->indepvars==NULL);
320 sys->indepvars = gl_create(10L);
321
322 IntegInitSymbols();
323
324 /* CONSOLE_DEBUG("About to visit..."); */
325 integrator_visit_system_vars(sys,&integrator_classify_indep_var);
326
327 /* CONSOLE_DEBUG("..."); */
328
329 result = integrator_check_indep_var(sys);
330 gl_free_and_destroy(sys->indepvars);
331 sys->indepvars = NULL;
332
333 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Returning result %d",result); */
334
335 return result;
336 }
337
338 /**
339 Analyse the system, either as DAE or as an ODE system, depending on the
340 solver engine selected.
341
342 @return 1 on success
343 */
344 int integrator_analyse(IntegratorSystem *sys){
345 switch(sys->engine){
346 case INTEG_LSODE: return integrator_analyse_ode(sys);
347 #ifdef ASC_WITH_IDA
348 case INTEG_IDA: return integrator_analyse_dae(sys);
349 #endif
350 case INTEG_UNKNOWN:
351 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No engine selected: can't analyse");
352 default:
353 ERROR_REPORTER_HERE(ASC_PROG_ERR
354 ,"The selected integration engine (%d) is not available"
355 ,sys->engine
356 );
357 }
358 return 0;
359 }
360
361 /**
362 To analyse a DAE we need to identify *ALL* variables in the system
363 Except for the highest-level derivatives of any present?
364 We also need to identify the independent variable (just one).
365
366 @TODO implement Pantelides algorithm in here?
367 @TODO prevent re-analysis without clearing out the data structures?
368 @return 1 on success
369 */
370 int integrator_analyse_dae(IntegratorSystem *sys){
371 struct Integ_var_t *info, *prev;
372 char *varname, *derivname;
373 struct var_variable **varlist;
374 int nvarlist;
375 int i, j;
376 int numstates;
377 int numy, nrels;
378 int yindex;
379 int maxderiv;
380
381 CONSOLE_DEBUG("Starting DAE analysis");
382 IntegInitSymbols();
383
384 assert(sys->indepvars==NULL);
385
386 sys->indepvars = gl_create(10L); /* t var info */
387 sys->dynvars = gl_create(200L); /* y and ydot var info */
388 sys->obslist = gl_create(100L); /* obs info */
389
390 if(sys->indepvars==NULL
391 || sys->dynvars==NULL
392 || sys->obslist==NULL
393 ){
394 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
395 return 0;
396 }
397
398 integrator_visit_system_vars(sys,&integrator_dae_classify_var);
399
400 CONSOLE_DEBUG("Found %lu observation variables:",gl_length(sys->obslist));
401 for(i=1; i<=gl_length(sys->obslist); ++i){
402 info = (struct Integ_var_t *)gl_fetch(sys->obslist, i);
403 varname = var_make_name(sys->system,info->i);
404 CONSOLE_DEBUG("observation[%d] = \"%s\"",i,varname);
405 ASC_FREE(varname);
406 }
407
408 /* CONSOLE_DEBUG("Checking found vars..."); */
409 if(gl_length(sys->dynvars)==0){
410 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No solver_var vars found to integrate (check 'ode_type'?).");
411 return 0;
412 }
413
414 CONSOLE_DEBUG("Found %lu vars.", gl_length(sys->dynvars));
415
416 maxderiv = 0;
417 numstates = 0;
418 for(i=1; i<=gl_length(sys->dynvars); ++i){
419 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
420 if(info->type==1 || info->type==0)numstates++;
421 if(maxderiv < info->type - 1)maxderiv = info->type - 1;
422 /* varname = var_make_name(sys->system,info->i);
423 CONSOLE_DEBUG("var[%d] = \"%s\": ode_index = %ld",i,varname,info->type);
424 ASC_FREE(varname); */
425 }
426 if(maxderiv == 0){
427 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No derivatives found (check 'ode_type' values for your vars).");
428 return 0;
429 }
430 if(maxderiv > 1){
431 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Higher-order derivatives found. You must provide a reduced order formulation for your model.");
432 return 0;
433 }
434
435 if(!integrator_check_indep_var(sys))return 0;
436
437 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
438 /* !!! NOTE THAT dynvars IS NOW REARRANGED !!! *sigh* bug hunting(!) */
439
440 CONSOLE_DEBUG("Variables rearranged to increasing ODE ID & TYPE (varindx = matrix order)");
441 for(i=1; i<=gl_length(sys->dynvars); ++i){
442 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
443 varname = var_make_name(sys->system,info->i);
444 CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld (varindx = %d)",i-1,varname,info->type,info->varindx);
445 ASC_FREE(varname);
446 }
447
448 /* link up variables with their derivatives */
449 prev = NULL;
450 for(i=1; i<=gl_length(sys->dynvars); ++i){ /* why does gl_list index with base 1??? */
451 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
452
453 if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){
454 varname = var_make_name(sys->system,info->i);
455 CONSOLE_DEBUG("Var \"%s\" is an algebraic variable",varname);
456 ASC_FREE(varname);
457 info->type = INTEG_STATE_VAR;
458 info->derivative_of = NULL;
459 }else{
460 if(prev==NULL || info->index != prev->index){
461 /* derivative, but without undifferentiated var present in model */
462 varname = var_make_name(sys->system,info->i);
463 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Derivative %d of \"%s\" is present without its un-differentiated equivalent"
464 , info->type-1
465 , varname
466 );
467 ASC_FREE(varname);
468 return 0;
469 }else if(info->type != prev->type + 1){
470 /* derivative, but missing the next-lower-order derivative */
471 derivname = var_make_name(sys->system,info->i);
472 varname = var_make_name(sys->system,prev->i);
473 ERROR_REPORTER_HERE(ASC_USER_ERROR
474 ,"Looking at \"%s\", expected a derivative (order %d) of \"%s\"."
475 ,varname
476 ,prev->type+1
477 ,derivname
478 );
479 ASC_FREE(varname);
480 ASC_FREE(derivname);
481 return 0;
482 }else{
483 /* variable with derivative */
484 varname = var_make_name(sys->system,prev->i);
485 derivname = var_make_name(sys->system,info->i);
486 CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
487 ASC_FREE(varname);
488 ASC_FREE(derivname);
489 info->derivative_of = prev;
490 }
491 }
492 prev = info;
493 }
494
495 /* record which vars have derivatives and which don't, and count 'states' */
496 numy = 0;
497 for(i=1; i<=gl_length(sys->dynvars); ++i){
498 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
499 if(info->derivative_of){
500 info->derivative_of->derivative = info;
501 }else{
502 numy++;
503 }
504 }
505
506 /* allocate storage for the 'y' and 'ydot' arrays */
507 sys->y = ASC_NEW_ARRAY(struct var_variable *,numy);
508 sys->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,numy);
509 sys->y_id = ASC_NEW_ARRAY(int, slv_get_num_master_vars(sys->system));
510
511 for(i=0; i<numy; ++i){
512 asc_assert(sys->ydot[i]==NULL);
513 }
514
515 /* now add variables and their derivatives to 'ydot' and 'y' */
516 yindex = 0;
517
518 CONSOLE_DEBUG("VARS IN THEIR MATRIX ORDER");
519 for(i=1; i<=gl_length(sys->dynvars); ++i){
520 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
521 if(info->derivative_of)continue;
522 if(info->derivative){
523 sys->y[yindex] = info->i;
524 assert(info->derivative);
525 sys->ydot[yindex] = info->derivative->i;
526 if(info->varindx >= 0){
527 sys->y_id[info->varindx - 1] = yindex;
528 CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
529 }
530 if(info->derivative->varindx >= 0){
531 sys->y_id[info->derivative->varindx - 1] = -1-yindex;
532 CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);
533 }
534 }else{
535 sys->y[yindex] = info ->i;
536 sys->ydot[yindex] = NULL;
537 if(info->varindx >= 0){
538 sys->y_id[info->varindx - 1] = yindex;
539 CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
540 }
541 }
542 yindex++;
543 }
544
545 nrels = slv_get_num_solvers_rels(sys->system);
546 if(numy != nrels){
547 ERROR_REPORTER_HERE(ASC_USER_ERROR
548 ,"System is not square: solver has %d rels, found %d system states"
549 ,nrels, numy
550 );
551 return 0;
552 }
553
554 CONSOLE_DEBUG("THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM",numy);
555
556 sys->n_y = numy;
557
558 if(!integrator_sort_obs_vars(sys))return 0;
559
560 return 1;
561 }
562
563 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
564 struct var_variable **vlist;
565 int i, vlen;
566
567 /* visit all the slv_system_t master var lists to collect vars */
568 /* find the vars mostly in this one */
569 vlist = slv_get_master_var_list(sys->system);
570 vlen = slv_get_num_master_vars(sys->system);
571 for (i=0;i<vlen;i++) {
572 (*visitfn)(sys, vlist[i], &i);
573 }
574
575 /*
576 CONSOLE_DEBUG("Checked %d vars",vlen);
577 integrator_print_var_stats(sys);
578 */
579
580 /* probably nothing here, but gotta check. */
581 vlist = slv_get_master_par_list(sys->system);
582 vlen = slv_get_num_master_pars(sys->system);
583 for (i=0;i<vlen;i++) {
584 (*visitfn)(sys, vlist[i], NULL);
585 }
586
587 /*
588 CONSOLE_DEBUG("Checked %d pars",vlen);
589 integrator_print_var_stats(sys);
590 */
591
592 /* might find t here */
593 vlist = slv_get_master_unattached_list(sys->system);
594 vlen = slv_get_num_master_unattached(sys->system);
595 for (i=0;i<vlen;i++) {
596 (*visitfn)(sys, vlist[i], NULL);
597 }
598
599 /* CONSOLE_DEBUG("Checked %d unattached",vlen); */
600 }
601 /**
602 @return 1 on success
603 */
604 int integrator_analyse_ode(IntegratorSystem *sys){
605 struct Integ_var_t *v1,*v2;
606 long half,i,len;
607 int happy=1;
608 char *varname1, *varname2;
609
610 CONSOLE_DEBUG("Starting ODE analysis");
611 IntegInitSymbols();
612
613 /* collect potential states and derivatives */
614 sys->indepvars = gl_create(10L); /* t var info */
615 sys->dynvars = gl_create(200L); /* y ydot var info */
616 sys->obslist = gl_create(100L); /* obs info */
617 if (sys->dynvars == NULL
618 || sys->obslist == NULL
619 || sys->indepvars == NULL
620 ){
621 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
622 return 0;
623 }
624
625 sys->nstates = sys->nderivs = 0;
626
627 integrator_visit_system_vars(sys,&integrator_ode_classify_var);
628
629 integrator_print_var_stats(sys);
630
631 if(!integrator_check_indep_var(sys))return 0;
632
633 /* check sanity of state and var lists */
634
635 len = gl_length(sys->dynvars);
636 half = len/2;
637 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
638
639 if (len % 2 || len == 0L || sys->nstates != sys->nderivs ) {
640 /* list length must be even for vars to pair off */
641 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
642 return 0;
643 }
644 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
645 if (gl_fetch(sys->dynvars,len)==NULL) {
646 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
647 return 0;
648 }
649 sys->states = gl_create(half); /* state vars Integ_var_t references */
650 sys->derivs = gl_create(half); /* derivative var atoms */
651 for (i=1;i < len; i+=2) {
652 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
653 v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
654 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
655 varname1 = var_make_name(sys->system,v1->i);
656 varname2 = var_make_name(sys->system,v2->i);
657
658 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
659 varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
660 varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
661 );
662 ASC_FREE(varname1);
663 ASC_FREE(varname2);
664 happy=0;
665 break;
666 } else {
667 gl_append_ptr(sys->states,(POINTER)v1);
668 gl_append_ptr(sys->derivs,(POINTER)v2->i);
669 }
670 }
671 if (!happy) {
672 return 0;
673 }
674 sys->n_y = half;
675 sys->y = ASC_NEW_ARRAY(struct var_variable *, half);
676 sys->y_id = ASC_NEW_ARRAY(long, half);
677 sys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
678 if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
679 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
680 return 0;
681 }
682 for (i = 1; i <= half; i++) {
683 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
684 sys->y[i-1] = v1->i;
685 sys->y_id[i-1] = v1->index;
686 sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
687 }
688
689 if(!integrator_sort_obs_vars(sys))return 0;
690
691 /* don't need the gl_lists now that we have arrays for everyone */
692 gl_destroy(sys->states);
693 gl_destroy(sys->derivs);
694 gl_free_and_destroy(sys->indepvars); /* we own the objects in indepvars */
695 gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
696 gl_free_and_destroy(sys->obslist); /* and obslist */
697 sys->states = NULL;
698 sys->derivs = NULL;
699 sys->indepvars = NULL;
700 sys->dynvars = NULL;
701 sys->obslist = NULL;
702
703 /* analysis completed OK */
704 return 1;
705 }
706
707 /**
708 Reindex observations. Sort if the user mostly numbered. Take natural order
709 if user just booleaned.
710
711 @return 1 on success
712 */
713 static int integrator_sort_obs_vars(IntegratorSystem *sys){
714 int half, i, len = 0;
715 struct Integ_var_t *v2;
716
717 half = sys->n_y;
718 len = gl_length(sys->obslist);
719 /* we shouldn't be seeing NULL here ever except if malloc fail. */
720 if (len > 1L) {
721 half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;
722 /* half != 0 now because we didn't collect 0 indexed vars */
723 for (i=2; i <= len; i++) {
724 if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {
725 /* change seen. sort and go on */
726 gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);
727 break;
728 }
729 }
730 }
731 for (i = half = 1; i <= len; i++) {
732 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
733 if (v2==NULL) {
734 /* we shouldn't be seeing NULL here ever except if malloc fail. */
735 gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */
736 } else {
737 Integ_SetObsId(v2->i,half);
738 v2->index = half++;
739 }
740 }
741
742 /* obslist now uniquely indexed, no nulls */
743 /* make into arrays */
744 half = gl_length(sys->obslist);
745 sys->obs = ASC_NEW_ARRAY(struct var_variable *,half);
746 sys->obs_id = ASC_NEW_ARRAY(long, half);
747 if ( sys->obs==NULL || sys->obs_id==NULL) {
748 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
749 return 0;
750 }
751 sys->n_obs = half;
752 for (i = 1; i <= half; i++) {
753 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
754 sys->obs[i-1] = v2->i;
755 sys->obs_id[i-1] = v2->index;
756 }
757
758 return 1;
759 }
760
761 static void integrator_print_var_stats(IntegratorSystem *sys){
762 int v = gl_length(sys->dynvars);
763 int i = gl_length(sys->indepvars);
764 CONSOLE_DEBUG("Currently %d vars, %d indep",v,i);
765 }
766
767 /**
768 Check sanity of the independent variable.
769
770 @return 1 on success
771 */
772 static int integrator_check_indep_var(IntegratorSystem *sys){
773 int len, i;
774 struct Integ_var_t *info;
775 char *varname;
776
777 /* check the sanity of the independent variable */
778 len = gl_length(sys->indepvars);
779 if (!len) {
780 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
781 return 0;
782 }
783 if (len > 1) {
784 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
785 FPRINTF(ASCERR,"Excess %ld independent variables found:",
786 len);
787 for(i=1; i <=len;i++) {
788 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
789 if(info==NULL)continue;
790
791 varname = var_make_name(sys->system,info->i);
792 FPRINTF(ASCERR," %s",varname);
793 ASC_FREE(varname);
794 }
795 FPRINTF(ASCERR , "\nSet the \"%s\" flag on all but one of these to %s >= 0.\n"
796 , SCP(STATEFLAG),SCP(STATEFLAG)
797 );
798 error_reporter_end_flush();
799 return 0;
800 }else{
801 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
802 sys->x = info->i;
803 }
804 return 1;
805 }
806
807 /*------------------------------------------------------------------------------
808 CLASSIFICATION OF VARIABLES (for ANALYSIS step)
809 */
810
811 #define INTEG_ADD_TO_LIST(info,TYPE,INDEX,VAR,VARINDX,LIST) \
812 info = ASC_NEW(struct Integ_var_t); \
813 if(info==NULL){ \
814 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory (INTEG_VAR_NEW)"); \
815 return; \
816 } \
817 info->type=TYPE; \
818 info->index=INDEX; \
819 info->i=VAR; \
820 info->derivative=NULL; \
821 info->derivative_of=NULL; \
822 if(VARINDX==NULL){ \
823 info->varindx = -1; \
824 }else{ \
825 info->varindx = *VARINDX; \
826 } \
827 gl_append_ptr(LIST,(void *)info); \
828 info = NULL
829
830 /**
831 In a DAE, it's either the (single) independent variable, or it's a
832 variable in the model.
833
834 I'm not sure what we should be doing with variables that are already
835 present as derivatives of other variables, I guess those ones need to be
836 removed from the list in a second pass?
837 */
838 void integrator_dae_classify_var(IntegratorSystem *sys
839 , struct var_variable *var, const int *varindx
840 ){
841 struct Integ_var_t *info;
842 long type,index;
843
844 /* filter for recognition of solver_vars */
845 var_filter_t vfilt;
846 vfilt.matchbits = VAR_SVAR;
847 vfilt.matchvalue = VAR_SVAR;
848
849 assert(var != NULL && var_instance(var)!=NULL );
850
851 if( var_apply_filter(var,&vfilt) ) {
852 if(!var_active(var)){
853 CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");
854 return;
855 }
856
857 /* only non-fixed variables are accepted */
858 if(!var_fixed(var)){
859 /* get the ode_type and ode_id of this solver_var */
860 type = DynamicVarInfo(var,&index);
861
862 if(type==INTEG_OTHER_VAR){
863 /* if the var's type is -1, it's independent */
864 INTEG_ADD_TO_LIST(info,INTEG_OTHER_VAR,0,var,varindx,sys->indepvars);
865 }else{
866 if(type < 0)type=0;
867 /* any other type of var is in the DAE system, at least for now */
868 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
869 }
870 }
871 #if 0
872 else{
873 /* fixed variable, only include it if ode_type == 1 */
874 type = DynamicVarInfo(var,&index);
875 if(type==INTEG_STATE_VAR){
876 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
877 }
878 }
879 #endif
880
881 /* if the var's obs_id > 0, add it to the observation list */
882 if(ObservationVar(var,&index) != NULL && index > 0L) {
883 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->obslist);
884 }
885 }
886 }
887
888 /**
889 Inspect a specific variable and work out what type it is (what 'ode_type' it
890 has) and what other variable(s) it corresponds to (ie dydt corresponds to
891 y as a derivative).
892
893 @TODO add ability to create new variables for 'missing' derivative vars?
894 */
895 void integrator_ode_classify_var(IntegratorSystem *sys, struct var_variable *var
896 , const int *varindx
897 ){
898 struct Integ_var_t *info;
899 long type,index;
900
901 var_filter_t vfilt;
902 vfilt.matchbits = VAR_SVAR;
903 vfilt.matchvalue = VAR_SVAR;
904
905 assert(var != NULL && var_instance(var)!=NULL );
906
907 if( var_apply_filter(var,&vfilt) ) {
908 /* it's a solver var: what type of variable? */
909 type = DynamicVarInfo(var,&index);
910
911 if(type==INTEG_ALGEBRAIC_VAR){
912 /* no action required */
913 }else if(type==INTEG_OTHER_VAR){
914 /* i.e. independent var */
915 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
916 }else if(type>=INTEG_STATE_VAR){
917 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
918 if(type == 1){
919 sys->nstates++;
920 }else if(type == 2){ /* what about higher-order derivatives? -- JP */
921 sys->nderivs++;
922 }else{
923 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Higher-order (>=2) derivatives are not supported in ODEs.");
924 } }
925
926 if(ObservationVar(var,&index) != NULL && index > 0L) {
927 INTEG_ADD_TO_LIST(info,0L,index,var,varindx,sys->obslist);
928 }
929 }
930 }
931
932 /**
933 Look at a variable and determine if it's the independent variable or not.
934 This is just for the purpose of the integrator_find_indep_var function,
935 which is a utility function provided for use by the GUI.
936 */
937 void integrator_classify_indep_var(IntegratorSystem *sys
938 , struct var_variable *var, const int *varindx
939 ){
940 struct Integ_var_t *info;
941 long type,index;
942
943 var_filter_t vfilt;
944 vfilt.matchbits = VAR_SVAR;
945 vfilt.matchvalue = VAR_SVAR;
946
947 /* CONSOLE_DEBUG("..."); */
948
949 assert(var != NULL && var_instance(var)!=NULL );
950
951 if( var_apply_filter(var,&vfilt) ) {
952 type = DynamicVarInfo(var,&index);
953
954 if(type==INTEG_OTHER_VAR){
955 /* i.e. independent var */
956 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
957 }
958 }
959 }
960
961 /**
962 Look at a variable, and if it is an 'ODE variable' (it has a child instance
963 named 'ode_type') return its type, which will be either:
964 - INTEG_OTHER_VAR (if 'ode_type' is -1)
965 - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
966 - INTEG_STATE_VAR (if 'ode_type' is 1)
967 - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
968
969 If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
970 there.
971 */
972 static long DynamicVarInfo(struct var_variable *v,long *index){
973 struct Instance *c, *d, *i;
974 i = var_instance(v);
975 assert(i!=NULL);
976 assert(STATEFLAG!=NULL);
977 assert(STATEINDEX!=NULL);
978 c = ChildByChar(i,STATEFLAG);
979 d = ChildByChar(i,STATEINDEX);
980 /* lazy evaluation is important in the following if */
981 if(c == NULL
982 || d == NULL
983 || InstanceKind(c) != INTEGER_INST
984 || InstanceKind(d) != INTEGER_INST
985 || !AtomAssigned(c)
986 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
987 ){
988 return INTEG_ALGEBRAIC_VAR;
989 }
990 if (index != NULL) {
991 *index = GetIntegerAtomValue(d);
992 }
993 return GetIntegerAtomValue(c);
994 }
995
996 /**
997 Looks at the given variable checks if it is an 'observation variable'. This
998 means that it has its 'obs_id' child instance set to a non-zero value.
999
1000 If the variable is an observation variable, its index value ('obs_id') is
1001 stuff into *index (provided index!=NULL), and the pointer to the original
1002 instance is rtruend.
1003
1004 If it's not an observation variable, we return NULL and *index is untouched.
1005 */
1006 static struct var_variable *ObservationVar(struct var_variable *v, long *index){
1007 struct Instance *c,*i;
1008 i = var_instance(v);
1009 assert(i!=NULL);
1010 c = ChildByChar(i,OBSINDEX);
1011 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1012 return NULL;
1013 }
1014 if (index != NULL) {
1015 *index = GetIntegerAtomValue(c);
1016 }
1017 return v;
1018 }
1019
1020 /*------------------------------------------------------------------------------
1021 RUNNING THE SOLVER
1022 */
1023
1024 /*
1025 Make the call to the actual integrator we've selected, for the range of
1026 time values specified. The sys contains all the specifics.
1027
1028 Return 1 on success
1029 */
1030 int integrator_solve(IntegratorSystem *sys, long i0, long i1){
1031
1032 long nstep;
1033 unsigned long start_index=0, finish_index=0;
1034 assert(sys!=NULL);
1035
1036 nstep = integrator_getnsamples(sys)-1;
1037 /* check for at least 2 steps and dimensionality of x vs steps here */
1038
1039 if (i0<0 || i1 <0) {
1040 /* dude, there's no way we're writing interactive stuff here... */
1041 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
1042 return 0;
1043 } else {
1044 start_index=i0;
1045 finish_index =i1;
1046 if (start_index >= (unsigned long)nstep) {
1047 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
1048 ,start_index,nstep
1049 );
1050 return 0;
1051 }
1052 if (finish_index > (unsigned long)nstep) {
1053 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
1054 ,finish_index,nstep
1055 );
1056 return 0;
1057 }
1058 }
1059
1060 if(finish_index <= start_index) {
1061 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
1062 ,start_index,finish_index
1063 );
1064 return 0;
1065 }
1066
1067 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1068
1069 /* now go and run the integrator */
1070 switch (sys->engine) {
1071 case INTEG_LSODE:
1072 return integrator_lsode_solve(sys, start_index, finish_index);
1073 break;
1074 #ifdef ASC_WITH_IDA
1075 case INTEG_IDA:
1076 return integrator_ida_solve(sys,start_index, finish_index);
1077 break;
1078 #endif
1079 default:
1080 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator (invalid, or not implemented yet)");
1081 return 0;
1082 }
1083 }
1084
1085 /*---------------------------------------------------------------
1086 HANDLING THE LIST OF TIMESTEMPS
1087 */
1088
1089 #define GETTER_AND_SETTER(TYPE,NAME) \
1090 void integrator_set_##NAME(IntegratorSystem *sys, TYPE val){ \
1091 sys->NAME=val; \
1092 } \
1093 TYPE integrator_get_##NAME(IntegratorSystem *sys){ \
1094 return sys->NAME; \
1095 }
1096
1097 GETTER_AND_SETTER(SampleList *,samples) /*;*/
1098 GETTER_AND_SETTER(double,maxstep) /*;*/
1099 GETTER_AND_SETTER(double,minstep) /*;*/
1100 GETTER_AND_SETTER(double,stepzero) /*;*/
1101 GETTER_AND_SETTER(int,maxsubsteps) /*;*/
1102 #undef GETTER_AND_SETTER
1103
1104 long integrator_getnsamples(IntegratorSystem *sys){
1105 assert(sys!=NULL);
1106 assert(sys->samples!=NULL);
1107 return samplelist_length(sys->samples);
1108 }
1109
1110 double integrator_getsample(IntegratorSystem *sys, long i){
1111 assert(sys!=NULL);
1112 assert(sys->samples!=NULL);
1113 return samplelist_get(sys->samples,i);
1114 }
1115
1116 void integrator_setsample(IntegratorSystem *sys, long i,double xi){
1117 assert(sys!=NULL);
1118 assert(sys->samples!=NULL);
1119 samplelist_set(sys->samples,i,xi);
1120 }
1121
1122 const dim_type *integrator_getsampledim(IntegratorSystem *sys){
1123 assert(sys!=NULL);
1124 assert(sys->samples!=NULL);
1125 return samplelist_dim(sys->samples);
1126 }
1127
1128 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *sys){
1129 return sys->currentstep;
1130 }
1131
1132 /*------------------------------------------------------------------------------
1133 GET/SET VALUE OF THE INDEP VARIABLE
1134 */
1135
1136 /**
1137 Retrieve the value of the independent variable (time) from ASCEND
1138 and return it as a double.
1139 */
1140 double integrator_get_t(IntegratorSystem *sys){
1141 assert(sys->x!=NULL);
1142 return var_value(sys->x);
1143 }
1144
1145 /**
1146 Set the value of the independent variable (time) in ASCEND.
1147 */
1148 void integrator_set_t(IntegratorSystem *sys, double value){
1149 var_set_value(sys->x, value);
1150 /* CONSOLE_DEBUG("set_t = %g", value); */
1151 }
1152
1153 /*------------------------------------------------------------------------------
1154 PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1155 */
1156 /**
1157 Retrieve the current values of the derivatives of the y-variables
1158 and stick them in the/an array that the integrator will use.
1159
1160 If the pointer 'y' is NULL, the necessary space is allocated (and
1161 must be freed somewhere else).
1162 */
1163 double *integrator_get_y(IntegratorSystem *sys, double *y) {
1164 long i;
1165
1166 if (y==NULL) {
1167 y = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1168 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1169 }
1170
1171 for (i=0; i< sys->n_y; i++) {
1172 assert(sys->y[i]!=NULL);
1173 y[i] = var_value(sys->y[i]);
1174 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1175 }
1176 return y;
1177 }
1178
1179 /**
1180 Take the values of the derivatives from the array that the integrator
1181 uses, and use them to update the values of the corresponding variables
1182 in ASCEND.
1183 */
1184 void integrator_set_y(IntegratorSystem *sys, double *y) {
1185 long i;
1186 #ifndef NDEBUG
1187 char *varname;
1188 #endif
1189
1190 for (i=0; i < sys->n_y; i++) {
1191 assert(sys->y[i]!=NULL);
1192 var_set_value(sys->y[i],y[i]);
1193 #ifndef NDEBUG
1194 varname = var_make_name(sys->system, sys->y[i]);
1195 /* CONSOLE_DEBUG("y[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, y[i]); */
1196 ASC_FREE(varname);
1197 #endif
1198 }
1199 }
1200
1201 /**
1202 Send the values of the derivatives of the 'y' variables to the solver.
1203 Allocate space for an array if necessary.
1204
1205 Any element in sys->ydot that is NULL will be passed over (the value
1206 won't be modified in dydx).
1207 */
1208 double *integrator_get_ydot(IntegratorSystem *sys, double *dydx) {
1209 long i;
1210
1211 if (dydx==NULL) {
1212 dydx = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1213 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1214 }
1215
1216 for (i=0; i < sys->n_y; i++) {
1217 if(sys->ydot[i]!=NULL){
1218 dydx[i] = var_value(sys->ydot[i]);
1219 }
1220 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1221 }
1222 return dydx;
1223 }
1224
1225 void integrator_set_ydot(IntegratorSystem *sys, double *dydx) {
1226 long i;
1227 #ifndef NDEBUG
1228 /* char *varname; */
1229 #endif
1230 for (i=0; i < sys->n_y; i++) {
1231 if(sys->ydot[i]!=NULL){
1232 var_set_value(sys->ydot[i],dydx[i]);
1233 #ifndef NDEBUG
1234 /* varname = var_make_name(sys->system, sys->ydot[i]);
1235 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1236 ASC_FREE(varname); */
1237 #endif
1238 }
1239 #ifndef NDEBUG
1240 /*else{
1241 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1242 }*/
1243 #endif
1244 }
1245 }
1246
1247 /**
1248 Retrieve the values of 'ode_atol' properties of each of y-variables,
1249 for use in setting absolute error tolerances for the Integrator.
1250
1251 If the pointer 'atol' is NULL, the necessary space is allocated (and
1252 must be freed somewhere else).
1253 */
1254 double *integrator_get_atol(IntegratorSystem *sys, double *atol){
1255 long i;
1256 char *varname;
1257
1258 if (atol==NULL) {
1259 atol = ASC_NEW_ARRAY_CLEAR(double, sys->n_y);
1260 }
1261
1262 for (i=0; i< sys->n_y; i++) {
1263 assert(sys->y[i]!=NULL);
1264 atol[i] = var_odeatol(sys->y[i]);
1265 assert(atol[i]!=-1);
1266 varname = var_make_name(sys->system,sys->y[i]);
1267 CONSOLE_DEBUG("%s.ode_atol = %8.2e",varname,atol[i]);
1268 ASC_FREE(varname);
1269 }
1270 return atol;
1271 }
1272
1273 /*-------------------------------------------------------------
1274 RETRIEVING OBSERVATION DATA
1275 */
1276
1277 /**
1278 This function takes the inst in the solver and returns the vector of
1279 observation variables that are located in the submodel d.obs array.
1280 */
1281 double *integrator_get_observations(IntegratorSystem *sys, double *obsi) {
1282 long i;
1283
1284 if (obsi==NULL) {
1285 obsi = ASC_NEW_ARRAY_CLEAR(double, sys->n_obs+1);
1286 }
1287
1288 /* C obsi[0] <==> ascend d.obs[1] */
1289
1290 for (i=0; i < sys->n_obs; i++) {
1291 obsi[i] = var_value(sys->obs[i]);
1292 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1293 }
1294 return obsi;
1295 }
1296
1297 struct var_variable *integrator_get_observed_var(IntegratorSystem *sys, const long i){
1298 assert(i>=0);
1299 assert(i<sys->n_obs);
1300 return sys->obs[i];
1301 }
1302
1303 /**
1304 @NOTE Although this shouldn't be required for implementation of solver
1305 engines, this is useful for GUI reporting of integration results.
1306 */
1307 struct var_variable *integrator_get_independent_var(IntegratorSystem *sys){
1308 return sys->x;
1309 }
1310
1311
1312 /*----------------------------------------------------
1313 Build an analytic jacobian for solving the state system
1314
1315 This necessarily ugly piece of code attempts to create a unique
1316 list of relations that explicitly contain the variables in the
1317 given input list. The utility of this information is that we know
1318 exactly which relations must be differentiated, to fill in the
1319 df/dy matrix. If the problem has very few derivative terms, this will
1320 be of great savings. If the problem arose from the discretization of
1321 a pde, then this will be not so useful. The decision wether to use
1322 this function or to simply differentiate the entire relations list
1323 must be done before calling this function.
1324
1325 Final Note: the callee owns the array, but not the array elements.
1326 */
1327 #define AVG_NUM_INCIDENT 4
1328
1329
1330 /**
1331 This function helps to arrange the observation variables in a sensible order.
1332 The 'obs_id' child instance of v, if present, is assigned the value of the
1333 given parameter 'index'.
1334 */
1335 static void Integ_SetObsId(struct var_variable *v, long index){
1336 struct Instance *c, *i;
1337 i = var_instance(v);
1338 assert(i!=NULL);
1339 c = ChildByChar(i,OBSINDEX);
1340 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1341 return;
1342 }
1343 SetIntegerAtomValue(c,index,0);
1344 }
1345
1346 /**
1347 Compares observation structs. NULLs should end up at far end.
1348 */
1349 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1350 if(v1 == NULL)return 1;
1351 if(v2 == NULL)return -1;
1352 if(v1->index > v2->index)return 1;
1353 if(v1->index == v2->index)return 0;
1354 return -1;
1355 }
1356
1357 /**
1358 Compares dynamic vars structs. NULLs should end up at far end.
1359 List should be sorted primarily by index and then by type, in order
1360 of increasing value of both.
1361 */
1362 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1363 if(v1 == NULL)return 1;
1364 if(v2 == NULL)return -1;
1365 if(v1->index > v2->index)return 1;
1366 if(v1->index != v2->index)return -1;
1367 if(v1->type > v2->type)return 1;
1368 return -1;
1369 }
1370 /*----------------------------
1371 Output handling to the GUI/interface.
1372 */
1373
1374 int integrator_set_reporter(IntegratorSystem *sys
1375 , IntegratorReporter *reporter
1376 ){
1377 assert(sys!=NULL);
1378 sys->reporter = reporter;
1379 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1380 return 1;
1381 }
1382
1383 int integrator_output_init(IntegratorSystem *sys){
1384 assert(sys!=NULL);
1385 assert(sys->reporter!=NULL);
1386 if(sys->reporter->init!=NULL){
1387 /* call the specified output function */
1388 return (*(sys->reporter->init))(sys);
1389 }
1390 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1391 return 1;
1392 }
1393
1394 int integrator_output_write(IntegratorSystem *sys){
1395 static int reported_already=0;
1396 assert(sys!=NULL);
1397 if(sys->reporter->write!=NULL){
1398 return (*(sys->reporter->write))(sys);
1399 }
1400 if(!reported_already){
1401 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1402 reported_already=1;
1403 }
1404 return 1;
1405 }
1406
1407 int integrator_output_write_obs(IntegratorSystem *sys){
1408 static int reported_already=0;
1409 assert(sys!=NULL);
1410 if(sys->reporter->write_obs!=NULL){
1411 return (*(sys->reporter->write_obs))(sys);
1412 }
1413 if(!reported_already){
1414 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1415 reported_already=1;
1416 }
1417 return 1;
1418 }
1419
1420 int integrator_output_close(IntegratorSystem *sys){
1421 assert(sys!=NULL);
1422 if(sys->reporter->close!=NULL){
1423 return (*(sys->reporter->close))(sys);
1424 }
1425 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1426 return 1;
1427 }
1428
1429 /**
1430 Decode status codes from the integrator, and output them via FPRINTF.
1431 */
1432 int integrator_checkstatus(slv_status_t status) {
1433 if (status.converged) {
1434 return 1;
1435 }
1436 if (status.diverged) {
1437 FPRINTF(stderr, "The derivative system did not converge. Integration will terminate.");
1438 return 0;
1439 }
1440 if (status.inconsistent) {
1441 FPRINTF(stderr, "A numerically inconsistent state was discovered while "
1442 "calculating derivatives. Integration will terminate.");
1443 return 0;
1444 }
1445 if (status.time_limit_exceeded) {
1446 FPRINTF(stderr, "The time limit was exceeded while calculating "
1447 "derivatives. Integration will terminate.");
1448 return 0;
1449 }
1450 if (status.iteration_limit_exceeded) {
1451 FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1452 "derivatives. Integration will terminate.");
1453 return 0;
1454 }
1455 if (status.panic) {
1456 FPRINTF(stderr, "The user patience limit was exceeded while "
1457 "calculating derivatives. Integration will terminate.");
1458 return 0;
1459 }
1460 return 0;
1461 }

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