new version Rmd for maturation.

exploration/NEA_calibration_offline/maturationPL.Rmd

+---
+date: "`r Sys.Date()`"
+author: "Your Name"
+title: "officedown template"
+output: 
+  officedown::rdocx_document:
+    mapstyles:
+      Normal: ['First Paragraph']
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE, fig.cap = TRUE)
+library(officedown)
+library(officer)
+
+library(tidyverse)
+library(tictoc) # see timeR
+
+
+
+```
+
+```{r}
+rm(list = ls())
+source("../GR3D_Rdescription/GR3Dfunction.R")
+
+load('SOS.rdata')
+```
+
+```{r}
+source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
+
+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(.)))
+
+
+growParStich <- Stich2020_sel %>% filter(catchment == "semelparous") %>% 
+  select(Linf, K,  L0_theo) %>% 
+  mutate(tempMinGrow = 1.6, tempOptGrow = 5, tempMaxGrow = 27.9, 
+         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") %>% 
+  # 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)
+
+```
+
+```{r local function}
+
+# 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, 
+                                                L0 = growPar$lengthAtHatching, 
+                                                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$age
+  
+  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, 
+                                               L0 = growPar$lengthAtHatching, 
+                                               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)
+}
+
+
+vonBertalanffyWithRandomVector = function(L, L0, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
+  if (sigma == 0) {
+    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
+    increment = exp(mu)
+  }
+  else {
+    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
+    increment = exp(randomVector * sigma + mu)
+  }
+  L = L + increment * tempEffect
+  return(L)
+}
+
+
+computeOgive = function(lengthTrajectories, lengthAtMaturity) {
+  ogive <- lengthTrajectories %>% 
+    distinct(age) %>% arrange(age) %>% 
+    left_join(lengthTrajectories %>% filter(L > lengthAtMaturity) %>% 
+                group_by(ind) %>% 
+                summarise(age = min(age), .groups = 'drop') %>% 
+                group_by(age) %>% 
+                summarise(mature = n(), .groups = 'drop'), 
+              by = 'age') %>%     
+    replace_na(list(mature = 0)) %>% 
+    mutate(immature = sum(mature) - cumsum(mature),
+           p = if_else(mature + immature == 0, 0,  mature/(mature + immature))) %>% 
+    select(age, p)
+  
+  return(ogive)
+} 
+
+
+```
+
+```{r test}
+
+tic()
+computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growParStich) %>% 
+  filter(season  == 'spring')
+toc()
+
+tic()
+computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 100000, growParStich) %>% 
+  filter(season  == 'spring')
+toc()
+
+computeMultipleLengthTrajectories(temperaturePattern %>%  filter(age < 10) , 
+                                  Nind = 10000, growPar) %>% 
+  filter(season  == 'spring') %>% 
+  computeOgive(lengthAtMaturity = 35) 
+
+```
+
+```{r ogive variability}
+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()
+```
+
+```{r growth curve envelope}
+# envelope of growth curves
+computeMultipleLengthTrajectories(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)) 
+
+
+# verify that all L < Linf
+computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growParStich) %>% 
+  filter(L>growParStich$Linf)
+```
+
+```{r optimisation}
+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) %>% 
+  filter(season  == 'spring')
+
+# calibrate only lengthAtMaturity
+objFn_A = function(par, ogiveObs, yearlyLengthTrajectories){
+  SSE = ogiveObs %>%
+    inner_join(computeOgive(yearlyLengthTrajectories,
+                            lengthAtMaturity = par[1]),
+               by = 'age', suffix = c('.obs','.pred')) %>%
+    mutate(squareError = (p.obs - p.pred)^2) %>%
+    summarise(SSE = sum(squareError)) %>%
+    unlist(use.names = FALSE)
+  return(SSE)
+  # epsilon = 1e-6 # to avoid log(0)
+  #   minusSumLogL = ogiveObs %>% 
+  #   inner_join(computeOgive(yearlyLengthTrajectories,
+  #                           lengthAtMaturity = par[1]), 
+  #              by = 'age', suffix = c('.obs','.pred')) %>% 
+  #   summarise(minusSumLogL = - sum(p.obs * log(p.pred + epsilon) + (1 - p.obs) * log(1 - p.pred + epsilon))) %>% 
+  #   unlist(use.names = FALSE)
+  # return(minusSumLogL)
+}
+
+
+objFn_A(par = c(lengthAtMaturity = 40), 
+        ogiveObs, 
+        yearlyLengthTrajectories = lengthTrajectoriesInSpring)
+
+res_A <- optim(par = c(lengthAtMaturity = 35),
+               fn = objFn_A,
+               ogiveObs = ogiveObs,
+               yearlyLengthTrajectories = lengthTrajectoriesInSpring,
+               method = 'Brent',
+               lower = 0,
+               upper = growPar$Linf)
+
+
+computeOgive(yearlyLengthTrajectories = lengthTrajectoriesInSpring,
+             lengthAtMaturity = res_A$par) %>% 
+  inner_join(ogiveObs, by = 'age', suffix = c(".pred", ".obs"))
+
+# ----------------------------------------------------------------------------------------------------
+# calibrate  lengthAtMaturity and Kopt
+objFn_B = function(par, growPar, ogiveObs,  temperaturePattern, Nind = 100000, seasonSelect = 'spring', sigmaDeltaLVonBert = .2){
+  growPar$K = par[2]
+  growPar$sigmaDeltaLVonBert = sigmaDeltaLVonBert
+  
+  ogivePredict <- computeMultipleLengthTrajectories(temperaturePattern , Nind = Nind, growPar) %>% 
+    filter(season  == seasonSelect) %>% 
+    computeOgive(lengthAtMaturity = par[1])
+  
+  SSE = ogiveObs %>% 
+    inner_join(ogivePredict, 
+               by = 'age', suffix = c('.obs','.pred')) %>% 
+    mutate(squareError = (p.obs - p.pred)^2) %>% 
+    summarise(SSE = sum(squareError)) %>% 
+    unlist(use.names = FALSE)
+  return(SSE)
+}
+
+objFn_B(par = c(lengthAtMaturity = 35, Kopt = 0.49 ),
+        growPar = growPar, 
+        ogiveObs = ogiveObs,
+        temperaturePattern = temperaturePattern, 
+        Nind = 100000, 
+        sigmaDeltaLVonBert = .2)
+
+res_B <- optim(par = c(lengthAtMaturity = 35, Kopt = 0.49 ),
+               fn = objFn_B,
+               growPar = growPar,
+               ogiveObs = ogiveObs,
+               temperaturePattern = temperaturePattern, 
+               Nind = 100000, 
+               sigmaDeltaLVonBert = .2 , 
+               control = list(trace = 1,
+                              reltol = 1e-5))
+
+save(res_B, file = 'res_B.Rdata')
+
+computeMultipleLengthTrajectories(temperaturePattern, 
+                                  Nind = 100000, 
+                                  growPar %>% mutate(K = res_B$par[2])) %>% 
+  filter(season  == 'spring') %>% 
+  computeOgive(lengthAtMaturity = res_B$par[1]) %>% 
+  inner_join(ogiveObs,  by = 'age', suffix = c('.pred', '.obs')) 
+
+computeMultipleLengthTrajectories(temperaturePattern, 
+                                  Nind = 10000, 
+                                  growPar %>% mutate(K = res_B$par[2])) %>% 
+  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 = 44)
+```
+
+```{r calibration K and lengthAtMaturity by gender with ogive}
+ #par = res_D$par
+
+objFn_D = function(par, fixedPar, ogivesObs, parGrowRef,  temperaturePattern, 
+                    Nind = 1000, seasonSelected = 'spring',
+                    RNGseed =1) {
+  
+  fullPar <- enframe(c(par, fixedPar)) %>%  pivot_wider() 
+  
+  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)  
+  
+  # SSE for female ogive
+  ogivePredictFemale  <-  computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , 
+                                                                          Nind = Nind, 
+                                                                          growPar = parFemale, 
+                                                                          RNGseed = RNGseed) %>%
+    filter(season == seasonSelected) %>% 
+    computeOgive(lengthAtMaturity = parFemale$lengthAtMaturity) %>% 
+        # add a penalty when p(max(age)) < 1
+    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')) %>% 
+    mutate(squareError = (p.obs - p.pred)^2) %>% 
+    summarise(SSE = sum(squareError)) %>% 
+    unlist(use.names = FALSE)
+  
+  # SSE for male ogive
+  ogivePredictMale <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , 
+                                                                      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)
+  
+  ### SSE for growth curve
+    femaleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , 
+                                                  Nind = 1, 
+                                                  growPar = parFemale %>% mutate(sigmaDeltaLVonBert=0) , 
+                                                  RNGseed = RNGseed) %>% 
+    filter(season == seasonSelected) %>% 
+    select(age, Lfemale = L)
+    
+  maleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , 
+                                                  Nind = 1, 
+                                                  growPar = parMale %>% mutate(sigmaDeltaLVonBert=0) , 
+                                                  RNGseed = RNGseed) %>% 
+    filter(season == seasonSelected) %>% 
+    select(age, Lmale = L)
+  
+  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 ) %>%
+    summarise(SSE = sum(SE)) %>%  unlist(use.names = FALSE)
+  
+
+  return((SSEmale + SSEfemale)*10 + SCE1growth)
+}
+
+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))
+
+
+#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 )
+
+
+
+# objFn_D(par = par,
+#         fixedPar = fixedPar, 
+#         ogivesObs = ogivesObs,
+#         parGrowRef = parGrowRef,
+#         temperaturePattern = temperaturePattern, 
+#         Nind = 5000)
+     
+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))
+
+```