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

Contents of /trunk/pygtk/integrator.cpp

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

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