if (FALSE) {
  # Run this part of the code only once
  install.packages("embryogrowth")
  
  install.packages("https://hebergement.universite-paris-saclay.fr/marcgirondot/CRAN/HelpersMG.tar.gz", repos=NULL, type="source")
  install.packages("https://hebergement.universite-paris-saclay.fr/marcgirondot/CRAN/embryogrowth.tar.gz", repos=NULL, type="source")
}

library(embryogrowth)
data(nest)

# Data from Girondot, M. & Kaska, Y. 2014. A model to predict the thermal reaction norm for the embryo growth rate from field data. J. Therm. Biol. (IF 2.361), 45, 96-102.
head(nest)


formated <- FormatNests(data=nest, previous=NULL, col.Time="Time", 
                        hatchling.metric.mean=39.33, 
                        hatchling.metric.sd=1.92)

#######################
# Gompertz model of SCL
#######################

# See nest video to understand the origin of these parameters
pfixed <- c(rK=1.208968)
M0 = 0.3470893 

# These values are just a first set of parameters
x4 <- c('DHA' = 109.31113503282113, 'DHH' = 617.80695919563857, 
        'T12H' = 306.38890489505093, 'Rho25' = 229.37265815800225)

plotR(parameters = x4, curve="ML")

resultNest_4p_SSM <- searchR(parameters=x4, fixed.parameters=pfixed, 
                             temperatures=formated, 
                             integral=integral.Gompertz, M0=M0)

pMCMC <- TRN_MHmcmc_p(resultNest_4p_SSM, accept=TRUE)

# Take care, it can be very long; several days
resultNest_mcmc_4p_SSM <- GRTRN_MHmcmc(result=resultNest_4p_SSM , 
                                       adaptive = TRUE          , 
                                       WAIC=TRUE                , 
                                       parametersMCMC=pMCMC     , 
                                       n.iter=10000             , 
                                       n.chains = 1             ,  
                                       n.adapt = 100            , 
                                       thin=10                  , 
                                       trace=100               )

plotR(resultNest_4p_SSM, resultmcmc = resultNest_mcmc_4p_SSM, 
      curve = "MCMC quantiles", ylim=c(0, 8), show.hist = TRUE, 
      ylimH = c(0, 1), atH=c(0, 0.1, 0.2))

plot(x=resultNest_4p_SSM, resultmcmc = resultNest_mcmc_4p_SSM, 
     xlim=c(0,70), ylimT=c(22, 32), GTRN.CI="MCMC", 
     replicate.CI = 100, las=1, 
     ylimS=c(0,45), series=1, 
     embryo.stages="Caretta caretta.SCL")

inests<- info.nests(x=resultNest_4p_SSM, 
                            resultmcmc = resultNest_mcmc_4p_SSM, 
                            GTRN.CI="MCMC", 
                            replicate.CI = 100, 
                            out="summary", 
                            embryo.stages="Caretta caretta.SCL")
inests <- inests$summary[, !is.na(inests$summary[1, ])]


####################################
# linear model of progress of growth
####################################

pfixed <- NULL
M0 = 0

x4 <- c('DHA' = 64.868697530424186, 'DHH' = 673.18292743646771, 
        'T12H' = 400.90952554047749, 'Rho25' = 82.217237723502123)

resultNest_4p_SSM_Linear <- searchR(parameters=x4, fixed.parameters=pfixed, 
                                    temperatures=formated, integral=integral.linear, M0=M0, 
                                    hatchling.metric=c(Mean=39.33, SD=1.92)/39.33)

plotR(resultNest_4p_SSM_Linear, ylim=c(0, 2), scaleY= 100000, curve = "ML")

pMCMC <- TRN_MHmcmc_p(resultNest_4p_SSM_Linear, accept=TRUE)

# Take care, it can be very long; several days on old computer
resultNest_mcmc_4p_SSM_Linear <- GRTRN_MHmcmc(result=resultNest_4p_SSM , 
                                       adaptive = TRUE          , 
                                       WAIC=TRUE                , 
                                       parametersMCMC=pMCMC     , 
                                       n.iter=10000             , 
                                       n.chains = 1             ,  
                                       n.adapt = 100            , 
                                       thin=10                  , 
                                       trace=100               )

plotR(resultNest_4p_SSM_Linear, resultmcmc = resultNest_mcmc_4p_SSM_Linear, 
      curve = "MCMC quantiles", ylim=c(0, 2), show.hist = TRUE, 
      ylimH = c(0, 1), atH=c(0, 0.1, 0.2))

plot(x=resultNest_4p_SSM_Linear, resultmcmc = resultNest_mcmc_4p_SSM_Linear, 
     xlim=c(0,70), ylimT=c(22, 32), GTRN.CI="MCMC", 
     replicate.CI = 100, las=1, 
     ylimS=c(0,1), series=1, 
     embryo.stages="Generic.ProportionDevelopment")

inests_linear <- info.nests(x=resultNest_4p_SSM_Linear, 
                            resultmcmc = resultNest_mcmc_4p_SSM_Linear, 
                            GTRN.CI="MCMC", 
                            replicate.CI = 100, 
                            out="summary", 
                            embryo.stages="Generic.ProportionDevelopment")
inests_linear <- inests_linear$summary[, !is.na(inests_linear$summary[1, ])]


###################################
# Comparison of TSP for both models
###################################

plot(inests$TSP.TimeWeighted.temperature.mean, 
     inests_linear$TSP.TimeWeighted.temperature.mean, 
     xlab = "Average temperature within TSP for Gompertz model in °C", 
     ylab = "Average temperature within TSP for linear model in °C", las=1, bty="n", 
     xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

plot(inests$TSP.GrowthWeighted.temperature.mean, 
     inests_linear$TSP.GrowthWeighted.temperature.mean, 
     xlab = "Average temperature within TSP for Gompertz model in °C", 
     ylab = "Average temperature within TSP for linear model in °C", las=1, bty="n", 
     xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

plot(inests$TSP.begin.mean/24/60, 
     inests_linear$TSP.begin.mean/24/60, 
     xlab = "Begin TSP for Gompertz model in min", 
     ylab = "Begin TSP for linear model in min", las=1, bty="n", 
     xlim = c(15, 30), ylim = c(15, 30)
     )
segments(x0=15, y0=15, x1=30, y1=30, lty=2, col="red")

plot(inests$TSP.end.mean/24/60, 
     inests_linear$TSP.end.mean/24/60, 
     xlab = "End TSP for Gompertz model in min", 
     ylab = "End TSP for linear model in min", las=1, bty="n", 
     xlim = c(30, 50), ylim = c(30, 50)
)
segments(x0=30, y0=30, x1=50, y1=50, lty=2, col="red")


plot(inests$TSP.length.mean/24/60, 
     inests_linear$TSP.length.mean/24/60, 
     xlab = "TSP length for Gompertz model in days", 
     ylab = "TSP length for linear model in days", las=1, bty="n", 
     xlim = c(14, 23), ylim = c(14, 23)
)
segments(x0=14, y0=14, x1=23, y1=23, lty=2, col="red")