/[ascend]/trunk/pygtk/integrator.cpp
ViewVC logotype

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 921 - (show annotations) (download) (as text)
Mon Nov 6 07:49:06 2006 UTC (13 years, 1 month ago) by johnpye
File MIME type: text/x-c++src
File size: 4637 byte(s)
Fixing up some malloc/free problems with  integrator.c, removing some debug output from slv3.
Working on fixing IDA for systems with inactive variables (ongoing).
1 #include "integrator.h"
2 #include "integratorreporter.h"
3 #include <stdexcept>
4 using namespace std;
5
6 /**
7 'creating' an integrator in the context of the GUI just means an object
8 we can store the parameters that will be later sent to the underlying
9 C-code API.
10 */
11 Integrator::Integrator(Simulation &simulation)
12 : simulation(simulation)
13 {
14 // create the C-level object
15 this->blsys = integrator_new(simulation.getSystem(),simulation.getModel().getInternalType());
16
17 samplelist = NULL;
18
19 // set default steps
20 setMinSubStep(0);
21 setMaxSubStep(0);
22 setMaxSubSteps(0);
23 setInitialSubStep(0);
24 }
25
26 Integrator::~Integrator(){
27 CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
28 CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
29 integrator_free(blsys);
30 CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
31 samplelist_free(samplelist);
32 }
33
34
35 void
36 Integrator::setReporter(IntegratorReporterCxx *reporter){
37 this->blsys->clientdata = reporter;
38 integrator_set_reporter(blsys,reporter->getInternalType());
39 CONSOLE_DEBUG("REPORTER HAS BEEN SET");
40 (*(this->blsys->reporter->init))(blsys);
41 CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
42 }
43
44 double
45 Integrator::getCurrentTime(){
46 return integrator_get_t(blsys);
47 }
48
49 long
50 Integrator::getCurrentStep(){
51 return integrator_getcurrentstep(blsys);
52 }
53
54 long
55 Integrator::getNumSteps(){
56 return integrator_getnsamples(blsys);
57 }
58
59 int
60 Integrator::findIndependentVar(){
61 return integrator_find_indep_var(blsys);
62 }
63
64 int
65 Integrator::analyse(){
66
67 int res;
68 res = integrator_analyse(blsys);
69
70 if(!res){
71 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Failed system analysis");
72 return 0;
73 }
74
75 return 1;
76 }
77
78 /**
79 @TODO what about root detection?
80 */
81 int
82 Integrator::solve(){
83
84 // check the integration limits
85 // trigger of the solution process
86 // report errors?
87
88 assert(samplelist!=NULL);
89 assert(samplelist->ns>0);
90 assert(blsys->reporter!=NULL);
91 assert(blsys->clientdata!=NULL);
92
93 int res;
94 res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
95
96 // communicate solver variable status back to the instance tree via 'interface_ptr'
97 simulation.processVarStatus();
98
99 if(!res){
100 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Failed integration");
101 return 0;
102 }
103
104 return 1;
105 }
106
107 int
108 Integrator::setEngine(IntegratorEngine engine){
109 return integrator_set_engine(this->blsys, engine);
110 }
111
112 int
113 Integrator::setEngine(int engine){
114 return integrator_set_engine(this->blsys, (IntegratorEngine)engine);
115 }
116
117 /**
118 Ideally this list would be dynamically generated based on what solvers
119 are available or are in memory.
120 */
121 map<int,string>
122 Integrator::getEngines() const{
123 map<int,string> m;
124 #ifdef ASC_WITH_LSODE
125 m.insert(pair<int,string>(INTEG_LSODE,"LSODE"));
126 #endif
127 #ifdef ASC_WITH_IDA
128 m.insert(pair<int,string>(INTEG_IDA,"IDA"));
129 #endif
130 return m;
131 }
132
133 string
134 Integrator::getEngineName() const{
135 map<int,string> m=getEngines();
136 map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
137 if(f==m.end()){
138 throw runtime_error("No engine selected");
139 }
140 return f->second;
141 }
142
143 void
144 Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
145 if(samplelist!=NULL){
146 ASC_FREE(samplelist);
147 }
148 const dim_type *d = units.getDimensions().getInternalType();
149 samplelist = samplelist_new(num+1, d);
150 double val = start;
151 double inc = (end-start)/(num);
152 for(unsigned long i=0;i<=num;++i){
153 samplelist_set(samplelist,i,val);
154 val += inc;
155 }
156 integrator_set_samples(blsys,samplelist);
157 }
158
159 vector<double>
160 Integrator::getCurrentObservations(){
161 double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
162 integrator_get_observations(blsys,d);
163 vector<double> v=vector<double>(d,d+getNumObservedVars());
164 // do I need to free d?
165 // can I do this in such a way as I avoid all this memory-copying?
166 return v;
167 }
168
169 Variable
170 Integrator::getObservedVariable(const long &i){
171 var_variable *v = integrator_get_observed_var(blsys,i);
172 return Variable(&simulation,v);
173 }
174
175 Variable
176 Integrator::getIndependentVariable(){
177 var_variable *v = integrator_get_independent_var(blsys);
178 if(v==NULL){
179 throw runtime_error("independent variable is null");
180 }
181 return Variable(&simulation,v);
182 }
183
184 int
185 Integrator::getNumVars(){
186 return blsys->n_y;
187 }
188
189 int
190 Integrator::getNumObservedVars(){
191 return blsys->n_obs;
192 }
193
194 void
195 Integrator::setMinSubStep(double n){
196 integrator_set_minstep(blsys,n);
197 }
198
199 void
200 Integrator::setMaxSubStep(double n){
201 integrator_set_maxstep(blsys,n);
202 }
203
204 void
205 Integrator::setInitialSubStep(double n){
206 integrator_set_stepzero(blsys,n);
207 }
208
209 void
210 Integrator::setMaxSubSteps(int n){
211 integrator_set_maxsubsteps(blsys,n);
212 }
213
214 IntegratorSystem *
215 Integrator::getInternalType(){
216 return blsys;
217 }

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