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

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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