From bcc08e41429fdaef0dffe76bce9cc47012ab4730 Mon Sep 17 00:00:00 2001
From: "patrick.lambert" <patrick.mh.lambert@inrae.fr>
Date: Tue, 20 Apr 2021 23:05:59 +0200
Subject: [PATCH] new calibration of grow

---
 exploration/GR3D_Rdescription/GR3Dfunction.R  | 108 ++-
 .../NEA_calibration_offline/maturationPL.Rmd  | 628 +++++++++++++++---
 2 files changed, 659 insertions(+), 77 deletions(-)

diff --git a/exploration/GR3D_Rdescription/GR3Dfunction.R b/exploration/GR3D_Rdescription/GR3Dfunction.R
index 9707741..2903bd6 100644
--- a/exploration/GR3D_Rdescription/GR3Dfunction.R
+++ b/exploration/GR3D_Rdescription/GR3Dfunction.R
@@ -54,7 +54,7 @@ vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma,
   return(L)
 }
 
-vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sigma, tempEffect ){
+vonBertalanffyWithNextIncrement = function(L,  Linf, K, timeStepDuration, sigma, tempEffect ){
   if (sigma == 0) {
     mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
     increment = exp(mu)
@@ -67,6 +67,112 @@ vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sig
   return(L)
 }
 
+# generate a cohort of length trajectories
+computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar){
+  ages = temperaturePattern$age
+  
+  res <- expand_grid(ind = seq.int(Nind),  age = ages) %>% 
+    inner_join(temperaturePattern, by = 'age') %>% 
+    arrange(age, ind) %>% 
+    mutate(temperatureEffect = temperatureEffect(temperature, 
+                                                 growPar$tempMinGrow, 
+                                                 growPar$tempOptGrow, 
+                                                 growPar$tempMaxGrow), 
+           L = if_else(age == 0, growPar$lengthAtHatching, 0))
+  
+  for (i in 2:length(ages)) {
+    previousAge  = ages[i - 1]
+    currentAge = ages[i]
+    tempEffect = res %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE)
+    previousL <- res %>% filter(age == previousAge) %>% select(L) %>% unlist(use.names = FALSE)
+    
+    currentL <- vonBertalanffyWithNextIncrement(L = previousL, 
+                                                Linf = growPar$Linf, 
+                                                K = growPar$kOpt,
+                                                timeStepDuration = currentAge - previousAge, 
+                                                sigma = growPar$sigmaDeltaLVonBert, 
+                                                tempEffect = tempEffect)
+    res = res %>% mutate(L =replace(L, age == ages[i], currentL) ) 
+  }
+  return(res)
+}
+
+computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern, 
+                                                           Nind = 10, 
+                                                           growPar, 
+                                                           RNGseed =1){
+  set.seed(RNGseed)
+  ages = temperaturePattern %>% 
+    distinct(age) %>% 
+    unlist(use.names = FALSE)
+  
+  res <- expand_grid(ind = seq.int(Nind),  age = ages) %>% 
+    mutate(random = rnorm(Nind * length(ages)))  %>% 
+    inner_join(temperaturePattern, by = 'age') %>% 
+    arrange(age, ind) %>% 
+    mutate(temperatureEffect = temperatureEffect(temperature, 
+                                                 growPar$tempMinGrow, 
+                                                 growPar$tempOptGrow, 
+                                                 growPar$tempMaxGrow), 
+           L = if_else(age == 0, growPar$lengthAtHatching, 0))
+  
+  for (i in 2:length(ages)) {
+    previousAge  = ages[i - 1]
+    currentAge = ages[i]
+    tempEffect <-  res %>% filter(age == currentAge) %>% 
+      select(temperatureEffect) %>% unlist(use.names = FALSE)
+    
+    previousL <- res %>% filter(age == previousAge) %>% 
+      select(L) %>% unlist(use.names = FALSE)
+    rnd <- res %>% filter(age == currentAge) %>% select(random) %>% unlist(use.names = FALSE)
+    currentL <- vonBertalanffyWithRandomVector(L = previousL, 
+                                               Linf = growPar$Linf, 
+                                               K = growPar$kOpt,
+                                               timeStepDuration = currentAge - previousAge, 
+                                               randomVector = rnd,
+                                               sigma = growPar$sigmaDeltaLVonBert, 
+                                               tempEffect = tempEffect)
+    res = res %>% mutate(L = replace(L, age == ages[i], currentL) ) 
+  }
+  return(res)
+}
+
+#computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 10, 
+#growPar = growParUnisex,
+# RNGseed = 1)
+
+vonBertalanffyWithRandomVector = function(L, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
+  if (sigma == 0) {
+    mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))),-Inf)
+    #mu =  log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration)))
+    increment = exp(mu)
+  }
+  else {
+    mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2, -Inf)
+    #      mu = log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2
+    increment = exp(randomVector * sigma + mu)
+  }
+  L = pmin(Linf, L + increment) 
+  return(L)
+}
+
+
+computeOgive = function(lengthTrajectories, lengthAtMaturity) {
+  ogive <-
+    lengthTrajectories %>% 
+    group_by(age) %>% 
+    summarise(nTotal=n()) %>% 
+    left_join(lengthTrajectories %>%  group_by(age) %>% 
+                filter(L > lengthAtMaturity) %>% 
+                summarise(mature = n(), .groups = 'drop'), 
+              by = 'age') %>% 
+    replace_na(list(mature = 0)) %>% 
+    mutate(p =  mature/nTotal) %>% 
+    select(age, p)
+  
+  return(ogive)
+} 
+
 # -------------------------------------------------------
 # Dispersal      
 # see see (Chapman et al., 2007; Holloway et al., 2016)
