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

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 942 - (show annotations) (download) (as text)
Sat Nov 25 05:26:47 2006 UTC (14 years, 2 months ago) by johnpye
File MIME type: text/x-c++src
File size: 6600 byte(s)
Incorporated 'SolverParameters' functionality into 'Integrator', both at C level and C++/Python.
This shouldn't break Tcl/Tk as default parameters will be set and used invisibly.
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 void
102 Integrator::solve(){
103
104 // check the integration limits
105 // trigger of the solution process
106 // report errors?
107
108 assert(samplelist!=NULL);
109 assert(samplelist->ns>0);
110 assert(blsys->reporter!=NULL);
111 assert(blsys->clientdata!=NULL);
112
113 int res;
114 res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
115
116 // communicate solver variable status back to the instance tree via 'interface_ptr'
117 simulation.processVarStatus();
118
119 if(!res){
120 throw runtime_error("Failed integration");
121 }
122 }
123
124 void
125 Integrator::setEngine(IntegratorEngine engine){
126 int res = integrator_set_engine(this->blsys, engine);
127 if(!res)return;
128 if(res==1)throw range_error("Unknown integrator");
129 if(res==2)throw range_error("Invalid integrator");
130 stringstream ss;
131 ss << "Unknown error in setEngine (res = " << res << ")";
132 throw runtime_error(ss.str());
133 }
134
135 void
136 Integrator::setEngine(int engine){
137 setEngine((IntegratorEngine)engine);
138 }
139
140 void
141 Integrator::setEngine(const string &name){
142 CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
143 IntegratorEngine engine = INTEG_UNKNOWN;
144 #ifdef ASC_WITH_LSODE
145 if(name=="LSODE")engine = INTEG_LSODE;
146 #endif
147 #ifdef ASC_WITH_IDA
148 if(name=="IDA")engine = INTEG_IDA;
149 #endif
150
151 setEngine(engine);
152 }
153
154 /**
155 Ideally this list would be dynamically generated based on what solvers
156 are available or are in memory.
157 */
158 map<int,string>
159 Integrator::getEngines(){
160 map<int,string> m;
161 #ifdef ASC_WITH_LSODE
162 m.insert(pair<int,string>(INTEG_LSODE,"LSODE"));
163 #endif
164 #ifdef ASC_WITH_IDA
165 m.insert(pair<int,string>(INTEG_IDA,"IDA"));
166 #endif
167 return m;
168 }
169
170 string
171 Integrator::getEngineName() const{
172 map<int,string> m=getEngines();
173 map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
174 if(f==m.end()){
175 throw runtime_error("No engine selected");
176 }
177 return f->second;
178 }
179
180 /**
181 @TODO what about conversion factors? Is an allowance being made?
182 */
183 void
184 Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
185 if(samplelist!=NULL){
186 ASC_FREE(samplelist);
187 }
188 const dim_type *d = units.getDimensions().getInternalType();
189 samplelist = samplelist_new(num+1, d);
190 double val = start;
191 double inc = (end-start)/(num);
192 for(unsigned long i=0;i<=num;++i){
193 samplelist_set(samplelist,i,val);
194 val += inc;
195 }
196 integrator_set_samples(blsys,samplelist);
197 }
198
199 /**
200 @TODO what about conversion factors? Is an allowance being made?
201 */
202 void
203 Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
204 if(samplelist!=NULL){
205 ASC_FREE(samplelist);
206 }
207 const dim_type *d = units.getDimensions().getInternalType();
208
209 if(start<=0)throw runtime_error("starting timestep needs to be > 0");
210 if(end<=0)throw runtime_error("end timestep needs to be > 0");
211 if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
212
213 samplelist = samplelist_new(num+1, d);
214 double val = start;
215 double inc = exp((log(end)-log(start))/num);
216 for(unsigned long i=0;i<=num;++i){
217 samplelist_set(samplelist,i,val);
218 CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
219 val *= inc;
220 }
221 integrator_set_samples(blsys,samplelist);
222 }
223
224 vector<double>
225 Integrator::getCurrentObservations(){
226 double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
227 integrator_get_observations(blsys,d);
228 vector<double> v=vector<double>(d,d+getNumObservedVars());
229 // do I need to free d?
230 // can I do this in such a way as I avoid all this memory-copying?
231 return v;
232 }
233
234 Variable
235 Integrator::getObservedVariable(const long &i){
236 var_variable *v = integrator_get_observed_var(blsys,i);
237 return Variable(&simulation,v);
238 }
239
240 Variable
241 Integrator::getIndependentVariable(){
242 var_variable *v = integrator_get_independent_var(blsys);
243 if(v==NULL){
244 throw runtime_error("independent variable is null");
245 }
246 return Variable(&simulation,v);
247 }
248
249 int
250 Integrator::getNumVars(){
251 return blsys->n_y;
252 }
253
254 int
255 Integrator::getNumObservedVars(){
256 return blsys->n_obs;
257 }
258
259 void
260 Integrator::setMinSubStep(double n){
261 integrator_set_minstep(blsys,n);
262 }
263
264 void
265 Integrator::setMaxSubStep(double n){
266 integrator_set_maxstep(blsys,n);
267 }
268
269 void
270 Integrator::setInitialSubStep(double n){
271 integrator_set_stepzero(blsys,n);
272 }
273
274 void
275 Integrator::setMaxSubSteps(int n){
276 integrator_set_maxsubsteps(blsys,n);
277 }
278
279 IntegratorSystem *
280 Integrator::getInternalType(){
281 return blsys;
282 }

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