trunk/pygtk/integrator.cpp

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

/**
    '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.
*/
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(){
    CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
    CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
    integrator_free(blsys);
    CONSOLE_DEBUG("DESTROYING samplelist at %p",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);
}

int 
Integrator::findIndependentVar(){
    return integrator_find_indep_var(blsys);
}

int
Integrator::analyse(){

    int res;
    res = integrator_analyse(blsys);

    if(!res){
        ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Failed system analysis");
        return 0;
    }

    return 1;
}

/**
    @TODO what about root detection?

    Integrate the function for the timesteps specified.

    This method will throw a runtime_error if integration fails.
*/
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);

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

    if(!res){
        throw runtime_error("Failed integration");
    }
}

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

    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;
#ifdef ASC_WITH_LSODE
    m.insert(pair<int,string>(INTEG_LSODE,"LSODE"));
#endif
#ifdef ASC_WITH_IDA
    m.insert(pair<int,string>(INTEG_IDA,"IDA"));
#endif
    return m;
}

string
Integrator::getEngineName() 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	/**
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	CONSOLE_DEBUG("DESTROYING Integrator (C++) at %p",this);
31	CONSOLE_DEBUG("DESTROYING IntegratorSystem at %p",blsys);
32	integrator_free(blsys);
33	CONSOLE_DEBUG("DESTROYING samplelist at %p",samplelist);
34	samplelist_free(samplelist);
35	}
36
37	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
45	void
46	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	Integrator::setReporter(IntegratorReporterCxx *reporter){
53	this->blsys->clientdata = reporter;
54	integrator_set_reporter(blsys,reporter->getInternalType());
55	CONSOLE_DEBUG("REPORTER HAS BEEN SET");
56	(*(this->blsys->reporter->init))(blsys);
57	CONSOLE_DEBUG("DONE TESTING OUTPUT_INIT");
58	}
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
97	Integrate the function for the timesteps specified.
98
99	This method will throw a runtime_error if integration fails.
100	*/
101	void
102	Integrator::solve(){
103
104	// check the integration limits
105	// trigger of the solution process
106	// report errors?
107
108	assert(samplelist!=NULL);
109	assert(samplelist->ns>0);
110	assert(blsys->reporter!=NULL);
111	assert(blsys->clientdata!=NULL);
112
113	int res;
114	res = integrator_solve(blsys, 0, samplelist_length(samplelist)-1);
115
116	// communicate solver variable status back to the instance tree via 'interface_ptr'
117	simulation.processVarStatus();
118
119	if(!res){
120	throw runtime_error("Failed integration");
121	}
122	}
123
124	void
125	Integrator::setEngine(IntegratorEngine engine){
126	int res = integrator_set_engine(this->blsys, engine);
127	if(!res)return;
128	if(res==1)throw range_error("Unknown integrator");
129	if(res==2)throw range_error("Invalid integrator");
130	stringstream ss;
131	ss << "Unknown error in setEngine (res = " << res << ")";
132	throw runtime_error(ss.str());
133	}
134
135	void
136	Integrator::setEngine(int engine){
137	setEngine((IntegratorEngine)engine);
138	}
139
140	void
141	Integrator::setEngine(const string &name){
142	CONSOLE_DEBUG("Setting integration engine to '%s'",name.c_str());
143	IntegratorEngine engine = INTEG_UNKNOWN;
144	#ifdef ASC_WITH_LSODE
145	if(name=="LSODE")engine = INTEG_LSODE;
146	#endif
147	#ifdef ASC_WITH_IDA
148	if(name=="IDA")engine = INTEG_IDA;
149	#endif
150
151	setEngine(engine);
152	}
153
154	/**
155	Ideally this list would be dynamically generated based on what solvers
156	are available or are in memory.
157	*/
158	map<int,string>
159	Integrator::getEngines(){
160	map<int,string> m;
161	#ifdef ASC_WITH_LSODE
162	m.insert(pair<int,string>(INTEG_LSODE,"LSODE"));
163	#endif
164	#ifdef ASC_WITH_IDA
165	m.insert(pair<int,string>(INTEG_IDA,"IDA"));
166	#endif
167	return m;
168	}
169
170	string
171	Integrator::getEngineName() const{
172	map<int,string> m=getEngines();
173	map<int,string>::iterator f = m.find(integrator_get_engine(blsys));
174	if(f==m.end()){
175	throw runtime_error("No engine selected");
176	}
177	return f->second;
178	}
179
180	/**
181	@TODO what about conversion factors? Is an allowance being made?
182	*/
183	void
184	Integrator::setLinearTimesteps(UnitsM units, double start, double end, unsigned long num){
185	if(samplelist!=NULL){
186	ASC_FREE(samplelist);
187	}
188	const dim_type *d = units.getDimensions().getInternalType();
189	samplelist = samplelist_new(num+1, d);
190	double val = start;
191	double inc = (end-start)/(num);
192	for(unsigned long i=0;i<=num;++i){
193	samplelist_set(samplelist,i,val);
194	val += inc;
195	}
196	integrator_set_samples(blsys,samplelist);
197	}
198
199	/**
200	@TODO what about conversion factors? Is an allowance being made?
201	*/
202	void
203	Integrator::setLogTimesteps(UnitsM units, double start, double end, unsigned long num){
204	if(samplelist!=NULL){
205	ASC_FREE(samplelist);
206	}
207	const dim_type *d = units.getDimensions().getInternalType();
208
209	if(start<=0)throw runtime_error("starting timestep needs to be > 0");
210	if(end<=0)throw runtime_error("end timestep needs to be > 0");
211	if(end <= start)throw runtime_error("end timestep needs to be > starting timestep");
212
213	samplelist = samplelist_new(num+1, d);
214	double val = start;
215	double inc = exp((log(end)-log(start))/num);
216	for(unsigned long i=0;i<=num;++i){
217	samplelist_set(samplelist,i,val);
218	CONSOLE_DEBUG("samplelist[%lu] = %f",i,val);
219	val *= inc;
220	}
221	integrator_set_samples(blsys,samplelist);
222	}
223
224	vector<double>
225	Integrator::getCurrentObservations(){
226	double *d = ASC_NEW_ARRAY(double,getNumObservedVars());
227	integrator_get_observations(blsys,d);
228	vector<double> v=vector<double>(d,d+getNumObservedVars());
229	// do I need to free d?
230	// can I do this in such a way as I avoid all this memory-copying?
231	return v;
232	}
233
234	Variable
235	Integrator::getObservedVariable(const long &i){
236	var_variable *v = integrator_get_observed_var(blsys,i);
237	return Variable(&simulation,v);
238	}
239
240	Variable
241	Integrator::getIndependentVariable(){
242	var_variable *v = integrator_get_independent_var(blsys);
243	if(v==NULL){
244	throw runtime_error("independent variable is null");
245	}
246	return Variable(&simulation,v);
247	}
248
249	int
250	Integrator::getNumVars(){
251	return blsys->n_y;
252	}
253
254	int
255	Integrator::getNumObservedVars(){
256	return blsys->n_obs;
257	}
258
259	void
260	Integrator::setMinSubStep(double n){
261	integrator_set_minstep(blsys,n);
262	}
263
264	void
265	Integrator::setMaxSubStep(double n){
266	integrator_set_maxstep(blsys,n);
267	}
268
269	void
270	Integrator::setInitialSubStep(double n){
271	integrator_set_stepzero(blsys,n);
272	}
273
274	void
275	Integrator::setMaxSubSteps(int n){
276	integrator_set_maxsubsteps(blsys,n);
277	}
278
279	IntegratorSystem *
280	Integrator::getInternalType(){
281	return blsys;
282	}