################################################################################
# Course SMACH - January 2014
# a classical SEIR model. From original publication of Zadoks (1971)
# not all features are specifyied in this version
# version $Date:: 2014-04-30 1#$ by $Author: brun $(v0: 2013-10-9)
################################################################################
zakoks.original.model = function (nlpd=4*10,nipd=1*10,dmfr=16,SITE0 = 5*10^9,weather, sdate = 145, ldate = 145+50, XLAT0=1){
    weather$TEEQ_L = 1/(nlpd)*pmax(0,(weather$TMIN+weather$TMAX)/2)
    weather$TEEQ_I = 1/(nipd)*pmax(0,(weather$TMIN+weather$TMAX)/2)
    
    XVAC <- rep(NA, ldate)
    XLAT <- rep(NA, ldate)
    XINF <- rep(NA, ldate)
    XCTR <- rep(NA, ldate)
    XTO1 <- rep(NA, ldate)
    # Latent : Boxcar train
    BOXL=matrix(0,nrow=ldate, ncol=2, dimnames =list(NULL, c("XLAT","sumTEEQ")))
    # Infectious : Boxcar train
    BOXI=matrix(0,nrow=ldate, ncol=2, dimnames =list(NULL, c("XINF","sumTEEQ")))

    XVAC[sdate] <- SITE0-XLAT0
    XLAT[sdate] <-  XLAT0
    XINF[sdate] <- 0
    XCTR[sdate] <- 0

    BOXL[sdate,"sumTEEQ"] = weather$TEEQ_L[sdate]
    BOXL[sdate,"XLAT"] = XLAT[sdate]
    #BOXI[sdate,"sumTEEQ"] = 0 #weather$TEEQ_I[sdate]
    #BOXI[sdate,"XINF"] = XINF[sdate]

    for (day in sdate:(ldate - 1)) {
        #if (day==inf_start){BOXI[day,"XINF"]=1, BOXI[sdate,"sumTEEQ"] = 1}
        # correction factor : feed back from total occupied sites
        cofr<-max(XVAC[day]/SITE0, 0)
        # rocc: rate of occupation : nb of sites Vacant=>Latent
        rocc = cofr * dmfr * sum(BOXI[,"XINF"],na.rm=TRUE)

        # rapp: rate of apparition : nb of sites Latent=>Infectant
        #print(BOXL[sdate:ldate,])
        pos_rapp = which(BOXL[,"sumTEEQ"]>=1&BOXL[,"XLAT"]>0)
        #print(pos_rapp)
        if(length(pos_rapp)>0){
          rapp=sum(BOXL[pos_rapp,"XLAT"])
          BOXL[pos_rapp,c("sumTEEQ","XLAT")]=c(NA,NA)
          }
          else {
          rapp=0
        }

        # rrem: rate of removal : nb of sites Infectant=>removed
        # rapp: rate of apparition : nb of sites Latent=>Infectant
        #print(BOXI[sdate:ldate,])
        pos_rrem = which(BOXI[,"sumTEEQ"]>=1&BOXI[,"XINF"]>0)
        #print(pos_rrem)
        if(length(pos_rrem)>0){
          rrem=sum(BOXI[pos_rrem,"XINF"])
          BOXI[pos_rrem,c("sumTEEQ","XINF")]=c(NA,NA)
          }
          else {
          rrem=0
        }

        XVAC[day+1] <- XVAC[day]-rocc
        XLAT[day + 1] <-XLAT[day] +rocc - rapp
        XINF[day+1] <- XINF[day] + rapp - rrem
        XCTR[day+1] <- XCTR[day] + rrem
        BOXL[1:(day+1),"sumTEEQ"] = BOXL[1:(day+1),"sumTEEQ"] + weather$TEEQ_L[day]
        BOXL[day+1,"XLAT"] = +rocc
        BOXI[1:(day+1),"sumTEEQ"] = BOXI[1:(day+1),"sumTEEQ"] + weather$TEEQ_I[day]
        BOXI[day+1,"XINF"] = +rapp

    }
    
    XTO1 = XLAT+XINF+XCTR
    XSEV = XINF+XCTR
    severity=XSEV/(XLAT+XINF+XCTR+XVAC)
    return(list(sim=data.frame(day = sdate:ldate, DACE = ((sdate:ldate)-sdate), XVAC = XVAC[sdate:ldate], XLAT = XLAT[sdate:ldate], XINF = XINF[sdate:ldate],XCTR = XCTR[sdate:ldate],XTO1=XTO1[sdate:ldate],  XSEV= XSEV[sdate:ldate],severity=severity[sdate:ldate]), param=c(nlpd=nlpd,nipd=nipd,dmfr=dmfr,SITE0 = SITE0)))
}

