/[ascend]/branches/georgy/ascend/integrator/integrator.c
ViewVC logotype

Contents of /branches/georgy/ascend/integrator/integrator.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3105 - (show annotations) (download) (as text)
Tue May 24 01:31:12 2016 UTC (8 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 43548 byte(s)
Creating new branch for [[User:Georgy]] for [[GSOC2016]], from trunk r3104.

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, see <http://www.gnu.org/licenses/>.
16 *//**
17 @file
18 Integrator API for ASCEND, for solving systems of ODEs and/or DAEs.
19 *//*
20 by John Pye, May 2006
21 based on parts of Integrators.c in the Tcl/Tk interface directory, heavily
22 modified to provide a non-GUI-specific API and modularised for multiple
23 integration engines.
24 */
25 #include <time.h>
26 #include <string.h>
27
28 #include <ascend/general/panic.h>
29 #include <ascend/general/ascMalloc.h>
30 #include <ascend/compiler/packages.h>
31 #include <ascend/compiler/link.h>
32
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, IntegratorSystem *sys);
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){
223 gl_destroy(L);
224 L = NULL;
225 }
226 init = 0;
227 return NULL;
228 }
229 if(!init){
230 L = gl_create(10);
231 init = 1; /* set init to 1 here now, since this will be called
232 recursively from the LoadArchiveLibrary call. */
233
234 /* CONSOLE_DEBUG("REGISTERING STANDARD SOLVER ENGINES"); */
235 for(i=0; defaultintegrators[i]!=NULL;++i){
236 error = package_load(defaultintegrators[i],NULL);
237 if(error){
238 ERROR_REPORTER_HERE(ASC_PROG_ERR
239 ,"Unable to register integrator '%s' (error %d)."
240 ,defaultintegrators[i],error
241 );
242 }else{
243 CONSOLE_DEBUG("Integrator '%s' registered OK",defaultintegrators[i]);
244 }
245 }
246 }
247 return L;
248 }
249
250 /**
251 Return gl_list of IntegratorInternals. C++ will use this to produce a
252 nice little list of integrator names that can be used in Python :-/
253 */
254 const struct gl_list_t *integrator_get_engines(){
255 return integrator_get_list(0);
256 }
257
258 struct gl_list_t *integrator_get_engines_growable(){
259 return integrator_get_list(0);
260 }
261
262 void integrator_free_engines(){
263 integrator_get_list(1);
264 }
265
266 /* return 0 on success */
267 int integrator_set_engine(IntegratorSystem *sys, const char *name){
268 struct gl_list_t *L = integrator_get_engines_growable();
269 int i;
270 const IntegratorInternals *I, *Ifound=NULL;
271 for(i=1; i <= gl_length(L); ++i){
272 I = gl_fetch(L,i);
273 if(strcmp(I->name,name)==0){
274 Ifound = I;
275 break;
276 }
277 }
278 if(Ifound){
279 /** @TODO tests for applicability of this engine... */
280
281 CONSOLE_DEBUG("Setting engine...");
282 if(Ifound->engine == sys->engine){
283 // already set...
284 return 0;
285 }
286
287 // clean up if old engine is there
288 if(sys->engine!=INTEG_UNKNOWN){
289 #ifdef DESTROY_DEBUG
290 CONSOLE_DEBUG("Freeing memory used by old integrator engine");
291 #endif
292 integrator_free_engine(sys);
293 #ifdef DESTROY_DEBUG
294 CONSOLE_DEBUG("done");
295 #endif
296 }
297
298 sys->engine = Ifound->engine;
299
300 sys->internals = Ifound;
301 }else{
302 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid engine name '%s'",name);
303 return 1;
304 }
305
306 asc_assert(sys->internals);
307 integrator_create_engine(sys);
308 return 0;
309 }
310
311 /**
312 @TODO rename this
313 */
314 const IntegratorInternals *integrator_get_engine(const IntegratorSystem *sys){
315 return sys->internals;
316 }
317
318 /**
319 Free any engine-specific data that was required for the solution of
320 this system. Note that this data is pointed to by sys->enginedata.
321
322 @TODO rename this, bad choice/confusing name.
323 */
324 void integrator_free_engine(IntegratorSystem *sys){
325 if(sys->engine==INTEG_UNKNOWN)return;
326 if(sys->enginedata){
327 if(sys->internals){
328 #ifdef DESTROY_DEBUG
329 CONSOLE_DEBUG("Running engine's freefn");
330 #endif
331 (sys->internals->freefn)(sys->enginedata);
332 #ifdef DESTROY_DEBUG
333 CONSOLE_DEBUG("Done with freefn");
334 #endif
335 sys->enginedata=NULL;
336 }else{
337 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to free engine data: no sys->internals");
338 }
339 }
340 #ifdef DESTROY_DEBUG
341 CONSOLE_DEBUG("Done with integrator_free_engine");
342 #endif
343 }
344
345 /**
346 Create enginedata memory if required for this solver. This doesn't include
347 allocating computation space, since we assume that this stage all we know
348 is that we want to use a specified integrator engine, not the full details
349 of the problem at hand. Allocating space inside enginedata should be done
350 during the solve stage (and freed inside integrator_free_engine)
351
352 @TODO rename this
353 */
354 void integrator_create_engine(IntegratorSystem *sys){
355 asc_assert(sys);
356 asc_assert(sys->engine!=INTEG_UNKNOWN);
357 asc_assert(sys->internals);
358 asc_assert(sys->internals->createfn);
359 asc_assert(sys->enginedata==NULL);
360 (sys->internals->createfn)(sys);
361 }
362
363 int integrator_register(const IntegratorInternals *integ){
364 /* get the current list of registered engines */
365 struct gl_list_t *L;
366 L = integrator_get_engines_growable();
367
368 CONSOLE_DEBUG("REGISTERING INTEGRATOR");
369 CONSOLE_DEBUG("There were %lu registered integrators", gl_length(integrator_get_list(0)));
370
371 int i;
372 IntegratorInternals *I;
373 for(i=1; i <= gl_length(L); ++i){
374 I = (IntegratorInternals *)gl_fetch(L,i);
375 if(strcmp(integ->name,I->name)==0){
376 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Integrator with name '%s' is already registered",integ->name);
377 return 0;
378 }
379 if(integ->engine == I->engine){
380 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Integrator with ID '%d' is already registered",integ->engine);
381 return 0;
382 }
383 }
384
385 CONSOLE_DEBUG("Adding engine '%s'",integ->name);
386
387 gl_append_ptr(L,(void *)integ);
388
389 CONSOLE_DEBUG("There are now %lu registered integrators", gl_length(integrator_get_list(0)));
390 return 0;
391 }
392
393 /*------------------------------------------------------------------------------
394 PARAMETERS
395 */
396
397 /**
398 Reset the parameters in this IntegratorSystem to the default ones for this
399 Integrator.
400
401 @return 0 on success, 1 on error
402 */
403 int integrator_params_default(IntegratorSystem *sys){
404 asc_assert(sys->engine!=INTEG_UNKNOWN);
405 asc_assert(sys->internals);
406 return (sys->internals->paramsdefaultfn)(sys);
407 }
408
409 int integrator_params_get(const IntegratorSystem *sys, slv_parameters_t *parameters){
410 asc_assert(sys!=NULL);
411 asc_assert(sys->params.num_parms > 0);
412 memcpy(parameters,&(sys->params),sizeof(slv_parameters_t));
413 return 0;
414 }
415
416 int integrator_params_set(IntegratorSystem *sys, const slv_parameters_t *parameters){
417 asc_assert(sys!=NULL);
418 asc_assert(parameters!=NULL);
419 memcpy(&(sys->params),parameters,sizeof(slv_parameters_t));
420 return 0;
421 }
422
423 /*------------------------------------------------------------------------------
424 ANALYSIS
425
426 Provide two modes in order to provide analysis suitable for solution of both
427 ODEs (the previous technique) and DAEs (new code). These share a common
428 "visit" method that needs to eventually be integrated with the code in
429 <solver/analyze.c>. For the moment, we're just hacking in to the compiler.
430 */
431
432 /**
433 Locate the independent variable.
434
435 @NOTE For the purpose of GUI design, this needs to work independent of the
436 integration engine being used. @ENDNOTE
437 */
438 int integrator_find_indep_var(IntegratorSystem *sys){
439 int result = 0;
440 #ifdef ANALYSE_DEBUG
441 char *varname;
442 #endif
443
444 /* if the indep var has been found, we don't look again (we assume the user won't fiddle with ode_type!) */
445 if(sys->x != NULL){
446 CONSOLE_DEBUG("sys->x already set");
447 return 0; /* success */
448 }
449
450 /* create a clear indepvars list */
451 if(sys->indepvars!=NULL){
452 ERROR_REPORTER_HERE(ASC_PROG_ERR,"indepvars should be NULL at this point");
453 return 1; /* error */
454 }
455 sys->indepvars = gl_create(10L);
456
457 IntegInitSymbols();
458
459 CONSOLE_DEBUG("Looking for independent var...");
460 integrator_visit_system_vars(sys,&integrator_classify_indep_var);
461 #ifdef ANALYSE_DEBUG
462 if(gl_length(sys->indepvars)){
463 CONSOLE_DEBUG("Found %lu indepvars",gl_length(sys->indepvars));
464 }else{
465 CONSOLE_DEBUG("NO INDEP VARS FOUND");
466 }
467 #endif
468
469 /* after visiting the instance tree, look at the candidates and return 0 on success */
470 result = integrator_check_indep_var(sys);
471
472 /* whatever happens, we clean up afterwards */
473 gl_free_and_destroy(sys->indepvars);
474 sys->indepvars = NULL;
475
476 #ifdef ANALYSE_DEBUG
477 asc_assert(sys->system);
478 if(!result){
479 asc_assert(sys->x);
480 varname = var_make_name(sys->system, sys->x);
481 CONSOLE_DEBUG("Indep var is '%s'",varname);
482 ASC_FREE(varname);
483 }else{
484 CONSOLE_DEBUG("No indep var was found");
485 }
486 #endif
487
488 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Returning result %d",result); */
489
490 return result;
491 }
492
493 /**
494 Perform whatever additional problem is required so that the system can be
495 integrated as a dynamical system with the IntegrationEngine chosen.
496
497 We can always assume that sys->system has had analyse_make_problem called on
498 it, so all the variable lists (etc) will be there already.
499
500 @return 0 on success
501 */
502 int integrator_analyse(IntegratorSystem *sys){
503 int res;
504
505 #ifdef ANALYSE_DEBUG
506 CONSOLE_DEBUG("Analysing integration system...");
507 #endif
508 asc_assert(sys);
509 if(sys->engine==INTEG_UNKNOWN){
510 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No engine selected: can't analyse");
511 return 1;
512 }
513 asc_assert(sys->engine!=INTEG_UNKNOWN);
514 asc_assert(sys->internals);
515
516 if(!sys->indepvars || !gl_length(sys->indepvars)){
517 if(integrator_find_indep_var(sys)){
518 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Independent variable problem: abandoning integration");
519 return 2;
520 }
521 #ifdef ANALYSE_DEBUG
522 else{
523 CONSOLE_DEBUG("got 0 from integrator_find_indep_var");
524 }
525 #endif
526 }
527
528 res = (sys->internals->analysefn)(sys);
529 #ifdef ANALYSE_DEBUG
530 CONSOLE_DEBUG("integrator_analyse returning %d",res);
531 #endif
532 return res;
533 }
534
535
536 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
537 struct var_variable **vlist;
538 int i, vlen;
539
540 /* visit all the slv_system_t master var lists to collect vars */
541 /* find the vars mostly in this one */
542 vlist = slv_get_solvers_var_list(sys->system);
543 vlen = slv_get_num_solvers_vars(sys->system);
544 for (i=0;i<vlen;i++) {
545 (*visitfn)(sys, vlist[i], &i);
546 }
547
548 /*
549 CONSOLE_DEBUG("Checked %d vars",vlen);
550 integrator_print_var_stats(sys);
551 */
552
553 /* probably nothing here, but gotta check. */
554 vlist = slv_get_master_par_list(sys->system);
555 vlen = slv_get_num_master_pars(sys->system);
556 for (i=0;i<vlen;i++) {
557 (*visitfn)(sys, vlist[i], NULL);
558 }
559
560 /*
561 CONSOLE_DEBUG("Checked %d pars",vlen);
562 integrator_print_var_stats(sys);
563 */
564
565 /* might find t here */
566 vlist = slv_get_master_unattached_list(sys->system);
567 vlen = slv_get_num_master_unattached(sys->system);
568 for (i=0;i<vlen;i++) {
569 (*visitfn)(sys, vlist[i], NULL);
570 }
571
572 /* CONSOLE_DEBUG("Checked %d unattached",vlen); */
573 }
574 /**
575 Analyse the ODE structure. We can assume that the independent variable was
576 already found.
577
578 @return 0 on success
579 */
580 int integrator_analyse_ode(IntegratorSystem *sys){
581 struct Integ_var_t *v1,*v2;
582 long half,i,len;
583 int happy=1;
584 char *varname1, *varname2;
585
586 asc_assert(sys->system!=NULL);
587
588 if(strcmp(slv_solver_name(slv_get_selected_solver(sys->system)),"QRSlv")!=0){
589 ERROR_REPORTER_HERE(ASC_PROG_ERR,"System must have solver 'QRSlv' assigned to it before integration");
590 return 2;
591 }
592 CONSOLE_DEBUG("Checked that NLA solver is set to '%s'",slv_solver_name(slv_get_selected_solver(sys->system)));
593
594 CONSOLE_DEBUG("Starting ODE analysis");
595 IntegInitSymbols();
596
597 /* collect potential states and derivatives */
598 sys->indepvars = gl_create(10L); /* t var info */
599 sys->dynvars = gl_create(200L); /* y ydot var info */
600 sys->obslist = gl_create(100L); /* obs info */
601 if (sys->dynvars == NULL
602 || sys->obslist == NULL
603 || sys->indepvars == NULL
604 ){
605 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
606 return 1;
607 }
608
609 sys->nstates = sys->nderivs = 0;
610
611 integrator_visit_system_vars(sys,&integrator_ode_classify_var);
612
613 integrator_print_var_stats(sys);
614
615 /* check sanity of state and var lists */
616
617 len = gl_length(sys->dynvars);
618 half = len/2;
619 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
620
621 if (len % 2 || len == 0L || sys->nstates != sys->nderivs ) {
622 /* list length must be even for vars to pair off */
623 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
624 return 3;
625 }
626 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
627 if (gl_fetch(sys->dynvars,len)==NULL) {
628 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
629 return 4;
630 }
631 sys->states = gl_create(half); /* state vars Integ_var_t references */
632 sys->derivs = gl_create(half); /* derivative var atoms */
633 for (i=1;i < len; i+=2) {
634 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
635 v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
636 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
637 varname1 = var_make_name(sys->system,v1->i);
638 varname2 = var_make_name(sys->system,v2->i);
639
640 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
641 varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
642 varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
643 );
644 ASC_FREE(varname1);
645 ASC_FREE(varname2);
646 happy=0;
647 break;
648 } else {
649 gl_append_ptr(sys->states,(POINTER)v1);
650 gl_append_ptr(sys->derivs,(POINTER)v2->i);
651 }
652 }
653 if (!happy) {
654 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem with ode_id and ode_type values");
655 return 5;
656 }
657 sys->n_y = half;
658 sys->y = ASC_NEW_ARRAY(struct var_variable *, half);
659 sys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
660 sys->y_id = ASC_NEW_ARRAY(int, half);
661 if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
662 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
663 return 6;
664 }
665 for (i = 1; i <= half; i++) {
666 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
667 sys->y[i-1] = v1->i;
668 sys->y_id[i-1] = v1->index;
669 sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
670 }
671
672 if(integrator_sort_obs_vars(sys)){
673 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error sorting observation variables");
674 return 7;
675 }
676
677 /* FIX all states */
678 for(i=0; i<sys->n_y; ++i){
679 if(!var_fixed(sys->y[i])){
680 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Fixing state %d",i);
681 var_set_fixed(sys->y[i], TRUE);
682 }
683 }
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 0;
699 }
700
701 /**
702 Reindex observations. Sort if the user mostly numbered. Take natural order
703 if user just booleaned.
704
705 @return 0 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(int, half);
741 if ( sys->obs==NULL || sys->obs_id==NULL) {
742 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
743 return 1;
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 0;
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 0 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 if(sys->x){
772 CONSOLE_DEBUG("Indep var already assigned");
773 return 1;
774 }
775
776 /* check the sanity of the independent variable */
777 len = gl_length(sys->indepvars);
778 if (!len) {
779 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
780 return 2;
781 }
782 if (len > 1) {
783 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
784 FPRINTF(ASCERR,"Excess %ld independent variables found:",
785 len);
786 for(i=1; i <=len;i++) {
787 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
788 if(info==NULL)continue;
789
790 varname = var_make_name(sys->system,info->i);
791 FPRINTF(ASCERR," %s",varname);
792 ASC_FREE(varname);
793 }
794 FPRINTF(ASCERR , "\nSet the '%s' flag on all but one of these to %s >= 0.\n"
795 , SCP(STATEFLAG),SCP(STATEFLAG)
796 );
797 error_reporter_end_flush();
798 return 3;
799 }else{
800 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
801 sys->x = info->i;
802 }
803 asc_assert(gl_length(sys->indepvars)==1);
804 asc_assert(sys->x);
805 return 0;
806 }
807
808 /*------------------------------------------------------------------------------
809 CLASSIFICATION OF VARIABLES (for ANALYSIS step)
810 */
811
812 #define INTEG_ADD_TO_LIST(info,TYPE,INDEX,VAR,VARINDX,LIST) \
813 info = ASC_NEW(struct Integ_var_t); \
814 if(info==NULL){ \
815 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory (INTEG_VAR_NEW)"); \
816 return; \
817 } \
818 info->type=TYPE; \
819 info->index=INDEX; \
820 info->i=VAR; \
821 info->derivative=NULL; \
822 info->derivative_of=NULL; \
823 if(VARINDX==NULL){ \
824 info->varindx = -1; \
825 }else{ \
826 info->varindx = *VARINDX; \
827 } \
828 gl_append_ptr(LIST,(void *)info); \
829 info = NULL
830
831 /**
832 In a DAE, it's either the (single) independent variable, or it's a
833 variable in the model.
834
835 I'm not sure what we should be doing with variables that are already
836 present as derivatives of other variables, I guess those ones need to be
837 removed from the list in a second pass?
838 */
839 void integrator_dae_classify_var(IntegratorSystem *sys
840 , struct var_variable *var, const int *varindx
841 ){
842 struct Integ_var_t *info;
843 long type,index;
844
845 /* filter for recognition of solver_vars */
846 var_filter_t vfilt;
847 vfilt.matchbits = VAR_SVAR;
848 vfilt.matchvalue = VAR_SVAR;
849
850 asc_assert(var != NULL && var_instance(var)!=NULL );
851
852 if( var_apply_filter(var,&vfilt) ) {
853 if(!var_active(var)){
854 CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");
855 return;
856 }
857
858 /* only non-fixed variables are accepted */
859 if(!var_fixed(var)){
860 /* get the ode_type and ode_id of this solver_var */
861 type = DynamicVarInfo(var,&index,sys);
862
863 if(type==INTEG_OTHER_VAR){
864 /* if the var's type is -1, it's independent */
865 INTEG_ADD_TO_LIST(info,INTEG_OTHER_VAR,0,var,varindx,sys->indepvars);
866 }else{
867 if(type < 0)type=0;
868 /* any other type of var is in the DAE system, at least for now */
869 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
870 }
871 }
872 #if 1
873 else{
874 /* fixed variable, only include it if ode_type == 1 */
875 type = DynamicVarInfo(var,&index,sys);
876 if(type==INTEG_STATE_VAR){
877 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
878 }
879 }
880 #endif
881
882 /* if the var's obs_id > 0, add it to the observation list */
883 if(ObservationVar(var,&index) != NULL && index > 0L) {
884 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->obslist);
885 }
886 }
887 }
888
889 /**
890 Inspect a specific variable and work out what type it is (what 'ode_type' it
891 has) and what other variable(s) it corresponds to (ie dydt corresponds to
892 y as a derivative).
893
894 @TODO add ability to create new variables for 'missing' derivative vars?
895 */
896 void integrator_ode_classify_var(IntegratorSystem *sys, struct var_variable *var
897 , const int *varindx
898 ){
899 struct Integ_var_t *info;
900 long type,index;
901
902 var_filter_t vfilt;
903 vfilt.matchbits = VAR_SVAR;
904 vfilt.matchvalue = VAR_SVAR;
905
906 asc_assert(var != NULL && var_instance(var)!=NULL );
907
908 if( var_apply_filter(var,&vfilt) ) {
909 /* it's a solver var: what type of variable? */
910 type = DynamicVarInfo(var,&index,sys);
911
912 if(type==INTEG_ALGEBRAIC_VAR){
913 /* no action required */
914 }else if(type==INTEG_OTHER_VAR){
915 /* i.e. independent var */
916 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
917 }else if(type>=INTEG_STATE_VAR){
918 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
919 if(type == 1){
920 sys->nstates++;
921 }else if(type == 2){ /* what about higher-order derivatives? -- JP */
922 sys->nderivs++;
923 }else{
924 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Higher-order (>=2) derivatives are not supported in ODEs.");
925 } }
926
927 if(ObservationVar(var,&index) != NULL && index > 0L) {
928 INTEG_ADD_TO_LIST(info,0L,index,var,varindx,sys->obslist);
929 }
930 }
931 }
932
933 /**
934 Look at a variable and determine if it's the independent variable or not.
935 This is just for the purpose of the integrator_find_indep_var function,
936 which is a utility function provided for use by the GUI.
937 */
938 void integrator_classify_indep_var(IntegratorSystem *sys
939 , struct var_variable *var, const int *varindx
940 ){
941 struct Integ_var_t *info;
942 long type,index;
943 var_filter_t vfilt;
944 #ifdef CLASSIFY_DEBUG
945 char *varname;
946 #endif
947
948 asc_assert(var != NULL && var_instance(var)!=NULL );
949 vfilt.matchbits = VAR_SVAR;
950 vfilt.matchvalue = VAR_SVAR;
951
952 #ifdef CLASSIFY_DEBUG
953 varname = var_make_name(sys->system,var);
954 #endif
955
956 if( var_apply_filter(var,&vfilt) ) {
957 type = DynamicVarInfo(var,&index,sys);
958
959 if(type==INTEG_OTHER_VAR){
960 /* i.e. independent var */
961 #ifdef CLASSIFY_DEBUG
962 CONSOLE_DEBUG("Var '%s' added to indepvars",varname);
963 #endif
964 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
965 #ifdef CLASSIFY_DEBUG
966 }else{
967 CONSOLE_DEBUG("Var '%s' not correct ode_type",varname);
968 #endif
969 }
970 #ifdef CLASSIFY_DEBUG
971 }else{
972 CONSOLE_DEBUG("Var '%s' failed filter",varname);
973 #endif
974 }
975
976 #ifdef CLASSIFY_DEBUG
977 ASC_FREE(varname);
978 #endif
979 }
980
981
982 /**
983 Look at a variable, and if it is an 'ODE variable' (it has a child instance
984 named 'ode_type') return its type, which will be either:
985 - INTEG_OTHER_VAR (if 'ode_type' is -1)
986 - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
987 - INTEG_STATE_VAR (if 'ode_type' is 1)
988 - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
989
990 If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
991 there.
992 */
993 static long DynamicVarInfo(struct var_variable *v,long *index, IntegratorSystem *sys){
994 struct Instance *c, *d, *i;
995 int type;
996
997 i = var_instance(v);
998
999 asc_assert(i!=NULL);
1000 asc_assert(STATEFLAG!=NULL);
1001 asc_assert(STATEINDEX!=NULL);
1002 c = ChildByChar(i,STATEFLAG);
1003 d = ChildByChar(i,STATEINDEX);
1004
1005
1006
1007 /* lazy evaluation is important in the following if */
1008 if(c == NULL
1009 || d == NULL
1010 || InstanceKind(c) != INTEGER_INST
1011 || InstanceKind(d) != INTEGER_INST
1012 || !AtomAssigned(c)
1013 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
1014 ){
1015 type = getOdeType(sys->instance,i);
1016 if(type != 0) {
1017 if(index !=NULL){
1018 *index = getOdeId(sys->instance,i);
1019 }
1020 return type;
1021 }
1022
1023 return INTEG_ALGEBRAIC_VAR;
1024 }
1025 if (index != NULL) {
1026 *index = GetIntegerAtomValue(d);
1027 }
1028 return GetIntegerAtomValue(c);
1029 }
1030
1031 /**
1032 Looks at the given variable checks if it is an 'observation variable'. This
1033 means that it has its 'obs_id' child instance set to a non-zero value.
1034
1035 If the variable is an observation variable, its index value ('obs_id') is
1036 stuff into *index (provided index!=NULL), and the pointer to the original
1037 instance is rtruend.
1038
1039 If it's not an observation variable, we return NULL and *index is untouched.
1040 */
1041 static struct var_variable *ObservationVar(struct var_variable *v, long *index){
1042 struct Instance *c,*i;
1043 i = var_instance(v);
1044 asc_assert(i!=NULL);
1045 c = ChildByChar(i,OBSINDEX);
1046 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1047 return NULL;
1048 }
1049 if (index != NULL) {
1050 *index = GetIntegerAtomValue(c);
1051 }
1052 return v;
1053 }
1054
1055 /*------------------------------------------------------------------------------
1056 RUNNING THE SOLVER
1057 */
1058
1059 /*
1060 Make the call to the actual integrator we've selected, for the range of
1061 time values specified. The sys contains all the specifics.
1062
1063 Return 0 on success
1064 */
1065 int integrator_solve(IntegratorSystem *sys, long i0, long i1){
1066
1067 long nstep;
1068 unsigned long start_index=0, finish_index=0;
1069 asc_assert(sys!=NULL);
1070
1071 asc_assert(sys->internals);
1072 asc_assert(sys->engine!=INTEG_UNKNOWN);
1073
1074 nstep = integrator_getnsamples(sys)-1;
1075 /* check for at least 2 steps and dimensionality of x vs steps here */
1076
1077 if (i0<0 || i1 <0) {
1078 /* removed completely inappropriate interactive code here */
1079 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
1080 return -1;
1081 } else {
1082 start_index=i0;
1083 finish_index =i1;
1084 if (start_index >= (unsigned long)nstep) {
1085 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
1086 ,start_index,nstep
1087 );
1088 return -2;
1089 }
1090 if (finish_index > (unsigned long)nstep) {
1091 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
1092 ,finish_index,nstep
1093 );
1094 return -3;
1095 }
1096 }
1097
1098 if(finish_index <= start_index) {
1099 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
1100 ,start_index,finish_index
1101 );
1102 return -4;
1103 }
1104
1105 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1106
1107 return (sys->internals->solvefn)(sys,start_index,finish_index);
1108 }
1109
1110 /*---------------------------------------------------------------
1111 HANDLING THE LIST OF TIMESTEMPS
1112 */
1113
1114 #define GETTER_AND_SETTER(TYPE,NAME) \
1115 void integrator_set_##NAME(IntegratorSystem *sys, TYPE val){ \
1116 sys->NAME=val; \
1117 } \
1118 TYPE integrator_get_##NAME(IntegratorSystem *sys){ \
1119 return sys->NAME; \
1120 }
1121
1122 GETTER_AND_SETTER(SampleList *,samples) /*;*/
1123 GETTER_AND_SETTER(double,maxstep) /*;*/
1124 GETTER_AND_SETTER(double,minstep) /*;*/
1125 GETTER_AND_SETTER(double,stepzero) /*;*/
1126 GETTER_AND_SETTER(int,maxsubsteps) /*;*/
1127 #undef GETTER_AND_SETTER
1128
1129 long integrator_getnsamples(IntegratorSystem *sys){
1130 asc_assert(sys!=NULL);
1131 asc_assert(sys->samples!=NULL);
1132 return samplelist_length(sys->samples);
1133 }
1134
1135 double integrator_getsample(IntegratorSystem *sys, long i){
1136 asc_assert(sys!=NULL);
1137 asc_assert(sys->samples!=NULL);
1138 return samplelist_get(sys->samples,i);
1139 }
1140
1141 void integrator_setsample(IntegratorSystem *sys, long i,double xi){
1142 asc_assert(sys!=NULL);
1143 asc_assert(sys->samples!=NULL);
1144 samplelist_set(sys->samples,i,xi);
1145 }
1146
1147 const dim_type *integrator_getsampledim(IntegratorSystem *sys){
1148 asc_assert(sys!=NULL);
1149 asc_assert(sys->samples!=NULL);
1150 return samplelist_dim(sys->samples);
1151 }
1152
1153 ASC_DLLSPEC long integrator_getcurrentstep(IntegratorSystem *sys){
1154 return sys->currentstep;
1155 }
1156
1157 /*------------------------------------------------------------------------------
1158 GET/SET VALUE OF THE INDEP VARIABLE
1159 */
1160
1161 /**
1162 Retrieve the value of the independent variable (time) from ASCEND
1163 and return it as a double.
1164 */
1165 double integrator_get_t(IntegratorSystem *sys){
1166 asc_assert(sys->x!=NULL);
1167 return var_value(sys->x);
1168 }
1169
1170 /**
1171 Set the value of the independent variable (time) in ASCEND.
1172 */
1173 void integrator_set_t(IntegratorSystem *sys, double value){
1174 var_set_value(sys->x, value);
1175 /* CONSOLE_DEBUG("set_t = %g", value); */
1176 }
1177
1178 /*------------------------------------------------------------------------------
1179 PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1180 */
1181 /**
1182 Retrieve the current values of the derivatives of the y-variables
1183 and stick them in the/an array that the integrator will use.
1184
1185 If the pointer 'y' is NULL, the necessary space is allocated (and
1186 must be freed somewhere else).
1187 */
1188 double *integrator_get_y(IntegratorSystem *sys, double *y) {
1189 long i;
1190
1191 if (y==NULL) {
1192 y = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1193 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1194 }
1195
1196 for (i=0; i< sys->n_y; i++) {
1197 asc_assert(sys->y[i]!=NULL);
1198 y[i] = var_value(sys->y[i]);
1199 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1200 }
1201 return y;
1202 }
1203
1204 /**
1205 Take the values of the differential variables from the array that the
1206 integrator uses, and use them to update the values of the corresponding
1207 variables in ASCEND.
1208 */
1209 void integrator_set_y(IntegratorSystem *sys, double *y) {
1210 long i;
1211 #ifdef SOLVE_DEBUG
1212 char *varname;
1213 #endif
1214
1215 for (i=0; i < sys->n_y; i++) {
1216 asc_assert(sys->y[i]!=NULL);
1217 var_set_value(sys->y[i],y[i]);
1218 #ifdef SOLVE_DEBUG
1219 varname = var_make_name(sys->system, sys->y[i]);
1220 CONSOLE_DEBUG("y[%ld] = %g --> '%s'", i+1, y[i], varname);
1221 ASC_FREE(varname);
1222 #endif
1223 }
1224 }
1225
1226 /**
1227 Send the values of the derivatives of the 'y' variables to the solver.
1228 Allocate space for an array if necessary.
1229
1230 Any element in sys->ydot that is NULL will be passed over (the value
1231 won't be modified in dydx).
1232 */
1233 double *integrator_get_ydot(IntegratorSystem *sys, double *dydx) {
1234 long i;
1235
1236 if (dydx==NULL) {
1237 dydx = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1238 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1239 }
1240
1241 for (i=0; i < sys->n_y; i++) {
1242 if(sys->ydot[i]!=NULL){
1243 dydx[i] = var_value(sys->ydot[i]);
1244 }
1245 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1246 }
1247 return dydx;
1248 }
1249
1250 void integrator_set_ydot(IntegratorSystem *sys, double *dydx) {
1251 long i;
1252 #ifdef SOLVE_DEBUG
1253 char *varname;
1254 #endif
1255 for (i=0; i < sys->n_y; i++) {
1256 if(sys->ydot[i]!=NULL){
1257 var_set_value(sys->ydot[i],dydx[i]);
1258 #ifdef SOLVE_DEBUG
1259 varname = var_make_name(sys->system, sys->ydot[i]);
1260 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1261 ASC_FREE(varname);
1262 }else{
1263 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1264 #endif
1265 }
1266 }
1267 }
1268
1269 /**
1270 Retrieve the values of 'ode_atol' properties of each of y-variables,
1271 for use in setting absolute error tolerances for the Integrator.
1272
1273 If the pointer 'atol' is NULL, the necessary space is allocated (and
1274 must be freed somewhere else).
1275 */
1276 double *integrator_get_atol(IntegratorSystem *sys, double *atol){
1277 long i;
1278 #ifdef ATOL_DEBUG
1279 char *varname;
1280 #endif
1281
1282 if (atol==NULL) {
1283 atol = ASC_NEW_ARRAY_CLEAR(double, sys->n_y);
1284 }
1285
1286 for (i=0; i< sys->n_y; i++) {
1287 asc_assert(sys->y[i]!=NULL);
1288 atol[i] = var_odeatol(sys->y[i]);
1289 asc_assert(atol[i]!=-1);
1290 #ifdef ATOL_DEBUG
1291 varname = var_make_name(sys->system,sys->y[i]);
1292 CONSOLE_DEBUG("%s.ode_atol = %8.2e",varname,atol[i]);
1293 ASC_FREE(varname);
1294 #endif
1295 }
1296 return atol;
1297 }
1298
1299 /*-------------------------------------------------------------
1300 RETRIEVING OBSERVATION DATA
1301 */
1302
1303 /**
1304 This function takes the inst in the solver and returns the vector of
1305 observation variables that are located in the submodel d.obs array.
1306 */
1307 double *integrator_get_observations(IntegratorSystem *sys, double *obsi) {
1308 long i;
1309
1310 if (obsi==NULL) {
1311 obsi = ASC_NEW_ARRAY_CLEAR(double, sys->n_obs+1);
1312 }
1313
1314 /* C obsi[0] <==> ascend d.obs[1] */
1315
1316 for (i=0; i < sys->n_obs; i++) {
1317 obsi[i] = var_value(sys->obs[i]);
1318 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1319 }
1320 return obsi;
1321 }
1322
1323 struct var_variable *integrator_get_observed_var(IntegratorSystem *sys, const long i){
1324 asc_assert(i>=0);
1325 asc_assert(i<sys->n_obs);
1326 return sys->obs[i];
1327 }
1328
1329 /**
1330 @NOTE Although this shouldn't be required for implementation of solver
1331 engines, this is useful for GUI reporting of integration results.
1332 */
1333 struct var_variable *integrator_get_independent_var(IntegratorSystem *sys){
1334 return sys->x;
1335 }
1336
1337 int integrator_write_matrix(const IntegratorSystem *sys, FILE *fp,const char *type){
1338 asc_assert(sys);
1339 asc_assert(sys->enginedata);
1340 asc_assert(sys->internals);
1341 asc_assert(sys->internals->name);
1342 if(sys->internals->writematrixfn){
1343 if(type!=NULL && strcmp(type,"")==0)type = NULL;
1344 return (sys->internals->writematrixfn)(sys,fp,type);
1345 }else{
1346 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Integrator '%s' defines no write_matrix function.",sys->internals->name);
1347 return -1;
1348 }
1349 }
1350
1351 int integrator_debug(const IntegratorSystem *sys, FILE *fp){
1352 asc_assert(sys);
1353 asc_assert(sys->enginedata);
1354 asc_assert(sys->internals);
1355 asc_assert(sys->internals->name);
1356 if(sys->internals->debugfn){
1357 return (sys->internals->debugfn)(sys,fp);
1358 }else{
1359 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Integrator '%s' defines no write_matrix function.",sys->internals->name);
1360 return -1;
1361 }
1362 }
1363
1364 /*----------------------------------------------------
1365 Build an analytic jacobian for solving the state system
1366
1367 This necessarily ugly piece of code attempts to create a unique
1368 list of relations that explicitly contain the variables in the
1369 given input list. The utility of this information is that we know
1370 exactly which relations must be differentiated, to fill in the
1371 df/dy matrix. If the problem has very few derivative terms, this will
1372 be of great savings. If the problem arose from the discretization of
1373 a pde, then this will be not so useful. The decision wether to use
1374 this function or to simply differentiate the entire relations list
1375 must be done before calling this function.
1376
1377 Final Note: the callee owns the array, but not the array elements.
1378 */
1379 #define AVG_NUM_INCIDENT 4
1380
1381
1382 /**
1383 This function helps to arrange the observation variables in a sensible order.
1384 The 'obs_id' child instance of v, if present, is assigned the value of the
1385 given parameter 'index'.
1386 */
1387 static void Integ_SetObsId(struct var_variable *v, long index){
1388 struct Instance *c, *i;
1389 i = var_instance(v);
1390 asc_assert(i!=NULL);
1391 c = ChildByChar(i,OBSINDEX);
1392 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1393 return;
1394 }
1395 SetIntegerAtomValue(c,index,0);
1396 }
1397
1398 /**
1399 Compares observation structs. NULLs should end up at far end.
1400 */
1401 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1402 if(v1 == NULL)return 1;
1403 if(v2 == NULL)return -1;
1404 if(v1->index > v2->index)return 1;
1405 if(v1->index == v2->index)return 0;
1406 return -1;
1407 }
1408
1409 /**
1410 Compares dynamic vars structs. NULLs should end up at far end.
1411 List should be sorted primarily by index and then by type, in order
1412 of increasing value of both.
1413 */
1414 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1415 if(v1 == NULL)return 1;
1416 if(v2 == NULL)return -1;
1417 if(v1->index > v2->index)return 1;
1418 if(v1->index != v2->index)return -1;
1419 if(v1->type > v2->type)return 1;
1420 return -1;
1421 }
1422 /*----------------------------
1423 Output handling to the GUI/interface.
1424 */
1425
1426 int integrator_set_reporter(IntegratorSystem *sys
1427 , IntegratorReporter *reporter
1428 ){
1429 asc_assert(sys!=NULL);
1430 sys->reporter = reporter;
1431 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1432 return 1;
1433 }
1434
1435 int integrator_output_init(IntegratorSystem *sys){
1436 asc_assert(sys!=NULL);
1437 asc_assert(sys->reporter!=NULL);
1438 if(sys->reporter->init!=NULL){
1439 /* call the specified output function */
1440 return (*(sys->reporter->init))(sys);
1441 }
1442 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1443 return 1;
1444 }
1445
1446 int integrator_output_write(IntegratorSystem *sys){
1447 static int reported_already=0;
1448 asc_assert(sys!=NULL);
1449 if(sys->reporter->write!=NULL){
1450 return (*(sys->reporter->write))(sys);
1451 }
1452 if(!reported_already){
1453 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1454 reported_already=1;
1455 }
1456 return 1;
1457 }
1458
1459 int integrator_output_write_obs(IntegratorSystem *sys){
1460 static int reported_already=0;
1461 asc_assert(sys!=NULL);
1462 if(sys->reporter->write_obs!=NULL){
1463 return (*(sys->reporter->write_obs))(sys);
1464 }
1465 if(!reported_already){
1466 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1467 reported_already=1;
1468 }
1469 return 1;
1470 }
1471
1472 int integrator_output_close(IntegratorSystem *sys){
1473 asc_assert(sys!=NULL);
1474 if(sys->reporter->close!=NULL){
1475 return (*(sys->reporter->close))(sys);
1476 }
1477 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1478 return 1;
1479 }
1480
1481 /**
1482 Decode status codes from the integrator, and output them via FPRINTF.
1483
1484 @return 0 on status ok (converged), <0 on unrecognised state, >0 on recognised error state.
1485 */
1486 int integrator_checkstatus(slv_status_t status) {
1487 if(status.converged){
1488 return 0;
1489 }
1490
1491 if(status.diverged){
1492 ERROR_REPORTER_HERE(ASC_PROG_ERR,"The derivative system did not converge."
1493 " Integration will terminate."); return 1;
1494 }else if(status.inconsistent){
1495 ERROR_REPORTER_HERE(ASC_PROG_ERR,"A numerically inconsistent state was discovered while"
1496 " calculating derivatives. Integration will terminate."); return 1;
1497 }else if(status.time_limit_exceeded){
1498 ERROR_REPORTER_HERE(ASC_PROG_ERR,"The time limit was exceeded while calculating"
1499 " derivatives. Integration will terminate."); return 1;
1500 }else if(status.iteration_limit_exceeded){
1501 ERROR_REPORTER_HERE(ASC_PROG_ERR,"The iteration limit was exceeded while calculating"
1502 " derivatives. Integration will terminate."); return 1;
1503 }else if(status.panic){
1504 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration cancelled by user."); return 1;
1505 }
1506
1507 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised solver state (non converged)");
1508 return -1;
1509 }

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