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{
    if(integrator_write_matrix(this->blsys, fp)){
        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	#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	// #define DESTROY_DEBUG
10
11	/**
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
16	@TODO at present the integrator requires a slv_system_t. This is the wrong
17	way around.
18	*/
19	Integrator::Integrator(Simulation &simulation)
20	: simulation(simulation)
21	{
22	// 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	#ifdef DESTROY_DEBUG
36	CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
37	CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
38	#endif
39	integrator_free(blsys);
40	#ifdef DESTROY_DEBUG
41	CONSOLE_DEBUG("done");
42	CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
43	#endif
44	if(samplelist)samplelist_free(samplelist);
45	}
46
47	SolverParameters
48	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
55	void
56	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	Integrator::setReporter(IntegratorReporterCxx *reporter){
63	this->blsys->clientdata = reporter;
64	integrator_set_reporter(blsys,reporter->getInternalType());
65	//CONSOLE_DEBUG("REPORTER HAS BEEN SET");
66	(*(this->blsys->reporter->init))(blsys);
67	//CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
68	}
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	/**
86	Find the independent variable in the system, or throw an exception if not found.
87	*/
88	void
89	Integrator::findIndependentVar(){
90	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	}
98
99	void
100	Integrator::analyse(){
101
102	int res;
103	/*
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	res = integrator_analyse(blsys);
108
109	if(res){
110	stringstream ss; ss << "Failed system analysis (error " << res << ")";
111	throw runtime_error(ss.str());
112	}
113	}
114
115	/**
116	@TODO what about root detection?
117
118	Integrate the function for the timesteps specified.
119
120	Method will throw a runtime_error if integrator_solve returns error (non zero)
121
122	@TODO does simulation.processVarStatus work for integrators like IDA???
123	*/
124	void
125	Integrator::solve(){
126
127	// check the integration limits
128	// trigger of the solution process
129	// report errors?
130
131	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
139	if(res){
140	stringstream ss;
141	ss << "Failed integration (integrator_solve returned " << res << ")";
142	throw runtime_error(ss.str());
143	}
144
145	// communicate solver variable status back to the instance tree via 'interface_ptr'
146	simulation.processVarStatus();
147	}
148
149	void
150	Integrator::writeMatrix(FILE *fp) const{
151	if(integrator_write_matrix(this->blsys, fp)){
152	throw runtime_error("Failed to write matrix");
153	}
154	}
155
156	void
157	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	Integrator::setEngine(IntegratorEngine engine){
165	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	}
173
174	void
175	Integrator::setEngine(int engine){
176	setEngine((IntegratorEngine)engine);
177	}
178
179	void
180	Integrator::setEngine(const string &name){
181	CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
182	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	if(engine==INTEG_UNKNOWN){
190	throw runtime_error("Unkown integrator name");
191	}
192	setEngine(engine);
193	}
194
195	/**
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	Integrator::getEngines(){
201	map<int,string> m;
202	const IntegratorLookup *list = integrator_get_engines();
203	while(list->id != INTEG_UNKNOWN){
204	if(list->name==NULL)throw runtime_error("list->name is NULL");
205	m.insert(pair<int,string>(list->id,list->name));
206	++list;
207	}
208	return m;
209	}
210
211	string
212	Integrator::getName() const{
213	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	}
220
221	/**
222	@TODO what about conversion factors? Is an allowance being made?
223	*/
224	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	for(unsigned long i=0;i<=num;++i){
234	samplelist_set(samplelist,i,val);
235	val += inc;
236	}
237	integrator_set_samples(blsys,samplelist);
238	}
239
240	/**
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	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	return Variable(&simulation,v);
279	}
280
281	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	int
291	Integrator::getNumVars(){
292	return blsys->n_y;
293	}
294
295	int
296	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	}