################################################################################
graph_epid_s=function(out,typel="s",all=TRUE, param=TRUE){
# a tip to show value=0 on the log plot
sim=out$sim
plot(sim$day,sim$XSEV, type=typel,lty=1, lwd=3, xlab="Time after first infection (day)", ylim=c(0,10^30),ylab="X (nb of infectious sites)")
if(all){
  lines(sim$day,sim$XLAT, type=typel,lty=1, lwd=2, col="orange")
  lines(sim$day,sim$XINF, type=typel,lty=1,lwd=2, col="red")
  lines(sim$day,sim$XCTR, type=typel,lty=1,lwd=2, col="grey")
  }
if(param){
text(seq(30,100,length.out=length(out$param)),max(sim$XSEV),names(out$param),cex=0.75)
text(seq(30,100,length.out=length(out$param)),0.95*max(sim$XSEV),(out$param),cex=0.75)
}
par(ann=FALSE,new=TRUE)
plot.default(sim$day,sim$severity,axes=FALSE,type="l",ylim=c(0,1),ylab='',lty=3, lwd=3,col="darkgrey")
axis(side=4)
mtext("severity (%)", side=4,line=2)
mtext(paste("AUDPC = ",round(sum(sim$severity,1))),side=3,line=0,cex=0.8)

legend("topleft",legend=c("XSEV","XLAT","XINF","XCTR","severity"),lty=c(1,1,2,1,3),
lwd=c(2,2,2,2,3),col=c("black","orange","red","grey","darkgrey"),
cex=0.6, ncol = 2)
return(invisible(out))
}
################################################################################
graph_epid=function(out,typel="s",all=TRUE, param=TRUE){
# a tip to show value=0 on the log plot
sim=out$sim
sim[sim==0]=0.5
plot(sim$DACE,sim$XSEV,log = "y", type=typel,lty=1, lwd=3,ylim=c(0.1,10^14),yaxt = "n", xlab="Time after first infection (day)", ylab="X (nb of infectious sites)")
axis(2, at=c(0.5,1,100,10^4,10^6,10^8,10^10,10^12), label=c("0","1","100","10^4","10^6","10^8","10^10","10^12"))
if(all){
  lines(sim$DACE,sim$XCTR, type=typel,lty=1, lwd=1, col="darkgrey")
  lines(sim$DACE,sim$XINF, type=typel,lty=2,lwd=2)
  }
if(param){
text(seq(5,30,length.out=length(out$param)),10^14,names(out$param),cex=0.75)
text(seq(5,30,length.out=length(out$param)),10^13.5,(out$param),cex=0.75)
}
}
################################################################################

library(ZeBook)
#sdate = as.numeric(format(strptime("1997-05-25", "%Y-%m-%d", tz=""), "%j"))
weather=subset(weather_FranceWest, WEYR==1997 & idsite==39)

################################################################################
# reproduce simulation from Zadoks 1971 (Phytopathology 61:441-598, 19 May 1971)
weather$TMIN = 10
weather$TMAX = 10
#dmfr : daily multiplication factor - number of effective spores produced by lesion
#nlpd : latent period (in degree.day) - latency for Puccinia triticina : ~10 days at 20°C
#nipd : infectious period (in degree.day) - for Puccinia triticina : ~20 days at 20°C

