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)

PieauDorrizi1981 <- structure(list(Stage = c(8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26), 
                                   Mass.min = c(NA, NA, NA, NA, NA, NA, 75, 110, 170, 210, 270, 340, 400, 650, 1100, 1600, 2500, 3500, 4000), 
                                   Mass.max = c(NA, NA, NA, NA, NA, NA, 85, 120, 180, 230, 300, 380, 500, 700, 1300, 1800, 3000, 4000, 6000), 
                                   Mass.mean = c(NA, NA, NA, NA, NA, NA, 80, 115, 175, 220, 285, 360, 450, 700, 1200, 1700, 2750, 3750, 5000), 
                                   SCL.min = c(NA, NA, NA, NA, NA, NA, NA, 5, 6.5, 7, 8, 9.5, 10, 12, 14, 16, 20, 23, 24), 
                                   SCL.max = c(NA, NA, NA, NA, NA, NA, NA, 6, 7, 8, 9, 10, 11, 13, 16, 17, 22, 24, 27), 
                                   SCL.mean = c(NA, NA, NA, NA, NA, NA, NA, 5.5, 6.75, 7.5, 8.5, 9.75, 10.5, 12.5, 15, 16.5, 21, 23.5, 25.5), 
                                   Age.25.min = c(6, 7, 9, 12, 14, 17, 22, 24, 27, 29, 31, 33, 35, 39, 44, 48, 55, 65, 72), 
                                   Age.25.max = c(7, 8, 10, 14, 16, 20, 23, 25, 28, 30, 32, 34, 37, 40, 46, 50, 60, 70, 82), 
                                   Age.25.mean = c(6.5, 7.5, 9.5, 13, 15, 18.5, 22.5, 24.5, 27.5, 29.5, 31.5, 33.5, 36, 39.5, 45, 49, 57.5, 67.5, 77), 
                                   Age.30.min = c(4, 4, 5, 7, 8, 10, 11, 12, 15, 16, 18, 19, 20, 23, 27, 30, 35, 43, 53 ), 
                                   Age.30.max = c(4, 5, 6, 8, 9, 11, 12, 13, 16, 17, 19, 20, 22, 24, 28, 31, 40, 48, 60), 
                                   Age.30.mean = c(4, 4.5, 5.5, 7.5, 8.5, 10.5, 11.5, 12.5, 15.5, 16.5, 18.5, 19.5, 21, 23.5, 27.5, 30.5, 37.5, 45.5, 56.5)), 
                              row.names = c(NA, 19L), class = "data.frame")

# From Pieau, C. & Dorizzi, M. 1981. Determination of temperature sensitive stages for sexual differentiation of the gonads in embryos of the turtle, Emys orbicularis. J. Morphol. (IF 1.563), 170, 373-382.

par(mar=c(4, 4, 1, 1)+0.4)

plot(0, 0, bty="n", las=1, type="n", xlim=c(0, 80), 
     xlab="Incubation days", 
     ylab="Straight carapace length (mm)", ylim=c(0, 30))

par(new=TRUE)
with(PieauDorrizi1981, 
     plot_errbar(x=Age.25.mean, y=SCL.mean, 
                 y.plus = SCL.max,
                 y.minus = SCL.min, 
                 x.minus = Age.25.min, 
                 x.plus = Age.25.max, 
                 errbar.tick = 1/100,
                 bty="n", las=1, type="p", xlim=c(0, 80), 
                 xlab="", axes=FALSE, 
                 ylab="", ylim=c(0, 30))
)

par(new=TRUE)
with(PieauDorrizi1981, 
     plot_errbar(x=Age.30.mean, y=SCL.mean, 
                 y.plus = SCL.max,
                 y.minus = SCL.min, 
                 x.minus = Age.30.min, 
                 x.plus = Age.30.max, 
                 errbar.tick = 1/100,
                 bty="n", las=1, type="p", xlim=c(0, 80), 
                 xlab="", 
                 ylab="", ylim=c(0, 30), axes=FALSE, pch=19, add=TRUE)
)

text(x=42, y=28, labels="30 °C")
text(x=70, y=20, labels="25 °C")

##########################################################
# Fit of the Gompertz model using Gaussian bivariate model
##########################################################

Gompertz <- function(data, par) {
  temperatures <- data[, "temperatures"]
  K <- par["K"]
  rT195 <- par["rT195"]
  rT25 <- par["rT25"]
  rT285 <- par["rT285"]
  rT30 <- par["rT30"]
  rT18 <- par["rT18"]
  rT35 <- par["rT35"]
  y0 <- par["y0"]
  rT <- sapply(temperatures, 
               function(Temp) {switch(paste0("X", Temp), 
                                      X19.5=rT195, X18=rT18, X25=rT25, 
                                      X28.5=rT285, X30=rT30, X35=rT35)})
  y <- abs(K)*exp(log(abs(y0)/abs(K))*exp(-rT*data[, "age.mean"]))
  return(unname(y))
}


library(mvtnorm)

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

# Note that max-min is supposed to be 4.SD
dataSCL_ec <- data.frame(temperatures=c(rep(25, 12), rep(30, 12)), SCL.mean=rep(dataSCL[, "SCL.mean"], 2), 
                         SCL.SD=rep((dataSCL[, "SCL.max"]-dataSCL[, "SCL.min"])/4, 2), 
                         age.mean=c(dataSCL[, "Age.25.mean"], dataSCL[, "Age.30.mean"]), 
                         age.SD=c((dataSCL[, "Age.25.max"]-dataSCL[, "Age.25.min"])/4, 
                                  (dataSCL[, "Age.30.max"]-dataSCL[, "Age.30.min"])/4))

