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)
Lepidochelys.olivacea <- subset(DatabaseTSD, Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0))

# Here the intersexes are included as males

#################

Lepidochelys.olivacea.WA <- subset(Lepidochelys.olivacea, subset=RMU.2023 == "West Atlantic")

males <- Lepidochelys.olivacea.WA$Males + ifelse(is.na(Lepidochelys.olivacea.WA$Intersexes), 0, Lepidochelys.olivacea.WA$Intersexes)
females <- Lepidochelys.olivacea.WA$Females
temperatures <- Lepidochelys.olivacea.WA$Incubation.temperature.set


tsdF.WA <- tsd(males = males                                                       , 
               females = females                                                   , 
               temperatures = temperatures                                         , 
               equation = "flexit**")

priors.WA = tsd_MHmcmc_p(result = tsdF.WA, 
                         default = c("dnorm", "dunif", "dunif"), 
                         accept = TRUE)

tsdF_mcmc.WA = tsd_MHmcmc(result = tsdF.WA, parametersMCMC = priors.WA, 
                          n.iter = 20000, thin = 10, n.adapt=1000)
plot(tsdF.WA, resultmcmc = tsdF_mcmc.WA, replicate.CI = 2000, xlim=c(24, 36))

tsdC.WA <- tsd(males = males                                                       , 
               females = females                                                   , 
               temperatures = temperatures                                         , 
               equation = "stairs"                                                 )
plot(tsdC.WA)

priors.WA = tsd_MHmcmc_p(result = tsdC.WA, 
                         default = c("dnorm", "dnorm", "dunif"), 
                         accept = TRUE)


tsdC_mcmc.WA = tsd_MHmcmc(result = tsdC.WA, parametersMCMC = priors.WA, 
                          n.iter = 20000, thin = 10, n.adapt=1000)

plot(tsdC_mcmc.WA, parameters = "TRTL", what="MarkovChain", main="TRTL")
plot(tsdC_mcmc.WA, parameters = "TRTH", what="MarkovChain", main="TRTH")
plot(tsdC_mcmc.WA, parameters = "SRTRT", what="MarkovChain", main="SRTRT")

as.parameters(tsdC_mcmc.WA, index = "quantile", probs = c(0.025, 0.5, 0.975))
as.parameters(tsdC_mcmc.WA, index = "quantile", probs = c(0.25, 0.5, 0.75))
plot(tsdC.WA, resultmcmc = tsdC_mcmc.WA, replicate.CI = 2000, xlim=c(24, 36))

library(loo)

tsdC_mcmc.WA.loo = loo(tsdC_mcmc.WA$WAIC)
tsdF_mcmc.WA.loo = loo(tsdF_mcmc.WA$WAIC)

loo_compare(list(Stairs.moldel=tsdC_mcmc.WA.loo, 
                 Fexit.model=tsdF_mcmc.WA.loo))

ic = c(Stairs=tsdC_mcmc.WA.loo$estimates["looic", "Estimate"], 
       Flexit=tsdF_mcmc.WA.loo$estimates["looic", "Estimate"])
exp(ic) / sum(exp(ic))

################
Lepidochelys.olivacea.EP <- subset(Lepidochelys.olivacea, subset=RMU.2023 == "East Pacific")

males <- Lepidochelys.olivacea.EP$Males + ifelse(is.na(Lepidochelys.olivacea.EP$Intersexes), 0, Lepidochelys.olivacea.EP$Intersexes)
females <- Lepidochelys.olivacea.EP$Females
temperatures <- Lepidochelys.olivacea.EP$Incubation.temperature.set
N <- males + females

tsdF.EP <- tsd(males = males                                                      , 
              females = females                                                   , 
              temperatures = temperatures                                         , 
              equation = "flexit**"                                               )

priors.EP = tsd_MHmcmc_p(result = tsdF.EP, 
                      default = c("dnorm", "dunif", "dunif"), 
                      accept = TRUE)

tsdF_mcmc.EP = tsd_MHmcmc(result = tsdF.EP, parametersMCMC = priors.EP, 
                          n.iter = 20000, thin = 10, n.adapt=1000)
plot(tsdF.EP, resultmcmc = tsdF_mcmc.EP, replicate.CI = 2000, xlim=c(24, 36))

