# Simulation de l'exercice propagation d'épidémie dans la classe.
# auteur VL.
rm(list=ls())

library(deSolve)
library(nls2)    # pas obligatoire (on peut utiliser nls au lieu de nls2)

# ------------------ DONNEES
  # Chemin du répertoire contenant les fichiers de données :
path = "C:/Users/quent/Desktop/cholera"
setwd(path)                      #Change le répertoire par défaut

d <- read.table(paste(path, "/data1_avec_t.txt", sep=""),
                  header=TRUE, sep="\t", na.strings="NA", dec=".")
d
plot(d$t, d$I1, xlab="temps", ylab="Nb infectés")

# ------------------ MODELE
  # modélisation par compartiments :
  # V=vecteur d'états, p=vecteur de paramètres
deriv <- function(t,y,p)
{
  with(as.list(c(p, y)), {
    N1 = S1 + I1 + D1 + R1
    dS1 = - beta*S1*I1/N1 - gamma*C/(k+C)*S1
    dI1 = beta*S1*I1/N1 + gamma*C/(k+C)*S1 - nu*I1
    dD1 = l*nu*I1
    dR1 = (1-l)*nu*I1
    dC = eta*gamma*I1/N-d*C
    list(c(dS1,dI1,dD1,dR1,dC))
  })
}

# ------------------ SIMULATION
  # test d'une simulation :
N = 550                          # nombre d'élèves
Tmax = max(d$t)                 # temps maximal
par <- list(beta = 0.1, gamma=5e-3, eta = 1, d = 0.00001, k=1e-3, nu=1/15, l=0.2) # vecteurs de paramètres
V <- c(S1=N-1, I1=1, D1=0, R1=0, C=1e-3)              # valeurs initiales du vecteur d'états
comparts <- ode(V, 1:Tmax, deriv, parms=par)
plot(comparts[,"time"], comparts[,"I1"],type="l",xlab="Temps",ylab="Ih(t)", lwd=4, col="red")
plot(comparts[,"time"], comparts[,"C"],type="l",xlab="Temps",ylab="Ih(t)", lwd=4, col="green")

# ------------------ ESTIMATION
# Fonction pour l'estimation des paramètres
# Retourne un vecteur I simulé, qui doit être de même taille que celui des données.
model_tofit <- function(beta, gamma, eta, times){
  V <- c(S1=N-1, I1=1, D1=0, R1=0, C=0)
  Fsim <- ode(V, times, deriv, parms=list("beta"=beta, "gamma"=gamma, "eta"=eta, "d"=0, "k"=1e-3, "nu"=1/15, "l"=0.2))
  return(Fsim[,"I1"])  # note : il faut que le vecteur times démarre à 1, sinon le temps initial est incorrect
                      # sinon mettre c(1, times) au lieu de times dans ode et enlever le 1er terme du vecteur renvoyé.
}


  # estimation des paramètres p :
obs <- d$I1[which(!is.na(d$I1))]
res<-nls2(obs~model_tofit(beta, gamma, eta, T), start=list(beta = 0.1, gamma=1e-3, eta = 0.5), data=list(T=d$t[which(!is.na(d$I1))]), trace=TRUE, nls.control(warnOnly=TRUE), algorithm = "port", lower=c(0,0,0,0,0))
  # data=list() sert pour garder des paramètrse constants
  # warn.only=TRUE permet de stocker le dernier résultat obtenu même si l'algo échoue
  # l'utilisation de l'algo "port" permet de fixer des bornes à l'espace de recherche des paramètres (ici on veut les garder positifs par exemple)
  # nls2 permet de relancer plusieurs fois l'optimisation en partant de différents points initiaux
  # (note : algos brute-force et grid-search ou random-search font seulement des échantillonages de la fonction objectif)

summary(res)
plot(d$t[which(!is.na(d$I1))], obs, xlab="Temps",ylab="I1(t)", col="red")
lines(d$t[which(!is.na(d$I1))], fitted(res),lwd=2, col="blue")

# # facultatif : prédiction et confiance
# confint(res)  # intervalles de confiance sur les estimées des paramètres
# conf<-predict(as.lm(res), interval="confidence")
# lines(d$t, conf[,"lwr"], col="green")
# lines(d$t, conf[,"upr"], col="green")
# conf<-predict(as.lm(res), interval="prediction")
# lines(d$t, conf[,"lwr"], col="dark green")
# lines(d$t, conf[,"upr"], col="dark green")

# ------------------ UTILISATION DES VALEURS ESTIMEES


model_pred <- function(beta, gamma, eta, times){
  V <- c(S1=N-1, I1=1, D1=0, R1=0, C=0)
  Fsim <- ode(V, times, deriv, parms=list("beta"=beta, "gamma"=gamma, "eta"=eta, "d"=0, "k"=1e-3, "nu"=1/15, "l"=0.2))
  return(Fsim)  # note : il faut que le vecteur times démarre à 1, sinon le temps initial est incorrect
  # sinon mettre c(1, times) au lieu de times dans ode et enlever le 1er terme du vecteur renvoyé.
}

beta_estim = coefficients(res)["beta"]
gamma_estim = coefficients(res)["gamma"]
eta_estim = coefficients(res)["gamma"]
  # prédiction :
newTmax=200
pred = model_pred(beta_estim, gamma_estim , eta_estim, 1:newTmax)
plot(1:newTmax, pred[,"I1"], type="l", xlab="Temps",ylab="I1(t)", col="red")
lines(1:newTmax, pred[,"D1"], xlab="Temps",ylab="I1(t)", col="black")