Hi Daniel: I corrected the typo, however, the problem is still there. It is weird because the P_diss function did not cause any problem. Regards! Dabing rm(list=ls()) N= 50 V2<- c(300,rep(3,29),rep(0.5,(N-30))) Si = 0.07 # initial solubility in jejunum unit is mg/ml Sf = 0.007 # final solubility in jejunum is mg/ml t0 = 20 # time for half phase transformation in jejunum in min p0 = 0.5 S <- rep(0.1,N) Sp <- rep(0.1,N) GImobility = 0.25 for (i in 2:N){ S[i] <- (Si-Sf)/(1+exp(p0*(((i-1)/GImobility)-t0)))+Sf Sp[i] <- (Si-Sf)/(1+exp(p0*(((i-1)/GImobility)-t0)))+Sf } z <- 50 # define a function for dissolution kinetics P_diss <- function(X,Y,V,S){ result<-pmin(z*(X^(1/3))*(S-Y/V),X) result[(X<=0)|(S Sp)] <- Y - Sp*V return(result) } test1 <- c(0,1,2,rep(0,(N-3))) test2 <- c(1,2,0,rep(0,(N-3))) P_diss(test1,test2,V2,S) P_precip(test1,test2,V2,Sp) On Tue, Jul 30, 2013 at 5:32 PM, Daniel Reed wrote: > 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@umich.edu > web: www.danielreed.org > > On Jul 30, 2013, at 6:11 PM, Dabing Chen 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 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 >> > >> return(result) > >> } > >> > >> Cheers, > >> Daniel > >> > >> _______________________________ > >> Daniel C. Reed, PhD. > >> Postdoctoral Fellow, > >> Dept. of Earth & Environmental Sciences, > >> University of Michigan, > >> Ann Arbor, MI, USA. > >> email: reeddc@umich.edu > >> web: www.danielreed.org > >> > >> On Jul 30, 2013, at 5:21 PM, Dabing Chen 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 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@umich.edu > >>>> web: www.danielreed.org > >>>> > >>>> On Jul 30, 2013, at 3:27 PM, Dabing Chen 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 > >>>>> 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@r-project.org [mailto: > >>>>>> r-sig-dynamic-models-bounces@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@r-project.org > >>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models > >>>>>> > >>>>>> _______________________________________________ > >>>>>> R-sig-dynamic-models mailing list > >>>>>> R-sig-dynamic-models@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@r-project.org > >>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models > >>>> > >>>> _______________________________________________ > >>>> R-sig-dynamic-models mailing list > >>>> R-sig-dynamic-models@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@r-project.org > >>> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models > >> > >> _______________________________________________ > >> R-sig-dynamic-models mailing list > >> R-sig-dynamic-models@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@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models > > _______________________________________________ > R-sig-dynamic-models mailing list > R-sig-dynamic-models@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models > [[alternative HTML version deleted]]