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

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1303 - (show annotations) (download) (as text)
Thu Mar 1 06:04:21 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-c++src
File size: 7776 byte(s)
silenced some debug 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
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
109 if(res){
110 stringstream ss; ss << "Failed system analysis (error " << res << ")";
111 throw runtime_error(ss.str());
112 }
113 }
114
115 /**
116 @TODO what about root detection?
117
118 Integrate the function for the timesteps specified.
119
120 Method will throw a runtime_error if integrator_solve returns error (non zero)
121
122 @TODO does simulation.processVarStatus work for integrators like IDA???
123 */
124 void
125 Integrator::solve(){
126
127 // check the integration limits
128 // trigger of the solution process
129 // report errors?
130
131 assert(samplelist!=NULL);
132 assert(samplelist->ns>0);
133 assert(blsys->reporter!=NULL);
134 assert(blsys->clientdata!=NULL);
135
136 int res;
137 res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
138
139 if(res){
140 stringstream ss;
141 ss << "Failed integration (integrator_solve returned " << res << ")";
142 throw runtime_error(ss.str());
143 }
144
145 // communicate solver variable status back to the instance tree via 'interface_ptr'
146 simulation.processVarStatus();
147 }
148
149 void
150 Integrator::writeMatrix(FILE *fp) const{
151 if(integrator_write_matrix(this->blsys, fp)){
152 throw runtime_error("Failed to write matrix");
153 }
154 }
155
156 void
157 Integrator::writeDebug(FILE *fp) const{
158 if(integrator_debug(this->blsys, fp)){
159 throw runtime_error("Failed to write debug output");
160 }
161 }
162
163 void
164 Integrator::setEngine(IntegratorEngine engine){
165 int res = integrator_set_engine(this->blsys, engine);
166 if(!res)return;
167 if(res==1)throw range_error("Unknown integrator");
168 if(res==2)throw range_error("Invalid integrator");
169 stringstream ss;
170 ss << "Unknown error in setEngine (res = " << res << ")";
171 throw runtime_error(ss.str());
172 }
173
174 void
175 Integrator::setEngine(int engine){
176 setEngine((IntegratorEngine)engine);
177 }
178
179 void
180 Integrator::setEngine(const string &name){
181 CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
182 IntegratorEngine engine = INTEG_UNKNOWN;
183 #ifdef ASC_WITH_LSODE
184 if(name=="LSODE")engine = INTEG_LSODE;
185 #endif
186 #ifdef ASC_WITH_IDA
187 if(name=="IDA")engine = INTEG_IDA;
188 #endif
189 if(engine==INTEG_UNKNOWN){
190 throw runtime_error("Unkown integrator name");
191 }
192 setEngine(engine);
193 }
194
195 /**
196 Ideally this list would be dynamically generated based on what solvers
197 are available or are in memory.
198 */
199 map<int,string>
200 Integrator::getEngines(){
201 map<int,string> m;
202 const IntegratorLookup *list = integrator_get_engines();
203 while(list->id != INTEG_UNKNOWN){
204 if(list->name==NULL)throw runtime_error("list->name is NULL");
205 m.insert(pair<int,string>(list->id,list->name));
206 ++list;
207 }
208 return m;
209 }
210
211 string
212 Integrator::getName() const{
213 map<int,string> m=getEngines();
214 map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
215 if(f==m.end()){
216 throw runtime_error("No engine selected");
217 }
218 return f->second;
219 }
220
221 /**
222 @TODO what about conversion factors? Is an allowance being made?
223 */
224 void
225 Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
226 if(samplelist!=NULL){
227 ASC_FREE(samplelist);
228 }
229 const dim_type *d = units.getDimensions().getInternalType();
230 samplelist = samplelist_new(num+1, d);
231 double val = start;
232 double inc = (end-start)/(num);
233 for(unsigned long i=0;i<=num;++i){
234 samplelist_set(samplelist,i,val);
235 val += inc;
236 }
237 integrator_set_samples(blsys,samplelist);
238 }
239
240 /**
241 @TODO what about conversion factors? Is an allowance being made?
242 */
243 void
244 Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
245 if(samplelist!=NULL){
246 ASC_FREE(samplelist);
247 }
248 const dim_type *d = units.getDimensions().getInternalType();
249
250 if(start<=0)throw runtime_error("starting timestep needs to be > 0");
251 if(end<=0)throw runtime_error("end timestep needs to be > 0");
252 if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
253
254 samplelist = samplelist_new(num+1, d);
255 double val = start;
256 double inc = exp((log(end)-log(start))/num);
257 for(unsigned long i=0;i<=num;++i){
258 samplelist_set(samplelist,i,val);
259 CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
260 val *= inc;
261 }
262 integrator_set_samples(blsys,samplelist);
263 }
264
265 vector<double>
266 Integrator::getCurrentObservations(){
267 double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
268 integrator_get_observations(blsys,d);
269 vector<double> v=vector<double>(d,d+getNumObservedVars());
270 // do I need to free d?
271 // can I do this in such a way as I avoid all this memory-copying?
272 return v;
273 }
274
275 Variable
276 Integrator::getObservedVariable(const long &i){
277 var_variable *v = integrator_get_observed_var(blsys,i);
278 return Variable(&simulation,v);
279 }
280
281 Variable
282 Integrator::getIndependentVariable(){
283 var_variable *v = integrator_get_independent_var(blsys);
284 if(v==NULL){
285 throw runtime_error("independent variable is null");
286 }
287 return Variable(&simulation,v);
288 }
289
290 int
291 Integrator::getNumVars(){
292 return blsys->n_y;
293 }
294
295 int
296 Integrator::getNumObservedVars(){
297 return blsys->n_obs;
298 }
299
300 void
301 Integrator::setMinSubStep(double n){
302 integrator_set_minstep(blsys,n);
303 }
304
305 void
306 Integrator::setMaxSubStep(double n){
307 integrator_set_maxstep(blsys,n);
308 }
309
310 void
311 Integrator::setInitialSubStep(double n){
312 integrator_set_stepzero(blsys,n);
313 }
314
315 void
316 Integrator::setMaxSubSteps(int n){
317 integrator_set_maxsubsteps(blsys,n);
318 }
319
320 IntegratorSystem *
321 Integrator::getInternalType(){
322 return blsys;
323 }

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