TRTL <- min(temperatures[females != 0])
TRTL <- (max(temperatures[temperatures < TRTL]) + TRTL) / 2
TRTH <- max(temperatures[males != 0])
TRTH <- (min(temperatures[temperatures > TRTH]) + TRTH) / 2
SRTRT <- which((temperatures > TRTL) & (temperatures < TRTH))
SR <- logit(sum(((males[SRTRT]/N[SRTRT])*N[SRTRT]))/sum(N[SRTRT]))

srt <- SR

SR <- -0.7634089103614673
SR <- 0.102857385439708

-sum(dbinom(prob=rep(invlogit(SR), length(SRTRT)), 
           size = N[SRTRT], x=males[SRTRT], log=TRUE))
  
getFromNamespace(".modelTSD", ns="embryogrowth")(par=c(TRTL=TRTL, TRTH=TRTH, SRTRT=SR), 
                                                 temperatures=temperatures, 
                                                 equation=equation)
getFromNamespace(".tsd_fit", ns="embryogrowth")(par=c(TRTL=TRTL, TRTH=TRTH, SRTRT=SR), 
                                                fixed.parameters = NULL, 
                                                males=males, N=N, 
                                                temperatures=temperatures, 
                                                equation=equation)


tsdC.EP <- tsd(males = males                                                       , 
               females = females                                                   , 
               temperatures = temperatures                                         , 
               equation = "stairs"                                                 )
plot(tsdC.EP)

priors.EP <- tsd_MHmcmc_p(result = tsdC.EP, 
                                   default = c("dnorm", "dnorm", "dunif"), 
                                   accept = TRUE)


tsdC_mcmc.EP <- tsd_MHmcmc(result = tsdC.EP, parametersMCMC = priors.EP, 
                          n.iter = 20000, thin = 10, n.adapt=1000)

plot(tsdC_mcmc.EP, parameters = "TRTL", what="MarkovChain", main="TRTL")
plot(tsdC_mcmc.EP, parameters = "TRTH", what="MarkovChain", main="TRTH")
plot(tsdC_mcmc.EP, parameters = "SRTRT", what="MarkovChain", main="SRTRT")

as.parameters(tsdC_mcmc.EP, index = "quantile", probs = c(0.025, 0.5, 0.975))
as.parameters(tsdC_mcmc.EP, index = "quantile", probs = c(0.25, 0.5, 0.75))
as.parameters(tsdC_mcmc.EP)

plot(tsdC.EP, resultmcmc = tsdC_mcmc.EP, replicate.CI = 2000, xlim=c(24, 36))

library(loo)

tsdC_mcmc.EP.loo = loo(tsdC_mcmc.EP$WAIC)
tsdF_mcmc.EP.loo = loo(tsdF_mcmc.EP$WAIC)

loo_compare(list(Stairs.moldel=tsdC_mcmc.EP.loo, 
                 Fexit.model=tsdF_mcmc.EP.loo))

ic = c(Stairs=tsdC_mcmc.EP.loo$estimates["looic", "Estimate"], 
       Flexit=tsdF_mcmc.EP.loo$estimates["looic", "Estimate"])
exp(ic) / sum(exp(ic))


################
Lepidochelys.olivacea.NEI <- subset(Lepidochelys.olivacea, subset=RMU.2023 == "Northeast Indian")

males <- Lepidochelys.olivacea.NEI$Males + ifelse(is.na(Lepidochelys.olivacea.NEI$Intersexes), 0, Lepidochelys.olivacea.NEI$Intersexes)
females <- Lepidochelys.olivacea.NEI$Females
temperatures <- Lepidochelys.olivacea.NEI$Incubation.temperature.set

tsdF.NEI <- tsd(males = males                                                       , 
                females = females                                                   , 
                temperatures = temperatures                                         , 
                equation = "flexit**"                                               )

priors.NEI <- tsd_MHmcmc_p(result = tsdF.NEI, 
                           default = c("dnorm", "dunif", "dunif"), 
                           accept = TRUE)

tsdF_mcmc.NEI = tsd_MHmcmc(result = tsdF.NEI, parametersMCMC = priors.NEI, 
                           n.iter = 20000, thin = 10, n.adapt=1000)
plot(tsdF.NEI, resultmcmc = tsdF_mcmc.NEI, replicate.CI = 2000, xlim=c(24, 36))

tsdC.NEI <- tsd(males = males                                                       , 
                females = females                                                   , 
                temperatures = temperatures                                         , 
                equation = "stairs"                                                 )