diff --git a/exploration/NEA_calibration_offline/maturationPL.Rmd b/exploration/NEA_calibration_offline/maturationPL.Rmd
index 87ca4b7..9253bb8 100644
--- a/exploration/NEA_calibration_offline/maturationPL.Rmd
+++ b/exploration/NEA_calibration_offline/maturationPL.Rmd
@@ -8,55 +8,61 @@ output:
       Normal: ['First Paragraph']
 ---
 
-```{r setup, include=FALSE}
+```{r setup - library , include=FALSE}
 knitr::opts_chunk$set(echo = TRUE, fig.cap = TRUE)
+
 library(officedown)
 library(officer)
-
 library(tidyverse)
 library(tictoc) # see timeR
-
+library(data.table)
+library(reader)
 
 
 ```
 
-```{r}
+```{r data}
 rm(list = ls())
 source("../GR3D_Rdescription/GR3Dfunction.R")
+source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
 
+#load("growth_optimal.RData")
 load('SOS.rdata')
+#load('res_optim.rdata')
+
+resA_par = read_csv2("res_A.csv", skip =1, col_names=F) %>% 
+  deframe()
+
+
 ```
 
-```{r}
-source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
+```{r load growth parameters and tempertaure pattern }
 
-growParXML = data.frame(t(unlist(fishXML[[1]][["processes"]][["processesEachStep"]][["species.Grow"]]))) %>% 
-  select(-synchronisationMode ) %>% 
-  mutate(linfVonBertForMale = fishXML[[1]][["linfVonBertForMale"]], 
-         linfVonBertForFemale = fishXML[[1]][["linfVonBertForFemale"]],
-         lengthAtHatching = fishXML[[1]][["lengthAtHatching"]],
-         lFirstMaturityForMale = fishXML[[1]][["lFirstMaturityForMale"]],
-         lFirstMaturityForFemale = fishXML[[1]][["lFirstMaturityForFemale"]]) %>% 
-  mutate_all(~as.numeric(as.character(.)))
+regional_metapop = c("semelparous", "southern iteroparous", "northern iteroparous")
+growParXML = growPar
 
+#Growth parameters from Stich et al, 2020
+#Temperature relative values from growth optimization 
 
-growParStich <- Stich2020_sel %>% filter(catchment == "semelparous") %>% 
-  select(Linf, K,  L0_theo) %>% 
+growParStich <- Stich2020_sel %>% 
+  filter(catchment %in% regional_metapop) %>% 
+  select(catchment, Linf, K,  L0_theo) %>% 
+  rename ("metapop" = "catchment") %>% 
   mutate(tempMinGrow = 1.6, tempOptGrow = 5, tempMaxGrow = 27.9, 
-         Linf =Linf/10, lengthAtHatching= L0_theo/10, kOpt = K,
+         Linf = Linf/10, lengthAtHatching= L0_theo/10, kOpt = K,
          sigmaDeltaLVonBert = 0.4)
 
 temperaturePattern <- growthInBasin %>% 
   group_by(metapop, age, season) %>% 
   summarise(temperature = mean(temperature_RIO), .groups = 'drop') %>% 
-  filter(metapop == "semelparous") %>% 
+  #filter(metapop == "semelparous") %>% 
   # mutate(season = case_when(
   #   (age - floor(age)) == 0.00 ~ 'spring',
   #   (age - floor(age)) == 0.25 ~ 'summer',
   #   (age - floor(age)) == 0.50 ~ 'fall',
   #   (age - floor(age)) == 0.75 ~ 'winter'
   # )) %>% 
