trunk/pygtk/integrator.cpp

#include "integrator.h"
#include "integratorreporter.h"
#include "solverparameters.h"
#include <stdexcept>
#include <sstream>
#include <cmath>
using namespace std;

// #define DESTROY_DEBUG

/**
    'creating' an integrator in the context of the GUI just means an object
    we can store the parameters that will be later sent to the underlying
    C-code API.

    @TODO at present the integrator requires a slv_system_t. This is the wrong
    way around.
*/
Integrator::Integrator(Simulation &simulation)
        : simulation(simulation)
{
    // create the C-level object
    this->blsys = integrator_new(simulation.getSystem(),simulation.getModel().getInternalType());

    samplelist = NULL;

    // set default steps
    setMinSubStep(0);
    setMaxSubStep(0);
    setMaxSubSteps(0);
    setInitialSubStep(0);
}

Integrator::~Integrator(){
#ifdef DESTROY_DEBUG
    CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
    CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
#endif
    integrator_free(blsys);
#ifdef DESTROY_DEBUG
    CONSOLE_DEBUG("done");
    CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
#endif
    if(samplelist)samplelist_free(samplelist);
}

SolverParameters
Integrator::getParameters() const{
    SolverParameters params;
    int res = integrator_params_get(blsys,&(params.getInternalType() ) );
    if(res)throw runtime_error("Failed to get integrator parameters");
    return params;
}

void
Integrator::setParameters(const SolverParameters &params){
    int res = integrator_params_set(blsys,&(params.getInternalTypeConst() ) );
    if(res)throw runtime_error("Failed to set integrator parameters");
}

void
Integrator::setReporter(IntegratorReporterCxx *reporter){
    this->blsys->clientdata = reporter;
    integrator_set_reporter(blsys,reporter->getInternalType());
    //CONSOLE_DEBUG("REPORTER HAS BEEN SET");
    (*(this->blsys->reporter->init))(blsys);
    //CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
}

double
Integrator::getCurrentTime(){
    return integrator_get_t(blsys);
}

long
Integrator::getCurrentStep(){
    return integrator_getcurrentstep(blsys);
}

long
Integrator::getNumSteps(){
    return integrator_getnsamples(blsys);
}

/**
    Find the independent variable in the system, or throw an exception if not found.
*/
void
Integrator::findIndependentVar(){
    int res = integrator_find_indep_var(blsys);

    if(res){
        stringstream ss;
        ss << "Independent variable not found (" << res << ")";
        throw runtime_error(ss.str());
    }
}

void
Integrator::analyse(){

    int res;
    /*
        Note, we never need to call analyze_make_system in any of the Integrator
        code, as it gets called by Simulation::build.
    */
    res = integrator_analyse(blsys);

    if(res){
        stringstream ss; ss << "Failed system analysis (error " << res << ")";
        throw runtime_error(ss.str());
    }
}

/**
    @TODO what about root detection?

    Integrate the function for the timesteps specified.

    Method will throw a runtime_error if integrator_solve returns error (non zero)

    @TODO does simulation.processVarStatus work for integrators like IDA???
*/
void
Integrator::solve(){

    // check the integration limits
    // trigger of the solution process
    // report errors?

    assert(samplelist!=NULL);
    assert(samplelist->ns>0);
    assert(blsys->reporter!=NULL);
    assert(blsys->clientdata!=NULL);

    int res;
    res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);

    if(res){
        stringstream ss;
        ss << "Failed integration (integrator_solve returned " << res << ")";
        throw runtime_error(ss.str());
    }

    // communicate solver variable status back to the instance tree via 'interface_ptr'
    simulation.processVarStatus();
}

void
Integrator::writeMatrix(FILE *fp,const char *type) const{
    if(integrator_write_matrix(this->blsys, fp, type)){
        throw runtime_error("Failed to write matrix");
    }
}

void
Integrator::writeDebug(FILE *fp) const{
    if(integrator_debug(this->blsys, fp)){
        throw runtime_error("Failed to write debug output");
    }
}

void
Integrator::setEngine(IntegratorEngine engine){
    int res = integrator_set_engine(this->blsys, engine);
    if(!res)return;
    if(res==1)throw range_error("Unknown integrator");
    if(res==2)throw range_error("Invalid integrator");
    stringstream ss;
    ss << "Unknown error in setEngine (res = " << res << ")";
    throw runtime_error(ss.str());
}

void
Integrator::setEngine(int engine){
    setEngine((IntegratorEngine)engine);
}