plot(tsdC.NEI)

priors.NEI <- tsd_MHmcmc_p(result = tsdC.NEI, 
                           default = c("dnorm", "dnorm", "dunif"), 
                           accept = TRUE)

tsdC_mcmc.NEI = tsd_MHmcmc(result = tsdC.NEI, parametersMCMC = priors.NEI, 
                           n.iter = 20000, thin = 10, n.adapt=1000)

plot(tsdC_mcmc.NEI, parameters = "TRTL", what="MarkovChain", main="TRTL")
plot(tsdC_mcmc.NEI, parameters = "TRTH", what="MarkovChain", main="TRTH")
plot(tsdC_mcmc.NEI, parameters = "SRTRT", what="MarkovChain", main="SRTRT")

as.parameters(tsdC_mcmc.NEI, index = "quantile", probs = c(0.025, 0.5, 0.975))
as.parameters(tsdC_mcmc.NEI, index = "quantile", probs = c(0.25, 0.5, 0.75))
plot(tsdC.NEI, resultmcmc = tsdC_mcmc.NEI, replicate.CI = 2000, xlim=c(24, 36))

library(loo)

tsdC_mcmc.NEI.loo = loo(tsdC_mcmc.NEI$WAIC)
tsdF_mcmc.NEI.loo = loo(tsdF_mcmc.NEI$WAIC)

loo_compare(list(Stairs.moldel=tsdC_mcmc.NEI.loo, 
                 Fexit.model=tsdF_mcmc.NEI.loo))

ic = c(Stairs=tsdC_mcmc.NEI.loo$estimates["looic", "Estimate"], 
       Flexit=tsdF_mcmc.NEI.loo$estimates["looic", "Estimate"])
exp(ic) / sum(exp(ic))


######################

males <- Lepidochelys.olivacea$Males + ifelse(is.na(Lepidochelys.olivacea$Intersexes), 0, Lepidochelys.olivacea$Intersexes)
females <- Lepidochelys.olivacea$Females
temperatures <- Lepidochelys.olivacea$Incubation.temperature.set

tsdF <- tsd(males = males                                                       , 
            females = females                                                   , 
            temperatures = temperatures                                         , 
            equation = "flexit**"                                               )

priors <- tsd_MHmcmc_p(result = tsdF, default = c("dnorm", "dunif", "dunif"), accept = TRUE)


tsdF_mcmc = tsd_MHmcmc(result = tsdF, parametersMCMC = priors, 
                       n.iter = 20000, thin = 10, n.adapt=1000)
plot(tsdF, resultmcmc = tsdF_mcmc, replicate.CI = 2000, xlim=c(24, 36))

tsdC <- tsd(males = males                                                       , 
            females = females                                                   , 
            temperatures = temperatures                                         , 
            equation = "stairs"                                                 )
plot(tsdC)

priors <-tsd_MHmcmc_p(result = tsdC, 
                                     default = c("dnorm", "dnorm", "dunif"), 
                                     accept = TRUE)


tsdC_mcmc = tsd_MHmcmc(result = tsdC, parametersMCMC = priors, 
                       n.iter = 20000, thin = 10, n.adapt=1000)

plot(tsdC_mcmc, parameters = "TRTL", what="MarkovChain", main="TRTL")
plot(tsdC_mcmc, parameters = "TRTH", what="MarkovChain", main="TRTH")
plot(tsdC_mcmc, parameters = "SRTRT", what="MarkovChain", main="SRTRT")

as.parameters(tsdC_mcmc, index = "quantile", probs = c(0.025, 0.5, 0.975))
as.parameters(tsdC_mcmc, index = "quantile", probs = c(0.25, 0.5, 0.75))
plot(tsdC, resultmcmc = tsdC_mcmc, replicate.CI = 2000, xlim=c(24, 36))

library(loo)

tsdC_mcmc.loo = loo(tsdC_mcmc$WAIC)
tsdF_mcmc.loo = loo(tsdF_mcmc$WAIC)

loo_compare(list(Stairs.moldel=tsdC_mcmc.loo, 
                 Fexit.model=tsdF_mcmc.loo))

ic = c(Stairs=tsdC_mcmc.loo$estimates["looic", "Estimate"], 
       Flexit=tsdF_mcmc.loo$estimates["looic", "Estimate"])
exp(ic) / sum(exp(ic))

# Note here that