parIni <- c('K' = 28.857952895393229, 
            'rT25' = 0.042344844024783335, 
            'rT30' = 0.071489670058642554, 
            'y0' = 0.32244561415883893)

Gompertz(dataSCL_ec, par=parIni)
ML.Gompertz.2D(data=dataSCL_ec, par=parIni)

fitsize.SCL.2D <- optim(parIni, ML.Gompertz.2D, data=dataSCL_ec)

cor(dataSCL_ec$SCL.mean, Gompertz(dataSCL_ec, par=fitsize.SCL.2D$par))^2

Bayes.Gompertz.2D <- function(data, x) {
  x <- abs(x)
  y <- Gompertz(data, x)
  L <- NULL
  for (i in seq_along(y)) {
    sigma <- matrix(c(data[i, "SCL.SD"]^2, 0, 0, data[i, "age.SD"]^2), nrow=2, byrow=TRUE, 
                    dimnames=list(c("SCL_SD_mm", "Days_SD"), c("SCL_SD_mm", "Days_SD")))
    L <- c(L, dmvnorm(x=matrix(c(SCL_SD_mm=data[i, "SCL.mean"], Days_SD=data[i, "age.mean"]), ncol=2), 
                      mean= matrix(c(SCL_SD_mm=y[i], Days_SD=data[i, "age.mean"]), ncol=2), 
                      sigma=sigma, log=TRUE))
  }
  L_ec <- -sum(L)
  attributes(L_ec) <- list(WAIC=unname(L))
  return(L_ec)
}

pMCMC <- structure(list(Density = c("dunif", "dunif", "dunif", "dunif"), 
                        Prior1 = c(0, 0, 0, 0), Prior2 = c(90, 1, 1, 2), 
                        SDProp = c(1, 1, 1, 1), 
                        Min = c(0, 0, 0, 0), Max = c(90, 1, 1, 2), 
                        Init = fitsize.SCL.2D$par), 
                   .Names = c("Density", "Prior1", "Prior2", "SDProp", "Min", "Max", "Init"), 
                   row.names = c("K", "rT25", "rT30", "y0"), class = "data.frame")
pMCMC$SDProp <- c('K' = 0.68231795867409029, 
                  'rT25' = 0.00075281517848764133, 
                  'rT30' = 0.0013980853314770495, 
                  'y0' = 0.030885898552990329)
mcmc_run.2D <- MHalgoGen(n.iter=50000, parameters=pMCMC, data=dataSCL_ec, 
                         likelihood=Bayes.Gompertz.2D, n.chains=1, n.adapt=100, thin=10, trace=100, 
                         adaptive = TRUE, n.datapoints = nrow(dataSCL_ec), WAIC.out = TRUE)

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

d(as.parameters(mcmc_run.2D))

as.parameters(mcmc_run.2D, index = "quantile")

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

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

outsp25 <- t(apply(par, MARGIN = 1, FUN=function(x) Gompertz(data.frame(temperatures=25, age.mean=0:floor(max(subset(dataSCL_ec, subset = temperatures==25, select=c("temperatures", "age.mean"))))), par=x)))
outsp30 <- t(apply(par, MARGIN = 1, FUN=function(x) Gompertz(data.frame(temperatures=30, age.mean=0:floor(max(subset(dataSCL_ec, subset = temperatures==30, select=c("temperatures", "age.mean"))))), par=x)))


rangqtiles25 <- apply(outsp25, MARGIN=2, function(x) {quantile(x, probs=c(0.025, 0.5, 0.975))})
rangqtiles30 <- apply(outsp30, 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=dta25$age.mean, y=dta25$SCL.mean, 
            errbar.y = 2*dta25$SCL.SD, errbar.x = 2*dta25$age.SD, 
            bty="n", las=1, 
            ylim=c(0, 30), xlab="Days", ylab="SCL mm", 
            xlim=c(0, 80))
par(new=TRUE)
plot_errbar(x=dta30$age.mean, y=dta30$SCL.mean, 
            errbar.y = 2*dta30$SCL.SD, errbar.x = 2*dta30$age.SD, 
            bty="n", las=1, 
            ylim=c(0, 30), xlab="", ylab="", 
            xlim=c(0, 80), axes=FALSE)


library("plotrix")

draw.ellipse(x=dta25$age.mean, y=dta25$SCL.mean, 
             a=2*dta25$age.SD, 
             b=2*dta25$SCL.SD, border = NA, 
             col=rgb(0.5, 0.5, 0.5, 0.5))

draw.ellipse(x=dta30$age.mean, y=dta30$SCL.mean, 
             a=2*dta30$age.SD, 
             b=2*dta30$SCL.SD, border = NA, 
             col=rgb(0.5, 0.5, 0.5, 0.5))


lines((0:(dim(rangqtiles25)[2]-1)), rangqtiles25["2.5%", ], lty=2, lwd=2)
lines((0:(dim(rangqtiles25)[2]-1)), rangqtiles25["97.5%", ], lty=2, lwd=2)
lines((0:(dim(rangqtiles25)[2]-1)), rangqtiles25["50%", ], lty=1, lwd=2)

lines((0:(dim(rangqtiles30)[2]-1)), rangqtiles30["2.5%", ], lty=2, lwd=2)
lines((0:(dim(rangqtiles30)[2]-1)), rangqtiles30["97.5%", ], lty=2, lwd=2)
lines((0:(dim(rangqtiles30)[2]-1)), rangqtiles30["50%", ], lty=1, 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(dataSCL_ec$SCL.mean, Gompertz(dataSCL_ec, 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"]/max(dataSCL_ec$SCL.mean)), 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"]/max(dataSCL_ec$SCL.mean)), digits = 4)))))

title(expression(italic("Emys orbicularis")))