-  select(age, season,temperature)
+  select(metapop,age, season,temperature)
 
 ```
 
@@ -82,7 +88,6 @@ computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, grow
     previousL <- res %>% filter(age == previousAge) %>% select(L) %>% unlist(use.names = FALSE)
     
     currentL <- vonBertalanffyWithNextIncrement(L = previousL, 
-                                                L0 = growPar$lengthAtHatching, 
                                                 Linf = growPar$Linf, 
                                                 K = growPar$kOpt,
                                                 timeStepDuration = currentAge - previousAge, 
@@ -98,7 +103,9 @@ computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
                                                            growPar, 
                                                            RNGseed =1){
   set.seed(RNGseed)
-  ages = temperaturePattern$age
+  ages = temperaturePattern %>% 
+    distinct(age) %>% 
+    unlist(use.names = FALSE)
   
   res <- expand_grid(ind = seq.int(Nind),  age = ages) %>% 
     mutate(random = rnorm(Nind * length(ages)))  %>% 
@@ -120,7 +127,6 @@ computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
       select(L) %>% unlist(use.names = FALSE)
     rnd <- res %>% filter(age == currentAge) %>% select(random) %>% unlist(use.names = FALSE)
     currentL <- vonBertalanffyWithRandomVector(L = previousL, 
-                                               L0 = growPar$lengthAtHatching, 
                                                Linf = growPar$Linf, 
                                                K = growPar$kOpt,
                                                timeStepDuration = currentAge - previousAge, 
@@ -132,22 +138,44 @@ computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
   return(res)
 }
 
+#computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 10, 
+#growPar = growParUnisex,
+# RNGseed = 1)
 
-vonBertalanffyWithRandomVector = function(L, L0, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
+vonBertalanffyWithRandomVector = function(L, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
   if (sigma == 0) {
-    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
+    mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))),-Inf)
+    #mu =  log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration)))
     increment = exp(mu)
   }
   else {
-    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
+  mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2, -Inf)
+  #      mu = log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2
     increment = exp(randomVector * sigma + mu)
   }
-  L = L + increment * tempEffect
+  L = pmin(Linf, L + increment) 
   return(L)
 }
 
 
 computeOgive = function(lengthTrajectories, lengthAtMaturity) {
+  ogive <-
+    lengthTrajectories %>% 
+    group_by(age) %>% 
+    summarise(nTotal=n()) %>% 
+    left_join(lengthTrajectories %>%  group_by(age) %>% 
+                filter(L > lengthAtMaturity) %>% 
+                summarise(mature = n(), .groups = 'drop'), 
+              by = 'age') %>% 
+    replace_na(list(mature = 0)) %>% 
+    mutate(p =  mature/nTotal) %>% 
+    select(age, p)
+  
+  return(ogive)
+} 
+
+#Deprecated
+computeOgive2 = function(lengthTrajectories, lengthAtMaturity) {
   ogive <- lengthTrajectories %>% 
     distinct(age) %>% arrange(age) %>% 
     left_join(lengthTrajectories %>% filter(L > lengthAtMaturity) %>% 
@@ -164,10 +192,10 @@ computeOgive = function(lengthTrajectories, lengthAtMaturity) {
   return(ogive)
 } 
 
-
 ```
 
-```{r test}
+```{r test }
+#Do not run 
 
 tic()
 computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growParStich) %>% 
@@ -175,7 +203,9 @@ computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growParStich)
 toc()
 
 tic()
-computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 100000, growParStich) %>% 
+computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern %>% filter(metapop == 'northern iteroparous' ), 
+                                                Nind = 10000, 
+                                                growParStich  %>% filter(metapop == 'northern iteroparous' )) %>% 
   filter(season  == 'spring')
 toc()
 
@@ -187,12 +217,31 @@ computeMultipleLengthTrajectories(temperaturePattern %>%  filter(age < 10) ,
 ```
 
 ```{r ogive variability}
+#Do not run 
 tic(msg = 'ogive variability')
 replicate(10, computeMultipleLengthTrajectories(temperaturePattern %>%  filter(age < 10) , 
                                                 Nind = 100000, growPar) %>% 
             filter(season  == 'spring') %>% 
             computeOgive(lengthAtMaturity = 35) %>% select(p) %>% unlist(use.names = FALSE))
 toc()
+
+#----------------------------------------
+replicate(10, computeMultipleLengthTrajectories(temperaturePattern %>%  filter(age < 10) , 
+                                                Nind = 100000, growPar) %>% 
+            filter(season  == 'spring') %>% 
+            computeOgive(lengthAtMaturity = 35) %>% select(p) %>% unlist(use.names = FALSE))
+
+#plot results 
+growParFemale = growPar %>% pivot_longer(1:11,names_to = "name",values_to = 'par') %>% 
+  tempMinGrow
+
+
+
+par_maturity_optim= enframe(c(res_A$par,fixedPar)) %>% 
+  pivot_wider() %>% 
+  select(contains("Maturity"))
+
+
 ```
 
 ```{r growth curve envelope}
@@ -220,7 +269,8 @@ computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growParStic
   filter(L>growParStich$Linf)
 ```
 
-```{r optimisation}
+```{r optimisation - OBSOLETE}
+#obsolete (presque) au k ou!
 ogiveObs = data.frame(age = seq.int(9), p = c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1))
 
 lengthTrajectoriesInSpring = computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growPar) %>% 
@@ -260,7 +310,7 @@ res_A <- optim(par = c(lengthAtMaturity = 35),
                upper = growPar$Linf)
 
 
-computeOgive(yearlyLengthTrajectories = lengthTrajectoriesInSpring,
+computeOgive(lengthTrajectories = lengthTrajectoriesInSpring,
              lengthAtMaturity = res_A$par) %>% 
   inner_join(ogiveObs, by = 'age', suffix = c(".pred", ".obs"))
 
@@ -329,15 +379,20 @@ computeMultipleLengthTrajectories(temperaturePattern,
   geom_abline(slope =0, intercept = 44)
 ```
 
-```{r calibration K and lengthAtMaturity by gender with ogive}
- #par = res_D$par
+```{r calibration K and lengthAtMaturity by gender with ogive, warning = "FALSE"}
+#par = res_D$par
+
+#----------------------------------------------#
+#------------ Optim function ------------------#
+#-----------------------------------------------#
 
 objFn_D = function(par, fixedPar, ogivesObs, parGrowRef,  temperaturePattern, 
-                    Nind = 1000, seasonSelected = 'spring',
-                    RNGseed =1) {
+                   Nind = 1000, seasonSelected = 'spring',
+                   RNGseed =1) {
   
-  fullPar <- enframe(c(par, fixedPar)) %>%  pivot_wider() 
+  fullPar <- enframe(c(par, fixedPar)) %>%  pivot_wider() #full set of parameters
   
+  #extract parameters by sex
   parFemale = fullPar %>% select(tempMinGrow, tempOptGrow, tempMaxGrow, 
                                  lengthAtHatching, 
                                  Linf = linfVonBert, 
@@ -359,7 +414,7 @@ objFn_D = function(par, fixedPar, ogivesObs, parGrowRef,  temperaturePattern,
                                                                           RNGseed = RNGseed) %>%
     filter(season == seasonSelected) %>% 
     computeOgive(lengthAtMaturity = parFemale$lengthAtMaturity) %>% 
-        # add a penalty when p(max(age)) < 1
+    # add a penalty when p(max(age)) < 1 - avoid immature individuals at maximuma age
     mutate(p = if_else(age == max(age) & p < 1, 10*(p-1)^2, p))
   
   SSEfemale = ogivesObs %>% select(age, p = p_female) %>% 
@@ -387,69 +442,490 @@ objFn_D = function(par, fixedPar, ogivesObs, parGrowRef,  temperaturePattern,
     unlist(use.names = FALSE)
   
   ### SSE for growth curve
-    femaleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , 
-                                                  Nind = 1, 
-                                                  growPar = parFemale %>% mutate(sigmaDeltaLVonBert=0) , 
-                                                  RNGseed = RNGseed) %>% 
+  #growth curve for female by using temperature pattern with sigma = 0 (mean_curve)
+  femaleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , 
+                                                                 Nind = 1, 
+                                                                 growPar = parFemale %>% mutate(sigmaDeltaLVonBert=0) , 
+                                                                 RNGseed = RNGseed) %>% 
     filter(season == seasonSelected) %>% 
-    select(age, Lfemale = L)
-    
+    select(age, Lfemale = L)#rename L by Lfemale 
+  
   maleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , 
-                                                  Nind = 1, 
-                                                  growPar = parMale %>% mutate(sigmaDeltaLVonBert=0) , 
-                                                  RNGseed = RNGseed) %>% 
+                                                               Nind = 1, 
+                                                               growPar = parMale %>% mutate(sigmaDeltaLVonBert=0) , 
+                                                               RNGseed = RNGseed) %>% 
+    filter(season == seasonSelected) %>% #1 value for a season 
+    select(age, Lmale = L)
+  
+  #SCE1 for growth curve from Stich et al, and growth curve by sex  
+  SCE1growth <- 
+    femaleCurve %>% 
+    inner_join(maleCurve, by = 'age') %>% 
+    mutate(Lref = vonBertalanffyGrowth(age, parGrowRef$L0, parGrowRef$Linf, parGrowRef$K),
+           L2genders = (Lfemale + Lmale)/2,
+           #SE = ((Lfemale + Lmale)/2 - Lref)^2 ) %>%
+           SE = (((Lfemale + Lmale)/2 - Lref)/Lref)^2) %>% #TODO TESTER 
+    summarise(SSE = sum(SE)) %>%  unlist(use.names = FALSE)
+  
+  SSEgrowAndMat = (SSEmale + SSEfemale)*10 + SCE1growth
+  
+  return(SSEgrowAndMat)# 10 for weigthed the sumSquare compared to the sum squared for gives so that both sum squared varied in the same interval
+}
+
+
+```
+
+```{r Camille - calibration by gender and metapop with ogive, warnin = "FALSE" }
+
+#----------------------------------------------#
+#------------ Optim function ------------------#
+#-----------------------------------------------#
+
+#Compute the SSE for one metapopulation 
+SSEforOneMetapop = function(fullPar, ogivesObs, parGrowRef,  temperaturePatternforOneMetapop, 
+                            Nind = 1000, seasonSelected = 'spring',
+                            RNGseed =1) {
+  
+  #tic(msg="par")
+  #extract parameters by sex
+  parFemale = fullPar %>% select(tempMinGrow, tempOptGrow, tempMaxGrow, 
+                                 lengthAtHatching, 
+                                 Linf = linfVonBert, 
+                                 kOpt = kOptForFemale,
+                                 sigmaDeltaLVonBert,
+                                 lengthAtMaturity = lFirstMaturityForFemale) 
+  
+  parMale = fullPar %>% select(tempMinGrow, tempOptGrow, tempMaxGrow, 
+                               lengthAtHatching, 
+                               Linf = linfVonBert, 
+                               kOpt = kOptForMale,
+                               sigmaDeltaLVonBert, 
+                               lengthAtMaturity = lFirstMaturityForMale)  
+  #toc()
+  
+  #tic(msg="ogive") #---------------------------------------------
+  # SSE for female ogive
+  ogivePredictFemale  <-  computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePatternforOneMetapop, 
+                                                                          Nind = Nind, 
+                                                                          growPar = parFemale, 
+                                                                          RNGseed = RNGseed) %>%
     filter(season == seasonSelected) %>% 
+    computeOgive(lengthAtMaturity = parFemale$lengthAtMaturity) %>% 
+    # add a penalty when p(max(age)) < 1 - avoid immature individuals at maximuma age
+    mutate(p = if_else(age == max(age) & p < 1, 10*(p-1)^2, p))
+  
+  SSEfemale = ogivesObs %>% 
+    select(age, p = p_female) %>% 
+    inner_join(ogivePredictFemale, 
+               by = 'age', suffix = c('.obs','.pred')) %>% 
+    #add penality 
+    mutate(squareError = (p.obs - p.pred)^2) %>% 
+    summarise(SSE = sum(squareError),) %>% 
+    unlist(use.names = FALSE)
+  
+  # SSE for male ogive
+  ogivePredictMale <- computeMultipleLengthTrajectoriesWithRandomSeed(
+    temperaturePatternforOneMetapop , 
+    Nind = Nind, 
+    growPar = parMale, 
+    RNGseed = RNGseed) %>% 
+    filter(season  == seasonSelected) %>% 
+    computeOgive(lengthAtMaturity = parMale$lengthAtMaturity) %>% 
+    # add a penalty when p(max(age)) < 1
+    mutate(p = if_else(age == max(age) & p < 1, 10*(p-1)^2, p))
+  
+  SSEmale = ogivesObs %>% 
+    select(age, p = p_male) %>%  
+    inner_join(ogivePredictMale, 
+               by = 'age', suffix = c('.obs','.pred')) %>% 
+    mutate(squareError = (p.obs - p.pred)^2) %>% 
+    summarise(SSE = sum(squareError)) %>% 
+    unlist(use.names = FALSE)
+  #toc()
+  
+  #tic(msg = 'growth')# -----------------------------------------------
+  ### SSE for growth curve
+  #growth curve for female by using temperature pattern with sigma = 0 (mean_curve)
+  femaleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed( 
+    temperaturePatternforOneMetapop , 
+    Nind = 1, 
+    growPar = parFemale %>% mutate(sigmaDeltaLVonBert=0) , 
+    RNGseed = RNGseed) %>% 
+    filter(season == seasonSelected) %>% 
+    select(age, Lfemale = L)#rename L by Lfemale 
+  
+  maleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(
+    temperaturePatternforOneMetapop, 
+    Nind = 1, 
+    growPar = parMale %>% mutate(sigmaDeltaLVonBert=0) , 
+    RNGseed = RNGseed) %>% 
+    filter(season == seasonSelected) %>% #1 value for a season 
     select(age, Lmale = L)
+  #toc()
   
+  #tic('SSE')
+  #SCE1 for growth curve from Stich et al, and growth curve by sex  
   SCE1growth <- 
- femaleCurve %>% 
+    femaleCurve %>% 
     inner_join(maleCurve, by = 'age') %>% 
     mutate(Lref = vonBertalanffyGrowth(age, parGrowRef$L0, parGrowRef$Linf, parGrowRef$K),
-          L2genders = (Lfemale + Lmale)/2,
-           SE = ((Lfemale + Lmale)/2 - Lref)^2 ) %>%
+           L2genders = (Lfemale + Lmale)/2,
+           #SE = ((Lfemale + Lmale)/2 - Lref)^2 ) %>%
+           SE = (((Lfemale + Lmale)/2 - Lref)/Lref)^2) %>% #TODO TESTER 
     summarise(SSE = sum(SE)) %>%  unlist(use.names = FALSE)
   
+  SSEgrowthAndOgives = (SSEmale + SSEfemale)*10 + SCE1growth # 10 for weigthed the sumSquare compared to the sum squared for gives so that both sum squared varied in the same interval
+  #toc()
+  return(SSEgrowthAndOgives)
+}
 
