[R-sig-dyn-mod] R-sig-dynamic-models Digest, Vol 141, Issue 4
Granjon David
dgranjon at ymail.com
Fri Feb 2 14:14:13 CET 2018
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
> ****************************************************
More information about the R-sig-dynamic-models
mailing list