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

Annotation of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 966 - (hide annotations) (download) (as text)
Thu Dec 14 14:04:54 2006 UTC (15 years, 8 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 johnpye 669 #include "integrator.h"
2     #include "integratorreporter.h"
3 johnpye 942 #include "solverparameters.h"
4 johnpye 669 #include <stdexcept>
5 johnpye 941 #include <sstream>
6 johnpye 942 #include <cmath>
7 johnpye 669 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 johnpye 921 CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
31     CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
32 johnpye 669 integrator_free(blsys);
33 johnpye 921 CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
34 johnpye 669 samplelist_free(samplelist);
35     }
36    
37 johnpye 942 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 johnpye 669
45     void
46 johnpye 942 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 johnpye 669 Integrator::setReporter(IntegratorReporterCxx *reporter){
53     this->blsys->clientdata = reporter;
54     integrator_set_reporter(blsys,reporter->getInternalType());
55 johnpye 966 //CONSOLE_DEBUG("REPORTER HAS BEEN SET");
56 johnpye 669 (*(this->blsys->reporter->init))(blsys);
57 johnpye 966 //CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
58 johnpye 669 }
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 johnpye 942
97     Integrate the function for the timesteps specified.
98    
99     This method will throw a runtime_error if integration fails.
100 johnpye 953
101     @TODO does simulation.processVarStatus work for integrators like IDA???
102 johnpye 669 */
103 johnpye 942 void
104 johnpye 669 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 johnpye 915
118 johnpye 669 if(!res){
119 johnpye 942 throw runtime_error("Failed integration");
120 johnpye 669 }
121 johnpye 953
122     // communicate solver variable status back to the instance tree via 'interface_ptr'
123     simulation.processVarStatus();
124 johnpye 669 }
125    
126 johnpye 941 void
127 johnpye 669 Integrator::setEngine(IntegratorEngine engine){
128 johnpye 941 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 johnpye 669 }
136    
137 johnpye 941 void
138 johnpye 669 Integrator::setEngine(int engine){
139 johnpye 941 setEngine((IntegratorEngine)engine);
140 johnpye 669 }
141    
142 johnpye 941 void
143 johnpye 938 Integrator::setEngine(const string &name){
144 johnpye 941 CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
145 johnpye 938 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 johnpye 952 if(engine==INTEG_UNKNOWN){
153     throw runtime_error("Unkown integrator name");
154     }
155 johnpye 941 setEngine(engine);
156 johnpye 938 }
157    
158 johnpye 669 /**
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 johnpye 941 Integrator::getEngines(){
164 johnpye 669 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 johnpye 942 /**
185     @TODO what about conversion factors? Is an allowance being made?
186     */
187 johnpye 669 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 johnpye 746 for(unsigned long i=0;i<=num;++i){
197 johnpye 669 samplelist_set(samplelist,i,val);
198     val += inc;
199     }
200     integrator_set_samples(blsys,samplelist);
201     }
202    
203 johnpye 942 /**
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 johnpye 669 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 johnpye 709 return Variable(&simulation,v);
242 johnpye 669 }
243    
244 johnpye 854 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 johnpye 669 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