trunk/ascxx/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);

    /* we *don't* initialise a reporter here, cause we don't want to own it */
}

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){
    blsys->clientdata = reporter; /* this *is* necessary, as well as the following */
    integrator_set_reporter(blsys,reporter->getInternalType());
    //CONSOLE_DEBUG("REPORTER HAS BEEN SET");
    (*(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 return-code '%d' from integrator_analyse",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);

    if(blsys->reporter==NULL){
        throw runtime_error("No reporter has been assigned to the integrator");
    }

    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(const string &name){
    int res = integrator_set_engine(this->blsys, name.c_str());
    if(!res)return; // all OK
    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());
}

/**
    Ideally this list would be dynamically generated based on what solvers
    are available or are in memory.
*/
vector<string>
Integrator::getEngines(){
    vector<string> m;
    const gl_list_t *L = integrator_get_engines();
    for(unsigned long i=1; i<=gl_length(L); ++i){
        const IntegratorInternals *I = (const IntegratorInternals *)gl_fetch(L,i);
        m.push_back(I->name);
    }
    return m;
}

string
Integrator::getName() const{
    const IntegratorInternals *I = integrator_get_engine(blsys);
    if(I==NULL){
        throw runtime_error("No engine selected");
    }
    return I->name;
}

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