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);
    CONSOLE_DEBUG("Got %d",res);

    if(res){
        CONSOLE_DEBUG("...which is bad");
        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	#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	CONSOLE_DEBUG("Got %d",res);
109
110	if(res){
111	CONSOLE_DEBUG("...which is bad");
112	stringstream ss;
113	ss << "Failed system analysis (error " << res << ")";
114	throw runtime_error(ss.str());
115	}
116	}
117
118	/**
119	@TODO what about root detection?
120
121	Integrate the function for the timesteps specified.
122
123	Method will throw a runtime_error if integrator_solve returns error (non zero)
124
125	@TODO does simulation.processVarStatus work for integrators like IDA???
126	*/
127	void
128	Integrator::solve(){
129
130	// check the integration limits
131	// trigger of the solution process
132	// report errors?
133
134	assert(samplelist!=NULL);
135	assert(samplelist->ns>0);
136	assert(blsys->reporter!=NULL);
137	assert(blsys->clientdata!=NULL);
138
139	int res;
140	res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
141
142	if(res){
143	stringstream ss;
144	ss << "Failed integration (integrator_solve returned " << res << ")";
145	throw runtime_error(ss.str());
146	}
147
148	// communicate solver variable status back to the instance tree via 'interface_ptr'
149	simulation.processVarStatus();
150	}
151
152	void
153	Integrator::writeMatrix(FILE fp,const char type) const{
154	if(integrator_write_matrix(this->blsys, fp, type)){
155	throw runtime_error("Failed to write matrix");
156	}
157	}
158
159	void
160	Integrator::writeDebug(FILE *fp) const{
161	if(integrator_debug(this->blsys, fp)){
162	throw runtime_error("Failed to write debug output");
163	}
164	}
165
166	void
167	Integrator::setEngine(IntegratorEngine engine){
168	int res = integrator_set_engine(this->blsys, engine);
169	if(!res)return;
170	if(res==1)throw range_error("Unknown integrator");
171	if(res==2)throw range_error("Invalid integrator");
172	stringstream ss;
173	ss << "Unknown error in setEngine (res = " << res << ")";
174	throw runtime_error(ss.str());
175	}
176
177	void
178	Integrator::setEngine(int engine){
179	setEngine((IntegratorEngine)engine);
180	}
181
182	void
183	Integrator::setEngine(const string &name){
184	CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
185	IntegratorEngine engine = INTEG_UNKNOWN;
186	#ifdef ASC_WITH_LSODE
187	if(name=="LSODE")engine = INTEG_LSODE;
188	#endif
189	#ifdef ASC_WITH_IDA
190	if(name=="IDA")engine = INTEG_IDA;
191	#endif
192	if(engine==INTEG_UNKNOWN){
193	throw runtime_error("Unkown integrator name");
194	}
195	setEngine(engine);
196	}
197
198	/**
199	Ideally this list would be dynamically generated based on what solvers
200	are available or are in memory.
201	*/
202	map<int,string>
203	Integrator::getEngines(){
204	map<int,string> m;
205	const IntegratorLookup *list = integrator_get_engines();
206	while(list->id != INTEG_UNKNOWN){
207	if(list->name==NULL)throw runtime_error("list->name is NULL");
208	m.insert(pair<int,string>(list->id,list->name));
209	++list;
210	}
211	return m;
212	}
213
214	string
215	Integrator::getName() const{
216	map<int,string> m=getEngines();
217	map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
218	if(f==m.end()){
219	throw runtime_error("No engine selected");
220	}
221	return f->second;
222	}
223
224	/**
225	@TODO what about conversion factors? Is an allowance being made?
226	*/
227	void
228	Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
229	if(samplelist!=NULL){
230	ASC_FREE(samplelist);
231	}
232	const dim_type *d = units.getDimensions().getInternalType();
233	samplelist = samplelist_new(num+1, d);
234	double val = start;
235	double inc = (end-start)/(num);
236	for(unsigned long i=0;i<=num;++i){
237	samplelist_set(samplelist,i,val);
238	val += inc;
239	}
240	integrator_set_samples(blsys,samplelist);
241	}
242
243	/**
244	@TODO what about conversion factors? Is an allowance being made?
245	*/
246	void
247	Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
248	if(samplelist!=NULL){
249	ASC_FREE(samplelist);
250	}
251	const dim_type *d = units.getDimensions().getInternalType();
252
253	if(start<=0)throw runtime_error("starting timestep needs to be > 0");
254	if(end<=0)throw runtime_error("end timestep needs to be > 0");
255	if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
256
257	samplelist = samplelist_new(num+1, d);
258	double val = start;
259	double inc = exp((log(end)-log(start))/num);
260	for(unsigned long i=0;i<=num;++i){
261	samplelist_set(samplelist,i,val);
262	// CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
263	val *= inc;
264	}
265	integrator_set_samples(blsys,samplelist);
266	}
267
268	vector<double>
269	Integrator::getCurrentObservations(){
270	double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
271	integrator_get_observations(blsys,d);
272	vector<double> v=vector<double>(d,d+getNumObservedVars());
273	// do I need to free d?
274	// can I do this in such a way as I avoid all this memory-copying?
275	return v;
276	}
277
278	Variable
279	Integrator::getObservedVariable(const long &i){
280	var_variable *v = integrator_get_observed_var(blsys,i);
281	return Variable(&simulation,v);
282	}
283
284	Variable
285	Integrator::getIndependentVariable(){
286	var_variable *v = integrator_get_independent_var(blsys);
287	if(v==NULL){
288	throw runtime_error("independent variable is null");
289	}
290	return Variable(&simulation,v);
291	}
292
293	int
294	Integrator::getNumVars(){
295	return blsys->n_y;
296	}
297
298	int
299	Integrator::getNumObservedVars(){
300	return blsys->n_obs;
301	}
302
303	void
304	Integrator::setMinSubStep(double n){
305	integrator_set_minstep(blsys,n);
306	}
307
308	void
309	Integrator::setMaxSubStep(double n){
310	integrator_set_maxstep(blsys,n);
311	}
312
313	void
314	Integrator::setInitialSubStep(double n){
315	integrator_set_stepzero(blsys,n);
316	}
317
318	void
319	Integrator::setMaxSubSteps(int n){
320	integrator_set_maxsubsteps(blsys,n);
321	}
322
323	IntegratorSystem *
324	Integrator::getInternalType(){
325	return blsys;
326	}