[R-sig-dyn-mod] R-sig-dynamic-models Digest, Vol 141, Issue 4

[Granjon David]

Issue solved correcting the « obvious » typo: 
Time <- seq(0, 100, 0.001)

Instead of time <- c(0,100,0.001)

Best regards 

David

> Le 2 févr. 2018 à 12:00, r-sig-dynamic-models-request at r-project.org a écrit :
> 
> Send R-sig-dynamic-models mailing list submissions to
>    r-sig-dynamic-models at r-project.org
> 
> To subscribe or unsubscribe via the World Wide Web, visit
>    https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> or, via email, send a message with subject or body 'help' to
>    r-sig-dynamic-models-request at r-project.org
> 
> You can reach the person managing the list at
>    r-sig-dynamic-models-owner at r-project.org
> 
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of R-sig-dynamic-models digest..."
> 
> 
> Today's Topics:
> 
>   1. Compiled code in deSolve: fail with Lorenz model (Granjon David)
> 
> 
> ----------------------------------------------------------------------
> 
> Message: 1
> Date: Thu, 1 Feb 2018 14:37:07 +0100
> From: Granjon David <dgranjon at ymail.com>
> To: r-sig-dynamic-models at r-project.org
> Subject: [R-sig-dyn-mod] Compiled code in deSolve: fail with Lorenz
>    model
> Message-ID: <C460FC2F-6A38-4DD8-B8EA-0C6D5916632A at ymail.com>
> Content-Type: text/plain; charset="UTF-8"
> 
> Dear Members,
> 
> I was curious to see the speed difference between deSolve, odeintr, RxODE (and maybe others) by integrating 1000 times the very famous Lorenz model.
> 
> 1) Below is the implementation using odeintr:
> 
> library(odeintr)
> x <- 1:1000
> ## with cpp ##
> Lorenz.sys = '
> double test1, test2, test3;
> test1 = 10.0 * (x[1] - x[0]);
> test2 = 28.0 * x[0] - x[1] - x[0] * x[2];
> test3 = -8.0 / 3.0 * x[2] + x[0] * x[1];
> dxdt[0] = test1;
> dxdt[1] = test2;
> dxdt[2] = test3;
> ' # Lorenz.sys
> code <- compile_sys("lorenz", Lorenz.sys)
> time_odeintr <- unlist(
>  lapply(1:length(x), FUN = function(i) {
>    system.time({out <- lorenz(rep(1, 3), 100, 0.001)})[3]
>  })
> )
> 
> 2) Below is a similar code corresponding to RxODE:
> 
> # ## With RxODE ##
> library(RxODE)
> lorenz_RxODE <- "
>   d/dt(x1) = k1 * (x2 - x1);
>   d/dt(x2) = k2 * x1 - x2 - x1 * x3;
>   d/dt(x3) = k3 * x3 + x1 * x2;
> "
> mod1 <- RxODE(model = lorenz_RxODE, modName = "lorenz_RxODE", 
>              wd = getwd(), do.compile = TRUE)
> inits <- rep(1, 3)
> theta <- c(k1 = 10, k2 = 28, k3 = -8/3)
> ev1 <- eventTable()
> ev1$add.sampling(seq(0, 100, by = 0.001))
> 
> time_RxODE <- unlist(
>  lapply(1:length(x), FUN = function(i) {
>    system.time(out <- mod1$run(theta, events = ev1, inits = inits))[3]
>  })
> )
> 
> 3) Finally, below is the implementation using deSolve (and compiled code). The integration fails, and I have no idea how to fix it. On the other hand, using deSolve without compiled code work perfectly.
> 
> ## with Compiled deSolve in C ##
> library(deSolve)
> system("R CMD SHLIB lorenz_desolve.c")
> dyn.load(paste("lorenz_desolve", .Platform$dynlib.ext, sep = ""))
> parms <- c(k1 = 10, k2 = 28, k3 = -8/3) 
> state <- c(y1 = 1, y2 = 1, y3 = 1) 
> times <- c(0, 100, 0.001)
> 
> time_compiled_deSolve <- unlist(
>   lapply(1:length(x), FUN = function(i) {
>     system.time(out <- ode(y = state , times, func = "derivs", parms = parms,  
>                             dllname = "lorenz_desolve", initfunc = "initmod", 
>                            nout = 1, outnames = "Sum"))[3]
>   })
> )
> 
> The c code is here: Very important to call it ? lorenz_desolve.c ? 
> 
> /* file mymod.c */
> #include <R.h>
> static double parms[3];
> 
> #define k1 parms[0]
> #define k2 parms[1]
> #define k3 parms[2]
> 
> /* initializer  */
> void initmod(void (* odeparms)(int *, double *))
> {
>  int N=3;
>  odeparms(&N, parms);
> }
> /* Derivatives and 1 output variable */
> void derivs (int *neq, double *t, double *y, double *ydot,
>             double *yout, int *ip)
> {
>  if (ip[0] <1) error("nout should be at least 1");
>  ydot[0] = k1 * (y[1] - y[0]);
>  ydot[1] = k2 * y[0] - y[1] - y[0] * y[2];
>  ydot[2] = k3 * y[2] + y[0] * y[1];
> 
>  yout[0] = y[0] + y[1] + y[2];
> }
> 
> 
> Do you have any idea of why deSolve fails, while RxODE and odeintr manage to do the job without setting any particular numerical conditions?
> 
> Best regards
> 
> Granjon David
> PhD
> University of Z?rich
> Institute of Physiology
> Winterthurerstr. 190, Y23 J 78
> 8057 Z?rich
> +41 (0) 44 635 50 56
>    [[alternative HTML version deleted]]
> 
> 
> 
> ------------------------------
> 
> Subject: Digest Footer
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> 
> ------------------------------
> 
> End of R-sig-dynamic-models Digest, Vol 141, Issue 4
> ****************************************************