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<Y/V)]<-0

  return(result)
}

P_precip <- function(X,Y,V,Sp){

  result <- pmin(z*(X^(1/3))*(Sp-Y/V),X)
  result[(X<0)] <- -X
  result[(Y/V > 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 <reeddc@umich.edu> 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 <dabing.c@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@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@umich.edu
> >> web: www.danielreed.org
> >>
> >> On Jul 30, 2013, at 5:21 PM, Dabing Chen <dabing.c@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@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@umich.edu
> >>>> web: www.danielreed.org
> >>>>
> >>>> On Jul 30, 2013, at 3:27 PM, Dabing Chen <dabing.c@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@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@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]]