void
Integrator::setEngine(const string &name){
    CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
    IntegratorEngine engine = INTEG_UNKNOWN;
#ifdef ASC_WITH_LSODE
    if(name=="LSODE")engine = INTEG_LSODE;
#endif
#ifdef ASC_WITH_IDA
    if(name=="IDA")engine = INTEG_IDA;
#endif
    if(engine==INTEG_UNKNOWN){
        throw runtime_error("Unkown integrator name");
    }
    setEngine(engine);
}

/**
    Ideally this list would be dynamically generated based on what solvers
    are available or are in memory.
*/
map<int,string>
Integrator::getEngines(){
    map<int,string> m;
    const IntegratorLookup *list = integrator_get_engines();
    while(list->id != INTEG_UNKNOWN){
        if(list->name==NULL)throw runtime_error("list->name is NULL");
        m.insert(pair<int,string>(list->id,list->name));
        ++list;
    }
    return m;
}

string
Integrator::getName() const{
    map<int,string> m=getEngines();
    map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
    if(f==m.end()){
        throw runtime_error("No engine selected");
    }
    return f->second;
}

/**
    @TODO what about conversion factors? Is an allowance being made?
*/
void
Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
    if(samplelist!=NULL){
        ASC_FREE(samplelist);
    }
    const dim_type *d = units.getDimensions().getInternalType();
    samplelist = samplelist_new(num+1, d);
    double val = start;
    double inc = (end-start)/(num);
    for(unsigned long i=0;i<=num;++i){
        samplelist_set(samplelist,i,val);
        val += inc;
    }
    integrator_set_samples(blsys,samplelist);
}

/**
    @TODO what about conversion factors? Is an allowance being made?
*/
void
Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
    if(samplelist!=NULL){
        ASC_FREE(samplelist);
    }
    const dim_type *d = units.getDimensions().getInternalType();

    if(start<=0)throw runtime_error("starting timestep needs to be > 0");
    if(end<=0)throw runtime_error("end timestep needs to be > 0");
    if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");

    samplelist = samplelist_new(num+1, d);
    double val = start;
    double inc = exp((log(end)-log(start))/num);
    for(unsigned long i=0;i<=num;++i){
        samplelist_set(samplelist,i,val);
        CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
        val *= inc;
    }
    integrator_set_samples(blsys,samplelist);
}

vector<double>
Integrator::getCurrentObservations(){
    double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
    integrator_get_observations(blsys,d);
    vector<double> v=vector<double>(d,d+getNumObservedVars());
    // do I need to free d?
    // can I do this in such a way as I avoid all this memory-copying?
    return v;
}

Variable
Integrator::getObservedVariable(const long &i){
    var_variable *v = integrator_get_observed_var(blsys,i);
    return Variable(&simulation,v);
}

Variable
Integrator::getIndependentVariable(){
    var_variable *v = integrator_get_independent_var(blsys);
    if(v==NULL){
        throw runtime_error("independent variable is null");
    }
    return Variable(&simulation,v);
}

int
Integrator::getNumVars(){
    return blsys->n_y;
}

int
Integrator::getNumObservedVars(){
    return blsys->n_obs;
}

void
Integrator::setMinSubStep(double n){
    integrator_set_minstep(blsys,n);
}

void
Integrator::setMaxSubStep(double n){
    integrator_set_maxstep(blsys,n);
}

void
Integrator::setInitialSubStep(double n){
    integrator_set_stepzero(blsys,n);
}

void
Integrator::setMaxSubSteps(int n){
    integrator_set_maxsubsteps(blsys,n);
}

