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

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1387 - (show annotations) (download) (as text)
Sat Apr 7 14:43:31 2007 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-c++src
File size: 8100 byte(s)
Added plot support in Integrator output tabs.
Some other minor debugging for pylab integration and idaanalyse output.
1 #include "integrator.h"
2 #include "integratorreporter.h"
3 #include "solverparameters.h"
4 #include <stdexcept>
5 #include <sstream>
6 #include <cmath>
7 using namespace std;
8
9 // #define DESTROY_DEBUG
10
11 /**
12 'creating' an integrator in the context of the GUI just means an object
13 we can store the parameters that will be later sent to the underlying
14 C-code API.
15
16 @TODO at present the integrator requires a slv_system_t. This is the wrong
17 way around.
18 */
19 Integrator::Integrator(Simulation &simulation)
20 : simulation(simulation)
21 {
22 // create the C-level object
23 this->blsys = integrator_new(simulation.getSystem(),simulation.getModel().getInternalType());
24
25 samplelist = NULL;
26
27 // set default steps
28 setMinSubStep(0);
29 setMaxSubStep(0);
30 setMaxSubSteps(0);
31 setInitialSubStep(0);
32
33 /* we *don't* initialise a reporter here, cause we don't want to own it */
34 }
35
36 Integrator::~Integrator(){
37 #ifdef DESTROY_DEBUG
38 CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
39 CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
40 #endif
41 integrator_free(blsys);
42 #ifdef DESTROY_DEBUG
43 CONSOLE_DEBUG("done");
44 CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
45 #endif
46 if(samplelist)samplelist_free(samplelist);
47 }
48
49 SolverParameters
50 Integrator::getParameters() const{
51 SolverParameters params;
52 int res = integrator_params_get(blsys,&(params.getInternalType() ) );
53 if(res)throw runtime_error("Failed to get integrator parameters");
54 return params;
55 }
56
57 void
58 Integrator::setParameters(const SolverParameters &params){
59 int res = integrator_params_set(blsys,&(params.getInternalTypeConst() ) );
60 if(res)throw runtime_error("Failed to set integrator parameters");
61 }
62
63 void
64 Integrator::setReporter(IntegratorReporterCxx *reporter){
65 blsys->clientdata = reporter; /* this *is* necessary, as well as the following */
66 integrator_set_reporter(blsys,reporter->getInternalType());
67 //CONSOLE_DEBUG("REPORTER HAS BEEN SET");
68 (*(blsys->reporter->init))(blsys);
69 //CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
70 }
71
72 double
73 Integrator::getCurrentTime(){
74 return integrator_get_t(blsys);
75 }
76
77 long
78 Integrator::getCurrentStep(){
79 return integrator_getcurrentstep(blsys);
80 }
81
82 long
83 Integrator::getNumSteps(){
84 return integrator_getnsamples(blsys);
85 }
86
87 /**
88 Find the independent variable in the system, or throw an exception if not found.
89 */
90 void
91 Integrator::findIndependentVar(){
92 int res = integrator_find_indep_var(blsys);
93
94 if(res){
95 stringstream ss;
96 ss << "Independent variable not found (" << res << ")";
97 throw runtime_error(ss.str());
98 }
99 }
100
101 void
102 Integrator::analyse(){
103
104 int res;
105 /*
106 Note, we never need to call analyze_make_system in any of the Integrator
107 code, as it gets called by Simulation::build.
108 */
109 res = integrator_analyse(blsys);
110 CONSOLE_DEBUG("Got return-code '%d' from integrator_analyse",res);
111
112 if(res){
113 CONSOLE_DEBUG("...which is bad");
114 stringstream ss;
115 ss << "Failed system analysis (error " << res << ")";
116 throw runtime_error(ss.str());
117 }
118 }
119
120 /**
121 @TODO what about root detection?
122
123 Integrate the function for the timesteps specified.
124
125 Method will throw a runtime_error if integrator_solve returns error (non zero)
126
127 @TODO does simulation.processVarStatus work for integrators like IDA???
128 */
129 void
130 Integrator::solve(){
131
132 // check the integration limits
133 // trigger of the solution process
134 // report errors?
135
136 assert(samplelist!=NULL);
137 assert(samplelist->ns>0);
138
139 if(blsys->reporter==NULL){
140 throw runtime_error("No reporter has been assigned to the integrator");
141 }
142
143 assert(blsys->clientdata!=NULL);
144
145 int res;
146 res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
147
148 if(res){
149 stringstream ss;
150 ss << "Failed integration (integrator_solve returned " << res << ")";
151 throw runtime_error(ss.str());
152 }
153
154 // communicate solver variable status back to the instance tree via 'interface_ptr'
155 simulation.processVarStatus();
156 }
157
158 void
159 Integrator::writeMatrix(FILE *fp,const char *type) const{
160 if(integrator_write_matrix(this->blsys, fp, type)){
161 throw runtime_error("Failed to write matrix");
162 }
163 }
164
165 void
166 Integrator::writeDebug(FILE *fp) const{
167 if(integrator_debug(this->blsys, fp)){
168 throw runtime_error("Failed to write debug output");
169 }
170 }
171
172 void
173 Integrator::setEngine(IntegratorEngine engine){
174 int res = integrator_set_engine(this->blsys, engine);
175 if(!res)return;
176 if(res==1)throw range_error("Unknown integrator");
177 if(res==2)throw range_error("Invalid integrator");
178 stringstream ss;
179 ss << "Unknown error in setEngine (res = " << res << ")";
180 throw runtime_error(ss.str());
181 }
182
183 void
184 Integrator::setEngine(int engine){
185 setEngine((IntegratorEngine)engine);
186 }
187
188 void
189 Integrator::setEngine(const string &name){
190 CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
191 IntegratorEngine engine = INTEG_UNKNOWN;
192 #ifdef ASC_WITH_LSODE
193 if(name=="LSODE")engine = INTEG_LSODE;
194 #endif
195 #ifdef ASC_WITH_IDA
196 if(name=="IDA")engine = INTEG_IDA;
197 #endif
198 if(engine==INTEG_UNKNOWN){
199 throw runtime_error("Unkown integrator name");
200 }
201 setEngine(engine);
202 }
203
204 /**
205 Ideally this list would be dynamically generated based on what solvers
206 are available or are in memory.
207 */
208 map<int,string>
209 Integrator::getEngines(){
210 map<int,string> m;
211 const IntegratorLookup *list = integrator_get_engines();
212 while(list->id != INTEG_UNKNOWN){
213 if(list->name==NULL)throw runtime_error("list->name is NULL");
214 m.insert(pair<int,string>(list->id,list->name));
215 ++list;
216 }
217 return m;
218 }
219
220 string
221 Integrator::getName() const{
222 map<int,string> m=getEngines();
223 map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
224 if(f==m.end()){
225 throw runtime_error("No engine selected");
226 }
227 return f->second;
228 }
229
230 /**
231 @TODO what about conversion factors? Is an allowance being made?
232 */
233 void
234 Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
235 if(samplelist!=NULL){
236 ASC_FREE(samplelist);
237 }
238 const dim_type *d = units.getDimensions().getInternalType();
239 samplelist = samplelist_new(num+1, d);
240 double val = start;
241 double inc = (end-start)/(num);
242 for(unsigned long i=0;i<=num;++i){
243 samplelist_set(samplelist,i,val);
244 val += inc;
245 }
246 integrator_set_samples(blsys,samplelist);
247 }
248
249 /**
250 @TODO what about conversion factors? Is an allowance being made?
251 */
252 void
253 Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
254 if(samplelist!=NULL){
255 ASC_FREE(samplelist);
256 }
257 const dim_type *d = units.getDimensions().getInternalType();
258
259 if(start<=0)throw runtime_error("starting timestep needs to be > 0");
260 if(end<=0)throw runtime_error("end timestep needs to be > 0");
261 if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
262
263 samplelist = samplelist_new(num+1, d);
264 double val = start;
265 double inc = exp((log(end)-log(start))/num);
266 for(unsigned long i=0;i<=num;++i){
267 samplelist_set(samplelist,i,val);
268 // CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
269 val *= inc;
270 }
271 integrator_set_samples(blsys,samplelist);
272 }
273
274 vector<double>
275 Integrator::getCurrentObservations(){
276 double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
277 integrator_get_observations(blsys,d);
278 vector<double> v=vector<double>(d,d+getNumObservedVars());
279 // do I need to free d?
280 // can I do this in such a way as I avoid all this memory-copying?
281 return v;
282 }
283
284 Variable
285 Integrator::getObservedVariable(const long &i){
286 var_variable *v = integrator_get_observed_var(blsys,i);
287 return Variable(&simulation,v);
288 }
289
290 Variable
291 Integrator::getIndependentVariable(){
292 var_variable *v = integrator_get_independent_var(blsys);
293 if(v==NULL){
294 throw runtime_error("independent variable is null");
295 }
296 return Variable(&simulation,v);
297 }
298
299 int
300 Integrator::getNumVars(){
301 return blsys->n_y;
302 }
303
304 int
305 Integrator::getNumObservedVars(){
306 return blsys->n_obs;
307 }
308
309 void
310 Integrator::setMinSubStep(double n){
311 integrator_set_minstep(blsys,n);
312 }
313
314 void
315 Integrator::setMaxSubStep(double n){
316 integrator_set_maxstep(blsys,n);
317 }
318
319 void
320 Integrator::setInitialSubStep(double n){
321 integrator_set_stepzero(blsys,n);
322 }
323
324 void
325 Integrator::setMaxSubSteps(int n){
326 integrator_set_maxsubsteps(blsys,n);
327 }
328
329 IntegratorSystem *
330 Integrator::getInternalType(){
331 return blsys;
332 }

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