[R-sig-dyn-mod] question on ode.1D
Daniel Reed
reeddc at umich.edu
Tue Jul 30 21:32:56 CEST 2013
Hi Dabing:
I find that the "vode" method is often quite fast when using ode.1D. Another option is rewriting advModel in C/C++/FORTRAN. Karline et al. have a vignette explaining that option here: http://cran.r-project.org/web/packages/deSolve/vignettes/compiledCode.pdf
Also, each time ode.1D calls advModel it (re)defines the functions GI, P_precip, and P_diss. I wonder if it would speed up if you defined those functions prior to advModel, so that they were only defined once per model run. Finally, instead of using a for loop, it would be much quicker if P_precip and P_diss operated on whole arrays at a time, if possible.
Cheers,
Daniel
_______________________________
Daniel C. Reed, PhD.
Postdoctoral Fellow,
Dept. of Earth & Environmental Sciences,
University of Michigan,
Ann Arbor, MI, USA.
email: reeddc at umich.edu
web: www.danielreed.org
On Jul 30, 2013, at 3:27 PM, Dabing Chen <dabing.c at gmail.com> wrote:
> Hi Karline:
> Thanks a lot for the reply. I rewrote the code, removed the the two
> variables that do not belong to any boxes.
> The program is much more robust now. However, it took 3 min to finish.
> Is there a way to have it run faster?
>
> Regards!
> Dabing
>
>
>
>
>
> rm(list=ls())
> library (deSolve)
> library(ReacTran)
>
>
> N= 50
>
> ka <- c(0,rep(0.03,29),rep(0.004,(N-30)))
>
>
> Dose <- 100 # initial dose in mg
>
> Ke <- 0.074 # rate of elimination
> Vs <- 6.6E3 # volume of distribution
>
> V2<- c(300,rep(3,29),rep(0.5,(N-30)))
>
>
>
> S <- rep(0.1,N)
> Sp <- rep(0.1,N)
>
> GImobility = 0.25
>
> z <- 50
>
> advModel <- function(t, state, parms) {
>
> with (as.list(parms), {
>
> Intestine1 <- state[1:N] # Tablet
> Intestine2 <- state[(N+1):(2*N)] # granule
> Intestine3 <- state[(2*N+1):(3*N)] # precipitate
> Intestine4 <- state[(3*N+1):(4*N)] # solution
> Intestine5 <- state[(4*N+1):(5*N)] # absorbed
>
>
> GI <- function(X,emptying) {
> dX <- rep(0,N)
> dX[1] <- -X[1]*emptying
> advection <- advection.1D(C = X[2:N], flux.up = emptying*X[1], v
> = GImobility, dx=1,adv.method="super")
> dX[2:N] <- advection$dC
> return(dX)
> }
>
>
> # define a function for dissolution kinetics
> P_diss <- function(X,Y,V,S){
> # X the solid remaining at time t, X0, the inital dose, r, the
> particle size
> if (X <= 0) {
> result <- 0 # avoid numerical error
> } else if (S < Y/V){
> result <- 0
> } else {
> result <- min(z*(X^(1/3))*(S-Y/V),X)
> }
> return(result)
> }
>
> # define the precipitation rate function
>
> P_precip <- function(X,Y,V,Sp){
> # X the solid remaining at time t, assuming the precipitate has
> particle size of 1 um
> if (X < 0) {
> result <- -X # avoid numerical error
>
> } else if(Y/V > Sp) {
> result <- V*( Y/V-Sp)
> } else {
> result <- -min(z*(X^(1/3))*(Sp-Y/V),X)
> }
> return(result)
> }
>
>
> ddisinte <- 0.03 # rate of tablet disintegration
>
> disso <- rep(0,N)
> precip <- rep(0,N)
>
> for (i in 1:N){
> disso[i] <- P_diss(Intestine2[i],Intestine4[i],V2[i],S[i])
> precip[i] <- P_precip(Intestine3[i],Intestine4[i],V2[i],Sp[i])
> }
>
>
> emptying <-
> 0.01+0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))
> emptying1 <- 0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))
> # dissolution in Stomach
> dIntestine1 <- GI(Intestine1,emptying1) - ddisinte*Intestine1
> dIntestine2 <- GI(Intestine2,emptying) - disso +
> ddisinte*Intestine1
> dIntestine3 <- GI(Intestine3,emptying) + precip
> dIntestine4 <- GI(Intestine4,emptying) + disso - precip -
> ka*Intestine4
> dIntestine5 <- ka*Intestine4
>
>
> return (list(c(dIntestine1,dIntestine2,dIntestine3,
> dIntestine4,dIntestine5)))
> })
> }
>
> times <- seq(1,1000,by =1)
>
> state <- c(Dose, rep(0,(5*N-1)))
>
> systime <- system.time(out <- ode.1D(func = advModel, y=state,times =times,
> parms = 0,nspec=5))
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> On Tue, Jul 30, 2013 at 7:04 AM, Karline Soetaert
> <Karline.Soetaert at nioz.nl>wrote:
>
>> Dabeng,
>>
>> Your complex example does not work, so I cannot really test it.
>>
>> However, ode.1D only works efficiently if (1) all components are defined
>> in an equal number of compartments (boxes), and (2) if for one component
>> there is transition only between adjacent compartments, e.g. between boxes
>> i-1, and i. (3) Also, interactions between different components can only
>> occur if in the same box.
>> This is so because ode.1D was meant to efficiently solve 1-D
>> reactive-transport equations.
>>
>> If these conditions are not valid for your model, then ode.1D will be very
>> inefficient.
>>
>> If your model is really sparse, you could try to use lsodes.
>>
>> Hope it helps
>>
>> Karline
>>
>>
>> -----Original Message-----
>> From: r-sig-dynamic-models-bounces at r-project.org [mailto:
>> r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Dabing Chen
>> Sent: donderdag 25 juli 2013 21:36
>> To: R simulation
>> Subject: [R-sig-dyn-mod] question on ode.1D
>>
>> Hi All:
>>
>> I am having some problems with the program below. I was able to
>> run when ka = 0.04. However, when I change it to 0.01, the program cannot
>> finish.
>>
>> I also tried using different integration methods in ode.1D, such
>> as "vode, daspk", which did not work. radau runs almost 3 min, the program
>> runs better when I left it blank, which took about 1 min to run.
>>
>> How can I make it faster and more robust?
>> Thanks a lot in advance.
>>
>> Dabing
>>
>>
>> rm(list=ls())
>> library (deSolve)
>> library(ReacTran)
>>
>> ka <- 0.04
>>
>> N= 100
>>
>> Dose <- 100 # initial dose in mg
>>
>> V1 <- 300
>> V2<- c(rep(3,30),rep(0.5,(N-30)))
>>
>>
>> Si = 0.02 # initial solubility in jejunum unit is mg/ml Sf = 0.02 # final
>> solubility in jejunum is mg/ml
>> t0 = 20 # time for half phase transformation in jejunum in min
>> p0 = 0.5
>>
>> SI <- rep(0.02,N)
>> GImobility = N/60/8
>>
>> for (i in 1:N){
>> SI[i] <- (Si-Sf)/(1+exp(p0*((i/GImobility)-t0)))+Sf
>> }
>>
>> Weibulla <- 2 # weibull function for time factor Weibullb <-0.9 # weibull
>> fuction for shape factor Tlag <- 0 # lag time for disintegration criticT <-
>> (15)^(1/Weibullb)*Weibulla # critical T for disintegration
>>
>> z <- 50
>>
>> ka <- 0.04
>>
>> para <- c(SGi <- 0.8, # inital solubility in Stomach unit is mg/ml
>> SGf <- 0.8,# final solubility in Stomach unit is mg/ml
>> TG <- 20, #time for half transformation in Stomach
>> PG <- 0.2) # rate of half transformation in Stomach
>>
>> advModel <- function(t, state, parms) {
>>
>> with (as.list(parms), {
>> Stomach1 <- state[1] # undisintegrated
>> Stomach2 <- state[2] # granule
>> Stomach3 <- state[3] # Precipitate
>> Stomach4 <- state[4] # dissolved
>> Intestine1 <- state[5:(N+4)] # undisintegated
>> Intestine2 <- state[(N+5):(2*N+4)] # granule
>> Intestine3 <- state[(2*N+5):(3*N+4)] # precipitate
>> Intestine4 <- state[(3*N+5):(4*N+4)] # dissolved
>> T1 <- state[4*N+5] # total undisintegrated mass
>> T2 <- state[4*N+6] # total granule
>> T3 <- state[4*N+7] # total precipitate
>> T4 <- state[4*N+8] # total in solution
>>
>> absorption <- state[4*N+9]
>> plasma <- state[4*N+10]
>>
>> SG <- SGi
>>
>> ddisinte <- 0.03
>> # wettability factor
>>
>> # define a function for dissolution kinetics
>> P_diss <- function(X,Y,V,S){
>> # X the solid remaining at time t, X0, the inital dose, r, the
>> particle size
>> if (X <= 0) {
>> result <- 0 # avoid numerical error
>> } else if (S < Y/V){
>> result <- 0
>> } else {
>> result <- min(z*(X^(1/3))*(S-Y/V),X)
>> }
>> return(result)
>> }
>>
>> # define the precipitation rate function
>>
>> P_precip <- function(X,Y,V,S){
>> # X the solid remaining at time t, assuming the precipitate has
>> particle size of 1 um
>> if (X <= 0) {
>> result <- 0 # avoid numerical error
>> } else if (Y/V > S){
>>
>> result <- V * ( Y/V-S)
>> } else {
>> result <- -min(z*(X^(1/3))*(S-Y/V),X) # + (V*( Y/V-S) +
>> min(z*(X^(1/3))*(S-Y/V),X))/(1+exp(1000*(S-Y/V)))
>> }
>> return(result)
>> }
>>
>>
>> emptying <-
>> 0.005+0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))
>> emptying1 <- 0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))
>>
>>
>> # dissolution in Stomach
>> disso1 <- P_diss(Stomach2,Stomach4,V1,SG)
>> precip1 <- P_precip(Stomach3,Stomach4,V1,SG)
>>
>> dStomach1 <- -ddisinte*Stomach1 - emptying1*Stomach1
>> dStomach2 <- ddisinte*Stomach1 - disso1 - emptying*Stomach2
>> dStomach3 <- precip1 -emptying*Stomach3
>> dStomach4 <- disso1 - precip1 - emptying*Stomach4
>>
>> disso2 <- rep(0,N)
>> precip2 <- rep(0,N)
>>
>> for (i in 1:N){
>> disso2[i] <- P_diss(Intestine2[i],Intestine4[i],V2[i],SI[i])
>> precip2[i] <- P_precip(Intestine3[i],Intestine4[i],V2[i],SI[i])
>> }
>>
>> Tran1 <- advection.1D(C = Intestine1, flux.up =
>> emptying1*Stomach1, v = GImobility, dx=1,adv.method="super")
>> Tran2 <- advection.1D(C = Intestine2, flux.up =
>> emptying*Stomach2, v = GImobility, dx=1,adv.method="super")
>> Tran3 <- advection.1D(C = Intestine3, flux.up =
>> emptying*Stomach3, v = GImobility, dx=1,adv.method="super")
>> Tran4 <- advection.1D(C = Intestine4, flux.up =
>> emptying*Stomach4, v = GImobility, dx=1,adv.method="super")
>>
>>
>> dIntestine1 <- Tran1$dC - ddisinte*Intestine1
>> dIntestine2 <- Tran2$dC - disso2 + ddisinte*Intestine1
>> dIntestine3 <- Tran3$dC + precip2
>> dIntestine4 <- Tran4$dC + disso2 - precip2 - ka*Intestine4
>>
>> dT1 <- sum(dIntestine1)
>> dT2 <- sum(dIntestine2)
>> dT3 <- sum(dIntestine3)
>> dT4 <- sum(dIntestine4)
>>
>>
>> dabsorption <- sum(ka*Intestine4)
>> dplasma <- dabsorption -0.002*plasma
>>
>>
>> return (list(c(dStomach1,dStomach2,dStomach3,dStomach4,
>> dIntestine1,dIntestine2,dIntestine3,dIntestine4,
>> dT1,dT2,dT3,dT4,
>> dabsorption,dplasma)))
>> })
>> }
>>
>> times <- seq(1,1000,by =1)
>>
>> state <- c(Dose,rep(0,3),rep(0,4*N),rep(0,4),0,0)
>>
>> systime <- system.time(out <- ode.1D(func = advModel, y=state,times =times,
>> parms = para,nspec=1))
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-dynamic-models mailing list
>> R-sig-dynamic-models at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>>
>> _______________________________________________
>> R-sig-dynamic-models mailing list
>> R-sig-dynamic-models at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
More information about the R-sig-dynamic-models
mailing list