with a flexible temperature effect in growth calibration

exploration/GR3D_Rdescription/GR3Dfunction.R

@@ -3,22 +3,42 @@
 #----------------------------------------------
 # temperature effect in GR3D
 # from Rosso
-temperatureEffect = function(tempWater, Tmin, Topt, Tmax){
+
+temperatureEffect = function(tempWater, Tmin, Topt, Tmax, effectFunction = Rosso){
+  return(effectFunction(tempWater, Tmin, Topt, Tmax))
+}
+
+
+
+modifiedLehman = function(temperature, Tmin, Topt, Tmax){
+  # Svirezhev, Yu. M., Krysanova, V. P. & Voinov, A. A. (1984) Mathematical modelling of a fish pond ecosystem. 
+  #   Ecological Modelling 21, 315–337.
+  
+  # rate(Tmin or Tmax) = exp(-4.6) = 1 %
+  
+  rate = rep(0, length(temperature))
+  rate[temperature < Topt] = exp(-4.6 * ((Topt - temperature[temperature < Topt]) / (Topt - Tmin))^4)
+  rate[temperature >= Topt] = exp(-4.6 * ((temperature[temperature >= Topt] - Topt) / (Tmax - Topt))^4)
+  
+  return(rate)
+}
+
+Rosso = function(tempWater, Tmin, Topt, Tmax){
   response = (tempWater - Tmin) * (tempWater - Tmax) / ((tempWater - Tmin) * (tempWater - Tmax) - (tempWater - Topt)^2)
   
   response[tempWater <= Tmin | tempWater >= Tmax] = 0
   return(response)
 }
 
-thermalRange = function(pct = 0.8, Tmin, Topt, Tmax){
-  lower = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax) - pct,
+thermalRange = function(pct = 0.8, Tmin, Topt, Tmax,...){
+  lower = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax,...) - pct,
            interval = c(Tmin, Topt))$root
-  upper = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax) - pct,
+  upper = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax,...) - pct,
                   interval = c(Topt, Tmax))$root
   return(c(lower = lower,upper = upper))
 }
 
-#optimalThermalRange(Tmin = 3, Topt = 17, Tmax =27)
+# thermalRange(Tmin = 3, Topt = 17, Tmax = 27, effectFunction = modifiedLehman )
 
 #----------------------------------------------
 # Growth simulation
@@ -69,7 +89,7 @@ vonBertalanffyWithNextIncrement = function(L,  Linf, K, timeStepDuration, sigma,
 }
 
 # generate a cohort of length trajectories
-computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar){
+computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar, ...){
   ages = temperaturePattern$age
   
   res <- expand_grid(ind = seq.int(Nind),  age = ages) %>% 
@@ -78,7 +98,7 @@ computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, grow
     mutate(temperatureEffect = temperatureEffect(temperature, 
                                                  growPar$tempMinGrow, 
                                                  growPar$tempOptGrow, 
-                                                 growPar$tempMaxGrow), 
+                                                 growPar$tempMaxGrow, ...), 
            L = if_else(age == 0, growPar$lengthAtHatching, 0))
   
   for (i in 2:length(ages)) {
@@ -101,7 +121,7 @@ computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, grow
 computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern, 
                                                            Nind = 10, 
                                                            growPar, 
-                                                           RNGseed =1){
+                                                           RNGseed =1, ...){
   set.seed(RNGseed)
   ages = temperaturePattern %>% 
     distinct(age) %>% 
@@ -114,7 +134,7 @@ computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
     mutate(temperatureEffect = temperatureEffect(temperature, 
                                                  growPar$tempMinGrow, 
                                                  growPar$tempOptGrow, 
-                                                 growPar$tempMaxGrow), 
+                                                 growPar$tempMaxGrow, ...), 
            L = if_else(age == 0, growPar$lengthAtHatching, 0))
   
   for (i in 2:length(ages)) {

exploration/NEA_calibration_offline/maturationPL.Rmd

@@ -30,7 +30,7 @@ source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
 load('SOS.rdata')
 #load('res_optim.rdata')
 
-resA_par = read_csv2("res_A.csv", skip =1, col_names=F) %>% 
+resA_par = read_csv2("res_A.csv", skip = 1, col_names = F) %>% 
   deframe()
 
 
@@ -735,31 +735,32 @@ par_0 = c(tempOptGrow = 5, #from optimisation on growth curves
           kOptForMale = .21, #XML - Alosa alosa 
           lFirstMaturityForMale = 40, #XML - Alosa sapidissima
           lFirstMaturityForFemale = 45,#XML - Alosa sapidissima
-          linfVonBert = 76, #XML - Alosa alosa 
-          sigmaDeltaLVonBert = .2)  #XML - Alosa alosa 
+          linfVonBert = 76) #XML - Alosa alosa 
+
 
 
 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
