/[ascend]/trunk/ascend/integrator/integrator.c
ViewVC logotype

Contents of /trunk/ascend/integrator/integrator.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2322 - (show annotations) (download) (as text)
Wed Dec 15 06:12:36 2010 UTC (12 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 43199 byte(s)
Refactoring ascMalloc.h, mem.h to belong in 'general' with goal of having no references in 'general' to functions in 'utilities'.
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
30 #include <ascend/utilities/ascPanic.h>
31 #include <ascend/general/ascMalloc.h>
32 #include <ascend/compiler/packages.h>
33
34 #include <ascend/system/slv_common.h>
35 #include <ascend/system/slv_stdcalls.h>
36 #include <ascend/system/block.h>
37
38 #include <ascend/solver/solver.h>
39
40 #include "integrator.h"
41 #include "samplelist.h"
42
43 #define ANALYSE_DEBUG
44 /* #define SOLVE_DEBUG */
45 /* #define CLASSIFY_DEBUG */
46 /* #define DESTROY_DEBUG */
47 /* #define ATOL_DEBUG */
48
49 /*------------------------------------------------------------------------------
50 The following names are of solver_var children or attributes
51 * we support (at least temporarily) to determine who is a state and
52 * who matching derivative.
53 * These should be supported directly in a future solveratominst.
54 */
55
56 static symchar *g_symbols[3];
57
58 #define STATEFLAG g_symbols[0]
59 /*
60 Integer child. 0= algebraic, 1 = state, 2 = derivative, 3 = 2nd deriv etc
61 independent variable is -1.
62 */
63 #define INTEG_OTHER_VAR -1L
64 #define INTEG_ALGEBRAIC_VAR 0L
65 #define INTEG_STATE_VAR 1L
66
67 #define STATEINDEX g_symbols[1]
68 /* Integer child. all variables with the same STATEINDEX value are taken to
69 * be derivatives of the same state variable. We really need a compiler
70 * that maintains this info by backpointers, but oh well until then.
71 */
72 #define OBSINDEX g_symbols[2]
73 /* Integer child. All variables with OBSINDEX !=0 will be sent to the
74 IntegratorOutputWriteObsFn allowing output to a file, graph, console, etc.
75 */
76
77
78 /** Temporary catcher of dynamic variable and observation variable data */
79 struct Integ_var_t {
80 long index;
81 long type;
82 struct Integ_var_t *derivative;
83 struct Integ_var_t *derivative_of;
84 struct var_variable *i;
85 int varindx; /**< index into slv_get_master_vars_list, or -1 if not there */
86 int isstate;
87 };
88
89 /*------------------------------------------------------------------------------
90 forward declarations
91 */
92
93 /* abstractions of setup/teardown procedures for the specific solvers */
94 void integrator_create_engine(IntegratorSystem *sys);
95 void integrator_free_engine(IntegratorSystem *sys);
96
97 IntegratorAnalyseFn integrator_analyse_ode;
98
99 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
100 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
101 IntegratorVarVisitorFn integrator_ode_classify_var;
102 IntegratorVarVisitorFn integrator_dae_classify_var;
103 IntegratorVarVisitorFn integrator_classify_indep_var;
104 IntegratorVarVisitorFn integrator_dae_show_var;
105
106 static int integrator_sort_obs_vars(IntegratorSystem *sys);
107 static void integrator_print_var_stats(IntegratorSystem *sys);
108 static int integrator_check_indep_var(IntegratorSystem *sys);
109
110 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2);
111 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2);
112 static void Integ_SetObsId(struct var_variable *v, long oindex);
113
114 static long DynamicVarInfo(struct var_variable *v,long *vindex);
115 static struct var_variable *ObservationVar(struct var_variable *v, long *oindex);
116 static void IntegInitSymbols(void);
117
118 /*------------------------------------------------------------------------------
119 INSTANTIATION AND DESTRUCTION
120 */
121
122 /**
123 Create a new IntegratorSystem and assign a slv_system_t to it.
124 */
125 IntegratorSystem *integrator_new(slv_system_t slvsys, struct Instance *inst){
126 IntegratorSystem *sys;
127
128 if (slvsys == NULL) {
129 ERROR_REPORTER_HERE(ASC_PROG_ERR,"slvsys is NULL!");
130 return NULL;
131 }
132
133 sys = ASC_NEW_CLEAR(IntegratorSystem);
134 sys->system = slvsys;
135 sys->instance = inst;
136
137 sys->engine = INTEG_UNKNOWN;
138 sys->internals = NULL;
139
140 sys->states = NULL; sys->derivs = NULL;
141 sys->dynvars = NULL; sys->obslist = NULL; sys->indepvars = NULL;
142
143 sys->y_id = NULL;
144 sys->obs_id = NULL;
145 sys->y = NULL;
146 sys->ydot = NULL;
147 sys->obs = NULL;
148 sys->n_y = 0;
149 return sys;
150 }
151
152 /**
153 Carefully trash any data in the IntegratorSystem that we own,
154 then destroy the IntegratorSystem struct.
155
156 Note that the integrator doesn't own the samplelist.
157
158 @param sys will be destroyed and set to NULL.
159 */
160 void integrator_free(IntegratorSystem *sys){
161 if(sys==NULL)return;
162
163 integrator_free_engine(sys);
164
165 /* CONSOLE_DEBUG("Engine freed, destroying internal data"); */
166
167 if(sys->states != NULL)gl_destroy(sys->states);
168 if(sys->derivs != NULL)gl_destroy(sys->derivs);
169
170 if(sys->dynvars != NULL)gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
171 if(sys->obslist != NULL)gl_free_and_destroy(sys->obslist); /* and obslist */
172 if (sys->indepvars != NULL)gl_free_and_destroy(sys->indepvars); /* and indepvars */
173
174 if(sys->y_id != NULL)ASC_FREE(sys->y_id);
175 if(sys->obs_id != NULL)ASC_FREE(sys->obs_id);
176 if(sys->y != NULL)ASC_FREE(sys->y);
177 if(sys->ydot != NULL)ASC_FREE(sys->ydot);
178 if(sys->obs != NULL)ASC_FREE(sys->obs);
179
180 slv_destroy_parms(&(sys->params));
181
182 ASC_FREE(sys);
183 #ifdef DESTROY_DEBUG
184 CONSOLE_DEBUG("Destroyed IntegratorSystem");
185 #endif
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 /**
204 Local function that holds the list of available integrators. The value
205 returned is NOT owned by the called.
206
207 @param free_space if 0, call as normal. if 1, free the list and maybe do some
208 cleaning up etc. Should be called whenever a simulation is destroyed.
209 */
210 static struct gl_list_t *integrator_get_list(int free_space){
211 static int init = 0;
212 static struct gl_list_t *L;
213 int i, error;
214 /* standard integrators that we will register */
215 static char *defaultintegrators[] = {
216 "lsode"
217 ,"ida"
218 ,NULL
219 };
220
221 if(free_space){
222 if(init && L)ASC_FREE(L);
223 init = 0;
224 return NULL;
225 }
226 if(!init){
227 L = gl_create(10);
228 init = 1; /* set init to 1 here now, since this will be called
229 recursively from the LoadArchiveLibrary call. */
230
231 /* CONSOLE_DEBUG("REGISTERING STANDARD SOLVER ENGINES"); */
232 for(i=0; defaultintegrators[i]!=NULL;++i){
233 error = package_load(defaultintegrators[i],NULL);
234 if(error){
235 ERROR_REPORTER_HERE(ASC_PROG_ERR
236 ,"Unable to register integrator '%s' (error %d)."
237 ,defaultintegrators[i],error
238 );
239 }else{
240 CONSOLE_DEBUG("Integrator '%s' registered OK",defaultintegrators[i]);
241 }
242 }
243 }
244 return L;
245 }
246
247 /**
248 Return gl_list of IntegratorInternals. C++ will use this to produce a
249 nice little list of integrator names that can be used in Python :-/
250 */
251 const struct gl_list_t *integrator_get_engines(){
252 return integrator_get_list(0);
253 }
254
255 struct gl_list_t *integrator_get_engines_growable(){
256 return integrator_get_list(0);
257 }
258
259 /* return 0 on success */
260 int integrator_set_engine(IntegratorSystem *sys, const char *name){
261 struct gl_list_t *L = integrator_get_engines_growable();
262 int i;
263 const IntegratorInternals *I, *Ifound=NULL;
264 for(i=1; i <= gl_length(L); ++i){
265 I = gl_fetch(L,i);
266 if(strcmp(I->name,name)==0){
267 Ifound = I;
268 break;
269 }
270 }
271 if(Ifound){
272 /** @TODO tests for applicability of this engine... */
273
274 CONSOLE_DEBUG("Setting engine...");
275 if(Ifound->engine == sys->engine){
276 // already set...
277 return 0;
278 }
279
280 // clean up if old engine is there
281 if(sys->engine!=INTEG_UNKNOWN){
282 #ifdef DESTROY_DEBUG
283 CONSOLE_DEBUG("Freeing memory used by old integrator engine");
284 #endif
285 integrator_free_engine(sys);
286 #ifdef DESTROY_DEBUG
287 CONSOLE_DEBUG("done");
288 #endif
289 }
290
291 sys->engine = Ifound->engine;
292
293 sys->internals = Ifound;
294 }else{
295 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid engine name '%s'",name);
296 return 1;
297 }
298
299 asc_assert(sys->internals);
300 integrator_create_engine(sys);
301 return 0;
302 }
303
304 /**
305 @TODO rename this
306 */
307 const IntegratorInternals *integrator_get_engine(const IntegratorSystem *sys){
308 return sys->internals;
309 }
310
311 /**
312 Free any engine-specific data that was required for the solution of
313 this system. Note that this data is pointed to by sys->enginedata.
314
315 @TODO rename this
316 */
317 void integrator_free_engine(IntegratorSystem *sys){
318 if(sys->engine==INTEG_UNKNOWN)return;
319 if(sys->enginedata){
320 if(sys->internals){
321 #ifdef DESTROY_DEBUG
322 CONSOLE_DEBUG("Running engine's freefn");
323 #endif
324 (sys->internals->freefn)(sys->enginedata);
325 #ifdef DESTROY_DEBUG
326 CONSOLE_DEBUG("Done with freefn");
327 #endif
328 sys->enginedata=NULL;
329 }else{
330 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to free engine data: no sys->internals");
331 }
332 }
333 #ifdef DESTROY_DEBUG
334 CONSOLE_DEBUG("Done with integrator_free_engine");
335 #endif
336 }
337
338 /**
339 Create enginedata memory if required for this solver. This doesn't include
340 allocating computation space, since we assume that this stage all we know
341 is that we want to use a specified integrator engine, not the full details
342 of the problem at hand. Allocating space inside enginedata should be done
343 during the solve stage (and freed inside integrator_free_engine)
344
345 @TODO rename this
346 */
347 void integrator_create_engine(IntegratorSystem *sys){
348 asc_assert(sys);
349 asc_assert(sys->engine!=INTEG_UNKNOWN);
350 asc_assert(sys->internals);
351 asc_assert(sys->internals->createfn);
352 asc_assert(sys->enginedata==NULL);
353 (sys->internals->createfn)(sys);
354 }
355
356 int integrator_register(const IntegratorInternals *integ){
357 /* get the current list of registered engines */
358 struct gl_list_t *L;
359 L = integrator_get_engines_growable();
360
361 CONSOLE_DEBUG("REGISTERING INTEGRATOR");
362 CONSOLE_DEBUG("There were %lu registered integrators", gl_length(integrator_get_list(0)));
363
364 int i;
365 IntegratorInternals *I;
366 for(i=1; i <= gl_length(L); ++i){
367 I = (IntegratorInternals *)gl_fetch(L,i);
368 if(strcmp(integ->name,I->name)==0){
369 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Integrator with name '%s' is already registered",integ->name);
370 return 0;
371 }
372 if(integ->engine == I->engine){
373 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Integrator with ID '%d' is already registered",integ->engine);
374 return 0;
375 }
376 }
377
378 CONSOLE_DEBUG("Adding engine '%s'",integ->name);
379
380 gl_append_ptr(L,(void *)integ);
381
382 CONSOLE_DEBUG("There are now %lu registered integrators", gl_length(integrator_get_list(0)));
383 return 0;
384 }
385
386 /*------------------------------------------------------------------------------
387 PARAMETERS
388 */
389
390 /**
391 Reset the parameters in this IntegratorSystem to the default ones for this
392 Integrator.
393
394 @return 0 on success, 1 on error
395 */
396 int integrator_params_default(IntegratorSystem *sys){
397 asc_assert(sys->engine!=INTEG_UNKNOWN);
398 asc_assert(sys->internals);
399 return (sys->internals->paramsdefaultfn)(sys);
400 }
401
402 int integrator_params_get(const IntegratorSystem *sys, slv_parameters_t *parameters){
403 asc_assert(sys!=NULL);
404 asc_assert(sys->params.num_parms > 0);
405 memcpy(parameters,&(sys->params),sizeof(slv_parameters_t));
406 return 0;
407 }
408
409 int integrator_params_set(IntegratorSystem *sys, const slv_parameters_t *parameters){
410 asc_assert(sys!=NULL);
411 asc_assert(parameters!=NULL);
412 memcpy(&(sys->params),parameters,sizeof(slv_parameters_t));
413 return 0;
414 }
415
416 /*------------------------------------------------------------------------------
417 ANALYSIS
418
419 Provide two modes in order to provide analysis suitable for solution of both
420 ODEs (the previous technique) and DAEs (new code). These share a common
421 "visit" method that needs to eventually be integrated with the code in
422 <solver/analyze.c>. For the moment, we're just hacking in to the compiler.
423 */
424
425 /**
426 Locate the independent variable.
427
428 @NOTE For the purpose of GUI design, this needs to work independent of the
429 integration engine being used. @ENDNOTE
430 */
431 int integrator_find_indep_var(IntegratorSystem *sys){
432 int result = 0;
433 #ifdef ANALYSE_DEBUG
434 char *varname;
435 #endif
436
437 /* if the indep var has been found, we don't look again (we assume the user won't fiddle with ode_type!) */
438 if(sys->x != NULL){
439 CONSOLE_DEBUG("sys->x already set");
440 return 0; /* success */
441 }
442
443 /* create a clear indepvars list */
444 if(sys->indepvars!=NULL){
445 ERROR_REPORTER_HERE(ASC_PROG_ERR,"indepvars should be NULL at this point");
446 return 1; /* error */
447 }
448 sys->indepvars = gl_create(10L);
449
450 IntegInitSymbols();
451
452 CONSOLE_DEBUG("Looking for independent var...");
453 integrator_visit_system_vars(sys,&integrator_classify_indep_var);
454 #ifdef ANALYSE_DEBUG
455 if(gl_length(sys->indepvars)){
456 CONSOLE_DEBUG("Found %lu indepvars",gl_length(sys->indepvars));
457 }else{
458 CONSOLE_DEBUG("NO INDEP VARS FOUND");
459 }
460 #endif
461
462 /* after visiting the instance tree, look at the candidates and return 0 on success */
463 result = integrator_check_indep_var(sys);
464
465 /* whatever happens, we clean up afterwards */
466 gl_free_and_destroy(sys->indepvars);
467 sys->indepvars = NULL;
468
469 #ifdef ANALYSE_DEBUG
470 asc_assert(sys->system);
471 if(!result){
472 asc_assert(sys->x);
473 varname = var_make_name(sys->system, sys->x);
474 CONSOLE_DEBUG("Indep var is '%s'",varname);
475 ASC_FREE(varname);
476 }else{
477 CONSOLE_DEBUG("No indep var was found");
478 }
479 #endif
480
481 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Returning result %d",result); */
482
483 return result;
484 }
485
486 /**
487 Perform whatever additional problem is required so that the system can be
488 integrated as a dynamical system with the IntegrationEngine chosen.
489
490 We can always assume that sys->system has had analyse_make_problem called on
491 it, so all the variable lists (etc) will be there already.
492
493 @return 0 on success
494 */
495 int integrator_analyse(IntegratorSystem *sys){
496 int res;
497
498 #ifdef ANALYSE_DEBUG
499 CONSOLE_DEBUG("Analysing integration system...");
500 #endif
501 asc_assert(sys);
502 if(sys->engine==INTEG_UNKNOWN){
503 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No engine selected: can't analyse");
504 return 1;
505 }
506 asc_assert(sys->engine!=INTEG_UNKNOWN);
507 asc_assert(sys->internals);
508
509 if(!sys->indepvars || !gl_length(sys->indepvars)){
510 if(integrator_find_indep_var(sys)){
511 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Independent variable problem: abandoning integration");
512 return 2;
513 }
514 #ifdef ANALYSE_DEBUG
515 else{
516 CONSOLE_DEBUG("got 0 from integrator_find_indep_var");
517 }
518 #endif
519 }
520
521 res = (sys->internals->analysefn)(sys);
522 #ifdef ANALYSE_DEBUG
523 CONSOLE_DEBUG("integrator_analyse returning %d",res);
524 #endif
525 return res;
526 }
527
528
529 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
530 struct var_variable **vlist;
531 int i, vlen;
532
533 /* visit all the slv_system_t master var lists to collect vars */
534 /* find the vars mostly in this one */
535 vlist = slv_get_solvers_var_list(sys->system);
536 vlen = slv_get_num_solvers_vars(sys->system);
537 for (i=0;i<vlen;i++) {
538 (*visitfn)(sys, vlist[i], &i);
539 }
540
541 /*
542 CONSOLE_DEBUG("Checked %d vars",vlen);
543 integrator_print_var_stats(sys);
544 */
545
546 /* probably nothing here, but gotta check. */
547 vlist = slv_get_master_par_list(sys->system);
548 vlen = slv_get_num_master_pars(sys->system);
549 for (i=0;i<vlen;i++) {
550 (*visitfn)(sys, vlist[i], NULL);
551 }
552
553 /*
554 CONSOLE_DEBUG("Checked %d pars",vlen);
555 integrator_print_var_stats(sys);
556 */
557
558 /* might find t here */
559 vlist = slv_get_master_unattached_list(sys->system);
560 vlen = slv_get_num_master_unattached(sys->system);
561 for (i=0;i<vlen;i++) {
562 (*visitfn)(sys, vlist[i], NULL);
563 }
564
565 /* CONSOLE_DEBUG("Checked %d unattached",vlen); */
566 }
567 /**
568 Analyse the ODE structure. We can assume that the independent variable was
569 already found.
570
571 @return 0 on success
572 */
573 int integrator_analyse_ode(IntegratorSystem *sys){
574 struct Integ_var_t *v1,*v2;
575 long half,i,len;
576 int happy=1;
577 char *varname1, *varname2;
578
579 asc_assert(sys->system!=NULL);
580
581 if(strcmp(slv_solver_name(slv_get_selected_solver(sys->system)),"QRSlv")!=0){
582 ERROR_REPORTER_HERE(ASC_PROG_ERR,"System must have solver 'QRSlv' assigned to it before integration");
583 return 2;
584 }
585 CONSOLE_DEBUG("Checked that NLA solver is set to '%s'",slv_solver_name(slv_get_selected_solver(sys->system)));
586
587 CONSOLE_DEBUG("Starting ODE analysis");
588 IntegInitSymbols();
589
590 /* collect potential states and derivatives */
591 sys->indepvars = gl_create(10L); /* t var info */
592 sys->dynvars = gl_create(200L); /* y ydot var info */
593 sys->obslist = gl_create(100L); /* obs info */
594 if (sys->dynvars == NULL
595 || sys->obslist == NULL
596 || sys->indepvars == NULL
597 ){
598 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
599 return 1;
600 }
601
602 sys->nstates = sys->nderivs = 0;
603
604 integrator_visit_system_vars(sys,&integrator_ode_classify_var);
605
606 integrator_print_var_stats(sys);
607
608 /* check sanity of state and var lists */
609
610 len = gl_length(sys->dynvars);
611 half = len/2;
612 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
613
614 if (len % 2 || len == 0L || sys->nstates != sys->nderivs ) {
615 /* list length must be even for vars to pair off */
616 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
617 return 3;
618 }
619 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
620 if (gl_fetch(sys->dynvars,len)==NULL) {
621 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
622 return 4;
623 }
624 sys->states = gl_create(half); /* state vars Integ_var_t references */
625 sys->derivs = gl_create(half); /* derivative var atoms */
626 for (i=1;i < len; i+=2) {
627 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
628 v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
629 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
630 varname1 = var_make_name(sys->system,v1->i);
631 varname2 = var_make_name(sys->system,v2->i);
632
633 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
634 varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
635 varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
636 );
637 ASC_FREE(varname1);
638 ASC_FREE(varname2);
639 happy=0;
640 break;
641 } else {
642 gl_append_ptr(sys->states,(POINTER)v1);
643 gl_append_ptr(sys->derivs,(POINTER)v2->i);
644 }
645 }
646 if (!happy) {
647 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem with ode_id and ode_type values");
648 return 5;
649 }
650 sys->n_y = half;
651 sys->y = ASC_NEW_ARRAY(struct var_variable *, half);
652 sys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
653 sys->y_id = ASC_NEW_ARRAY(int, half);
654 if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
655 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
656 return 6;
657 }
658 for (i = 1; i <= half; i++) {
659 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
660 sys->y[i-1] = v1->i;
661 sys->y_id[i-1] = v1->index;
662 sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
663 }
664
665 if(integrator_sort_obs_vars(sys)){
666 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error sorting observation variables");
667 return 7;
668 }
669
670 /* FIX all states */
671 for(i=0; i<sys->n_y; ++i){
672 if(!var_fixed(sys->y[i])){
673 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Fixing state %d",i);
674 var_set_fixed(sys->y[i], TRUE);
675 }
676 }
677
678 /* don't need the gl_lists now that we have arrays for everyone */
679 gl_destroy(sys->states);
680 gl_destroy(sys->derivs);
681 gl_free_and_destroy(sys->indepvars); /* we own the objects in indepvars */
682 gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
683 gl_free_and_destroy(sys->obslist); /* and obslist */
684 sys->states = NULL;
685 sys->derivs = NULL;
686 sys->indepvars = NULL;
687 sys->dynvars = NULL;
688 sys->obslist = NULL;
689
690 /* analysis completed OK */
691 return 0;
692 }
693
694 /**
695 Reindex observations. Sort if the user mostly numbered. Take natural order
696 if user just booleaned.
697
698 @return 0 on success
699 */
700 static int integrator_sort_obs_vars(IntegratorSystem *sys){
701 int half, i, len = 0;
702 struct Integ_var_t *v2;
703
704 half = sys->n_y;
705 len = gl_length(sys->obslist);
706 /* we shouldn't be seeing NULL here ever except if malloc fail. */
707 if (len > 1L) {
708 half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;
709 /* half != 0 now because we didn't collect 0 indexed vars */
710 for (i=2; i <= len; i++) {
711 if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {
712 /* change seen. sort and go on */
713 gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);
714 break;
715 }
716 }
717 }
718 for (i = half = 1; i <= len; i++) {
719 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
720 if (v2==NULL) {
721 /* we shouldn't be seeing NULL here ever except if malloc fail. */
722 gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */
723 } else {
724 Integ_SetObsId(v2->i,half);
725 v2->index = half++;
726 }
727 }
728
729 /* obslist now uniquely indexed, no nulls */
730 /* make into arrays */
731 half = gl_length(sys->obslist);
732 sys->obs = ASC_NEW_ARRAY(struct var_variable *,half);
733 sys->obs_id = ASC_NEW_ARRAY(int, half);
734 if ( sys->obs==NULL || sys->obs_id==NULL) {
735 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
736 return 1;
737 }
738 sys->n_obs = half;
739 for (i = 1; i <= half; i++) {
740 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
741 sys->obs[i-1] = v2->i;
742 sys->obs_id[i-1] = v2->index;
743 }
744
745 return 0;
746 }
747
748 static void integrator_print_var_stats(IntegratorSystem *sys){
749 int v = gl_length(sys->dynvars);
750 int i = gl_length(sys->indepvars);
751 CONSOLE_DEBUG("Currently %d vars, %d indep",v,i);
752 }
753
754 /**
755 Check sanity of the independent variable.
756
757 @return 0 on success
758 */
759 static int integrator_check_indep_var(IntegratorSystem *sys){
760 int len, i;
761 struct Integ_var_t *info;
762 char *varname;
763
764 if(sys->x){
765 CONSOLE_DEBUG("Indep var already assigned");
766 return 1;
767 }
768
769 /* check the sanity of the independent variable */
770 len = gl_length(sys->indepvars);
771 if (!len) {
772 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
773 return 2;
774 }
775 if (len > 1) {
776 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
777 FPRINTF(ASCERR,"Excess %ld independent variables found:",
778 len);
779 for(i=1; i <=len;i++) {
780 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
781 if(info==NULL)continue;
782
783 varname = var_make_name(sys->system,info->i);
784 FPRINTF(ASCERR," %s",varname);
785 ASC_FREE(varname);
786 }
787 FPRINTF(ASCERR , "\nSet the '%s' flag on all but one of these to %s >= 0.\n"
788 , SCP(STATEFLAG),SCP(STATEFLAG)
789 );
790 error_reporter_end_flush();
791 return 3;
792 }else{
793 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
794 sys->x = info->i;
795 }
796 asc_assert(gl_length(sys->indepvars)==1);
797 asc_assert(sys->x);
798 return 0;
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 asc_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 1
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 asc_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 var_filter_t vfilt;
937 #ifdef CLASSIFY_DEBUG
938 char *varname;
939 #endif
940
941 asc_assert(var != NULL && var_instance(var)!=NULL );
942 vfilt.matchbits = VAR_SVAR;
943 vfilt.matchvalue = VAR_SVAR;
944
945 #ifdef CLASSIFY_DEBUG
946 varname = var_make_name(sys->system,var);
947 #endif
948
949 if( var_apply_filter(var,&vfilt) ) {
950 type = DynamicVarInfo(var,&index);
951
952 if(type==INTEG_OTHER_VAR){
953 /* i.e. independent var */
954 #ifdef CLASSIFY_DEBUG
955 CONSOLE_DEBUG("Var '%s' added to indepvars",varname);
956 #endif
957 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
958 #ifdef CLASSIFY_DEBUG
959 }else{
960 CONSOLE_DEBUG("Var '%s' not correct ode_type",varname);
961 #endif
962 }
963 #ifdef CLASSIFY_DEBUG
964 }else{
965 CONSOLE_DEBUG("Var '%s' failed filter",varname);
966 #endif
967 }
968
969 #ifdef CLASSIFY_DEBUG
970 ASC_FREE(varname);
971 #endif
972 }
973
974 /**
975 Look at a variable, and if it is an 'ODE variable' (it has a child instance
976 named 'ode_type') return its type, which will be either:
977 - INTEG_OTHER_VAR (if 'ode_type' is -1)
978 - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
979 - INTEG_STATE_VAR (if 'ode_type' is 1)
980 - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
981
982 If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
983 there.
984 */
985 static long DynamicVarInfo(struct var_variable *v,long *index){
986 struct Instance *c, *d, *i;
987 i = var_instance(v);
988 asc_assert(i!=NULL);
989 asc_assert(STATEFLAG!=NULL);
990 asc_assert(STATEINDEX!=NULL);
991 c = ChildByChar(i,STATEFLAG);
992 d = ChildByChar(i,STATEINDEX);
993 /* lazy evaluation is important in the following if */
994 if(c == NULL
995 || d == NULL
996 || InstanceKind(c) != INTEGER_INST
997 || InstanceKind(d) != INTEGER_INST
998 || !AtomAssigned(c)
999 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
1000 ){
1001 return INTEG_ALGEBRAIC_VAR;
1002 }
1003 if (index != NULL) {
1004 *index = GetIntegerAtomValue(d);
1005 }
1006 return GetIntegerAtomValue(c);
1007 }
1008
1009 /**
1010 Looks at the given variable checks if it is an 'observation variable'. This
1011 means that it has its 'obs_id' child instance set to a non-zero value.
1012
1013 If the variable is an observation variable, its index value ('obs_id') is
1014 stuff into *index (provided index!=NULL), and the pointer to the original
1015 instance is rtruend.
1016
1017 If it's not an observation variable, we return NULL and *index is untouched.
1018 */
1019 static struct var_variable *ObservationVar(struct var_variable *v, long *index){
1020 struct Instance *c,*i;
1021 i = var_instance(v);
1022 asc_assert(i!=NULL);
1023 c = ChildByChar(i,OBSINDEX);
1024 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1025 return NULL;
1026 }
1027 if (index != NULL) {
1028 *index = GetIntegerAtomValue(c);
1029 }
1030 return v;
1031 }
1032
1033 /*------------------------------------------------------------------------------
1034 RUNNING THE SOLVER
1035 */
1036
1037 /*
1038 Make the call to the actual integrator we've selected, for the range of
1039 time values specified. The sys contains all the specifics.
1040
1041 Return 0 on success
1042 */
1043 int integrator_solve(IntegratorSystem *sys, long i0, long i1){
1044
1045 long nstep;
1046 unsigned long start_index=0, finish_index=0;
1047 asc_assert(sys!=NULL);
1048
1049 asc_assert(sys->internals);
1050 asc_assert(sys->engine!=INTEG_UNKNOWN);
1051
1052 nstep = integrator_getnsamples(sys)-1;
1053 /* check for at least 2 steps and dimensionality of x vs steps here */
1054
1055 if (i0<0 || i1 <0) {
1056 /* removed completely inappropriate interactive code here */
1057 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
1058 return -1;
1059 } else {
1060 start_index=i0;
1061 finish_index =i1;
1062 if (start_index >= (unsigned long)nstep) {
1063 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
1064 ,start_index,nstep
1065 );
1066 return -2;
1067 }
1068 if (finish_index > (unsigned long)nstep) {
1069 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
1070 ,finish_index,nstep
1071 );
1072 return -3;
1073 }
1074 }
1075
1076 if(finish_index <= start_index) {
1077 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
1078 ,start_index,finish_index
1079 );
1080 return -4;
1081 }
1082
1083 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1084
1085 return (sys->internals->solvefn)(sys,start_index,finish_index);
1086 }
1087
1088 /*---------------------------------------------------------------
1089 HANDLING THE LIST OF TIMESTEMPS
1090 */
1091
1092 #define GETTER_AND_SETTER(TYPE,NAME) \
1093 void integrator_set_##NAME(IntegratorSystem *sys, TYPE val){ \
1094 sys->NAME=val; \
1095 } \
1096 TYPE integrator_get_##NAME(IntegratorSystem *sys){ \
1097 return sys->NAME; \
1098 }
1099
1100 GETTER_AND_SETTER(SampleList *,samples) /*;*/
1101 GETTER_AND_SETTER(double,maxstep) /*;*/
1102 GETTER_AND_SETTER(double,minstep) /*;*/
1103 GETTER_AND_SETTER(double,stepzero) /*;*/
1104 GETTER_AND_SETTER(int,maxsubsteps) /*;*/
1105 #undef GETTER_AND_SETTER
1106
1107 long integrator_getnsamples(IntegratorSystem *sys){
1108 asc_assert(sys!=NULL);
1109 asc_assert(sys->samples!=NULL);
1110 return samplelist_length(sys->samples);
1111 }
1112
1113 double integrator_getsample(IntegratorSystem *sys, long i){
1114 asc_assert(sys!=NULL);
1115 asc_assert(sys->samples!=NULL);
1116 return samplelist_get(sys->samples,i);
1117 }
1118
1119 void integrator_setsample(IntegratorSystem *sys, long i,double xi){
1120 asc_assert(sys!=NULL);
1121 asc_assert(sys->samples!=NULL);
1122 samplelist_set(sys->samples,i,xi);
1123 }
1124
1125 const dim_type *integrator_getsampledim(IntegratorSystem *sys){
1126 asc_assert(sys!=NULL);
1127 asc_assert(sys->samples!=NULL);
1128 return samplelist_dim(sys->samples);
1129 }
1130
1131 ASC_DLLSPEC long integrator_getcurrentstep(IntegratorSystem *sys){
1132 return sys->currentstep;
1133 }
1134
1135 /*------------------------------------------------------------------------------
1136 GET/SET VALUE OF THE INDEP VARIABLE
1137 */
1138
1139 /**
1140 Retrieve the value of the independent variable (time) from ASCEND
1141 and return it as a double.
1142 */
1143 double integrator_get_t(IntegratorSystem *sys){
1144 asc_assert(sys->x!=NULL);
1145 return var_value(sys->x);
1146 }
1147
1148 /**
1149 Set the value of the independent variable (time) in ASCEND.
1150 */
1151 void integrator_set_t(IntegratorSystem *sys, double value){
1152 var_set_value(sys->x, value);
1153 /* CONSOLE_DEBUG("set_t = %g", value); */
1154 }
1155
1156 /*------------------------------------------------------------------------------
1157 PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1158 */
1159 /**
1160 Retrieve the current values of the derivatives of the y-variables
1161 and stick them in the/an array that the integrator will use.
1162
1163 If the pointer 'y' is NULL, the necessary space is allocated (and
1164 must be freed somewhere else).
1165 */
1166 double *integrator_get_y(IntegratorSystem *sys, double *y) {
1167 long i;
1168
1169 if (y==NULL) {
1170 y = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1171 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1172 }
1173
1174 for (i=0; i< sys->n_y; i++) {
1175 asc_assert(sys->y[i]!=NULL);
1176 y[i] = var_value(sys->y[i]);
1177 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1178 }
1179 return y;
1180 }
1181
1182 /**
1183 Take the values of the differential variables from the array that the
1184 integrator uses, and use them to update the values of the corresponding
1185 variables in ASCEND.
1186 */
1187 void integrator_set_y(IntegratorSystem *sys, double *y) {
1188 long i;
1189 #ifdef SOLVE_DEBUG
1190 char *varname;
1191 #endif
1192
1193 for (i=0; i < sys->n_y; i++) {
1194 asc_assert(sys->y[i]!=NULL);
1195 var_set_value(sys->y[i],y[i]);
1196 #ifdef SOLVE_DEBUG
1197 varname = var_make_name(sys->system, sys->y[i]);
1198 CONSOLE_DEBUG("y[%ld] = %g --> '%s'", i+1, y[i], varname);
1199 ASC_FREE(varname);
1200 #endif
1201 }
1202 }
1203
1204 /**
1205 Send the values of the derivatives of the 'y' variables to the solver.
1206 Allocate space for an array if necessary.
1207
1208 Any element in sys->ydot that is NULL will be passed over (the value
1209 won't be modified in dydx).
1210 */
1211 double *integrator_get_ydot(IntegratorSystem *sys, double *dydx) {
1212 long i;
1213
1214 if (dydx==NULL) {
1215 dydx = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1216 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1217 }
1218
1219 for (i=0; i < sys->n_y; i++) {
1220 if(sys->ydot[i]!=NULL){
1221 dydx[i] = var_value(sys->ydot[i]);
1222 }
1223 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1224 }
1225 return dydx;
1226 }
1227
1228 void integrator_set_ydot(IntegratorSystem *sys, double *dydx) {
1229 long i;
1230 #ifdef SOLVE_DEBUG
1231 char *varname;
1232 #endif
1233 for (i=0; i < sys->n_y; i++) {
1234 if(sys->ydot[i]!=NULL){
1235 var_set_value(sys->ydot[i],dydx[i]);
1236 #ifdef SOLVE_DEBUG
1237 varname = var_make_name(sys->system, sys->ydot[i]);
1238 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1239 ASC_FREE(varname);
1240 }else{
1241 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1242 #endif
1243 }
1244 }
1245 }
1246
1247 /**
1248 Retrieve the values of 'ode_atol' properties of each of y-variables,
1249 for use in setting absolute error tolerances for the Integrator.
1250
1251 If the pointer 'atol' is NULL, the necessary space is allocated (and
1252 must be freed somewhere else).
1253 */
1254 double *integrator_get_atol(IntegratorSystem *sys, double *atol){
1255 long i;
1256 #ifdef ATOL_DEBUG
1257 char *varname;
1258 #endif
1259
1260 if (atol==NULL) {
1261 atol = ASC_NEW_ARRAY_CLEAR(double, sys->n_y);
1262 }
1263
1264 for (i=0; i< sys->n_y; i++) {
1265 asc_assert(sys->y[i]!=NULL);
1266 atol[i] = var_odeatol(sys->y[i]);
1267 asc_assert(atol[i]!=-1);
1268 #ifdef ATOL_DEBUG
1269 varname = var_make_name(sys->system,sys->y[i]);
1270 CONSOLE_DEBUG("%s.ode_atol = %8.2e",varname,atol[i]);
1271 ASC_FREE(varname);
1272 #endif
1273 }
1274 return atol;
1275 }
1276
1277 /*-------------------------------------------------------------
1278 RETRIEVING OBSERVATION DATA
1279 */
1280
1281 /**
1282 This function takes the inst in the solver and returns the vector of
1283 observation variables that are located in the submodel d.obs array.
1284 */
1285 double *integrator_get_observations(IntegratorSystem *sys, double *obsi) {
1286 long i;
1287
1288 if (obsi==NULL) {
1289 obsi = ASC_NEW_ARRAY_CLEAR(double, sys->n_obs+1);
1290 }
1291
1292 /* C obsi[0] <==> ascend d.obs[1] */
1293
1294 for (i=0; i < sys->n_obs; i++) {
1295 obsi[i] = var_value(sys->obs[i]);
1296 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1297 }
1298 return obsi;
1299 }
1300
1301 struct var_variable *integrator_get_observed_var(IntegratorSystem *sys, const long i){
1302 asc_assert(i>=0);
1303 asc_assert(i<sys->n_obs);
1304 return sys->obs[i];
1305 }
1306
1307 /**
1308 @NOTE Although this shouldn't be required for implementation of solver
1309 engines, this is useful for GUI reporting of integration results.
1310 */
1311 struct var_variable *integrator_get_independent_var(IntegratorSystem *sys){
1312 return sys->x;
1313 }
1314
1315 int integrator_write_matrix(const IntegratorSystem *sys, FILE *fp,const char *type){
1316 asc_assert(sys);
1317 asc_assert(sys->enginedata);
1318 asc_assert(sys->internals);
1319 asc_assert(sys->internals->name);
1320 if(sys->internals->writematrixfn){
1321 return (sys->internals->writematrixfn)(sys,fp,type);
1322 }else{
1323 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Integrator '%s' defines no write_matrix function.",sys->internals->name);
1324 return -1;
1325 }
1326 }
1327
1328 int integrator_debug(const IntegratorSystem *sys, FILE *fp){
1329 asc_assert(sys);
1330 asc_assert(sys->enginedata);
1331 asc_assert(sys->internals);
1332 asc_assert(sys->internals->name);
1333 if(sys->internals->debugfn){
1334 return (sys->internals->debugfn)(sys,fp);
1335 }else{
1336 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Integrator '%s' defines no write_matrix function.",sys->internals->name);
1337 return -1;
1338 }
1339 }
1340
1341 /*----------------------------------------------------
1342 Build an analytic jacobian for solving the state system
1343
1344 This necessarily ugly piece of code attempts to create a unique
1345 list of relations that explicitly contain the variables in the
1346 given input list. The utility of this information is that we know
1347 exactly which relations must be differentiated, to fill in the
1348 df/dy matrix. If the problem has very few derivative terms, this will
1349 be of great savings. If the problem arose from the discretization of
1350 a pde, then this will be not so useful. The decision wether to use
1351 this function or to simply differentiate the entire relations list
1352 must be done before calling this function.
1353
1354 Final Note: the callee owns the array, but not the array elements.
1355 */
1356 #define AVG_NUM_INCIDENT 4
1357
1358
1359 /**
1360 This function helps to arrange the observation variables in a sensible order.
1361 The 'obs_id' child instance of v, if present, is assigned the value of the
1362 given parameter 'index'.
1363 */
1364 static void Integ_SetObsId(struct var_variable *v, long index){
1365 struct Instance *c, *i;
1366 i = var_instance(v);
1367 asc_assert(i!=NULL);
1368 c = ChildByChar(i,OBSINDEX);
1369 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1370 return;
1371 }
1372 SetIntegerAtomValue(c,index,0);
1373 }
1374
1375 /**
1376 Compares observation structs. NULLs should end up at far end.
1377 */
1378 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1379 if(v1 == NULL)return 1;
1380 if(v2 == NULL)return -1;
1381 if(v1->index > v2->index)return 1;
1382 if(v1->index == v2->index)return 0;
1383 return -1;
1384 }
1385
1386 /**
1387 Compares dynamic vars structs. NULLs should end up at far end.
1388 List should be sorted primarily by index and then by type, in order
1389 of increasing value of both.
1390 */
1391 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1392 if(v1 == NULL)return 1;
1393 if(v2 == NULL)return -1;
1394 if(v1->index > v2->index)return 1;
1395 if(v1->index != v2->index)return -1;
1396 if(v1->type > v2->type)return 1;
1397 return -1;
1398 }
1399 /*----------------------------
1400 Output handling to the GUI/interface.
1401 */
1402
1403 int integrator_set_reporter(IntegratorSystem *sys
1404 , IntegratorReporter *reporter
1405 ){
1406 asc_assert(sys!=NULL);
1407 sys->reporter = reporter;
1408 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1409 return 1;
1410 }
1411
1412 int integrator_output_init(IntegratorSystem *sys){
1413 asc_assert(sys!=NULL);
1414 asc_assert(sys->reporter!=NULL);
1415 if(sys->reporter->init!=NULL){
1416 /* call the specified output function */
1417 return (*(sys->reporter->init))(sys);
1418 }
1419 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1420 return 1;
1421 }
1422
1423 int integrator_output_write(IntegratorSystem *sys){
1424 static int reported_already=0;
1425 asc_assert(sys!=NULL);
1426 if(sys->reporter->write!=NULL){
1427 return (*(sys->reporter->write))(sys);
1428 }
1429 if(!reported_already){
1430 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1431 reported_already=1;
1432 }
1433 return 1;
1434 }
1435
1436 int integrator_output_write_obs(IntegratorSystem *sys){
1437 static int reported_already=0;
1438 asc_assert(sys!=NULL);
1439 if(sys->reporter->write_obs!=NULL){
1440 return (*(sys->reporter->write_obs))(sys);
1441 }
1442 if(!reported_already){
1443 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1444 reported_already=1;
1445 }
1446 return 1;
1447 }
1448
1449 int integrator_output_close(IntegratorSystem *sys){
1450 asc_assert(sys!=NULL);
1451 if(sys->reporter->close!=NULL){
1452 return (*(sys->reporter->close))(sys);
1453 }
1454 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1455 return 1;
1456 }
1457
1458 /**
1459 Decode status codes from the integrator, and output them via FPRINTF.
1460
1461 @return 0 on status ok (converged), <0 on unrecognised state, >0 on recognised error state.
1462 */
1463 int integrator_checkstatus(slv_status_t status) {
1464 if(status.converged){
1465 return 0;
1466 }
1467
1468 if(status.diverged){
1469 ERROR_REPORTER_HERE(ASC_PROG_ERR,"The derivative system did not converge."
1470 " Integration will terminate."); return 1;
1471 }else if(status.inconsistent){
1472 ERROR_REPORTER_HERE(ASC_PROG_ERR,"A numerically inconsistent state was discovered while"
1473 " calculating derivatives. Integration will terminate."); return 1;
1474 }else if(status.time_limit_exceeded){
1475 ERROR_REPORTER_HERE(ASC_PROG_ERR,"The time limit was exceeded while calculating"
1476 " derivatives. Integration will terminate."); return 1;
1477 }else if(status.iteration_limit_exceeded){
1478 ERROR_REPORTER_HERE(ASC_PROG_ERR,"The iteration limit was exceeded while calculating"
1479 " derivatives. Integration will terminate."); return 1;
1480 }else if(status.panic){
1481 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration cancelled by user."); return 1;
1482 }
1483
1484 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised solver state (non converged)");
1485 return -1;
1486 }

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