[R-sig-dyn-mod] question on ode.1D
Daniel Reed
reeddc at umich.edu
Tue Jul 30 23:32:30 CEST 2013
Hi Dabing:
I think there's a typo: should 'S' in the second line of the function be 'Sp'? Does that fix your problem?
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 6:11 PM, Dabing Chen <dabing.c at gmail.com> wrote:
> Hi Daniel:
>
> I was doing the same thing with P-precip:
>
> P_precip <- function(X,Y,V,Sp){
> result <- pmin(z*(X^(1/3))*(S-Y/V),X)
> result[(X<0)]<- -X
> result[(Y/V>Sp)] <- V*(Y/V-Sp)
>
> return(result)
> }
>
> an error message says:
> Warning message:
> In result[(Y/V > Sp)] <- V * (Y/V - Sp) :
> number of items to replace is not a multiple of replacement length
>
> How should I rewrite the code?
>
> Thanks a lot.
>
> Dabing
>
>
> On Tue, Jul 30, 2013 at 4:39 PM, Daniel Reed <reeddc at umich.edu> wrote:
>
>> 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
>>
>> _______________________________________________
>> 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