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

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 966 - (show annotations) (download) (as text)
Thu Dec 14 14:04:54 2006 UTC (13 years, 10 months ago) by johnpye
File MIME type: text/x-c++src
File size: 6757 byte(s)
Cleaned up code comments in units.h
Fixed a bug with starting timesteps in LSODE (i think)
Added function to output a SampleList to the console (for debugging)
Removed debug output from base/generic/test/SConscript.
Removed 'custom' headers for UnitsM in ascpy.i
Made a sane default ctor for UnitsM
Added some more test cases to the Python unittest suite (still a problem running multiple solver-based tests in a single run)
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 int
76 Integrator::findIndependentVar(){
77 return integrator_find_indep_var(blsys);
78 }
79
80 int
81 Integrator::analyse(){
82
83 int res;
84 res = integrator_analyse(blsys);
85
86 if(!res){
87 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Failed system analysis");
88 return 0;
89 }
90
91 return 1;
92 }
93
94 /**
95 @TODO what about root detection?
96
97 Integrate the function for the timesteps specified.
98
99 This method will throw a runtime_error if integration fails.
100
101 @TODO does simulation.processVarStatus work for integrators like IDA???
102 */
103 void
104 Integrator::solve(){
105
106 // check the integration limits
107 // trigger of the solution process
108 // report errors?
109
110 assert(samplelist!=NULL);
111 assert(samplelist->ns>0);
112 assert(blsys->reporter!=NULL);
113 assert(blsys->clientdata!=NULL);
114
115 int res;
116 res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
117
118 if(!res){
119 throw runtime_error("Failed integration");
120 }
121
122 // communicate solver variable status back to the instance tree via 'interface_ptr'
123 simulation.processVarStatus();
124 }
125
126 void
127 Integrator::setEngine(IntegratorEngine engine){
128 int res = integrator_set_engine(this->blsys, engine);
129 if(!res)return;
130 if(res==1)throw range_error("Unknown integrator");
131 if(res==2)throw range_error("Invalid integrator");
132 stringstream ss;
133 ss << "Unknown error in setEngine (res = " << res << ")";
134 throw runtime_error(ss.str());
135 }
136
137 void
138 Integrator::setEngine(int engine){
139 setEngine((IntegratorEngine)engine);
140 }
141
142 void
143 Integrator::setEngine(const string &name){
144 CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
145 IntegratorEngine engine = INTEG_UNKNOWN;
146 #ifdef ASC_WITH_LSODE
147 if(name=="LSODE")engine = INTEG_LSODE;
148 #endif
149 #ifdef ASC_WITH_IDA
150 if(name=="IDA")engine = INTEG_IDA;
151 #endif
152 if(engine==INTEG_UNKNOWN){
153 throw runtime_error("Unkown integrator name");
154 }
155 setEngine(engine);
156 }
157
158 /**
159 Ideally this list would be dynamically generated based on what solvers
160 are available or are in memory.
161 */
162 map<int,string>
163 Integrator::getEngines(){
164 map<int,string> m;
165 #ifdef ASC_WITH_LSODE
166 m.insert(pair<int,string>(INTEG_LSODE,"LSODE"));
167 #endif
168 #ifdef ASC_WITH_IDA
169 m.insert(pair<int,string>(INTEG_IDA,"IDA"));
170 #endif
171 return m;
172 }
173
174 string
175 Integrator::getEngineName() const{
176 map<int,string> m=getEngines();
177 map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
178 if(f==m.end()){
179 throw runtime_error("No engine selected");
180 }
181 return f->second;
182 }
183
184 /**
185 @TODO what about conversion factors? Is an allowance being made?
186 */
187 void
188 Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
189 if(samplelist!=NULL){
190 ASC_FREE(samplelist);
191 }
192 const dim_type *d = units.getDimensions().getInternalType();
193 samplelist = samplelist_new(num+1, d);
194 double val = start;
195 double inc = (end-start)/(num);
196 for(unsigned long i=0;i<=num;++i){
197 samplelist_set(samplelist,i,val);
198 val += inc;
199 }
200 integrator_set_samples(blsys,samplelist);
201 }
202
203 /**
204 @TODO what about conversion factors? Is an allowance being made?
205 */
206 void
207 Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
208 if(samplelist!=NULL){
209 ASC_FREE(samplelist);
210 }
211 const dim_type *d = units.getDimensions().getInternalType();
212
213 if(start<=0)throw runtime_error("starting timestep needs to be > 0");
214 if(end<=0)throw runtime_error("end timestep needs to be > 0");
215 if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
216
217 samplelist = samplelist_new(num+1, d);
218 double val = start;
219 double inc = exp((log(end)-log(start))/num);
220 for(unsigned long i=0;i<=num;++i){
221 samplelist_set(samplelist,i,val);
222 CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
223 val *= inc;
224 }
225 integrator_set_samples(blsys,samplelist);
226 }
227
228 vector<double>
229 Integrator::getCurrentObservations(){
230 double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
231 integrator_get_observations(blsys,d);
232 vector<double> v=vector<double>(d,d+getNumObservedVars());
233 // do I need to free d?
234 // can I do this in such a way as I avoid all this memory-copying?
235 return v;
236 }
237
238 Variable
239 Integrator::getObservedVariable(const long &i){
240 var_variable *v = integrator_get_observed_var(blsys,i);
241 return Variable(&simulation,v);
242 }
243
244 Variable
245 Integrator::getIndependentVariable(){
246 var_variable *v = integrator_get_independent_var(blsys);
247 if(v==NULL){
248 throw runtime_error("independent variable is null");
249 }
250 return Variable(&simulation,v);
251 }
252
253 int
254 Integrator::getNumVars(){
255 return blsys->n_y;
256 }
257
258 int
259 Integrator::getNumObservedVars(){
260 return blsys->n_obs;
261 }
262
263 void
264 Integrator::setMinSubStep(double n){
265 integrator_set_minstep(blsys,n);
266 }
267
268 void
269 Integrator::setMaxSubStep(double n){
270 integrator_set_maxstep(blsys,n);
271 }
272
273 void
274 Integrator::setInitialSubStep(double n){
275 integrator_set_stepzero(blsys,n);
276 }
277
278 void
279 Integrator::setMaxSubSteps(int n){
280 integrator_set_maxsubsteps(blsys,n);
281 }
282
283 IntegratorSystem *
284 Integrator::getInternalType(){
285 return blsys;
286 }

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