+             lengthAtHatching = 2.8,# XML - Alosa sapidissima
+             sigmaDeltaLVonBert = .5)
 
 #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)
+res_Aquat <- 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 
@@ -938,8 +939,327 @@ res_D <- optim(par = par_2,
 
 ```
 
+
 ```{r}
-data.frame(L=rnorm(10)) %>% 
-  mutate(mu=ifelse(L<0, 0, L*10))
+# 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()
+
+
+
+parGrowRef = growParStich  %>% 
+  select(L0 = lengthAtHatching, K, Linf, metapop)
+
+```
+
+```{r}
+SSEGrowth = function(par, fixedPar, regional_metapop,  parGrowRef,  temperaturePattern, 
+                     seasonSelected = 'spring',
+                     RNGseed =1, ...) {
+  # for both genders
+  fullPar = enframe(c(par, fixedPar)) %>% pivot_wider() %>% 
+    rename(Linf = linfVonBert)
+  
+  # result
+  SSE_all_metapop = matrix(NA, length(regional_metapop), 1)
+  
+  i = 1 
+  for (metapopulation in regional_metapop) {
+    #tic(msg = 'growth')# -----------------------------------------------
+    ### SSE for growth curve
+    #growth curve for female by using temperature pattern with sigma = 0 (mean_curve)
+    growthCurve <- computeMultipleLengthTrajectoriesWithRandomSeed( 
+      temperaturePattern = temperaturePattern %>%  filter(metapop == metapopulation), 
+      Nind = 1, 
+      growPar = fullPar, 
+      RNGseed = RNGseed,
+      ...) %>% 
+      filter(season == seasonSelected) %>% 
+      select(age, L)
+    #toc()
+    
+    #SCE for growth curve from Stich et al
+    parGrow = parGrowRef %>% filter(metapop == metapopulation)
+    SSE_all_metapop[i] <- 
+      growthCurve %>% 
+      mutate(Lref = vonBertalanffyGrowth(age, parGrow$L0, parGrow$Linf, parGrow$K),
+             SE = (L - Lref)^2) %>% 
+      summarise(SSE = sum(SE)) %>%  unlist(use.names = FALSE)
+    
+    i = i + 1
+  }
+  return(sum(SSE_all_metapop))
+}
+
+
+SSEGrowth(par = c(lengthAtHatching = 3.0, kOpt = .5, linfVonBert = 60, tempOptGrow = 10), 
+          fixedPar = c(sigmaDeltaLVonBert = 0, tempMinGrow = 1.6, tempMaxGrow = 27.9) , 
+          regional_metapop,  
+          parGrowRef,  
+          temperaturePattern, 
+          seasonSelected = 'spring',
+          RNGseed = 1, 
+          effectFunction = modifiedLehman)
+
+#lower = c(0,0,0, -2, 1.6, 27.9)
+#upper = c(+Inf, +Inf, +Inf, 1.6, 27.9, 35)
+
+fixedPar = c(sigmaDeltaLVonBert = 0 )
+optimGrowth =  optim(par = c(lengthAtHatching = 3.0, kOpt = .5, linfVonBert = 60,  tempOptGrow = 10, tempMinGrow = 1.0, tempMaxGrow = 30),
+                     fn = SSEGrowth,
+                     # method = "L-BFGS-B", 
+                     control = list(trace = 1,
+                                    maxit = 1000),
+                     #hessian = TRUE,
+                     fixedPar = fixedPar, 
+                     # lower = lower,
+                     # upper = upper, 
+                     regional_metapop = regional_metapop,
+                     parGrowRef = parGrowRef,
+                     temperaturePattern = temperaturePattern, 
+                     seasonSelected = "spring",
+                     RNGseed = 1, 
+                     effectFunction = modifiedLehman)
+
+#set of optim parameters for growth
+fullParOptimGrowth = enframe(c(optimGrowth$par,fixedPar)) %>% 
+  pivot_wider()  %>% 
+  rename(Linf = linfVonBert)
+
+
+# --------------------------------------------- graphs
+growthOptimForAllMetapop = list()
+for (metapopulation in regional_metapop) {
+  
+  parGrowRefOneMetapop = parGrowRef %>% filter(metapop == metapopulation)
+  # simulated growth
+  simulatedOneMetapop <-  computeMultipleLengthTrajectoriesWithRandomSeed( 
+    temperaturePattern = temperaturePattern %>%  filter(metapop == metapopulation), 
+    Nind = 1, 
+    growPar = fullParOptimGrowth, 
+    RNGseed = 1,
+    effectFunction = modifiedLehman) %>% 
+    select(age, L) %>% mutate(metapop = metapopulation, type = "simulated")
+  
+  # observed growth
+  parGrowRefOneMetapop = parGrowRef %>% filter(metapop == metapopulation)
+  observedOneMetapop <- simulatedOneMetapop %>% select(age) %>% 
+    mutate(L = vonBertalanffyGrowth(age, L0 = parGrowRefOneMetapop$L0, Linf = parGrowRefOneMetapop$Linf, K = parGrowRefOneMetapop$K),
+           metapop = metapopulation, type = "Stich")
+  
+  growthOptimForAllMetapop[[metapopulation]] = bind_rows(simulatedOneMetapop, observedOneMetapop )
+}
+
+rbindlist(growthOptimForAllMetapop) %>% 
+  # filter(type == "Stich") %>% 
+  ggplot() + 
+  geom_path(aes(x = age, y = L, col = metapop, linetype = type), show.legend = TRUE) +
+  facet_grid(~metapop)
+
+tibble(temperature = seq(-5,50,1)) %>% 
+  mutate(tempEffect = modifiedLehman(temperature, Tmin =fullParOptim$tempMinGrow,
+                                     Topt = fullParOptim$tempOptGrow,
+                                     Tmax = fullParOptim$tempMaxGrow)) %>% 
+  ggplot() +
+  geom_path(aes(x = temperature, y = tempEffect)) +
+  geom_vline(xintercept = c( 1.6, 27.9), linetype = 2) +
+  geom_hline(yintercept = .8, linetype = 3)
+
+# =================================================
+# calibrate ogives 
+
+SSEOgives = function(par, fixedPar, regional_metapop,  temperaturePattern, 
+                     seasonSelected = 'spring',
+                     Nind = 1000,
+                     RNGseed = 1, ...) {
+  
+  #fixedPar already a tibble
+  par <- par %>% enframe() %>% pivot_wider() %>% 
+    bind_cols(fixedPar)
+  # result
+  SSE_all_metapop = matrix(NA, length(regional_metapop), 1)
+  
+  i = 1 
+  for (metapopulation in regional_metapop) {
+    # ------------------------------------------ female
+    parGrowthFemale =  par %>%
+      mutate(kOpt = kOpt + par$kappaKOpt,
+             sigmaDeltaLVonBert = par$sigmaDeltaLVonBert)
+    
+    ogivePredictFemale <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern %>% filter(metapop == metapopulation), 
+                                                                          Nind = Nind, 
+                                                                          growPar = parGrowthFemale, 
+                                                                          RNGseed = RNGseed, ...) %>% 
+      filter(season == seasonSelected) %>% 
+      computeOgive(lengthAtMaturity = par$lFirstMaturityForFemale) 
+    
+    SSEfemale = ogivesObs %>% 
+      select(age, p = p_female) %>% 
+      inner_join(ogivePredictFemale, 
+                 by = 'age', suffix = c('.obs','.pred')) %>% 
+      summarise(SSE = sum((p.obs - p.pred)^2),) %>% 
+      unlist(use.names = FALSE)
+    
+    # ------------------------------------- male
+    parGrowthMale =  parGrowthFemale %>%
+      mutate(kOpt = kOpt - par$kappaKOpt)
+    
+    ogivePredictMale <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern %>% filter(metapop == metapopulation), 
+                                                                        Nind = Nind, 
+                                                                        growPar = parGrowthMale, 
+                                                                        RNGseed = RNGseed, ...) %>% 
+      filter(season == seasonSelected) %>% 
+      computeOgive(lengthAtMaturity = par$lFirstMaturityForMale) 
+    
+    SSEmale = ogivesObs %>% 
+      select(age, p = p_male) %>% 
+      inner_join(ogivePredictMale, 
+                 by = 'age', suffix = c('.obs','.pred')) %>% 
+      summarise(SSE = sum((p.obs - p.pred)^2),) %>% 
+      unlist(use.names = FALSE)
+    
+    # ------------------------------------------ total SSE 
+    SSE_all_metapop[i] = SSEfemale + SSEmale
+    i = i + 1
+  }
+  return(sum(SSE_all_metapop))
+}
+
+# ------------------------------------ test
+tic()
+SSEOgives(ogivePar0, fixedPar = ogivefixedPar, regional_metapop,  temperaturePattern, 
+          seasonSelected = 'spring',
+          Nind = 5000,
+          RNGseed = 1, 
+          effectFunction = modifiedLehman)
+toc()
+
+# --------------------------------- calibration
+
+ogivefixedPar = fullParOptimGrowth %>% select(-sigmaDeltaLVonBert) %>% 
+  mutate(kappaKOpt = 0)
+ogivePar0 =  c(lFirstMaturityForFemale = 43, lFirstMaturityForMale = 38, sigmaDeltaLVonBert = 0.2 )
+
+optimOgives =  optim(par = ogivePar0,
+                     fn = SSEOgives,
+                     # method = "L-BFGS-B", 
+                     control = list(trace = 1,
+                                    maxit = 1000),
+                     #hessian = TRUE,
+                     fixedPar = ogivefixedPar, 
+                     # lower = lower,
+                     # upper = upper, 
+                     regional_metapop = regional_metapop,
+                     temperaturePattern = temperaturePattern, 
+                     seasonSelected = "spring",
+                     Nind = 5000,
+                     RNGseed = 1, 
+                     effectFunction = modifiedLehman)
+
+
+fullParOptim = optimOgives$par %>% enframe() %>% pivot_wider() %>% 
+  bind_cols(ogivefixedPar)
+
+# ----------------------------------------------------- graph
+# ------------------- growth per gender
+growthOptimForAllMetapop = list()
+for (metapopulation in regional_metapop) {
+  
+  parGrowRefOneMetapop = parGrowRef %>% filter(metapop == metapopulation)
+  # simulated growth
+  simulatedOneMetapopFemale <-  computeMultipleLengthTrajectoriesWithRandomSeed( 
+    temperaturePattern = temperaturePattern %>%  filter(metapop == metapopulation), 
+    Nind = 1, 
+    growPar = fullParOptim %>% mutate(kOpt = kOpt + kappaKOpt), 
+    RNGseed = 1,
+    effectFunction = modifiedLehman) %>% 
+    select(age, L) %>% mutate(metapop = metapopulation, gender = "female", type = "simulated") 
+  
+  simulatedOneMetapopMale <-  computeMultipleLengthTrajectoriesWithRandomSeed( 
+    temperaturePattern = temperaturePattern %>%  filter(metapop == metapopulation), 
+    Nind = 1, 
+    growPar = fullParOptim %>% mutate(kOpt = kOpt - kappaKOpt), 
+    RNGseed = 1,
+    effectFunction = modifiedLehman) %>% 
+    select(age, L) %>% mutate(metapop = metapopulation, gender = "male", type = "simulated") 
+  
+  # observed growth
+  parGrowRefOneMetapop = parGrowRef %>% filter(metapop == metapopulation)
+  observedOneMetapop <- simulatedOneMetapop %>% select(age) %>% 
+    mutate(L = vonBertalanffyGrowth(age, L0 = parGrowRefOneMetapop$L0, Linf = parGrowRefOneMetapop$Linf, K = parGrowRefOneMetapop$K),
+           metapop = metapopulation, gender = "both", type = "Stich")
+  
+  growthOptimForAllMetapop[[metapopulation]] = bind_rows(simulatedOneMetapopFemale, simulatedOneMetapopMale,  observedOneMetapop )
+}
+
+
+rbindlist(growthOptimForAllMetapop) %>% 
+  # filter(type == "Stich") %>% 
+  ggplot() + 
+  geom_path(aes(x = age, y = L, col = metapop, linetype = gender), show.legend = TRUE) +
+  facet_grid(~metapop)
+
+
+# ------------------- ogives per gender
+predictedOgivesForAllMetapop = list()
+for (metapopulation in regional_metapop) {
+  # ------------- ogive for female
+  parGrowthFemale =  fullParOptim %>%
+    mutate(kOpt  = kOpt + kappaKOpt)
+  
+  ogivePredictFemale <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern %>% filter(metapop == metapopulation), 
+                                                                        Nind = 5000, 
+                                                                        growPar = parGrowthFemale, 
+                                                                        RNGseed = 1, 
+                                                                        effectFunction = modifiedLehman) %>% 
+    filter(season == 'spring') %>% 
+    computeOgive(lengthAtMaturity = parGrowthFemale$lFirstMaturityForFemale) %>% 
+    mutate(metapop = metapopulation, gender = "female", type = "simulated") 
+
+  # ------------- ogive for male
+  parGrowthMale =  fullParOptim %>%
+    mutate(kOpt  = kOpt - kappaKOpt)
+  
+  ogivePredictMale <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern %>% filter(metapop == metapopulation), 
+                                                                      Nind = 5000, 
+                                                                      growPar = parGrowthMale, 
+                                                                      RNGseed = 1, 
+                                                                      effectFunction = modifiedLehman) %>% 
+    filter(season == 'spring') %>% 
+    computeOgive(lengthAtMaturity = parGrowthFemale$lFirstMaturityForMale) %>% 
+    mutate(metapop = metapopulation, gender = "male", type = "simulated") 
+  
+  # ------------observed ogives and save in output
+   predictedOgivesForAllMetapop[[metapopulation]] <- ogivesObs %>% pivot_longer(cols = c(p_female,p_male), names_to = 'gender', values_to = 'p') %>% 
+    mutate(gender = str_remove(gender, 'p_'),
+           type = 'observed') %>% 
+    bind_rows(ogivePredictFemale, ogivePredictMale)
+}
+
+rbindlist(predictedOgivesForAllMetapop) %>%
+    ggplot() + 
+  geom_line(aes(x = age, y = p, col = metapop, linetype = type), show.legend = FALSE) +
+  facet_grid(gender ~ metapop)
+  
+rbindlist(predictedOgivesForAllMetapop) %>%
+  filter(type == 'observed') %>% 
+  ggplot() +
+  geom_line(aes(x = age, y = p, col = gender), show.legend = TRUE) +
+    facet_grid(~metapop)
+
+
+ogivesObs %>% pivot_longer(cols = c(p_female,p_male), names_to = 'gender', values_to = 'p') %>% 
+    mutate(gender = str_remove(gender, 'p_')) %>% 
+  ggplot() +
+  geom_line(aes(x = age, y = p, col = gender), show.legend = TRUE) +
+    facet_grid(~metapop)
+
 ```
 
+