performance - Calculate more time points from ode() of the "deSolve" R package without increasing runtime -

performance - Calculate more time points from ode() of the "deSolve" R package without increasing runtime -

i wrote function: odesystem = function(t, states, parameters), contains ode-system , solved documented r packages "desolve" written karline soetaert, thomas petzoldt , r. woodrow setzer. documentation of bundle comprehensive , many examples. gives me confidence in programming , memory optimization skills.

however, when solving ode-system daily intervals instead of monthly intervals time takes calculate state-values specified moments increases tenfold. there might bit additional calculations reaching exact required moments in time, both cases same internal dynamic time steps should have been made. did not expect such big drop in runtime.

the phone call ode() in “desolve” looks this:

out <- as.data.frame(ode(states, t=times, func=odesystem), parms=parameters, method="ode45"))

i used 2 variants times times = seq(0, 100*365, by=365/12) # 100 years, 1 time point per month times = seq(0, 100*365, by=1) # 100 years, 1 time point per day

calling info points per month

user scheme elapsed 4.59 0.00 4.58

calling info points per month , cmpfun() on function containing odesystem

user scheme elapsed 4.39 0.00 4.38

calling info points per day

user scheme elapsed 44.41 0.00 44.46

calling info points per day , cmpfun() on function containing odesystem

user scheme elapsed 43.01 0.00 43.17

the runtime measured system.time() increases factor 10 when switching monthly intervals daily intervals. matters not improve much using cmpfun() on function containing ode-system.

(the output "out" assigned when function phone call ode() done. pre-assigning "out" yields no performance gain.)

question 1: looking reason why there decrease in runtime/performance.

(i expect in internals of desolve package.)

question 2: given reply question 1, how can improve runtime without resorting dynamic link libraries?

pre-assign memory become “out” might help (using knowledge on time steps in “times”), not known internal variable in ode() affect.

#### clear currrent lists memory rm(list=ls()) ### load libraries # library(rootsolve);library(ggplot2); library(base);library(desolve);library(stringr);library(compiler);library(data.table); #### constants dpy=365;durx1 = 40*dpy;rh = 1/durx1;durx4 = 365/12;rx4 = 1/durx4;durx6 = 365/12;rx6 = 1/durx6;durx2 = 80;rx2 = 1/durx2;durx3 = 31;rx3 = 1/durx3;durx7 = 20*365/12;rx7 = 1/durx7;durx5 = 29;rx5 = 1/durx5;durx8 = 200;rx8 = 1/durx8;fs = 0.013;fr = 8/100;fl = .03;fp = .03;ff = .05;x1zero = 1000;uddur = 365/12*5;rk = rx3*(1/uddur);fd1 = .05;fd2 = .05;durbt = 4;bt = 1/durbt;lx11 = 14;rf = 1/lx11;durx11 = 5;rx11 = 1/durx11;inix12 = 0;ph = 1;frac_im = 0;durx9 = dpy*5;ini_x2 = 1;sp = .90;fpx1 = 5;nf = fpx1*x1zero;rt1 = fd1*rx4;rt2 = fd2*rx6;px1 = 0*sp;px2 = 1/80*sp;px3 = .50*sp;px4 = .5*sp;px6 = .5*sp;px7 = 1/100*sp;px5 = px3;px9 = 0*sp;px8 = 1*sp;rx9 = 1/durx9; #### vector parameters parameters = c(rh, rx3, rx4, rx6, rx2, rx8, rx7, rx5, rk, rt1, rt2, bt, rf, nf, rx11, px1, px2, px3, px4, px6, px7, px5, px9, px8, rx9, x1zero) ### states contains initial conditions states = c( x1 =x1zero-1,x2=1,x3=0, x4=0, x5=0,x6=0, x7=0, x8=0, x9=0, x10=nf,x11=0,x12=0, x13 = 0) ### function ode scheme odesystem = function(t,states,parameters){ with(as.list(c(states,parameters)),{ ### functions x1part = (px2*x2 + px3*x3 + px4*x4 + px6*x6 + px7*x7 + px5*x5 + px9*x9 + px8*x8); prob1 = bt * x12 / x1zero; lf = bt * x1part/x1zero; advertisement = rk*(x3+x5+x4+x6)+rt1*x4+rt2*x6; ### fluxes j1 = prob1*x1; j2 = fs*rx2*x2; j3 = (1-fs)*rx2*x2; j4 = (1-fp)*rx3*x3 ; j5 = fp*rx3*x3; j6 = (1-ff)*rx4*x4; j7 = ff*rx4*x4; j8 = rx6*x6; j9 = fr*rx7*x7; j10 = rx5*x5; j11 = (1-fr)*(1-fl)*rx7*x7; j12 = (1-fr)*fl*rx7*x7; j13 = rx8*x8; j14 = rh*x3; j15 = rh*x1; j16 = rh*x2; j17 = rh*x4; j18 = rh*x6; j19 = rh*x5; j20 = rh*x8; j21 = rh*x7; j22 = rh*x9; j23 = rk*x3; j24 = rk*x4; j25 = rk*x6; j26 = rt2*x6; j27 = rh*x1zero; j28 = rt1*x4; j29 = ad; j30 = rk*x5; j31 = rf*x12; j32 = rf*x11; j33 = rf*x10; j34 = lf*x10; j35 = rx11*x11; j36 = rf*nf; j37 = rx9*x9; j38 = 0; j39 = 0; j40 = 0; j41 = 0; j42 = 0; j43 = 0; flux1=j4/x1zero*1e4*dpy; flux2=j12/x1zero*1e4*dpy; # rate of alter dx1 = - j1 - j15 + j27 + j29 + j37 dx2 = + j1 - j2 - j3 - j16 - j40 dx3 = + j2 - j4 - j5 - j14 - j23 - j41 dx4 = + j4 - j6 - j7 - j17 - j24 - j28 dx5 = + j9 - j10 - j19 - j30 - j43 dx6 = + j7 - j8 + j10 - j18 - j25 - j26 dx7 = + j5 + j6 + j8 - j9 - j11 - j12 - j21 dx8 = + j12 - j13 - j20 - j42 dx9 = + j3 + j11 + j13 - j22 - j37 + j40 + j41 + j42 + j43 dx10 = - j33 - j34 + j36 dx11 = - j32 + j34 - j35 dx12 = - j31 + j35 dx13 = + j38 - j39 # homecoming rate of alter list(c(dx1,dx2,dx3,dx4,dx5,dx6,dx7,dx8,dx9,dx10,dx11,dx12,dx13),flux1,flux2,prob1) }) } ## compiled version of ode scheme function cfodesystem=cmpfun(odesystem) #### time points calculated times = seq(0, 100*365,by=365/12) # 100 year, time points per month #times = seq(0, 100*365,by=1) # 100 year, time points per day ### calculations system.time(out <- as.data.frame(ode(states, t=times, func=odesystem, parms=parameters, method="ode45"))) #system.time(out <- as.data.frame(ode(states, t=times, func=cfodesystem, parms=parameters, method="ode45"))) ### longitudinal plots of each variable, flux1 , 2 , prob1 (i in seq(from=2, to=dim(out)[2], by=1) ) { tempdata <- out[c("time",names(out)[i])] tempdata$time= tempdata$time/365 templabel <-names(out)[i] plot(tempdata,col = "black","l",xlab="time (years)",ylab=templabel, xlim=c(0, max(tempdata$time)), ylim=c(0, signif(max(tempdata[2]),2))) }

so writing question, prompted me desolve internals , larn bit (and maybe speed own code up).

question 1

the ode function called number of times solve function (maybe less number of timepoints), called 1 time timepoint evaluate additional algebraic equations. if add together 30x timepoints, add together runtime less factor of 30, due things setup , teardown.

question 2

there few things can speed things without resorting c code (though first-class option)

use different solver (e.g. lsoda, on scheme around 5x faster ode45) when using lsoda, increment hmax allow adaptive timestepping more freedom integrate ahead (another speedup) rewrite code avoid utilize of with, instead using array accesses (possibly temporary named variables).

r performance ode