/[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 944 - (show annotations) (download) (as text)
Sat Nov 25 10:46:13 2006 UTC (14 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 42730 byte(s)
Implemented ATOLVECT, ATOL, RTOL parameters for the IDA integrator.
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
439 fprintf(stderr,"\n\n\nSORTED VARS\n");
440 for(i=1; i<=gl_length(sys->dynvars); ++i){
441 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
442 varname = var_make_name(sys->system,info->i);
443 CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld (varindx = %d)",i,varname,info->type,info->varindx);
444 ASC_FREE(varname);
445 }
446
447 /* link up variables with their derivatives */
448 prev = NULL;
449 for(i=1; i<=gl_length(sys->dynvars); ++i){ /* why does gl_list index with base 1??? */
450 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
451
452 if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){
453 varname = var_make_name(sys->system,info->i);
454 CONSOLE_DEBUG("Var \"%s\" is an algebraic variable",varname);
455 ASC_FREE(varname);
456 info->type = INTEG_STATE_VAR;
457 info->derivative_of = NULL;
458 }else{
459 if(prev==NULL || info->index != prev->index){
460 /* derivative, but without undifferentiated var present in model */
461 varname = var_make_name(sys->system,info->i);
462 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Derivative %d of \"%s\" is present without its un-differentiated equivalent"
463 , info->type-1
464 , varname
465 );
466 ASC_FREE(varname);
467 return 0;
468 }else if(info->type != prev->type + 1){
469 /* derivative, but missing the next-lower-order derivative */
470 derivname = var_make_name(sys->system,info->i);
471 varname = var_make_name(sys->system,prev->i);
472 ERROR_REPORTER_HERE(ASC_USER_ERROR
473 ,"Looking at \"%s\", expected a derivative (order %d) of \"%s\"."
474 ,varname
475 ,prev->type+1
476 ,derivname
477 );
478 ASC_FREE(varname);
479 ASC_FREE(derivname);
480 return 0;
481 }else{
482 /* variable with derivative */
483 varname = var_make_name(sys->system,prev->i);
484 derivname = var_make_name(sys->system,info->i);
485 CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
486 ASC_FREE(varname);
487 ASC_FREE(derivname);
488 info->derivative_of = prev;
489 }
490 }
491 prev = info;
492 }
493
494 /* record which vars have derivatives and which don't, and count 'states' */
495 numy = 0;
496 for(i=1; i<=gl_length(sys->dynvars); ++i){
497 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
498 if(info->derivative_of){
499 info->derivative_of->derivative = info;
500 }else{
501 numy++;
502 }
503 }
504
505 /* allocate storage for the 'y' and 'ydot' arrays */
506 sys->y = ASC_NEW_ARRAY(struct var_variable *,numy);
507 sys->ydot = ASC_NEW_ARRAY(struct var_variable *,numy);
508 sys->y_id = ASC_NEW_ARRAY(int, slv_get_num_master_vars(sys->system));
509
510 /* now add variables and their derivatives to 'ydot' and 'y' */
511 yindex = 0;
512
513 for(i=1; i<=gl_length(sys->dynvars); ++i){
514 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
515 if(info->derivative_of)continue;
516 if(info->derivative){
517 sys->y[yindex] = info->i;
518 assert(info->derivative);
519 sys->ydot[yindex] = info->derivative->i;
520 if(info->varindx >= 0){
521 sys->y_id[info->varindx] = yindex;
522 CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
523 }
524 if(info->derivative->varindx >= 0){
525 sys->y_id[info->derivative->varindx] = -1-yindex;
526 CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);
527 }
528 }else{
529 sys->y[yindex] = info ->i;
530 sys->ydot[yindex] = NULL;
531 if(info->varindx >= 0){
532 sys->y_id[info->varindx] = yindex;
533 CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
534 }
535 }
536 yindex++;
537 }
538
539 nrels = slv_get_num_solvers_rels(sys->system);
540 if(numy != nrels){
541 ERROR_REPORTER_HERE(ASC_USER_ERROR
542 ,"System is not square: solver has %d rels, found %d system states"
543 ,nrels, numy
544 );
545 return 0;
546 }
547
548 CONSOLE_DEBUG("THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM",numy);
549
550 sys->n_y = numy;
551
552 if(!integrator_sort_obs_vars(sys))return 0;
553
554 return 1;
555 }
556
557 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
558 struct var_variable **vlist;
559 int i, vlen;
560
561 /* visit all the slv_system_t master var lists to collect vars */
562 /* find the vars mostly in this one */
563 vlist = slv_get_master_var_list(sys->system);
564 vlen = slv_get_num_master_vars(sys->system);
565 for (i=0;i<vlen;i++) {
566 (*visitfn)(sys, vlist[i], &i);
567 }
568
569 /*
570 CONSOLE_DEBUG("Checked %d vars",vlen);
571 integrator_print_var_stats(sys);
572 */
573
574 /* probably nothing here, but gotta check. */
575 vlist = slv_get_master_par_list(sys->system);
576 vlen = slv_get_num_master_pars(sys->system);
577 for (i=0;i<vlen;i++) {
578 (*visitfn)(sys, vlist[i], NULL);
579 }
580
581 /*
582 CONSOLE_DEBUG("Checked %d pars",vlen);
583 integrator_print_var_stats(sys);
584 */
585
586 /* might find t here */
587 vlist = slv_get_master_unattached_list(sys->system);
588 vlen = slv_get_num_master_unattached(sys->system);
589 for (i=0;i<vlen;i++) {
590 (*visitfn)(sys, vlist[i], NULL);
591 }
592
593 /* CONSOLE_DEBUG("Checked %d unattached",vlen); */
594 }
595 /**
596 @return 1 on success
597 */
598 int integrator_analyse_ode(IntegratorSystem *sys){
599 struct Integ_var_t *v1,*v2;
600 long half,i,len;
601 int happy=1;
602 char *varname1, *varname2;
603
604 CONSOLE_DEBUG("Starting ODE analysis");
605 IntegInitSymbols();
606
607 /* collect potential states and derivatives */
608 sys->indepvars = gl_create(10L); /* t var info */
609 sys->dynvars = gl_create(200L); /* y ydot var info */
610 sys->obslist = gl_create(100L); /* obs info */
611 if (sys->dynvars == NULL
612 || sys->obslist == NULL
613 || sys->indepvars == NULL
614 ){
615 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
616 return 0;
617 }
618
619 sys->nstates = sys->nderivs = 0;
620
621 integrator_visit_system_vars(sys,&integrator_ode_classify_var);
622
623 integrator_print_var_stats(sys);
624
625 if(!integrator_check_indep_var(sys))return 0;
626
627 /* check sanity of state and var lists */
628
629 len = gl_length(sys->dynvars);
630 half = len/2;
631 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
632
633 if (len % 2 || len == 0L || sys->nstates != sys->nderivs ) {
634 /* list length must be even for vars to pair off */
635 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
636 return 0;
637 }
638 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
639 if (gl_fetch(sys->dynvars,len)==NULL) {
640 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
641 return 0;
642 }
643 sys->states = gl_create(half); /* state vars Integ_var_t references */
644 sys->derivs = gl_create(half); /* derivative var atoms */
645 for (i=1;i < len; i+=2) {
646 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
647 v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
648 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
649 varname1 = var_make_name(sys->system,v1->i);
650 varname2 = var_make_name(sys->system,v2->i);
651
652 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
653 varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
654 varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
655 );
656 ASC_FREE(varname1);
657 ASC_FREE(varname2);
658 happy=0;
659 break;
660 } else {
661 gl_append_ptr(sys->states,(POINTER)v1);
662 gl_append_ptr(sys->derivs,(POINTER)v2->i);
663 }
664 }
665 if (!happy) {
666 return 0;
667 }
668 sys->n_y = half;
669 sys->y = ASC_NEW_ARRAY(struct var_variable *, half);
670 sys->y_id = ASC_NEW_ARRAY(long, half);
671 sys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
672 if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
673 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
674 return 0;
675 }
676 for (i = 1; i <= half; i++) {
677 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
678 sys->y[i-1] = v1->i;
679 sys->y_id[i-1] = v1->index;
680 sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
681 }
682
683 if(!integrator_sort_obs_vars(sys))return 0;
684
685 /* don't need the gl_lists now that we have arrays for everyone */
686 gl_destroy(sys->states);
687 gl_destroy(sys->derivs);
688 gl_free_and_destroy(sys->indepvars); /* we own the objects in indepvars */
689 gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
690 gl_free_and_destroy(sys->obslist); /* and obslist */
691 sys->states = NULL;
692 sys->derivs = NULL;
693 sys->indepvars = NULL;
694 sys->dynvars = NULL;
695 sys->obslist = NULL;
696
697 /* analysis completed OK */
698 return 1;
699 }
700
701 /**
702 Reindex observations. Sort if the user mostly numbered. Take natural order
703 if user just booleaned.
704
705 @return 1 on success
706 */
707 static int integrator_sort_obs_vars(IntegratorSystem *sys){
708 int half, i, len = 0;
709 struct Integ_var_t *v2;
710
711 half = sys->n_y;
712 len = gl_length(sys->obslist);
713 /* we shouldn't be seeing NULL here ever except if malloc fail. */
714 if (len > 1L) {
715 half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;
716 /* half != 0 now because we didn't collect 0 indexed vars */
717 for (i=2; i <= len; i++) {
718 if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {
719 /* change seen. sort and go on */
720 gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);
721 break;
722 }
723 }
724 }
725 for (i = half = 1; i <= len; i++) {
726 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
727 if (v2==NULL) {
728 /* we shouldn't be seeing NULL here ever except if malloc fail. */
729 gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */
730 } else {
731 Integ_SetObsId(v2->i,half);
732 v2->index = half++;
733 }
734 }
735
736 /* obslist now uniquely indexed, no nulls */
737 /* make into arrays */
738 half = gl_length(sys->obslist);
739 sys->obs = ASC_NEW_ARRAY(struct var_variable *,half);
740 sys->obs_id = ASC_NEW_ARRAY(long, half);
741 if ( sys->obs==NULL || sys->obs_id==NULL) {
742 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
743 return 0;
744 }
745 sys->n_obs = half;
746 for (i = 1; i <= half; i++) {
747 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
748 sys->obs[i-1] = v2->i;
749 sys->obs_id[i-1] = v2->index;
750 }
751
752 return 1;
753 }
754
755 static void integrator_print_var_stats(IntegratorSystem *sys){
756 int v = gl_length(sys->dynvars);
757 int i = gl_length(sys->indepvars);
758 CONSOLE_DEBUG("Currently %d vars, %d indep",v,i);
759 }
760
761 /**
762 Check sanity of the independent variable.
763
764 @return 1 on success
765 */
766 static int integrator_check_indep_var(IntegratorSystem *sys){
767 int len, i;
768 struct Integ_var_t *info;
769 char *varname;
770
771 /* check the sanity of the independent variable */
772 len = gl_length(sys->indepvars);
773 if (!len) {
774 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
775 return 0;
776 }
777 if (len > 1) {
778 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
779 FPRINTF(ASCERR,"Excess %ld independent variables found:",
780 len);
781 for(i=1; i <=len;i++) {
782 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
783 if(info==NULL)continue;
784
785 varname = var_make_name(sys->system,info->i);
786 FPRINTF(ASCERR," %s",varname);
787 ASC_FREE(varname);
788 }
789 FPRINTF(ASCERR , "\nSet the \"%s\" flag on all but one of these to %s >= 0.\n"
790 , SCP(STATEFLAG),SCP(STATEFLAG)
791 );
792 error_reporter_end_flush();
793 return 0;
794 }else{
795 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
796 sys->x = info->i;
797 }
798 return 1;
799 }
800
801 /*------------------------------------------------------------------------------
802 CLASSIFICATION OF VARIABLES (for ANALYSIS step)
803 */
804
805 #define INTEG_ADD_TO_LIST(info,TYPE,INDEX,VAR,VARINDX,LIST) \
806 info = ASC_NEW(struct Integ_var_t); \
807 if(info==NULL){ \
808 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory (INTEG_VAR_NEW)"); \
809 return; \
810 } \
811 info->type=TYPE; \
812 info->index=INDEX; \
813 info->i=VAR; \
814 info->derivative=NULL; \
815 info->derivative_of=NULL; \
816 if(VARINDX==NULL){ \
817 info->varindx = -1; \
818 }else{ \
819 info->varindx = *VARINDX; \
820 } \
821 gl_append_ptr(LIST,(void *)info); \
822 info = NULL
823
824 /**
825 In a DAE, it's either the (single) independent variable, or it's a
826 variable in the model.
827
828 I'm not sure what we should be doing with variables that are already
829 present as derivatives of other variables, I guess those ones need to be
830 removed from the list in a second pass?
831 */
832 void integrator_dae_classify_var(IntegratorSystem *sys
833 , struct var_variable *var, const int *varindx
834 ){
835 struct Integ_var_t *info;
836 long type,index;
837
838 /* filter for recognition of solver_vars */
839 var_filter_t vfilt;
840 vfilt.matchbits = VAR_SVAR;
841 vfilt.matchvalue = VAR_SVAR;
842
843 assert(var != NULL && var_instance(var)!=NULL );
844
845 if( var_apply_filter(var,&vfilt) ) {
846 if(!var_active(var)){
847 CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");
848 return;
849 }
850
851 /* only non-fixed variables are accepted */
852 if(!var_fixed(var)){
853 /* get the ode_type and ode_id of this solver_var */
854 type = DynamicVarInfo(var,&index);
855
856 if(type==INTEG_OTHER_VAR){
857 /* if the var's type is -1, it's independent */
858 INTEG_ADD_TO_LIST(info,INTEG_OTHER_VAR,0,var,varindx,sys->indepvars);
859 }else{
860 if(type < 0)type=0;
861 /* any other type of var is in the DAE system, at least for now */
862 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
863 }
864 }
865 #if 0
866 else{
867 /* fixed variable, only include it if ode_type == 1 */
868 type = DynamicVarInfo(var,&index);
869 if(type==INTEG_STATE_VAR){
870 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
871 }
872 }
873 #endif
874
875 /* if the var's obs_id > 0, add it to the observation list */
876 if(ObservationVar(var,&index) != NULL && index > 0L) {
877 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->obslist);
878 }
879 }
880 }
881
882 /**
883 Inspect a specific variable and work out what type it is (what 'ode_type' it
884 has) and what other variable(s) it corresponds to (ie dydt corresponds to
885 y as a derivative).
886
887 @TODO add ability to create new variables for 'missing' derivative vars?
888 */
889 void integrator_ode_classify_var(IntegratorSystem *sys, struct var_variable *var
890 , const int *varindx
891 ){
892 struct Integ_var_t *info;
893 long type,index;
894
895 var_filter_t vfilt;
896 vfilt.matchbits = VAR_SVAR;
897 vfilt.matchvalue = VAR_SVAR;
898
899 assert(var != NULL && var_instance(var)!=NULL );
900
901 if( var_apply_filter(var,&vfilt) ) {
902 /* it's a solver var: what type of variable? */
903 type = DynamicVarInfo(var,&index);
904
905 if(type==INTEG_ALGEBRAIC_VAR){
906 /* no action required */
907 }else if(type==INTEG_OTHER_VAR){
908 /* i.e. independent var */
909 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
910 }else if(type>=INTEG_STATE_VAR){
911 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
912 if(type == 1){
913 sys->nstates++;
914 }else if(type == 2){ /* what about higher-order derivatives? -- JP */
915 sys->nderivs++;
916 }else{
917 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Higher-order (>=2) derivatives are not supported in ODEs.");
918 } }
919
920 if(ObservationVar(var,&index) != NULL && index > 0L) {
921 INTEG_ADD_TO_LIST(info,0L,index,var,varindx,sys->obslist);
922 }
923 }
924 }
925
926 /**
927 Look at a variable and determine if it's the independent variable or not.
928 This is just for the purpose of the integrator_find_indep_var function,
929 which is a utility function provided for use by the GUI.
930 */
931 void integrator_classify_indep_var(IntegratorSystem *sys
932 , struct var_variable *var, const int *varindx
933 ){
934 struct Integ_var_t *info;
935 long type,index;
936
937 var_filter_t vfilt;
938 vfilt.matchbits = VAR_SVAR;
939 vfilt.matchvalue = VAR_SVAR;
940
941 /* CONSOLE_DEBUG("..."); */
942
943 assert(var != NULL && var_instance(var)!=NULL );
944
945 if( var_apply_filter(var,&vfilt) ) {
946 type = DynamicVarInfo(var,&index);
947
948 if(type==INTEG_OTHER_VAR){
949 /* i.e. independent var */
950 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
951 }
952 }
953 }
954
955 /**
956 Look at a variable, and if it is an 'ODE variable' (it has a child instance
957 named 'ode_type') return its type, which will be either:
958 - INTEG_OTHER_VAR (if 'ode_type' is -1)
959 - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
960 - INTEG_STATE_VAR (if 'ode_type' is 1)
961 - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
962
963 If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
964 there.
965 */
966 static long DynamicVarInfo(struct var_variable *v,long *index){
967 struct Instance *c, *d, *i;
968 i = var_instance(v);
969 assert(i!=NULL);
970 assert(STATEFLAG!=NULL);
971 assert(STATEINDEX!=NULL);
972 c = ChildByChar(i,STATEFLAG);
973 d = ChildByChar(i,STATEINDEX);
974 /* lazy evaluation is important in the following if */
975 if(c == NULL
976 || d == NULL
977 || InstanceKind(c) != INTEGER_INST
978 || InstanceKind(d) != INTEGER_INST
979 || !AtomAssigned(c)
980 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
981 ){
982 return INTEG_ALGEBRAIC_VAR;
983 }
984 if (index != NULL) {
985 *index = GetIntegerAtomValue(d);
986 }
987 return GetIntegerAtomValue(c);
988 }
989
990 /**
991 Looks at the given variable checks if it is an 'observation variable'. This
992 means that it has its 'obs_id' child instance set to a non-zero value.
993
994 If the variable is an observation variable, its index value ('obs_id') is
995 stuff into *index (provided index!=NULL), and the pointer to the original
996 instance is rtruend.
997
998 If it's not an observation variable, we return NULL and *index is untouched.
999 */
1000 static struct var_variable *ObservationVar(struct var_variable *v, long *index){
1001 struct Instance *c,*i;
1002 i = var_instance(v);
1003 assert(i!=NULL);
1004 c = ChildByChar(i,OBSINDEX);
1005 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1006 return NULL;
1007 }
1008 if (index != NULL) {
1009 *index = GetIntegerAtomValue(c);
1010 }
1011 return v;
1012 }
1013
1014 /*------------------------------------------------------------------------------
1015 RUNNING THE SOLVER
1016 */
1017
1018 /*
1019 Make the call to the actual integrator we've selected, for the range of
1020 time values specified. The sys contains all the specifics.
1021
1022 Return 1 on success
1023 */
1024 int integrator_solve(IntegratorSystem *sys, long i0, long i1){
1025
1026 long nstep;
1027 unsigned long start_index=0, finish_index=0;
1028 assert(sys!=NULL);
1029
1030 nstep = integrator_getnsamples(sys)-1;
1031 /* check for at least 2 steps and dimensionality of x vs steps here */
1032
1033 if (i0<0 || i1 <0) {
1034 /* dude, there's no way we're writing interactive stuff here... */
1035 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
1036 return 0;
1037 } else {
1038 start_index=i0;
1039 finish_index =i1;
1040 if (start_index >= (unsigned long)nstep) {
1041 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
1042 ,start_index,nstep
1043 );
1044 return 0;
1045 }
1046 if (finish_index > (unsigned long)nstep) {
1047 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
1048 ,finish_index,nstep
1049 );
1050 return 0;
1051 }
1052 }
1053
1054 if(finish_index <= start_index) {
1055 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
1056 ,start_index,finish_index
1057 );
1058 return 0;
1059 }
1060
1061 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1062
1063 /* now go and run the integrator */
1064 switch (sys->engine) {
1065 case INTEG_LSODE:
1066 return integrator_lsode_solve(sys, start_index, finish_index);
1067 break;
1068 #ifdef ASC_WITH_IDA
1069 case INTEG_IDA:
1070 return integrator_ida_solve(sys,start_index, finish_index);
1071 break;
1072 #endif
1073 default:
1074 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator (invalid, or not implemented yet)");
1075 return 0;
1076 }
1077 }
1078
1079 /*---------------------------------------------------------------
1080 HANDLING THE LIST OF TIMESTEMPS
1081 */
1082
1083 #define GETTER_AND_SETTER(TYPE,NAME) \
1084 void integrator_set_##NAME(IntegratorSystem *sys, TYPE val){ \
1085 sys->NAME=val; \
1086 } \
1087 TYPE integrator_get_##NAME(IntegratorSystem *sys){ \
1088 return sys->NAME; \
1089 }
1090
1091 GETTER_AND_SETTER(SampleList *,samples) /*;*/
1092 GETTER_AND_SETTER(double,maxstep) /*;*/
1093 GETTER_AND_SETTER(double,minstep) /*;*/
1094 GETTER_AND_SETTER(double,stepzero) /*;*/
1095 GETTER_AND_SETTER(int,maxsubsteps) /*;*/
1096 #undef GETTER_AND_SETTER
1097
1098 long integrator_getnsamples(IntegratorSystem *sys){
1099 assert(sys!=NULL);
1100 assert(sys->samples!=NULL);
1101 return samplelist_length(sys->samples);
1102 }
1103
1104 double integrator_getsample(IntegratorSystem *sys, long i){
1105 assert(sys!=NULL);
1106 assert(sys->samples!=NULL);
1107 return samplelist_get(sys->samples,i);
1108 }
1109
1110 void integrator_setsample(IntegratorSystem *sys, long i,double xi){
1111 assert(sys!=NULL);
1112 assert(sys->samples!=NULL);
1113 samplelist_set(sys->samples,i,xi);
1114 }
1115
1116 const dim_type *integrator_getsampledim(IntegratorSystem *sys){
1117 assert(sys!=NULL);
1118 assert(sys->samples!=NULL);
1119 return samplelist_dim(sys->samples);
1120 }
1121
1122 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *sys){
1123 return sys->currentstep;
1124 }
1125
1126 /*------------------------------------------------------------------------------
1127 GET/SET VALUE OF THE INDEP VARIABLE
1128 */
1129
1130 /**
1131 Retrieve the value of the independent variable (time) from ASCEND
1132 and return it as a double.
1133 */
1134 double integrator_get_t(IntegratorSystem *sys){
1135 assert(sys->x!=NULL);
1136 return var_value(sys->x);
1137 }
1138
1139 /**
1140 Set the value of the independent variable (time) in ASCEND.
1141 */
1142 void integrator_set_t(IntegratorSystem *sys, double value){
1143 var_set_value(sys->x, value);
1144 /* CONSOLE_DEBUG("set_t = %g", value); */
1145 }
1146
1147 /*------------------------------------------------------------------------------
1148 PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1149 */
1150 /**
1151 Retrieve the current values of the derivatives of the y-variables
1152 and stick them in the/an array that the integrator will use.
1153
1154 If the pointer 'y' is NULL, the necessary space is allocated (and
1155 must be freed somewhere else).
1156 */
1157 double *integrator_get_y(IntegratorSystem *sys, double *y) {
1158 long i;
1159
1160 if (y==NULL) {
1161 y = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1162 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1163 }
1164
1165 for (i=0; i< sys->n_y; i++) {
1166 assert(sys->y[i]!=NULL);
1167 y[i] = var_value(sys->y[i]);
1168 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1169 }
1170 return y;
1171 }
1172
1173 /**
1174 Take the values of the derivatives from the array that the integrator
1175 uses, and use them to update the values of the corresponding variables
1176 in ASCEND.
1177 */
1178 void integrator_set_y(IntegratorSystem *sys, double *y) {
1179 long i;
1180 #ifndef NDEBUG
1181 char *varname;
1182 #endif
1183
1184 for (i=0; i < sys->n_y; i++) {
1185 assert(sys->y[i]!=NULL);
1186 var_set_value(sys->y[i],y[i]);
1187 #ifndef NDEBUG
1188 varname = var_make_name(sys->system, sys->y[i]);
1189 /* CONSOLE_DEBUG("y[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, y[i]); */
1190 ASC_FREE(varname);
1191 #endif
1192 }
1193 }
1194
1195 /**
1196 Send the values of the derivatives of the 'y' variables to the solver.
1197 Allocate space for an array if necessary.
1198
1199 Any element in sys->ydot that is NULL will be passed over (the value
1200 won't be modified in dydx).
1201 */
1202 double *integrator_get_ydot(IntegratorSystem *sys, double *dydx) {
1203 long i;
1204
1205 if (dydx==NULL) {
1206 dydx = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1207 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1208 }
1209
1210 for (i=0; i < sys->n_y; i++) {
1211 if(sys->ydot[i]!=NULL){
1212 dydx[i] = var_value(sys->ydot[i]);
1213 }
1214 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1215 }
1216 return dydx;
1217 }
1218
1219 void integrator_set_ydot(IntegratorSystem *sys, double *dydx) {
1220 long i;
1221 #ifndef NDEBUG
1222 /* char *varname; */
1223 #endif
1224 for (i=0; i < sys->n_y; i++) {
1225 if(sys->ydot[i]!=NULL){
1226 var_set_value(sys->ydot[i],dydx[i]);
1227 #ifndef NDEBUG
1228 /* varname = var_make_name(sys->system, sys->ydot[i]);
1229 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1230 ASC_FREE(varname); */
1231 #endif
1232 }
1233 #ifndef NDEBUG
1234 /*else{
1235 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1236 }*/
1237 #endif
1238 }
1239 }
1240
1241 /**
1242 Retrieve the values of 'ode_atol' properties of each of y-variables,
1243 for use in setting absolute error tolerances for the Integrator.
1244
1245 If the pointer 'atol' is NULL, the necessary space is allocated (and
1246 must be freed somewhere else).
1247 */
1248 double *integrator_get_atol(IntegratorSystem *sys, double *atol){
1249 long i;
1250 char *varname;
1251
1252 if (atol==NULL) {
1253 atol = ASC_NEW_ARRAY_CLEAR(double, sys->n_y);
1254 }
1255
1256 for (i=0; i< sys->n_y; i++) {
1257 assert(sys->y[i]!=NULL);
1258 atol[i] = var_odeatol(sys->y[i]);
1259 assert(atol[i]!=-1);
1260 varname = var_make_name(sys->system,sys->y[i]);
1261 CONSOLE_DEBUG("%s.ode_atol = %8.2e",varname,atol[i]);
1262 ASC_FREE(varname);
1263 }
1264 return atol;
1265 }
1266
1267 /*-------------------------------------------------------------
1268 RETRIEVING OBSERVATION DATA
1269 */
1270
1271 /**
1272 This function takes the inst in the solver and returns the vector of
1273 observation variables that are located in the submodel d.obs array.
1274 */
1275 double *integrator_get_observations(IntegratorSystem *sys, double *obsi) {
1276 long i;
1277
1278 if (obsi==NULL) {
1279 obsi = ASC_NEW_ARRAY_CLEAR(double, sys->n_obs+1);
1280 }
1281
1282 /* C obsi[0] <==> ascend d.obs[1] */
1283
1284 for (i=0; i < sys->n_obs; i++) {
1285 obsi[i] = var_value(sys->obs[i]);
1286 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1287 }
1288 return obsi;
1289 }
1290
1291 struct var_variable *integrator_get_observed_var(IntegratorSystem *sys, const long i){
1292 assert(i>=0);
1293 assert(i<sys->n_obs);
1294 return sys->obs[i];
1295 }
1296
1297 /**
1298 @NOTE Although this shouldn't be required for implementation of solver
1299 engines, this is useful for GUI reporting of integration results.
1300 */
1301 struct var_variable *integrator_get_independent_var(IntegratorSystem *sys){
1302 return sys->x;
1303 }
1304
1305
1306 /*----------------------------------------------------
1307 Build an analytic jacobian for solving the state system
1308
1309 This necessarily ugly piece of code attempts to create a unique
1310 list of relations that explicitly contain the variables in the
1311 given input list. The utility of this information is that we know
1312 exactly which relations must be differentiated, to fill in the
1313 df/dy matrix. If the problem has very few derivative terms, this will
1314 be of great savings. If the problem arose from the discretization of
1315 a pde, then this will be not so useful. The decision wether to use
1316 this function or to simply differentiate the entire relations list
1317 must be done before calling this function.
1318
1319 Final Note: the callee owns the array, but not the array elements.
1320 */
1321 #define AVG_NUM_INCIDENT 4
1322
1323
1324 /**
1325 This function helps to arrange the observation variables in a sensible order.
1326 The 'obs_id' child instance of v, if present, is assigned the value of the
1327 given parameter 'index'.
1328 */
1329 static void Integ_SetObsId(struct var_variable *v, long index){
1330 struct Instance *c, *i;
1331 i = var_instance(v);
1332 assert(i!=NULL);
1333 c = ChildByChar(i,OBSINDEX);
1334 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1335 return;
1336 }
1337 SetIntegerAtomValue(c,index,0);
1338 }
1339
1340 /**
1341 Compares observation structs. NULLs should end up at far end.
1342 */
1343 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1344 if(v1 == NULL)return 1;
1345 if(v2 == NULL)return -1;
1346 if(v1->index > v2->index)return 1;
1347 if(v1->index == v2->index)return 0;
1348 return -1;
1349 }
1350
1351 /**
1352 Compares dynamic vars structs. NULLs should end up at far end.
1353 List should be sorted primarily by index and then by type, in order
1354 of increasing value of both.
1355 */
1356 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1357 if(v1 == NULL)return 1;
1358 if(v2 == NULL)return -1;
1359 if(v1->index > v2->index)return 1;
1360 if(v1->index != v2->index)return -1;
1361 if(v1->type > v2->type)return 1;
1362 return -1;
1363 }
1364 /*----------------------------
1365 Output handling to the GUI/interface.
1366 */
1367
1368 int integrator_set_reporter(IntegratorSystem *sys
1369 , IntegratorReporter *reporter
1370 ){
1371 assert(sys!=NULL);
1372 sys->reporter = reporter;
1373 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1374 return 1;
1375 }
1376
1377 int integrator_output_init(IntegratorSystem *sys){
1378 assert(sys!=NULL);
1379 assert(sys->reporter!=NULL);
1380 if(sys->reporter->init!=NULL){
1381 /* call the specified output function */
1382 return (*(sys->reporter->init))(sys);
1383 }
1384 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1385 return 1;
1386 }
1387
1388 int integrator_output_write(IntegratorSystem *sys){
1389 static int reported_already=0;
1390 assert(sys!=NULL);
1391 if(sys->reporter->write!=NULL){
1392 return (*(sys->reporter->write))(sys);
1393 }
1394 if(!reported_already){
1395 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1396 reported_already=1;
1397 }
1398 return 1;
1399 }
1400
1401 int integrator_output_write_obs(IntegratorSystem *sys){
1402 static int reported_already=0;
1403 assert(sys!=NULL);
1404 if(sys->reporter->write_obs!=NULL){
1405 return (*(sys->reporter->write_obs))(sys);
1406 }
1407 if(!reported_already){
1408 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1409 reported_already=1;
1410 }
1411 return 1;
1412 }
1413
1414 int integrator_output_close(IntegratorSystem *sys){
1415 assert(sys!=NULL);
1416 if(sys->reporter->close!=NULL){
1417 return (*(sys->reporter->close))(sys);
1418 }
1419 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1420 return 1;
1421 }
1422
1423 /**
1424 Decode status codes from the integrator, and output them via FPRINTF.
1425 */
1426 int integrator_checkstatus(slv_status_t status) {
1427 if (status.converged) {
1428 return 1;
1429 }
1430 if (status.diverged) {
1431 FPRINTF(stderr, "The derivative system did not converge. Integration will terminate.");
1432 return 0;
1433 }
1434 if (status.inconsistent) {
1435 FPRINTF(stderr, "A numerically inconsistent state was discovered while "
1436 "calculating derivatives. Integration will terminate.");
1437 return 0;
1438 }
1439 if (status.time_limit_exceeded) {
1440 FPRINTF(stderr, "The time limit was exceeded while calculating "
1441 "derivatives. Integration will terminate.");
1442 return 0;
1443 }
1444 if (status.iteration_limit_exceeded) {
1445 FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1446 "derivatives. Integration will terminate.");
1447 return 0;
1448 }
1449 if (status.panic) {
1450 FPRINTF(stderr, "The user patience limit was exceeded while "
1451 "calculating derivatives. Integration will terminate.");
1452 return 0;
1453 }
1454 return 0;
1455 }

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