Commit 7706ae6e authored by Poulet Camille's avatar Poulet Camille
Browse files

Maturation and growth calibration for all metapop

parent e099100a
......@@ -8,11 +8,11 @@ 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
......@@ -20,28 +20,27 @@ library(tictoc) # see timeR
```
```{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')
```
```{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,
sigmaDeltaLVonBert = 0.4)
......@@ -49,14 +48,14 @@ growParStich <- Stich2020_sel %>% filter(catchment == "semelparous") %>%
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)
```
......@@ -98,7 +97,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))) %>%
......@@ -133,6 +134,9 @@ computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
}
vonBertalanffyWithRandomVector = function(L, L0, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
if (sigma == 0) {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
......@@ -167,7 +171,7 @@ computeOgive = function(lengthTrajectories, lengthAtMaturity) {
```
```{r test}
```{r test }
#Do not run
tic()
......@@ -222,7 +226,7 @@ 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))
......@@ -333,11 +337,15 @@ computeMultipleLengthTrajectories(temperaturePattern,
```
```{r calibration K and lengthAtMaturity by gender with ogive, warning = "FALSE"}
#par = res_D$par
#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() #full set of parameters
......@@ -363,7 +371,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 - avoid immature individuals at maximuma age
# 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) %>%
......@@ -392,74 +400,302 @@ objFn_D = function(par, fixedPar, ogivesObs, parGrowRef, temperaturePattern,
### SSE for growth curve
#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) %>%
femaleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern ,
Nind = 1,
growPar = parFemale %>% mutate(sigmaDeltaLVonBert=0) ,
RNGseed = RNGseed) %>%
filter(season == seasonSelected) %>%
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 %>%
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
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
}
return((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
```
```{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) {
#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)
# 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')) %>%
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)
### 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)
#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)
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
return(SSEgrowthAndOgives)
}
#TODO: ajuster les bonnes valeurs par rapport aux ogives de la literature
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))
#Fonction to be optimzed over
SSEforAllMetapop = function(par, fixedPar, regional_metapop,ogiveObs, 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)
}
#TODO verifier si ce sont bien les parametres de Stich ou les valeurs ajustés sur les paramètres de Stich expeté sur le Top ? A retravailler pour voir la manière dont on peu enchainer les processus.
#Quelle source de paramètre choisit-on ?
parGrowRef = growParStich %>% select(L0 = lengthAtHatching, K, Linf)
#Verification
SSEforAllMetapop(regional_metapop = regional_metapop,
par = par_0,
fixedPar = fixedPar,
ogiveObs = ogiveObserved,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
Nind = 1000,
seasonSelected = "spring",
RNGseed = 1)
#TODO: modifer les paramètres fixés et les paramètres à optimiser à mon bon vouloir ahah
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 )
# computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern = temperaturePattern %>% filter(metapop == "semelparous"),
# Nind = 10,
# growPar = parFemale,
# RNGseed = 1)
# objFn_D(par = par,
# fixedPar = fixedPar,
# ogivesObs = ogivesObs,
# parGrowRef = parGrowRef,
# temperaturePattern = temperaturePattern,
# Nind = 5000)
res_D <- optim(par = par,
fn = objFn_D,
```
```{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
ogivesObserved = 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
linfVonBert = 76, #XML - Alosa alosa
lFirstMaturityForMale = 40, #XML - Alosa sapidissima
lFirstMaturityForFemale = 45,#XML - Alosa sapidissima
sigmaDeltaLVonBert = .2) #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
#optimisation
res_A <- optim(par = par_0,
fn = SSEforAllMetapop,
method = "BFGS",
control = list(trace = 1,
maxit = 1000),
hessian = TRUE,
fixedPar = fixedPar,
ogivesObs = ogivesObs,
regional_metapop = regional_metapop,
ogivesObs = ogivesObserved,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
Nind = 5000,
RNGseed =1,
control = list(trace = 1,
maxit = 1000))
RNGseed =1)
#Case 2
par_1 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 76,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45,
sigmaDeltaLVonBert = .2,
lengthAtHatching = 2.8)
fixedPar_1 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9)
#Case 3
par_2 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 76,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45)
fixedPar_2 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9,
sigmaDeltaLVonBert = .2,
lengthAtHatching = 2.8)
```
```{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))
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment