diff --git a/exploration/NEA_calibration_offline/1-Calibration_Thermal_Curves.Rmd b/exploration/NEA_calibration_offline/1-Calibration_Thermal_Curves.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..9ca6ee700e0d0bbf939a322b86ccea45a11fba70
--- /dev/null
+++ b/exploration/NEA_calibration_offline/1-Calibration_Thermal_Curves.Rmd
@@ -0,0 +1,693 @@
+---
+title: "Calibration_Thermal_Response"
+author: "Camille POULET"
+date: "22/02/2021"
+output: word_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = FALSE, include = TRUE, 
+                      warning = FALSE,
+                      fig.cap = TRUE)
+
+rm(list = ls())
+```
+
+```{r library, include=FALSE}
+library(dplyr)
+library(readr)
+library(ggplot2)
+library(gam)
+library(mgcv)
+library(XML)
+library(cowplot)
+library(gridExtra)
+library(voxel)
+library(officedown)
+library(flextable)
+library(tidymv)
+
+#library(lmerTest)
+
+```
+
+```{r upload data, warning=FALSE}
+
+source("../GR3D_Rdescription/GR3Dfunction.R")
+source("../GR3D_Rdescription/GR3D_NEA_env_data.R")
+source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
+
+
+```
+
+After reproduce, spawners either died (semelparous) or survived reproduction (primo-spawners).
+
+#Survival pre-reproduction 
+
+It concerns only spawners during reproduction, and this phase has been shown highly dependent on water (Walburg and Nichols, 1967).
+For the American shad (Alosa sapidissima), spawning occurs between 8°C to 26°C, with most of individuals reproduced between 14-21°C (Limburg et al, 2003).
+Leggett and Whitney (1972) found a median temperature for spawning of 16-20°C for York, St Johns and Connecticut Rivers.
+
+## Thermal range using GAMs
+
+To adjust the temperature thresholds and enhance model predictive capacity at the northern edge of species distribution, parameter values for spawners survival before and after spawning were obtained using generalized additive models (GAM) (Hastie and Tibshirani 1990). Observations of presence/absence across the 1900-1950 period (pristine distribution) are predicted regarding the average temperature times series observed across the 98 watersheds composing the environment of GR3D-U.S.
+Either summer or spring temperature are used as environmental variables depending the season at which each process occurs. 
+
+In GR3D, spawning occurs in river at the beginning of spring, so the spawner survival during reproduction was linked to spring temperature time-series in river. 
+
+```{r build a df for presence/absence over the 1900-1950 period }
+
+#Join the preence/asbsence to river temperature 
+nea_presence_temp<- inner_join(tempInriver %>% 
+                                 filter(between(year,  1900, 1950)) %>% 
+                                 group_by(basin_id, basin_name) %>% 
+                                 mutate(spring_avg = mean(spring_river_temperature),
+                                        summer_avg = mean(summer_river_temperature),
+                                        fall_avg = mean(fall_river_temperature),
+                                        winter_avg = mean(winter_river_temperature)), 
+                               nea_presence,
+                               by = c("basin_id","basin_name"))
+
+```
+
+
+```{r}
+rectTemperature = function(Tmin, Topt, Tmax, pct = 0.8, col ='red'){
+  optRange = thermalRange(pct, Tmin , Topt, Tmax)
+  rect = data.frame(xmin = c(-Inf, Tmin, optRange["upper"], Tmax),
+                    xmax = c(Tmin, optRange["lower"], Tmax, Inf),
+                    ymin = -Inf,
+                    ymax = Inf, 
+                    alpha = c(.3,.6,.6,.3),
+                    col = c(col, 'grey', 'grey', col))
+  return(rect)
+}
+
+```
+
+
+```{r Observations of presence/absence according to temperature, include = TRUE, fig.cap= 'Temperature experienced by spawners during their river stage in spring'' }
+
+#plot of temperature ranges
+obs_spg <- nea_presence_temp %>% 
+  ggplot() +
+  geom_point(aes(spring_river_temperature,obs_1900_1950), col = "orange", alpha = 0.3) +
+  labs(y = "Observed presences", x = "Spring river temperature (°C)") +
+  theme_bw()
+
+#temperature range experienced by fish in river in spring in which presence or absence are recorded 
+nea_presence_temp%>% 
+  group_by(obs_1900_1950) %>% 
+  summarise(Tmin = min(spring_river_temperature), Tmax = max(spring_river_temperature), .groups = 'drop') %>% 
+  flextable() %>% 
+  autofit()
+
+```
+
+```{r Gam on presences absences according to spring river temperature, include = TRUE  }
+
+#gam  with mgcv
+gam_spg <-  gam(obs_1900_1950 ~ s(spring_river_temperature),
+                data = nea_presence_temp, 
+                family = binomial("logit"))
+#Outputs
+#statistics and p-values
+summary(gam_spg) 
+#check
+par(mfrow = c(2,2))
+gam.check(gam_spg)
+
+#-------- with mgcv ----------#
+#predictions 
+#temp_spg <- seq(0, max(nea_presence_temp$spring_river_temperature), length.out = 100)
+#log distribution
+#y_pred_spg <- predict.gam(gam_spg, data.frame(spring_river_temperature = temp_spg)) 
+#pred_spg <- predict_gam(gam_spg)
+#natural distribution
+#y_pred_natural_spg <- 1/(1 + exp(-y_pred_spg))
+
+#-------- with tidyr ----------#
+pred_spg <- predict_gam(gam_spg) %>% #predict_gam extract the predicted values from the gam buid with mgcv
+  mutate(natural_fit = 1/(1+exp(-fit))) %>% #tranform into natural scales
+  mutate(natural_fit_80_pct = natural_fit * 0.80) #80% of survival  
+
+tempRangeOpt <-pred_spg %>% 
+  mutate(natural_fit_dec = as.numeric(sprintf("%.3f", natural_fit))) %>% 
+  filter(natural_fit_dec >= 1) 
+
+```
+
+
+```{r Functions for survival computation, include = TRUE  }
+
+S = 100  #number of spawners arriving at spawning grounds 
+temperature = seq(8,30,0.1) #temperature time series 
+#add the actual thermal range 
+rectTempRange = rectTemperature(Tmin = survivePar$TminSurv, 
+                                Topt = survivePar$ToptSurv,
+                                Tmax = survivePar$TmaxSurv, col = 'blue')
+
+
+#------- Compute spawners survival in river before reproduction --------#
+SpawnerSurvivalPreRep <- spawnerSurvivalPreReproduction(nea_presence_temp$spring_river_temperature, survivePar$TminSurv,survivePar$ToptSurv, survivePar$TmaxSurv)
+
+#------- Plot results from GAM with ggplot------------------------------#
+#with ggplot
+gam_gp_spg <- plotGAM(gam_spg, smooth.cov = "spring_river_temperature")+
+  geom_rug(data = nea_presence_temp,
+           aes_string(y = "obs_1900_1950", x = "spring_river_temperature"), alpha = 0.2)+
+  theme_bw()
+
+#natural scales
+natural_pred <-pred_spg %>% 
+  ggplot(aes(spring_river_temperature))+
+  geom_line(aes(y = natural_fit))+
+  geom_line(aes(y = natural_fit_80_pct), col = 'red', alpha = 0.3)+
+  theme_bw()+
+  annotate("rect", xmin = rectTempRange$xmin, xmax = rectTempRange$xmax, ymin = rectTempRange$ymin, ymax = rectTempRange$ymax, alpha = rectTempRange$alpha, fill = rectTempRange$col) 
+
+
+#plot all plots in the same window with plot_grid 
+plot_grid(obs_spg, gam_gp_spg, natural_pred, labels=c("A","B","C"), ncol = 3, nrow = 1)
+
+
+
+# plot(y_pred_natural_spg ~ temp_spg, lwd =2, type="l", 
+#  xlab = "Spring river temperature (°C)", 
+# ylab = "Predicted Presences/Absences", 
+#  xlim = c(0, 30))
+# par(new = TRUE)
+# plot(SpawnerSurvivalPreRep ~ nea_presence_temp$spring_river_temperature, 
+# lty=2, 
+# lwd=2,
+# col = "red",
+# xlab =" ",
+# ylab = " ", 
+# xlim = c(0,30), 
+#  xaxt="n", 
+# yaxt = "n")
+
+```
+
+```{r optimization, include = TRUE }
+
+#Found Tmin and Tmax to have a survival of 80% accross these ranges
+
+#--------------------------------------#
+#--Initialisation of parameters values--#
+#---------------------------------------#
+
+par_0 =  c(tempMinSurv = survivePar$TminSurv,
+           tempOptSurv = survivePar$ToptSurv,
+           tempMaxSurv = survivePar$TmaxSurv)
+
+
+par_1 =  c(tempMinSurv = 4.15,
+           tempOptSurv = survivePar$ToptSurv,
+           tempMaxSurv = 25.25)
+
+
+par_ini = par_0[c('tempMinSurv', 'tempMaxSurv')]
+par_ini_1 = par_1[c('tempMinSurv', 'tempMaxSurv')]
+par_ini_2 = par_1[c(('tempMinSurv'))]
+par_ini_3 = par_0[c(('tempOptSurv'))]
+
+fixedPar =  c(tempOptSurv = 17)
+fixedPar_2 = c(tempOptSurv = 17, tempMaxSurv = 26)
+fixedPar_3 = c(tempMinSurv = 4.15, tempMaxSurv = 25.25)
+fixedPar_4 = c(tempMinSurv = 8, tempMaxSurv = 26)
+
+#--------------------------------------#
+#----------- Optimisation --------------#
+#---------------------------------------#
+
+spSurvival_df = nea_presence_temp 
+
+spSurvival_df %>% 
+  group_by(obs_1900_1950) %>% 
+  summarise(Tmin = min(spring_river_temperature),
+            Tmax = max(spring_river_temperature), .groups='drop') %>% 
+  #filter(obs_1900_1950 == 1) %>% 
+  
+  objectiveFunction = function(par, data, fixedPar){
+    
+    df <- data %>% 
+      group_by(basin_id, basin_name,year) %>% 
+      mutate(spSurvival = spawnerSurvivalPreReproduction(spring_river_temperature,
+                                                         Tmin = fixedPar['tempMinSurv'],
+                                                         Topt = par ['tempOptSurv'],
+                                                         Tmax =fixedPar['tempMaxSurv'])) 
+    tempRange <- df %>% 
+      group_by(obs_1900_1950) %>% 
+      summarise(Tmin = min(spring_river_temperature),
+                Tmax = max(spring_river_temperature), .groups='drop') %>% 
+      filter(obs_1900_1950 == 1)
+    
+    
+    df<-df %>% 
+      filter(obs_1900_1950 == 1) %>% 
+      group_by(basin_id,basin_name) %>% 
+      mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.80, fifelse (spring_river_temperature == tempRange$Tmax, 0.80,
+                                                                                                 spSurvival))) %>% 
+      select(basin_id,basin_name,spring_river_temperature,spSurvival,spSurvival_opt,obs_1900_1950) %>% 
+      arrange(spring_river_temperature,spSurvival)
+    #mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.8 * spSurvival, fifelse (spring_river_temperature == tempRange$Tmax, 0.8 * spSurvival,
+    #spSurvival)))
+    
+    
+    SSE = df %>% 
+      group_by(basin_id,basin_name) %>% 
+      filter(spSurvival != spSurvival_opt) %>% 
+      mutate(SSE = (spSurvival - spSurvival_opt)^2) %>% 
+      ungroup() %>% 
+      summarise(sum(SSE), .groups = 'drop') %>% 
+      unlist(use.names = FALSE)
+    
+    return(SSE)
+  }
+
+#objectiveFunction(par_0,spSurvival_df)
+res <- optim(par = par_1, 
+             fn = objectiveFunction, 
+             data = spSurvival_df)
+#method = "L-BFGS-B",
+#lower = c(-10,14,26),
+#upper = c(10,21,30))
+
+
+res_3 <- optim(par = par_ini_3, 
+               fn = objectiveFunction, 
+               data = spSurvival_df,
+               fixedPar =fixedPar_3,
+               #method = "L-BFGS-B",
+               #lower = c(8),
+               upper = c(21))
+
+res_4 <- optim(par = par_ini_3, 
+               fn = objectiveFunction, 
+               data = spSurvival_df,
+               fixedPar =fixedPar_4)
+#method = "L-BFGS-B",
+#lower = c(8),
+
+#--------------------------------------#
+#----------- Plots curves --------------#
+#---------------------------------------#
+
+spSurvival_df <- spSurvival_df %>%  
+  group_by(basin_id, basin_name,year) %>% 
+  mutate(spSurvivalOpt = spawnerSurvivalPreReproduction(spring_river_temperature,
+                                                        Tmin = res$par['tempMinSurv'],
+                                                        Topt = res$par['tempOptSurv'],
+                                                        Tmax = res$par['tempMaxSurv']),
+         spSurvival = spawnerSurvivalPreReproduction(spring_river_temperature,
+                                                     Tmin = survivePar$TminSurv,
+                                                     Topt = survivePar$ToptSurv,
+                                                     Tmax = survivePar$TmaxSurv)) %>% 
+  select(basin_id,basin_name,year,spring_river_temperature,obs_1900_1950,spSurvival, spSurvivalOpt)
+
+
+spSurvival_df %>% 
+  ggplot(aes(spring_river_temperature)) +
+  geom_line(aes(y= spSurvival), col = "orange", alpha = 0.3) +
+  geom_line(aes(y= spSurvivalOpt), col = "red", alpha = 0.3)+
+  theme_bw()
+
+#------------------------
+
+
+```
+
+
+#Juvenile survival during the reproduction 
+
+After hatching, juveniles spend their first summer in freshwater before emigrating seaward, in the inshore basin to spend fall. 
+Considering that spawner reproduce in spring and a larvae survival 14dph, we used the spring riverine temperature for thermal range optimisation. 
+
+
+TO COMPLETE BY THERMAL RANGE 
+
+```{r optimisation with gam }
+
+S = 100  #number of spawners arriving at spawning grounds 
+temperature = seq(8,30,0.1) #temperature time series 
+#add the actual thermal range 
+rectTempRangeLarv = rectTemperature(Tmin = reproducePar$TminRep, 
+                                    Topt = reproducePar$ToptRep,
+                                    Tmax = reproducePar$TmaxRep, col = 'blue')
+
+#natural scales
+natural_pred_larv <-pred_spg %>% 
+  ggplot(aes(spring_river_temperature))+
+  geom_line(aes(y = natural_fit))+
+  geom_line(aes(y = natural_fit_80_pct), col = 'red', alpha = 0.3)+
+  theme_bw()+
+  annotate("rect", xmin = rectTempRangeLarv$xmin, xmax = rectTempRangeLarv$xmax, ymin = rectTempRangeLarv$ymin, ymax = rectTempRangeLarv$ymax, alpha = rectTempRangeLarv$alpha, fill = rectTempRangeLarv$col) 
+
+
+#plot all plots in the same window with plot_grid 
+plot_grid(obs_spg, gam_gp_spg, natural_pred_larv, labels=c("A","B","C"), ncol = 3, nrow = 1)
+
+
+
+
+```
+
+```{r optimization, include = TRUE }
+
+#Found Tmin and Tmax to have a survival of 80% accross these ranges
+
+#Juvenile and recruits survival
+RecruitsSurvival <- stockRecruitementRelationship(Triver = nea_presence_temp$spring_river_temperature,
+                                                  Tmin = reproducePar$TminRep,
+                                                  Topt = reproducePar$ToptRep,
+                                                  Tmax = reproducePar$TmaxRep,
+                                                  survivalStock = spawnerSurvivalPreReproduction(nea_presence_temp$spring_river_temperature, survivePar$TminSurv,survivePar$ToptSurv, survivePar$TmaxSurv))
+
+#--------------------------------------#
+#--Initialisation of parameters values--#
+#---------------------------------------#
+
+par_0 =  c(tempMinRep = reproducePar$TminRep,
+           tempOptRep = reproducePar$ToptRep,
+           tempMaxRep = reproducePar$TmaxRep)
+
+
+par_1 =  c(tempMinRep = 4.15,
+           tempOptRep = survivePar$ToptSurv,
+           tempMaxRep = 25.25)
+
+
+par_ini = par_0[c('tempMinSurv', 'tempMaxSurv')]
+par_ini_1 = par_1[c('tempMinSurv', 'tempMaxSurv')]
+par_ini_2 = par_1[c(('tempMinSurv'))]
+par_ini_3 = par_0[c(('tempOptSurv'))]
+
+fixedPar =  c(tempOptSurv = 17)
+fixedPar_2 = c(tempOptSurv = 17, tempMaxSurv = 26)
+fixedPar_3 = c(tempMinSurv = 4.15, tempMaxSurv = 25.25)
+fixedPar_4 = c(tempMinSurv = 8, tempMaxSurv = 26)
+
+#--------------------------------------#
+#----------- Optimisation --------------#
+#--------------------------------------#
+
+spSurvival_df = nea_presence_temp 
+
+spSurvival_df %>% 
+  group_by(obs_1900_1950) %>% 
+  summarise(Tmin = min(spring_river_temperature),
+            Tmax = max(spring_river_temperature), .groups='drop') %>% 
+  #filter(obs_1900_1950 == 1) %>% 
+  
+objectiveFunctionLarv = function(par, data, fixedPar){
+    
+    df <- data %>% 
+      group_by(basin_id, basin_name,year) %>% 
+      mutate(spLarvSurvival = stockRecruitementRelationship(spring_river_temperature,
+                                                            Tmin = par['tempMinRep'],
+                                                            Topt = par ['tempOptRep'],
+                                                            Tmax = par['tempMaxRep'],
+                                                            survivalStock = spawnerSurvivalPreReproduction(spring_river_temperature, survivePar$TminSurv,survivePar$ToptSurv, survivePar$TmaxSurv))) 
+
+                              
+tempRange <- df %>% 
+  group_by(obs_1900_1950) %>% 
+  summarise(Tmin = min(spring_river_temperature),
+            Tmax = max(spring_river_temperature), .groups='drop') %>%
+  filter(obs_1900_1950 == 1)
+                              
+                              
+df<-df %>% 
+  filter(obs_1900_1950 == 1) %>% 
+  group_by(basin_id,basin_name) %>% 
+  mutate(spLarvSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.80, fifelse (spring_river_temperature == tempRange$Tmax, 0.80,
+                                                                                                                           spLarvSurvival))) %>% 
+                                select(basin_id,basin_name,spring_river_temperature,spLarvSurvival,spLarvSurvival_opt,obs_1900_1950) %>% 
+  arrange(spring_river_temperature,spLarvSurvival)
+#mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.8 * spSurvival, fifelse (spring_river_temperature == tempRange$Tmax, 0.8 * spSurvival,
+                              #spSurvival)))
+                              
+                              
+SSE = df %>% 
+  group_by(basin_id,basin_name) %>% 
+  filter(spLarvSurvival != spLarvSurvival_opt) %>% 
+  mutate(SSE = (spLarvSurvival - spLarvSurvival_opt)^2) %>% 
+  ungroup() %>% 
+  summarise(sum(SSE), .groups = 'drop') %>% 
+  unlist(use.names = FALSE)
+                              
+return(SSE)
+  }
+
+#objectiveFunction(par_0,spSurvival_df)
+res <- optim(par = par_1, 
+             fn = objectiveFunctionLarv, 
+             data = spSurvival_df)
+#method = "L-BFGS-B",
+#lower = c(-10,14,26),
+#upper = c(10,21,30))
+
+
+res_3 <- optim(par = par_ini_3, 
+               fn = objectiveFunction, 
+               data = spSurvival_df,
+               fixedPar =fixedPar_3,
+               #method = "L-BFGS-B",
+               #lower = c(8),
+               upper = c(21))
+
+res_4 <- optim(par = par_ini_3, 
+               fn = objectiveFunction, 
+               data = spSurvival_df,
+               fixedPar =fixedPar_4)
+#method = "L-BFGS-B",
+#lower = c(8),
+
+#--------------------------------------#
+#----------- Plots curves --------------#
+#---------------------------------------#
+
+spSurvival_df <- spSurvival_df %>%  
+  group_by(basin_id, basin_name,year) %>% 
+  mutate(spSurvivalOpt = spawnerSurvivalPreReproduction(spring_river_temperature,
+                                                        Tmin = res$par['tempMinSurv'],
+                                                        Topt = res$par['tempOptSurv'],
+                                                        Tmax = res$par['tempMaxSurv']),
+         spSurvival = spawnerSurvivalPreReproduction(spring_river_temperature,
+                                                     Tmin = survivePar$TminSurv,
+                                                     Topt = survivePar$ToptSurv,
+                                                     Tmax = survivePar$TmaxSurv)) %>% 
+  select(basin_id,basin_name,year,spring_river_temperature,obs_1900_1950,spSurvival, spSurvivalOpt)
+
+
+spSurvival_df %>% 
+  ggplot(aes(spring_river_temperature)) +
+  geom_line(aes(y= spSurvival), col = "orange", alpha = 0.3) +
+  geom_line(aes(y= spSurvivalOpt), col = "red", alpha = 0.3)+
+  theme_bw()
+
+#------------------------
+
+
+```
+
+
+
+
+
+
+```{r Observations of presence/absence according to temperature }
+
+#In spring
+obs_spg <- nea_presence_temp %>% 
+  ggplot() +
+  geom_point(aes(spring_river_temperature,obs_1900_1950), col = "orange", alpha = 0.3) +
+  labs(y = "Observed presences", x = "Spring river temperature (°C)") +
+  theme_classic()
+
+#In summer 
+obs_sum <- nea_presence_temp%>% 
+  ggplot() +
+  geom_point(aes(summer_river_temperature,obs_1900_1950), col = "red", alpha = 0.1) +
+  labs(y = "", x = "Summer river temperature (°C)") +
+  theme_classic()
+
+#plot_grid(spp, sup, labels = c("A", "B"), ncol = 2, nrow = 1)
+#grid.arrange(spp, sup, ncol = 2, nrow = 1)
+
+```
+
+```{r Gam on presences absences according to river temperature, include = TRUE  }
+
+#Spring 
+gam_spg <-  gam(obs_1900_1950 ~ s(spring_river_temperature),
+                data = nea_presence_temp, 
+                family = binomial("logit"))
+#Summer
+gam_sum <-  gam(obs_1900_1950 ~ s(summer_river_temperature),
+                data = nea_presence_temp, 
+                family = binomial("logit"))
+
+#Outputs
+#statistics and p-values
+summary(gam_spg) 
+summary(gam_sum)
+
+#check
+par(mfrow = c(2,2))
+gam.check(gam_sum)
+
+par(mfrow = c(2,2))
+gam.check(gam_spg)
+
+#predictions 
+temp_spg <- seq(0, max(nea_presence_temp$spring_river_temperature), length.out = 100)
+temp_sum <- seq(0, max(nea_presence_temp$summer_river_temperature), length.out = 100)
+
+#log distribution
+y_pred_spg <- predict.gam(gam_spg, data.frame(spring_river_temperature = temp_spg)) 
+y_pred_sum <- predict.gam(gam_sum, data.frame(summer_river_temperature = temp_sum)) 
+
+#natural distribution
+y_pred_natural_sum <- 1/(1 + exp(-y_pred_spg))
+y_pred_natural_spg <- 1/(1 + exp(-y_pred_sum))
+
+```
+
+```{r Functions for survival computation, include = TRUE  }
+
+S = 100  #number of spawners arriving at spawning grounds 
+temperature = seq(8,30,0.1) #temperature time series 
+
+#Spawners survival in river before reproduction 
+SpawnerSurvivalPreRep <- spawnerSurvivalPreReproduction(nea_presence_temp$spring_river_temperature, survivePar$TminSurv,survivePar$ToptSurv, survivePar$TmaxSurv)
+
+#Juvenile and recruits survival
+RecruitsSurvival <- stockRecruitementRelationship(Triver = nea_presence_temp$spring_river_temperature,
+                                                  Tmin = reproducePar$TminRep,
+                                                  Topt = reproducePar$ToptRep,
+                                                  Tmax = reproducePar$TmaxRep,
+                                                  survivalStock = S)
+
+#Spawners survival in river after reproduction
+TrefSummer <- nea_presence_temp %>% 
+  ungroup() %>% 
+  dplyr::summarise(mean(summer_river_temperature)) %>% 
+  unlist()
+
+
+SpawnerSurvivalPostRep <- spawnerSurvivalPostReproduction(nea_presence_temp$summer_river_temperature, Tref  = TrefSummer,
+                                                          coeffa = 1,
+                                                          coeffb = 0.1)
+
+```
+
+```{r plots results, include = TRUE}
+
+#----------Plot results from GAM----------#
+gam_gp_spg <- plotGAM(gam_spg, smooth.cov = "spring_river_temperature")+
+  geom_rug(data = nea_presence_temp,
+           aes_string(y = "obs_presence", x = "spring_river_temperature"), alpha = 0.2)
+
+#plot.gam(gam_spg, ylim = c(-100,100), xlab = "Spring river temperature (°C)")
+
+gam_gp_sum <- plotGAM(gam_sum, smooth.cov = "summer_river_temperature")+
+  geom_rug(data = nea_presence_temp,
+           aes_string(y = "obs_presence", x = "summer_river_temperature"), alpha = 0.2)
+
+#plot.gam(gam_sum, ylim = c(-100,100), xlab = "Summer river temperature (°C)")
+
+#----------Plot predictions from GAM ----------#
+
+plot(y_pred_natural_spg ~ temp_spg, lwd =2, type="l", 
+     xlab = "Spring river temperature (°C)", 
+     ylab = "Predicted Presences/Absences", 
+     xlim = c(0, 30))
+par(new = TRUE)
+plot(SpawnerSurvivalPreRep ~ nea_presence_temp$spring_river_temperature, 
+     lty=2, 
+     lwd=2,
+     col = "red",
+     xlab =" ",
+     ylab = " ", 
+     xlim = c(0,30), 
+     xaxt="n", 
+     yaxt = "n")
+par(new = TRUE)
+plot(RecruitsSurvival ~ nea_presence_temp$spring_river_temperature, 
+     lty = 2, 
+     lwd = 2,
+     col = "orange",
+     xlab =" ",
+     ylab = " ", 
+     xlim = c(0,30), 
+     xaxt="n", 
+     yaxt = "n")
+
+
+#For summer 
+plot(y_pred_natural_sum ~ temp_sum, lwd =2, type="l", 
+     xlab = "Spring river temperature (°C)", 
+     ylab = "Predicted Presences/Absences", 
+     xlim = c(0, 30))
+par(new = TRUE)
+plot(SpawnerSurvivalPostRep ~ nea_presence_temp$summer_river_temperature, 
+     lty=2, 
+     lwd=2,
+     col = "red",
+     xlab =" ",
+     ylab = " ", 
+     xlim = c(0,30), 
+     xaxt="n", 
+     yaxt = "n")
+```
+
+#Calibration using an optimization 
+
+
+```{r optimization, include = TRUE }
+
+FinalFunction <- function (temp, RiverSurvivalPreReproduction, StockRecruitementRelationship, RiverSurvivalPostReproduction_2){
+  
+  ###--------------------- River Survival -----------------
+  #Define the number of spawners which survived during the upstream migration and reach the spawning grounds 
+  SpRiver = RiverSurvivalPreReproduction(Tmin, Topt, Tmax, S)
+  StockAfterSurv = S * SpRiver
+  ###--------------------- SR relationship -----------------
+  #Define the number of recruits produced by the surival spawners 
+  Rj = StockRecruitementRelationship(Tmin, Topt, Tmax, StockAfterSurv)
+  
+  ###--------------------- River survival after reproduce -----------------
+  
+  SpRiverPostSpawn = RiverSurvivalPostReproduction_2(temp,TrefSim, coeffa ,coeffb)
+  
+  StockAfterSurvPostSpawn = StockAfterSurv * SpRiverPostSpawn
+  
+  ###--------------------- Sea Survival -----------------
+  #Define the number of juveniles which survived after 5 years at sea. 
+  
+  #Probabilité de survie d'un individu en mer 
+  SpSea = exp (-Zsea * PeriodAtSea) 
+  
+  StockAfterSurvSea = Rj * SpSea #Will begin new spawners 
+  
+  return(StockAfterSurvSea)
+  
+}
+
+ObjectiveFunction = function(temperature, RiverSurvivalPreReproduction,StockRecruitementRelationship, RiverSurvivalPostReproduction_2){
+  squareDistance <- (FinalFunction(temperature, RiverSurvivalPreReproduction, StockRecruitementRelationship, RiverSurvivalPostReproduction_2) - S)^2
+  return(squareDistance)
+}
+
+optimise(ObjectiveFunction, c(10,20), RiverSurvivalPreReproduction(8,17,26.,S), StockRecruitementRelationship(10,17,30,S))
+
+
+
+```
+
+
+
+
+