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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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