I am trying to build a generic stochastic differential equation solver coded by a class de_solver which takes some set of differential equations given by a model class. This model is fed to the solver through a class sde which works interface between model and solver, recasting equations of the original model in two subsets of deterministic + stochastic equations.
There are several problems with my code at the moment which I believe have to deal with inheritance and type conversion. In particular I cannot manage to pass the member functions where model equations are specified to the integrator. First say I have a model class of this type:
class model{
...
public:
size_t N_eq;
model(...); // Assign N_eq and parameters and other stuff
...
int equations_det(double t, const double y_[], double dy_[]); // Define deterministic RHS of equations
static int equations_gsl_wrapper(double t, const double y_[], double dy_[],void *params); // A wrapper that make equations_det suitable for a GSL solver
void equations_stoch(double t,double dt,const double y_[], double dy_[],); // The stochastic RHS of the model equations
static void ode_s_wrapper(double t,const double y_[], double dy_[], void *params); // A wrapper to pass equations_det to the stochastic integrator (in the class `de_solver`).
int simulate(t,tfin,tstep); // The actual simulator which will invoke the 'de_solver'
};
The static specification follows from the need to use GSL integrators for the deterministic part of the model, as pointed out in this post.
Then, the interfacing class sde is such that:
class sde{
...
public:
size_t N_eq;
sde(size_t N_,
int (*dtr)(double, const double*, double*, void*),
void (*stc)(double, double, const double*, double*, void*));
int (*deterministic)(double t, const double* y_[], double* dy_[], void * params);
void (*stochastic)(double t,const double y_[], double dy_[], void *params);
//Then again, akin to the model class, use some `wrappers`
static int deterministic_wrapper(double t, const double y_[], double dy_[], void * params);
static void stochastic_wrapper(double t, const double y_[], double *dy_, void *params);
};
The idea is to have inside the sde class, member functions that are inherited from the whatever model is given. As for the two 'wrappers,' I will clarify on the reason why I am introducing them shortly.
Finally the de_solver class is such that:
class de_solver{
sde *sys;
public:
de_solver(sde *system); // Will initialize solver with the system put in the `sde` form
...
void integrate(void *params, double *ts, double **sol);
};
model is declared/defined in separate files (model.H and model.CPP) with respect to sde and de_solver (solvers.h and solvers.cpp).
The idea is to have a simulate member function in the model class such that
int model::simulate(double t, double tfin, double dt){
// Prepare solver
// 1. Create the `sde` object from model `sys`
sde recast_sys(NEQ, model::deterministic_wrapper, model::stochastic_wrapper);
// 2. Instantiate solver with recast system
de_solver integrator(&recast_sys);
// Run simulation
double *ts = ... // Output time instants
double **sol = ... // Output solution
void *params_base = static_cast<void*>(std::addressof(this));
integrator.integrate(params_base);
return 1; // In practice there is some error check on this return condition (omitted here for brevity)
}
In summary sys invokes the integrator which works on deterministic and stochastic parts of the model equation as provided through recast_sys. Because the deterministic part of the integrator relies on GSL solvers, I am using the additional parameter argument to pass the pointer to the solver to the actually class member function. In this fashion, inside the integrator.integrate member function I have (see above mentioned post)
de_solver::integrate(void *params_base, ...){
...
// I allocate an array of two void pointers: the first to the `model` class (assumed to be passed by `params_base`), and the second to the `sde` class
void **params = (void**)calloc(2,sizeof(void*));
params[0] = params_base;
params[1] = reinterpret_cast<void *>(std::addressof(sys)); // the recast system as private member of the sde class
gsl_odeiv2_driver * d;
gsl_odeiv2_system system = {sys->deterministic_wrapper, nullptr, sys->NEQ, params};
d = gsl_odeiv2_driver_alloc_y_new (&system, gsl_odeiv2_step_bsimp, opts.dt, opts.atol, opts.rtol);
...
}
int sde::deterministic_wrapper(double t, const double y_[], double dy_[], void * params){
assert(params);
return(static_cast<sde*>(params[1])->deterministic(t,y_,dy_,params)); // This will issue an error: ‘void*’ is not a pointer-to-object type
}
int model::equations_gsl_wrapper(double t, const double y_[], double dy_[], void * params){
assert(params);
return(static_cast<model*>(params[0])->ode_gsl(t,y_,dy_)); // This will issue an error: ‘void*’ is not a pointer-to-object type
}
Allocating an array of two pointers to voids was taken from this post. However, it seems that the once used in the wrapper, it produces an error maybe because arithmetic on void arrays is not clear (as pointed out here)?
At the moment I cannot compile my code for the errors reported above. Also for some reason, the compiler tells me that the this pointer in the model.simulate member function generates
error: use of deleted function ‘const _Tp* std::addressof(const _Tp&&) [with _Tp = model*]’
I suspect that I am messing around with static and non-static member functions and I am not passing them correctly. Any input would be appreciated.
I have not fully grasped your entire problem, but the answers seem straightforward enough:
error: use of deleted function ‘const _Tp* std::addressof(const _Tp&&) [with _Tp = model*]’
The compiler tells you that you are attempting to take the address of the pointer to your instance, because you did std::addressof(this). this already is the pointer you are looking for (a model*), cast it to void* directly.
return(static_cast<model*>(params[0])->ode_gsl(t,y_,dy_)); // This will issue an error: ‘void*’ is not a pointer-to-object type
When you do params[0] you already dereferenced the pointer once (the inbuilt x[y] is exactly equivalent to *(x + y)) and would end up with a plain void, which the compiler does not like. Your array of void pointers has another level of indirection:
int deterministic_wrapper(double t, /* ... */ void * params) {
void** voidPtrs = static_cast<void**>(params);
return static_cast<sde*>(voidPtrs[1])->deterministic(t, /* ... */ params);
}
See here.
The better solution would be to just create a struct that holds pointers to your classes (your array of void pointers is basically that, just close enough to three star programming that it had you make the above mistake):
struct myParams
{
model* _model;
sde* _sde;
};
// ...
int deterministic_wrapper(double t, /* ... */ void * params) {
return static_cast<myParams*>(params)->_sde->deterministic(t, /* ... */ params);
}
Note that no matter what, you will have to clean up that memory. If you allocated that array you need to free it, if you new'd the struct you need to delete it, and if you create it on the stack you need to ensure it lives long enough for all these functions to finish.