IntegratorSystem *
Integrator::getInternalType(){
    return blsys;
}
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	johnpye	1303	// #define DESTROY_DEBUG
10
11	johnpye	669	/**
12			'creating' an integrator in the context of the GUI just means an object
13			we can store the parameters that will be later sent to the underlying
14			C-code API.
15	johnpye	1199
16			@TODO at present the integrator requires a slv_system_t. This is the wrong
17			way around.
18	johnpye	669	*/
19			Integrator::Integrator(Simulation &simulation)
20			: simulation(simulation)
21	johnpye	1087	{
22	johnpye	669	// create the C-level object
23			this->blsys = integrator_new(simulation.getSystem(),simulation.getModel().getInternalType());
24
25			samplelist = NULL;
26
27			// set default steps
28			setMinSubStep(0);
29			setMaxSubStep(0);
30			setMaxSubSteps(0);
31			setInitialSubStep(0);
32			}
33
34			Integrator::~Integrator(){
35	johnpye	1303	#ifdef DESTROY_DEBUG
36	johnpye	921	CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
37			CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
38	johnpye	1303	#endif
39	johnpye	669	integrator_free(blsys);
40	johnpye	1303	#ifdef DESTROY_DEBUG
41	johnpye	1221	CONSOLE_DEBUG("done");
42	johnpye	921	CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
43	johnpye	1303	#endif
44	johnpye	1201	if(samplelist)samplelist_free(samplelist);
45	johnpye	669	}
46
47	johnpye	1087	SolverParameters
48	johnpye	942	Integrator::getParameters() const{
49			SolverParameters params;
50			int res = integrator_params_get(blsys,&(params.getInternalType() ) );
51			if(res)throw runtime_error("Failed to get integrator parameters");
52			return params;
53			}
54	johnpye	669
55			void
56	johnpye	942	Integrator::setParameters(const SolverParameters &params){
57			int res = integrator_params_set(blsys,&(params.getInternalTypeConst() ) );
58			if(res)throw runtime_error("Failed to set integrator parameters");
59			}
60
61			void
62	johnpye	669	Integrator::setReporter(IntegratorReporterCxx *reporter){
63			this->blsys->clientdata = reporter;
64			integrator_set_reporter(blsys,reporter->getInternalType());
65	johnpye	1087	//CONSOLE_DEBUG("REPORTER HAS BEEN SET");
66	johnpye	669	(*(this->blsys->reporter->init))(blsys);
67	johnpye	966	//CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
68	johnpye	669	}
69
70			double
71			Integrator::getCurrentTime(){
72			return integrator_get_t(blsys);
73			}
74
75			long
76			Integrator::getCurrentStep(){
77			return integrator_getcurrentstep(blsys);
78			}
79
80			long
81			Integrator::getNumSteps(){
82			return integrator_getnsamples(blsys);
83			}
84
85	johnpye	1087	/**
86			Find the independent variable in the system, or throw an exception if not found.
87			*/
88			void
89	johnpye	669	Integrator::findIndependentVar(){
90	johnpye	1087	int res = integrator_find_indep_var(blsys);
91
92			if(res){
93			stringstream ss;
94			ss << "Independent variable not found (" << res << ")";
95			throw runtime_error(ss.str());
96			}
97	johnpye	669	}
98
99	johnpye	979	void
100	johnpye	669	Integrator::analyse(){
101
102			int res;
103	johnpye	1069	/*
104			Note, we never need to call analyze_make_system in any of the Integrator
105			code, as it gets called by Simulation::build.
106			*/
107	johnpye	669	res = integrator_analyse(blsys);
108
109	johnpye	1069	if(res){
110	johnpye	1097	stringstream ss; ss << "Failed system analysis (error " << res << ")";
111			throw runtime_error(ss.str());
112	johnpye	669	}
113			}
114
115			/**
116			@TODO what about root detection?
117	johnpye	942
118			Integrate the function for the timesteps specified.
119
120	johnpye	1049	Method will throw a runtime_error if integrator_solve returns error (non zero)
121	johnpye	953
122			@TODO does simulation.processVarStatus work for integrators like IDA???
123	johnpye	669	*/
124	johnpye	942	void
125	johnpye	669	Integrator::solve(){
126
127			// check the integration limits
128			// trigger of the solution process
129			// report errors?
130	johnpye	1087
131	johnpye	669	assert(samplelist!=NULL);
132			assert(samplelist->ns>0);
133			assert(blsys->reporter!=NULL);
134			assert(blsys->clientdata!=NULL);
135
136			int res;
137			res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
138	johnpye	915
139	johnpye	1049	if(res){
140			stringstream ss;
141			ss << "Failed integration (integrator_solve returned " << res << ")";
142			throw runtime_error(ss.str());
143	johnpye	669	}
144	johnpye	953
145			// communicate solver variable status back to the instance tree via 'interface_ptr'
146			simulation.processVarStatus();
147	johnpye	669	}
148
149	johnpye	941	void
150	johnpye	1306	Integrator::writeMatrix(FILE fp,const char type) const{
151			if(integrator_write_matrix(this->blsys, fp, type)){
152	johnpye	1129	throw runtime_error("Failed to write matrix");
153			}
154			}
155
156			void
157	johnpye	1196	Integrator::writeDebug(FILE *fp) const{
158			if(integrator_debug(this->blsys, fp)){
159			throw runtime_error("Failed to write debug output");
160			}
161			}
162
163			void
164	johnpye	669	Integrator::setEngine(IntegratorEngine engine){
165	johnpye	941	int res = integrator_set_engine(this->blsys, engine);
166			if(!res)return;
167			if(res==1)throw range_error("Unknown integrator");
168			if(res==2)throw range_error("Invalid integrator");
169			stringstream ss;
170			ss << "Unknown error in setEngine (res = " << res << ")";
171			throw runtime_error(ss.str());
172	johnpye	669	}
173
174	johnpye	941	void
175	johnpye	669	Integrator::setEngine(int engine){
176	johnpye	941	setEngine((IntegratorEngine)engine);
177	johnpye	669	}
178
179	johnpye	941	void
180	johnpye	938	Integrator::setEngine(const string &name){
181	johnpye	941	CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
182	johnpye	938	IntegratorEngine engine = INTEG_UNKNOWN;
183			#ifdef ASC_WITH_LSODE
184			if(name=="LSODE")engine = INTEG_LSODE;
185			#endif
186			#ifdef ASC_WITH_IDA
187			if(name=="IDA")engine = INTEG_IDA;
188			#endif
189	johnpye	952	if(engine==INTEG_UNKNOWN){
190			throw runtime_error("Unkown integrator name");
191			}
192	johnpye	941	setEngine(engine);
193	johnpye	938	}
194
195	johnpye	669	/**
196			Ideally this list would be dynamically generated based on what solvers
197			are available or are in memory.
198			*/
199			map<int,string>
200	johnpye	941	Integrator::getEngines(){
201	johnpye	669	map<int,string> m;
202	johnpye	977	const IntegratorLookup *list = integrator_get_engines();
203	johnpye	969	while(list->id != INTEG_UNKNOWN){
204	johnpye	1071	if(list->name==NULL)throw runtime_error("list->name is NULL");
205	johnpye	969	m.insert(pair<int,string>(list->id,list->name));
206			++list;
207			}
208	johnpye	669	return m;
209			}
210
211			string
212	johnpye	1140	Integrator::getName() const{
213	johnpye	669	map<int,string> m=getEngines();
214			map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
215			if(f==m.end()){
216			throw runtime_error("No engine selected");
217			}
218			return f->second;
219	johnpye	1087	}
220	johnpye	669
221	johnpye	942	/**
222			@TODO what about conversion factors? Is an allowance being made?
223			*/
224	johnpye	669	void
225			Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
226			if(samplelist!=NULL){
227			ASC_FREE(samplelist);
228			}
229			const dim_type *d = units.getDimensions().getInternalType();
230			samplelist = samplelist_new(num+1, d);
231			double val = start;
232			double inc = (end-start)/(num);
233	johnpye	746	for(unsigned long i=0;i<=num;++i){
234	johnpye	669	samplelist_set(samplelist,i,val);
235			val += inc;
236			}
237			integrator_set_samples(blsys,samplelist);
238			}
239
240	johnpye	942	/**
241			@TODO what about conversion factors? Is an allowance being made?
242			*/
243			void
244			Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
245			if(samplelist!=NULL){
246			ASC_FREE(samplelist);
247			}
248			const dim_type *d = units.getDimensions().getInternalType();
249
250			if(start<=0)throw runtime_error("starting timestep needs to be > 0");
251			if(end<=0)throw runtime_error("end timestep needs to be > 0");
252			if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
253
254			samplelist = samplelist_new(num+1, d);
255			double val = start;
256			double inc = exp((log(end)-log(start))/num);
257			for(unsigned long i=0;i<=num;++i){
258			samplelist_set(samplelist,i,val);
259			CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
260			val *= inc;
261			}
262			integrator_set_samples(blsys,samplelist);
263			}
264
265	johnpye	669	vector<double>
266			Integrator::getCurrentObservations(){
267			double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
268			integrator_get_observations(blsys,d);
269			vector<double> v=vector<double>(d,d+getNumObservedVars());
270			// do I need to free d?
271			// can I do this in such a way as I avoid all this memory-copying?
272			return v;
273			}
274
275			Variable
276			Integrator::getObservedVariable(const long &i){
277			var_variable *v = integrator_get_observed_var(blsys,i);
278	johnpye	709	return Variable(&simulation,v);
279	johnpye	669	}
280
281	johnpye	854	Variable
282			Integrator::getIndependentVariable(){
283			var_variable *v = integrator_get_independent_var(blsys);
284			if(v==NULL){
285			throw runtime_error("independent variable is null");
286			}
287			return Variable(&simulation,v);
288			}
289
290	johnpye	669	int
291			Integrator::getNumVars(){
292			return blsys->n_y;
293			}
294
295	johnpye	1087	int
296	johnpye	669	Integrator::getNumObservedVars(){
297			return blsys->n_obs;
298			}
299
300			void
301			Integrator::setMinSubStep(double n){
302			integrator_set_minstep(blsys,n);
303			}
304
305			void
306			Integrator::setMaxSubStep(double n){
307			integrator_set_maxstep(blsys,n);
308			}
309
310			void
311			Integrator::setInitialSubStep(double n){
312			integrator_set_stepzero(blsys,n);
313			}
314
315			void
316			Integrator::setMaxSubSteps(int n){
317			integrator_set_maxsubsteps(blsys,n);
318			}
319
320			IntegratorSystem *
321			Integrator::getInternalType(){
322			return blsys;
323			}