################################################################################
# Course SMACH - January 2014
# my first dynamics models
# version $Date:: 2014-04-30 1#$ by $Author: brun $(v0: 2013-10-9)
################################################################################
#1a) Logistic model, an exact solution exists.
# x : fraction of disease host material
# r: proportionality factor, named apparent infection rate
logistic.model = function (r, X0, Time){

    X = X0 / ((1-X0 )*exp(-r*Time)+X0)
    return(data.frame(Time = Time, X = X))
}

sim1a = logistic.model(0.1, 0.001, 0:100)
dev.new()
plot(sim1a, ylab="X, fraction of diseased host material",type="l", lty=2, main="logistic model")
legend("topleft", legend=c("exact solution"), lty=c(2))

#1b) Logistic model, integration of Euler with a time step dt.
# x : fraction of disease host material
# r: proportionality factor, named apparent infection rate
logistic.dt.model = function (r, X0, duration,dt){
    # K : number of time steps
    K = duration/dt+1
    # 1 state variable, as 1 vector initialized to NA
    X = rep(NA, K)
    # Initialize state variable
    X[1] = X0
    # time variable
    Time = 0
    k = 1
    while (Time < duration){
        # compute rate of change
        dX = r * X[k] * (1 - X[k]) * dt
        # update state variable
        X[k + 1] = X[k] + dX

    # Update time
    Time = Time + dt
    k = k + 1
    }
    return(data.frame(Time = ((1:K)-1)*dt, X = X[1:K]))
}

sim1b = logistic.dt.model(0.1, 0.001, 10, 1)
dev.new()
plot(sim1b, ylab="X, fraction of diseased host material",type="l", main="logistic model")
lines(sim1a, lty=2)
legend("topleft", legend=c("exact solution", "numerical solution"), lty=c(2,1))

################################################################################
#2) Van der Planck (1963) logistic model, integration of Euler with a time step dt.
# x : fraction of disease host material
# Rc: 
# p: latent period (day)
# i: infectious period (day)
# When p=0 and i->inf (i> duration), this model give same result as logistic model.
vanderplanck.model = function (Rc, p, i, X0, duration,dt){
    # K : number of time steps
    K = duration/dt+1
    # 1 state variable, as 1 vector initialized to NA
    X = rep(NA, K)
    # Initialize state variable
    X[1] = X0
    # time variable
    Time = 0
    k = 1
    while (Time < duration){
        # compute rate of change
        if ((Time-p)<0) {X_kt_p = 0} else {X_kt_p = X[(Time-p)/dt + 1] }
        if ((Time-i-p)<0) {X_kt_i_p = 0} else {X_kt_i_p = X[(Time-i-p)/dt + 1] }   

        dX = Rc * (X_kt_p - X_kt_i_p) * (1 - X[k]) * dt
        # update state variable
        X[k + 1] = X[k] + dX
        # Update time
        Time = Time + dt
        k = k + 1
        }
    return(data.frame(Time = ((1:K)-1)*dt, X = X[1:K]))
}

sim2a = vanderplanck.model(1.5, p=10, i=5, X0=0.001, 50, 1)
sim2b = vanderplanck.model(1.5, p=10, i=5, X0=0.001, 50, 0.01)

dev.new()
plot(sim2a, ylab="X, fraction of diseased host material",type="l", main="logistic model with latency and infectious period")
lines(sim2b, lty=2)
legend("topleft", legend=c("dt=1", "dt=0.01"), lty=c(1,2))

################################################################################
# application to Puccinia triticina
# some data from Bancal et al. (2008)
setwd("D:/data/FormationSmach/data/leaf_rust")
rust.data = read.table("leaf_rust_wheat.data.csv", sep=";", header=TRUE, dec=".", stringsAsFactors=FALSE)
# select data for low Nitrogen only
sdata=subset(rust.data, Nitrogen=="low" & Year==1997)
sdata$day_anth=sdata$day-as.numeric(format(strptime(sdata$date_anth, "%Y-%m-%d", tz=""),"%j"))

dev.new()
plot(sdata$day_anth, sdata$DS/100, xlim=c(0,50), ylim=c(0,1), ylab="severity (-)", pch=1, main="data for leaf_rust and logistic models")
#latent period, depend on cultivar : p = 7 to 13 days
# infectious period : i = 10 to 40 (in controlled condition)

#
sim3a = logistic.model(r=0.3, X0=0.001, Time=0:50)
lines(sim3a, lw=2,lty=1)

# after several try for i, p, X0 and Rc
sim3b = vanderplanck.model(Rc=5, p=9, i=15, X0=0.001, duration=50, dt=0.01)
lines(sim3b, lw=2,lty=2)

legend("topleft", legend=c("data", "logistic", "van der planck logistic"), lty=c(NA,1,2), pch=c(1,NA, NA))

################################################################################
# end of file