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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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