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.

    @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){
        throw runtime_error("Failed integration");
    }

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

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