Hi Dabing, the problem is that in P_precip "result[(Y/V > Sp)]" and "Y - Sp*V" have different lengths - something like that should help: P_precip <- function(X,Y,V,Sp){ result <- pmin(z*(X^(1/3))*(Sp-Y/V),X) # the two conditions cond1 <- X<0 cond2 <- (Y/V > Sp) # indexing with both conditions result[cond1] <- -X[cond1] result[cond2] <- (Y - Sp*V)[cond2] return(result) } Cheers, Bernhard On 31.07.2013 15:52, Dabing Chen wrote: > 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 > 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 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]] > > > > _______________________________________________ > 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]]