-  return((SSEmale + SSEfemale)*10 + SCE1growth)
+
+#Fonction to be optimzed over 
+SSEforAllMetapop = function(par, fixedPar, regional_metapop,ogivesObs, parGrowRef, temperaturePattern, Nind, seasonSelected, RNGseed = 1){
+  
+  fullPar = enframe(c(par, fixedPar)) %>%  pivot_wider()
+  
+  SSE_all_metapop = matrix(NA, length(regional_metapop), 1)
+  
+  i = 1 
+  
+  for (metapopulation in regional_metapop){
+    
+    SSE_all_metapop[i] <- SSEforOneMetapop(
+      fullPar = fullPar, 
+      ogivesObs = ogivesObs %>%  filter(metapop == metapopulation),
+      parGrowRef = parGrowRef %>%  filter(metapop == metapopulation),
+      temperaturePatternforOneMetapop = temperaturePattern %>%  filter(metapop == metapopulation),
+      Nind = Nind, 
+      seasonSelected = seasonSelected,
+      RNGseed = RNGseed)
+    i = i + 1 
+    
+  }
+  
+  SSE = sum(SSE_all_metapop)
+  
+  return(SSE)
+  
 }
 
-ogivesObs = data.frame(age = seq.int(9), p_female = c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1),
-                       p_male =  c(0, 0, 0.01, 0.15, 0.45, 0.75, 0.92, 1, 1))
+#Verification
+tic()
+SSEforAllMetapop(par = par_0,
+                 fixedPar = fixedPar,
+                 regional_metapop = regional_metapop,
+                 ogivesObs = ogivesObs,
+                 parGrowRef = parGrowRef,
+                 temperaturePattern =  temperaturePattern,
+                 Nind = 1000,
+                 seasonSelected = "spring",
+                 RNGseed = 1)
+toc()
+tic()
+metapopulation = "semelparous"
+SSEforOneMetapop(fullPar = enframe(c(par_0, fixedPar)) %>%  pivot_wider(), 
+                 ogivesObs = ogivesObs %>%  filter(metapop == metapopulation), 
+                 parGrowRef = parGrowRef %>%  filter(metapop == metapopulation), 
+                 temperaturePatternforOneMetapop = temperaturePattern %>% filter(metapop == metapopulation),
+                 Nind = 1000, 
+                 seasonSelected = 'spring',
+                 RNGseed = 1)
+toc()
+
+fixedPar0 = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020  
+              tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020 
+              lengthAtHatching = 2.8)
+#  avirer
+parFemale = enframe(c(resA_par, fixedPar0)) %>%  pivot_wider() %>% select(tempMinGrow, tempOptGrow, tempMaxGrow, 
+                                                                          lengthAtHatching, 
+                                                                          Linf = linfVonBert, 
+                                                                          kOpt = kOptForFemale,
+                                                                          sigmaDeltaLVonBert,
+                                                                          lengthAtMaturity = lFirstMaturityForFemale) 
+
+
+#================================ ICI =======================================
+rm(lengthTrajectories)
+lengthTrajectories = computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern = temperaturePattern %>% filter(metapop == "semelparous"),
+                                                                     Nind = 1000,
+                                                                     growPar = parFemale,
+                                                                     RNGseed = 1) %>% 
+  arrange(ind,age) %>% 
+  filter(season  == "spring") 
+
+
+lengthTrajectories %>% 
+  ggplot(aes(x = age, y= L, color = as.factor(ind))) +
+geom_path( show.legend = FALSE)
+
+
+lengthTrajectories %>%  
+  group_by(age) %>% 
+  summarise(L = mean(L)) %>% 
+  print(n=Inf)
+
+lengthTrajectories %>%  
+  computeOgive(lengthAtMaturity = parFemale$lengthAtMaturity) %>% 
+  print(n=Inf)
+
+
+
+
+vonBertalanffyWithRandomVector(L = c(2.8, 2.8), 
+                               Linf = parFemale$Linf, 
+                               timeStepDuration = 0.25, 
+                               K = parFemale$kOpt, 
+                               randomVector =  c(0.184, -0.253),
+                               sigma = parFemale$sigmaDeltaLVonBert, 
+                               tempEffect = c(0.0639, 0.0639))
+```
+
+```{r set parameters values for calibration, warning = "FALSE"}
+
+#----------------------------------------------#
+#------------ Starting parameters values ------------------#
+#-----------------------------------------------#
+
+# 
+# ogivesObs = data.frame(age = seq.int(9), p_female = c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1),
+#                        p_male =  c(0, 0, 0.07, 0.29, 0.60, 0.86, 1, 1, 1))
+
+# for 3 metapopulations
+ogivesObs = expand_grid(age = seq.int(9),  metapop = regional_metapop) %>% 
+  arrange(metapop, age) %>% 
+  group_by(metapop) %>% 
+  mutate(p_female = fifelse(
+    metapop == "semelparous", c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1), fifelse (metapop == "southern iteroparous", c(0, 0, 0,0.04,0.27,0.64,0.81,0.9,1), c(0,0,0,0.04,0.69,0.69,0.9,1,1)))) %>% 
+  mutate(p_male = fifelse(
+    metapop == "semelparous", c(0, 0, 0.07, 0.29, 0.60, 0.86, 1, 1, 1), fifelse (metapop == "southern iteroparous", c(0,0,0,0.12,0.5,0.79,0.8,0.84,1), c(0,0,0.03,0.21,0.82,1,1,1,1)))) %>% 
+  ungroup()
+
+
+#----------Growth parameters------------#
+#----------------------------------------
+
+#With parameters values from Stich et al and Topt from optimization on growParUnisex
+#metapopSelected = "semelparous"
+
+parGrowRef = growParStich  %>% 
+  select(L0 = lengthAtHatching, K, Linf, metapop)
+
+#With parameters values from optimization 
+#parGrowRef = growParOptim %>% select(L0 = lengthAtHatching, K, Linf)
+
+#----------Parameters to be optimized over -----------#
+#----------------------------------------------------#
+#Case 1
+#parameters from optimisation and XML file 
+par_0 = c(tempOptGrow = 5, #from optimisation on growth curves 
+          kOptForFemale = 0.32, #XML - Alosa alosa 
+          kOptForMale = .21, #XML - Alosa alosa 
+          lFirstMaturityForMale = 40, #XML - Alosa sapidissima
+          lFirstMaturityForFemale = 45)#XML - Alosa sapidissima
+
+
+fixedPar = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020  
+             tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020 
+             lengthAtHatching = 2.8,# XML - Alosa sapidissima
+             linfVonBert = 76, #XML - Alosa alosa 
+             sigmaDeltaLVonBert = .2)  #XML - Alosa alosa 
+#optimisation
+
+# XML - Alosa sapidissia)
+res_Abis <- optim(par = par_0,
+                  fn = SSEforAllMetapop,
+                  #method = "L-BFGS-B", 
+                  control = list(trace = 1,
+                                 maxit = 1000),
+                  #hessian = TRUE,
+                  fixedPar = fixedPar, 
+                  regional_metapop = regional_metapop,
+                  ogivesObs = ogivesObs,
+                  parGrowRef = parGrowRef,
+                  temperaturePattern = temperaturePattern, 
+                  seasonSelected = "spring",
+                  Nind = 5000,
+                  RNGseed = 1)
+
+
+#plot results 
+par_maturity_optim= enframe(c(res_A$par,fixedPar)) %>% 
+  pivot_wider() %>% 
+  select(contains("Maturity"))
+
+par_grow_optim_female = enframe(c(res_A$par,fixedPar)) %>% 
+  pivot_wider() %>% 
+  select(-c('kOptForMale','lFirstMaturityForMale')) 
+
+
+#growth curves
+growthCurves_female = computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , Nind = 10000, 
+                                                                      par_grow_optim_female %>% 
+                                                                        pivot_wider(), RNGseed = 1)
+
+
+computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , Nind = 1000, growParStich) %>% 
+  group_by(age) %>% 
+  summarise(Lmin = min(L), 
+            L025 = quantile(L, .25), 
+            L25 = quantile(L, .25), 
+            Lmed = median(L), 
+            L75 = quantile(L, .75), 
+            L975 = quantile(L, .975), 
+            Lmax = max(L)) %>% 
+  ggplot(aes(x = age)) +
+  geom_ribbon(aes(ymin = Lmin, ymax = Lmax), fill = "gray30", alpha = .2) + 
+  geom_ribbon(aes(ymin = L025, ymax = L975), fill = "gray30", alpha = .2) + 
+  geom_ribbon(aes(ymin = L25, ymax = L75), fill = "gray30", alpha = .2) + 
+  geom_line(aes(y = Lmed)) +
+  labs(x = 'Age (y)', y = 'Length (cm)')  +
+  geom_abline(slope = 0, intercept = c(40.3, 42)) 
+
+
 
 
-#TODO verifier si ce sont bien les parametres de Stich
-parGrowRef = growParStich %>% select(L0 = lengthAtHatching, K, Linf)
 
