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)

kaska99.SCL <- subset(stages, subset=(Species == "Caretta caretta"), select=c("Stage", "SCL_Mean_mm", "SCL_SD_mm", "Days_Begin", "Days_End"))

kaska99.SCL[kaska99.SCL$Stage==31, "Days_Begin"] <- 51
kaska99.SCL[kaska99.SCL$Stage==31, "Days_End"] <- 62
kaska99.SCL <- na.omit(kaska99.SCL)

kaska99.SCL <- cbind(kaska99.SCL, 
                     Days_Mean=(kaska99.SCL[, "Days_Begin"]+kaska99.SCL[, "Days_End"])/2)
kaska99.SCL <- cbind(kaska99.SCL, 
                     Days_SD=(kaska99.SCL[, "Days_End"]-kaska99.SCL[, "Days_Begin"])/4)


Gompertz <- function(x, par) {
  K <- par["K"]
  rT <- par["rT"]
  y0 <- par["y0"]
  y <- abs(K)*exp(log(abs(y0)/abs(K))*exp(-rT*x))
  return(y)
}

library(mvtnorm)

ML.Gompertz.2D <- function(x, par) {
  par <- abs(par)
  y <- Gompertz(x, par)
  L <- 0
  for (i in seq_along(y)) {
    sigma <- matrix(c(kaska99.SCL$SCL_SD_mm[i]^2, 0, 0, kaska99.SCL$Days_SD[i]^2), nrow=2, byrow=TRUE, 
                    dimnames=list(c("SCL_SD_mm", "Days_SD"), c("SCL_SD_mm", "Days_SD")))
    L <- L -dmvnorm(x=c(SCL_SD_mm=kaska99.SCL$SCL_Mean_mm[i], Days_SD=kaska99.SCL$Days_Mean[i]), 
                    mean= c(SCL_SD_mm=y[i], Days_SD=kaska99.SCL$Days_Mean[i]), sigma=sigma, log=TRUE)
  }
  return(L)
}

parIni <- structure(c(48.66977358, 0.06178453, 0.38640902), 
                    .Names = c("K", "rT", "y0"))

fitsize.SCL.2D <- optim(parIni, ML.Gompertz.2D, x=kaska99.SCL[, "Days_Mean"])

cor(kaska99.SCL$SCL_Mean_mm, Gompertz(kaska99.SCL[, "Days_Mean"], par=fitsize.SCL.2D$par))^2

Bayes.Gompertz.2D <- function(data, x) {
  x <- abs(x)
  y <- Gompertz(data, x)
  L <- 0
  for (i in seq_along(y)) {
    sigma <- matrix(c(kaska99.SCL$SCL_SD_mm[i]^2, 0, 0, kaska99.SCL$Days_SD[i]^2), nrow=2, byrow=TRUE, 
                    dimnames=list(c("SCL_SD_mm", "Days_SD"), c("SCL_SD_mm", "Days_SD")))
    L <- L -dmvnorm(x=c(SCL_SD_mm=kaska99.SCL$SCL_Mean_mm[i], Days_SD=kaska99.SCL$Days_Mean[i]), 
                    mean= c(SCL_SD_mm=y[i], Days_SD=kaska99.SCL$Days_Mean[i]), sigma=sigma, log=TRUE)
  }
  return(L)
}

pMCMC <- structure(list(Density = c("dunif", "dunif", "dunif"), 
                        Prior1 = c(0, 0, 0), Prior2 = c(90, 1, 2), 
                        SDProp = c(1, 1, 1), 
                        Min = c(0, 0, 0), Max = c(90, 1, 2), 
                        Init = fitsize.SCL.2D$par), 
                   .Names = c("Density", "Prior1", "Prior2", "SDProp", "Min", "Max", "Init"), 
                   row.names = c("K", "rT", "y0"), class = "data.frame")

mcmc_run.2D <- MHalgoGen(n.iter=50000, parameters=pMCMC, data=kaska99.SCL[, "Days_Mean"], 
                         likelihood=Bayes.Gompertz.2D, n.chains=1, n.adapt=100, thin=10, 
                         trace=100, 
                         adaptive = TRUE)