# unlimited number of Site
SITE0= 10^30
sdate = 1
ldate = 150
out=zakoks.original.model(nlpd=16*10,nipd=30*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate )
################################################################################
# Latent period varies from 8-14 days in the field with temperatures of 10-25 C.
# With temperatures of 5-10 C latent periods of several weeks are common.
# typical infectious period for leaf rust about 30 days
par(mfrow=c(3,1), mar=c(4,4,2,4))
graph_epid_s(zakoks.original.model(nlpd=8*10,nipd=30*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate, XLAT0=SITE0*1e-14),typel="l")
graph_epid_s(zakoks.original.model(nlpd=16*10,nipd=30*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate, XLAT0=SITE0*1e-14),typel="l")
graph_epid_s(zakoks.original.model(nlpd=24*10,nipd=30*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate, XLAT0=SITE0*1e-14),typel="l")

################################################################################
# Figure 7.
par(mfrow=c(3,1), mar=c(3,3,0,0))
graph_epid(zakoks.original.model(nlpd=4*10,nipd=1*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 ))
graph_epid(zakoks.original.model(nlpd=4*10,nipd=2*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 ),typel="l")
graph_epid(zakoks.original.model(nlpd=4*10,nipd=4*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 ),typel="l")

################################################################################
# Figure 8.
dev.new()
par(mfcol=c(2,2), mar=c(3,3,1,0))
#A)Various latent period
graph_epid(zakoks.original.model(nlpd=4*10,nipd=10*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 ),typel="l", all=FALSE, param=FALSE)
title("Various latent period",cex.main=0.75)
for (nlpd in c(4,8,12,16)*10){
    sim=zakoks.original.model(nlpd=nlpd,nipd=10*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 )$sim
    lines(sim$DACE,sim$XSEV,lty=1, lwd=3)
    text(53, sim$XSEV[sim$DACE==50], nlpd,cex=0.8, col="darkgrey")
    }
#B)Various infectious period
graph_epid(zakoks.original.model(nlpd=4*10,nipd=1*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 ),typel="s", all=FALSE, param=FALSE)
title("Various infectious period",cex.main=0.75)
for (nipd in c(2,3,4,16)*10){
    sim=zakoks.original.model(nlpd=4*10,nipd=nipd,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 )$sim
    lines(sim$DACE,sim$XSEV,lty=1, lwd=3)
    text(50, sim$XSEV[sim$DACE==50], nipd,cex=0.8, col="darkgrey")
}

#C)Various infectious period on  cumulative total of removal
graph_epid(zakoks.original.model(nlpd=4*10,nipd=10*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 ),typel="l", all=FALSE, param=FALSE)
title("Various infectious period on XCTR",cex.main=0.75)
for (nipd in c(4,8,12,16)*10){
    sim=zakoks.original.model(nlpd=4*10,nipd=nipd,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 )$sim
    lines(sim$DACE,sim$XCTR,lty=1, lwd=1, col="grey")
    text(50, sim$XCTR[sim$DACE==50], nipd,cex=0.8, col="darkgrey")
}
#D)Various multiplication factor
graph_epid(zakoks.original.model(nlpd=4*10,nipd=8*10,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 ),typel="l", all=FALSE, param=FALSE)
title("Various infectious period",cex.main=0.75)
for (dmfr in c(4,8,12,16)){
    sim=zakoks.original.model(nlpd=4*10,nipd=8*10,dmfr=dmfr,SITE0 = SITE0, weather, sdate = sdate, ldate = sdate+55 )$sim
    lines(sim$DACE,sim$XSEV,lty=1, lwd=3)
    text(50, sim$XSEV[sim$DACE==50], dmfr,cex=0.8, col="darkgrey")
}



# Figure 9B
#SITE0= 5*10^9 # from Zadoks 1971


# end of file