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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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