Skip to content
GitLab
Commit bcc08e41 authored by patrick.lambert's avatar patrick.lambert
Browse files

new calibration of grow

parent e3b5329a
exploration_GR3D_process SpawnerRunAnalysis develop
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 659 additions and 77 deletions
+659 -77
exploration/GR3D_Rdescription/GR3Dfunction.R
+ 107
− 1
View file @ bcc08e41
......@@ -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)
......
exploration/NEA_calibration_offline/maturationPL.Rmd
+ 552
− 76
View file @ bcc08e41
......@@ -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))
```
Supports Markdown
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