-par = c(tempOptGrow = 5, 
-        kOptForFemale = 0.32, 
-        kOptForMale = .21, 
-        linfVonBert = 76, 
-        lFirstMaturityForMale = 40, 
-        lFirstMaturityForFemale = 45,
-        sigmaDeltaLVonBert = .2)
 
-fixedPar  = c(tempMinGrow = 1.6,  
-              tempMaxGrow = 27.9,
-              lengthAtHatching = parGrowRef$L0 )
 
 
+#growth curve according to metapopulations 
+growthCurves %>% 
+  ggplot(aes(x = age)) + 
+  geom_line(aes(y = L, col = metapop), linetype = "dashed") +
+  geom_line(aes(y = LStich, col = metapop))+
+  labs( x = 'Age (year)', y = 'Fork length (cm)', color = "Metapopulation (Stich)") +
+  #scale_x_continuous(limits = c(0, 10)) +
+  geom_hline(yintercept = 45, color = "pink", linetype =1)+ #female
+  geom_hline(yintercept = 40, color = "lightblue", linetype = 1)+
+  geom_vline(aes(age), xintercept = 8, linetype = 1)+
+  facet_wrap(.~ metapop, ncol =1)
 
-# objFn_D(par = par,
-#         fixedPar = fixedPar, 
-#         ogivesObs = ogivesObs,
-#         parGrowRef = parGrowRef,
-#         temperaturePattern = temperaturePattern, 
-#         Nind = 5000)
-     
-res_D <- optim(par = par,
-               fn = objFn_D,
-               fixedPar = fixedPar, 
+#temperature effect 
+temp <- seq(0,45,0.1)
+
+growthCurves %>% 
+  ggplot(aes(x = temperature)) + 
+  geom_line(aes(y = temperatureEffect)) +
+  geom_jitter(aes(y = temperatureEffect, col = metapop))
+
+
+#Case 2
+
+fixedPar_1 = c(tempMinGrow = 1.6,  
+               tempMaxGrow = 27.9,
+               lengthAtHatching = 5.87)
+
+
+res_B <- optim(par = par_0,
+               fn = SSEforAllMetapop,
+               #method = "BFGS", 
+               control = list(trace = 1,
+                              maxit = 1000),
+               #hessian = TRUE,
+               fixedPar = fixedPar_1, 
+               regional_metapop = regional_metapop,
                ogivesObs = ogivesObs,
                parGrowRef = parGrowRef,
                temperaturePattern = temperaturePattern, 
+               seasonSelected = "spring",
                Nind = 5000,
-               RNGseed =1,
+               RNGseed =1)
+
+
+save(res_A, res_B, res_C, res_D,file = "res_optim.Rdata")
+
+
+#Case 3
+par_1 = c(tempOptGrow = 5, 
+          kOptForFemale = 0.32, 
+          kOptForMale = .21, 
+          linfVonBert = 76, 
+          lFirstMaturityForMale = 40, 
+          lFirstMaturityForFemale = 45,
+          lengthAtHatching = 2.8,
+          sigmaDeltaLVonBert = .2)
+
+fixedPar_2 = c(tempMinGrow = 1.6,  
+               tempMaxGrow = 27.9)
+
+
+res_C <- optim(par = par_1,
+               fn = SSEforAllMetapop,
+               #method = "BFGS", 
                control = list(trace = 1,
-                              maxit = 1000))
+                              maxit = 1000),
+               #hessian = TRUE,
+               fixedPar = fixedPar_2, 
+               regional_metapop = regional_metapop,
+               ogivesObs = ogivesObs,
+               parGrowRef = parGrowRef,
+               temperaturePattern = temperaturePattern, 
+               seasonSelected = "spring",
+               Nind = 5000,
+               RNGseed =1)
+
+
+par_2 = c(tempOptGrow = 5, 
+          kOptForFemale = 0.32, 
+          kOptForMale = .21, 
+          linfVonBert = 47.5, 
+          lFirstMaturityForMale = 40, 
+          lFirstMaturityForFemale = 45,
+          lengthAtHatching = 5.87)
+
+fixedPar_3 = c(tempMinGrow = 1.6,  
+               tempMaxGrow = 27.9,
+               sigmaDeltaLVonBert = .2)
+
+res_D <- optim(par = par_2,
+               fn = SSEforAllMetapop,
+               #method = "BFGS", 
+               control = list(trace = 1,
+                              maxit = 1000),
+               #hessian = TRUE,
+               fixedPar = fixedPar_3, 
+               regional_metapop = regional_metapop,
+               ogivesObs = ogivesObs,
+               parGrowRef = parGrowRef,
+               temperaturePattern = temperaturePattern, 
+               seasonSelected = "spring",
+               Nind = 5000,
+               RNGseed =1)
+
+
+
+
+
+
+
+
+
+
 
 ```
+
+```{r optimization, warning = "FALSE"}
+
+#----------------------------------------------#
+#------------ Optmization----------------------#
+#----------------------------------------------#
+
+# res_D <- optim(par = par,
+#                fn = objFn_D,
+#                fixedPar = fixedPar, 
+#                ogivesObs = ogivesObs,
+#                parGrowRef = parGrowRef,
+#                temperaturePattern = temperaturePattern, 
+#                Nind = 5000,
+#                RNGseed =1,
+#                control = list(trace = 1,
+#                               maxit = 1000))
+
+
+```
+
+```{r}
+data.frame(L=rnorm(10)) %>% 
+  mutate(mu=ifelse(L<0, 0, L*10))
+```
+
-- 
GitLab