[R-sig-dyn-mod] question on ode.1D
Daniel Reed
reeddc at umich.edu
Tue Jul 30 22:39:55 CEST 2013
Hi Dabing:
I'm not sure you'd need sapply or lappy. Perhaps something like this would work… (NB: untested code)
P_diss <- function(X,Y,V,S){
result<-pmin(z*(X^(1/3))*(S-Y/V),X)
result[(X<=0)|(S<Y/V)]<-0
return(result)
}
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 5:21 PM, Dabing Chen <dabing.c at gmail.com> wrote:
> Hi Daniel:
>
> Thanks a lot for the tip.
> With vode, the program can finish in 40 sec, instead of 4 min
> sometimes.
> I am having difficulties in apply the P_precip and P_diss on the
> matrix. As you can see, in the function, there is if (Y/V > S) criteria,
> they always give me error message.
> I was trying to use Sapply or Lapply. Which one should do the work?
> Regards!
>
> Dabing
>
>
>
>
>
>
> On Tue, Jul 30, 2013 at 3:32 PM, Daniel Reed <reeddc at umich.edu> wrote:
>
>> 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
>>
>> _______________________________________________
>> 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