model <- function(t, state, parms) {
  state <- ifelse(state<0, 0, state)
  with(as.list(c(state,parms)), {
    f <- v*R/(R+k)
    tlag <- t - lambda
    if (tlag < 0) lags <- rep(0,5) # no initial infection
    else lags <- lagvalue(tlag)    # returns 5 lags 
    dtR <- - f*e*(B0+(1-s)*B1)
    dtB0 <- f*B0 - delta*B0*P
    dtB1 <- f*(1-s)*B1
    dtM <- delta*B0*P - delta*lags[2]*lags[5] 
    dtP <- b*delta*lags[2]*lags[5] - delta*P*(B0+B1)
    return(list(c(dtR, dtB0, dtB1, dtM, dtP)))  
  }) 
}  

odeModel <- function(t, state, parms) {
  state <- ifelse(state<0, 0, state)
  with(as.list(c(state,parms)), {
    f <- v*R/(R+k)
    dtR <- - f*e*(B0+(1-s)*B1)
    dtB0 <- f*B0 - delta*B0*P
    dtB1 <- (1-s)*f*B1
    dtM <- delta*B0*P - M/lambda
    dtP <- b*M/lambda - delta*P*(B0+B1)
    return(list(c(dtR, dtB0, dtB1, dtM, dtP)))  
  }) 
} 

p <- c(b=80,delta=5e-8,e=5e-7,lambda=0.4,v=1.4,k=1,s=0.001)
s <- c(R=350,B0=5e6,B1=500,M=0,P=5e6)
tweak <- "nsol$B<-nsol[,\"B0\"]+nsol[,\"B1\"]+nsol[,\"M\"]"
run(7,0.1,ymin=1e3,ymax=1e10,log="y",delay=TRUE,tweak=tweak)

fig2B0 <- read.csv("data/LevinPG13Fig2B0.csv",header=TRUE)
fig2B  <- read.csv("data/LevinPG13Fig2.csv",header=TRUE)
timePlot(fig2B0,ymin=1e3,ymax=1e10,log="y")
timePlot(fig2B,ymin=1e3,ymax=1e10,log="y")

free <- c("v")
s <- c(R=350,B0=5189702,B1=0,M=0,P=0)
f0 <- fit(fig2B0,tstep=0.1,ymin=1e3,ymax=1e10,log="y",fun=log1p,lower=0,free=free,tweak=tweak,delay=TRUE)

p[free] <- f0$par
#p["s"] <- 0.25
s <- c(R=350,B0=4382322,B1=1000,M=0,P=5.58e6)
free <- c("B1","b","delta","lambda","s")
f1 <- fit(fig2B,tstep=0.1,ymin=1e3,ymax=1e10,log="y",fun=log1p,free=free,tweak=tweak,lower=0,delay=TRUE)

# Original parameters now for the ODE model:

p <- c(b=80,delta=5e-8,e=5e-7,lambda=0.4,v=1.4,k=1,s=0.001)
s <- c(R=350,B0=5189702,B1=0,M=0,P=0)
run(7,0.1,ymin=1e3,ymax=1e10,log="y",tweak=tweak)

free <- c("v")
f0 <- fit(fig2B0,odes=odeModel,tstep=0.1,ymin=1e3,ymax=1e10,log="y",fun=log1p,lower=0,free=free,tweak=tweak)

p[free] <- f0$par
s <- c(R=350,B0=4382322,B1=1000,M=0,P=5.58e6)
free <- c("B1","b","delta","lambda","s")
f1 <- fit(fig2B,odes=odeModel,tstep=0.1,ymin=1e3,ymax=1e10,log="y",fun=log1p,free=free,tweak=tweak,lower=0)