
#====================================================
#       Project: PhenoBR
#       Authors: Tiphaine Mace, Eliel Gonzalez Garcia, Gyorgy Kover, Dominique Hazard and Masoomeh Taghipoor    
#       Date: December 2020
#====================================================

rm(list=ls(all=TRUE)) 
dev.off() 

    library(lattice)
    library(plotrix)
    library(deSolve)
    library(FactoMineR)
    library(Hmisc)
    library("corrplot")

#======================================
load("MydataBCS.Rdata") #original data with BCSm and initial valu of BCS
load("fichierfinalGroups.Rdata")
    #1. pardata: final results with appartenance to different groups of parities
    #2. Gi : different groups from 1 to 6 function of parities
    #3. char2: file with time related parameters 
    #4. param.cor: another representation of final results that was used for correlation analysis

#==========================================================
#======== simulation for one or several animals ============* 
#==========================================================

## the perturbation ODE function - 

BCSModel2 <- function(Time, State, Pars) {
  kb1 = Pars[1]
  kb2 = Pars[2]
  kb3 = Pars[3]
  kp1 = Pars[4]
  kp2 = Pars[5]
  kp3 = Pars[6]
  if (is.na(tb1)==TRUE) { #absence of 1st pertrubation
    tb1=0
    te1=0
    }
  if (is.na(tb2)==TRUE) { #absence of 2nd pertrubation
    tb2=tb2N
    te2=0
  } 
  if (is.na(tb3)==TRUE) { #absence of 3d pertrubation
    kp3=0
    tb3=tb3N
    te3=0
    kb3=0
  }
  #
  with(as.list(c(State, Pars)), {
    #
    dBCS1     =  kb1*(BCSm -BCS1)*Time.pos.bis(Time, tb1,tb2)-
      kp1*(Pmax - P1)*Time.pos.bis(Time, tb1,te1)
    dBCS2 =  kb2*(BCSm -BCS2)*Time.pos.bis(Time, tb2,tb3)-
      kp2*(Pmax - P2)*Time.pos.bis(Time, tb2,te2)
    dBCS3 =  kb3*(BCSm -BCS3)*Time.pos.ter(Time, tb3)-
      kp3*(Pmax - P3)*Time.pos.bis(Time, tb3,te3)    
    dP1    =   kp1*(Pmax - P1)*Time.pos.bis(Time, tb1,te1)
    dP2    =   kp2*(Pmax - P2)*Time.pos.bis(Time, tb2,te2)
    dP3    =   kp3*(Pmax - P3)*Time.pos.bis(Time, tb3,te3)    
    dBCS = kb1*(BCSm -BCS1)*Time.pos.bis(Time, tb1,tb2)- #total variation
      kp1*(Pmax - P1)*Time.pos.bis(Time, tb1,te1)+
      kb2*(BCSm -BCS2)*Time.pos.bis(Time, tb2,tb3)-
      kp2*(Pmax - P2)*Time.pos.bis(Time, tb2,te2)+
      kb3*(BCSm -BCS3)*Time.pos.ter(Time, tb3)-
      kp3*(Pmax - P3)*Time.pos.bis(Time, tb3,te3)
    
    return(list(c( dBCS1,dBCS2,dBCS3,dP1,dP2,dP3,dBCS)))
  })
}

Time.pos.bis = function(Time, t1,t2){ 
  ifelse(Time-t1>=0 & Time-t2<=0,1,0)
}

Time.pos.ter = function(Time, t1){
  ifelse(Time-t1>=0,1,0)
}

#days of measures for p1, p2 and p3 
  d.p1 = c(0,40,80,160,190,230,250,310,330)
  d.p2 = (2-1)*360+d.p1
  d.p3 = (3-1)*360+d.p1
  
  tb2N=d.p1[8]+150
  tb3N=d.p2[8]+150

# Plot
  
par(mfrow=c(2,2))
ID=c(param.cor$ID[param.cor$Group=="G1"][20], #select IDs as a function of appartenance to a group
     param.cor$ID[param.cor$Group=="G13"][2],
     param.cor$ID[param.cor$Group=="G23"][20],
    param.cor$ID[param.cor$Group=="G123"][20])

for(i in 1:length(ID)){
  print(i)
  Pmax = 2
  BCSm= param.cor$BCSm[param.cor$ID==ID[i]]
  BCSopt=BCSm-0.5
  tb1=param.cor$tb1[param.cor$ID==ID[i]]
  tb2=param.cor$tb2[param.cor$ID==ID[i]]
  tb3=param.cor$tb3[param.cor$ID==ID[i]]
  te1=param.cor$te1[param.cor$ID==ID[i]]
  te2=param.cor$te2[param.cor$ID==ID[i]]
  te3=param.cor$te3[param.cor$ID==ID[i]]
  
# solve the model
  
  char.ID = char2[char2$Ewe.ID%in%ID[i],]
  
  if (char.ID$bcsinit[1]>0){
    yinit <- c(BCS1=char.ID$bcsinit[1],BCS2=char.ID$bcsinit[2],
               BCS3=char.ID$bcsinit[3], P1=0, P2=0,P3=0,BCS=char.ID$bcsinit[1])
  }else if (char.ID$bcsinit[2]>0 && char.ID$bcsinit[1]==0){
    yinit <- c(BCS1=char.ID$bcsinit[1],BCS2=char.ID$bcsinit[2],
               BCS3=char.ID$bcsinit[3], P1=0, P2=0,P3=0,BCS=char.ID$bcsinit[2])
  }else if (char.ID$bcsinit[3]>0 && char.ID$bcsinit[1]==0 && char.ID$bcsinit[2]==0){
    yinit <- c(BCS1=char.ID$bcsinit[1],BCS2=char.ID$bcsinit[2],
               BCS3=char.ID$bcsinit[3], P1=0, P2=0,P3=0,BCS=char.ID$bcsinit[3])
  }
  
    times <- seq(from = 0, to = max(MydataBCS$day[MydataBCS$Ewe==ID[i]]), by = 1)
  
  ##  
  ##----------------------------- 
  par=param.cor[param.cor$ID==ID[i],2:7]
  par[is.na(par)==TRUE]=0
 
  yout <- ode(y= yinit, times = times, func = BCSModel2,
              parms = par,
              atol = 1e-10, rtol = 1e-10,method="rk4")
 
  plot(yout[ , 1], yout[ ,8], type = "l", xlab = "days of age", ylab = "BCS",
          lwd = 2, ylim = c(1,5),cex.lab=1.5,xlim=c(1,1100))
  points(MydataBCS$day[MydataBCS$Ewe==ID[i]] , MydataBCS$BCS[MydataBCS$Ewe==ID[i]], col="darkgrey")
   abline(v=c(tb1, te1,tb2,te2,tb3,te3),lty=2, col="grey")
   text(c(tb1, te1,tb2,te2,tb3,te3),
        c(1.5,1.5,1.5,1.5,1.5,1.5),c("tb1", "te1","tb2","te2","tb3","te3"))
   abline(h=BCSm,lty=2, col="grey")
  text(15,c(BCSm),c("BCSm"))
  
}

