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

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 941 - (show annotations) (download) (as text)
Fri Nov 24 10:46:32 2006 UTC (13 years, 11 months ago) by johnpye
File MIME type: text/x-c++src
File size: 5166 byte(s)
Changed integrator_set_engine to return 0 on success.
Fixed Integrator::setEngine to throw range_error / IndexError on invalid selection.
Test suite contains testIDA that works now (more tests yet to come)
1 #include "integrator.h"
2 #include "integratorreporter.h"
3 #include <stdexcept>
4 #include <sstream>
5 using namespace std;
6
7 /**
8 'creating' an integrator in the context of the GUI just means an object
9 we can store the parameters that will be later sent to the underlying
10 C-code API.
11 */
12 Integrator::Integrator(Simulation &simulation)
13 : simulation(simulation)
14 {
15 // create the C-level object
16 this->blsys = integrator_new(simulation.getSystem(),simulation.getModel().getInternalType());
17
18 samplelist = NULL;
19
20 // set default steps
21 setMinSubStep(0);
22 setMaxSubStep(0);
23 setMaxSubSteps(0);
24 setInitialSubStep(0);
25 }
26
27 Integrator::~Integrator(){
28 CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
29 CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
30 integrator_free(blsys);
31 CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
32 samplelist_free(samplelist);
33 }
34
35
36 void
37 Integrator::setReporter(IntegratorReporterCxx *reporter){
38 this->blsys->clientdata = reporter;
39 integrator_set_reporter(blsys,reporter->getInternalType());
40 CONSOLE_DEBUG("REPORTER HAS BEEN SET");
41 (*(this->blsys->reporter->init))(blsys);
42 CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
43 }
44
45 double
46 Integrator::getCurrentTime(){
47 return integrator_get_t(blsys);
48 }
49
50 long
51 Integrator::getCurrentStep(){
52 return integrator_getcurrentstep(blsys);
53 }
54
55 long
56 Integrator::getNumSteps(){
57 return integrator_getnsamples(blsys);
58 }
59
60 int
61 Integrator::findIndependentVar(){
62 return integrator_find_indep_var(blsys);
63 }
64
65 int
66 Integrator::analyse(){
67
68 int res;
69 res = integrator_analyse(blsys);
70
71 if(!res){
72 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Failed system analysis");
73 return 0;
74 }
75
76 return 1;
77 }
78
79 /**
80 @TODO what about root detection?
81 */
82 int
83 Integrator::solve(){
84
85 // check the integration limits
86 // trigger of the solution process
87 // report errors?
88
89 assert(samplelist!=NULL);
90 assert(samplelist->ns>0);
91 assert(blsys->reporter!=NULL);
92 assert(blsys->clientdata!=NULL);
93
94 int res;
95 res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
96
97 // communicate solver variable status back to the instance tree via 'interface_ptr'
98 simulation.processVarStatus();
99
100 if(!res){
101 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Failed integration");
102 return 0;
103 }
104
105 return 1;
106 }
107
108 void
109 Integrator::setEngine(IntegratorEngine engine){
110 int res = integrator_set_engine(this->blsys, engine);
111 if(!res)return;
112 if(res==1)throw range_error("Unknown integrator");
113 if(res==2)throw range_error("Invalid integrator");
114 stringstream ss;
115 ss << "Unknown error in setEngine (res = " << res << ")";
116 throw runtime_error(ss.str());
117 }
118
119 void
120 Integrator::setEngine(int engine){
121 setEngine((IntegratorEngine)engine);
122 }
123
124 void
125 Integrator::setEngine(const string &name){
126 CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
127 IntegratorEngine engine = INTEG_UNKNOWN;
128 #ifdef ASC_WITH_LSODE
129 if(name=="LSODE")engine = INTEG_LSODE;
130 #endif
131 #ifdef ASC_WITH_IDA
132 if(name=="IDA")engine = INTEG_IDA;
133 #endif
134
135 setEngine(engine);
136 }
137
138 /**
139 Ideally this list would be dynamically generated based on what solvers
140 are available or are in memory.
141 */
142 map<int,string>
143 Integrator::getEngines(){
144 map<int,string> m;
145 #ifdef ASC_WITH_LSODE
146 m.insert(pair<int,string>(INTEG_LSODE,"LSODE"));
147 #endif
148 #ifdef ASC_WITH_IDA
149 m.insert(pair<int,string>(INTEG_IDA,"IDA"));
150 #endif
151 return m;
152 }
153
154 string
155 Integrator::getEngineName() const{
156 map<int,string> m=getEngines();
157 map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
158 if(f==m.end()){
159 throw runtime_error("No engine selected");
160 }
161 return f->second;
162 }
163
164 void
165 Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
166 if(samplelist!=NULL){
167 ASC_FREE(samplelist);
168 }
169 const dim_type *d = units.getDimensions().getInternalType();
170 samplelist = samplelist_new(num+1, d);
171 double val = start;
172 double inc = (end-start)/(num);
173 for(unsigned long i=0;i<=num;++i){
174 samplelist_set(samplelist,i,val);
175 val += inc;
176 }
177 integrator_set_samples(blsys,samplelist);
178 }
179
180 vector<double>
181 Integrator::getCurrentObservations(){
182 double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
183 integrator_get_observations(blsys,d);
184 vector<double> v=vector<double>(d,d+getNumObservedVars());
185 // do I need to free d?
186 // can I do this in such a way as I avoid all this memory-copying?
187 return v;
188 }
189
190 Variable
191 Integrator::getObservedVariable(const long &i){
192 var_variable *v = integrator_get_observed_var(blsys,i);
193 return Variable(&simulation,v);
194 }
195
196 Variable
197 Integrator::getIndependentVariable(){
198 var_variable *v = integrator_get_independent_var(blsys);
199 if(v==NULL){
200 throw runtime_error("independent variable is null");
201 }
202 return Variable(&simulation,v);
203 }
204
205 int
206 Integrator::getNumVars(){
207 return blsys->n_y;
208 }
209
210 int
211 Integrator::getNumObservedVars(){
212 return blsys->n_obs;
213 }
214
215 void
216 Integrator::setMinSubStep(double n){
217 integrator_set_minstep(blsys,n);
218 }
219
220 void
221 Integrator::setMaxSubStep(double n){
222 integrator_set_maxstep(blsys,n);
223 }
224
225 void
226 Integrator::setInitialSubStep(double n){
227 integrator_set_stepzero(blsys,n);
228 }
229
230 void
231 Integrator::setMaxSubSteps(int n){
232 integrator_set_maxsubsteps(blsys,n);
233 }
234
235 IntegratorSystem *
236 Integrator::getInternalType(){
237 return blsys;
238 }

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