[R-sig-dyn-mod] deSolve with compiled code and varying size of parms

[Thomas Petzoldt]

Hi,

Am 30.03.2017 um 00:41 schrieb Tim Keitt:
> Unlike more
> integrated solvers, that framework puts more onus on the user to not supply
> invalid parameters, or perhaps it would be better to say that the user has
> to avoid certain corner cases.

Tim made a good point here. It is not easy to prevent all corner cases 
that potentially break R, if compiled code is allowed, both in odeintr 
(based on ODEINT), deSolve (based partly on ODEPACK) and other packages 
that were built on top of them.

In the meantime, I had a look again in the deSolve sources. The 
initialisation routine is a little bit tricky to enable communication 
between the user-dll and deSolve.dll. It would be possible, in 
principle, to relax the "confusion over the length ..." error message or 
to enable passing more complex parameter constructions other than double 
precision vectors, i.e. a complete SEXP that could then be 
de-constructed by the user.

I assume that the decision about the length check was made with the 
intent to make the interface easy to use, and to request a minimum of 
self-check from the user (here, thinking about the parameter length 
him/herself) -- so I would leave it as is.

Besides this intrinsically fixed match between the C and R vector, there 
are several methods that can make parameter passing more flexible, i.e. 
using C vectors directly instead of #define macros, or employ UNIONs 
between numbered vector arguments and named arguments.

One can then pass a bigger vector with some reserve (e.g. Nmax 100), and 
then to fill only the parameters that are needed.

A small example is shown below, another example with UNIONs can be 
posted on request.

Thomas


-------------- next part --------------
/* file lorenz.c */
#include <R.h>

#define Nmax 3

double p[Nmax];


/* initializer  */
void initmod(void (* odeparms)(int *, double *))
{
    int N = Nmax;
    odeparms(&N, p);
}

/* Derivatives */
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip)
{
    ydot[0] = p[0] * y[0] + y[1] * y[2];
    ydot[1] = p[1] * (y[1] - y[2]);
    ydot[2] = - y[0] * y[1] + p[2] * y[1] - y[2];
}
-------------- next part --------------
library(deSolve)
library(scatterplot3d)

system("R CMD SHLIB lorenzc_pvar.c")
dyn.load("lorenzc_pvar.dll")

parameters <- c(a = -8/3, b = -10, c =  28)
state <- c(X = 1, Y = 1, Z = 1)
times <- seq(0, 100, by = 0.01)

out <- ode(state, times, func = "derivs", parms = parameters,
           dllname = "lorenzc_pvar", initfunc = "initmod")
plot(out)
scatterplot3d(out[,-1], type="l")

dyn.unload("lorenzc_pvar.dll")