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

Contents of /trunk/pygtk/integrator.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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