plot(mcmc_run.2D, parameters="K", what="MarkovChain")
plot(mcmc_run.2D, parameters="rT", what="MarkovChain")
plot(mcmc_run.2D, parameters="y0", what="MarkovChain")

1-rejectionRate(as.mcmc(mcmc_run.2D))

par <- mcmc_run.2D$resultMCMC[[1]]

outsp <- t(apply(par, MARGIN = 1, FUN=function(x) Gompertz(0:70, par=x)))

rangqtiles <- apply(outsp, MARGIN=2, function(x) {quantile(x, probs=c(0.025, 0.5, 0.975))})

# pdf(file = "Figure Table.pdf", width=12, height = 7, pointsize = 18)
par(mar=c(4, 4, 2, 1))

plot_errbar(x=kaska99.SCL[, "Days_Mean"], y=kaska99.SCL[, "SCL_Mean_mm"], 
            errbar.y = 2*kaska99.SCL[, "SCL_SD_mm"], bty="n", las=1, 
            ylim=c(0, 50),
            xlim=c(0, 70), pch=3, 
            xlab="Days of incubation", ylab="SCL in mm", 
            x.plus = kaska99.SCL[, "Days_End"], 
            x.minus = kaska99.SCL[, "Days_Begin"])

library("plotrix")

draw.ellipse(x=kaska99.SCL[, "Days_Mean"], y=kaska99.SCL[, "SCL_Mean_mm"], 
             a=(kaska99.SCL[, "Days_End"]-kaska99.SCL[, "Days_Begin"])/2, 
             b=kaska99.SCL[, "SCL_SD_mm"]*2, border = NA, 
             col=rgb(0.5, 0.5, 0.5, 0.5))




lines((0:70)[1:61], rangqtiles["2.5%", ][1:61], lty=2, lwd=2)
lines((0:70)[1:61], rangqtiles["97.5%", ][1:61], lty=2, lwd=2)
lines((0:70)[1:61], rangqtiles["50%", ][1:61], lty=3, lwd=2)

# dev.off()

text(x=ScalePreviousPlot(x=0.2, y=0.8)$x, y=ScalePreviousPlot(x=0.2, y=0.8)$y, pos=4, labels=as.expression(bquote(r^2 ~ "=" ~ .(format(cor(kaska99.SCL$SCL_Mean_mm, Gompertz(kaska99.SCL[, "Days_Mean"], par=fitsize.SCL.2D$par))^2, digits = 3, scientific = FALSE, big.mark = " ")) )))

text(x=ScalePreviousPlot(x=0.7, y=0.1)$x, y=ScalePreviousPlot(x=0.7, y=0.1)$y, pos=4, labels=as.expression(bquote("K =" ~ .(format(x = unname(fitsize.SCL.2D$par["K"]), digits = 4)))))
text(x=ScalePreviousPlot(x=0.7, y=0.15)$x, y=ScalePreviousPlot(x=0.7, y=0.15)$y, pos=4, labels=as.expression(bquote(r ~ "" [K] ~ "=" ~ .(format(x = unname(fitsize.SCL.2D$par["K"]/39.33), digits = 4)))))

text(x=ScalePreviousPlot(x=0.7, y=0.20)$x, y=ScalePreviousPlot(x=0.7, y=0.20)$y, pos=4, 
     labels=as.expression(bquote(y * "" [0] ~ "=" ~ .(format(x = unname(fitsize.SCL.2D$par["y0"]), digits = 4)))))

text(x=ScalePreviousPlot(x=0.7, y=0.25)$x, y=ScalePreviousPlot(x=0.7, y=0.25)$y, pos=4, 
     labels=as.expression(bquote(r ~ "" [y] * "" [0] ~ "=" ~ .(format(x = unname(fitsize.SCL.2D$par["y0"]/39.33), digits = 4)))))

title(expression(italic("Caretta caretta")))