Commit 87b7bb9e authored by patrick.lambert's avatar patrick.lambert
Browse files

with a new computeOgive

parent 10525849
......@@ -39,6 +39,7 @@ vonBertalanffyInverse = function(L, L0, Linf, K){
# von Bertalanffy increment
# pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche
vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){
tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26)
L = matrix(nrow = nStep + 1)
......@@ -156,8 +157,24 @@ vonBertalanffyWithRandomVector = function(L, Linf, K, timeStepDuration, randomVe
return(L)
}
computeOgive = function(lengthTrajectories, lengthAtMaturity){
ogive <- lengthTrajectories %>%
group_by(age) %>%
summarise(nTotal = n()) %>%
left_join(lengthTrajectories %>% group_by(age) %>%
filter(L >= lengthAtMaturity) %>%
summarise(taller = n(), .groups = 'drop'),
by = 'age') %>%
replace_na(list(taller = 0)) %>%
mutate(mature = c(0,diff(taller)),
immature = nTotal - cumsum(mature),
p = if_else(mature + immature > 0 , mature / (mature + immature), 1)) %>%
select(age, p)
return(ogive)
}
computeOgive = function(lengthTrajectories, lengthAtMaturity) {
# deprecated
computeOgive3 = function(lengthTrajectories, lengthAtMaturity) {
ogive <-
lengthTrajectories %>%
group_by(age) %>%
......
......@@ -67,130 +67,131 @@ temperaturePattern <- growthInBasin %>%
```
```{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,
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)
}
#Deprecated
computeOgive2 = 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)
}
# TRENSFERED IN GR3D functiond
# # 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)
# }
#
# #Deprecated
# computeOgive2 = 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)
# }
```
......@@ -632,52 +633,64 @@ SSEforOneMetapop(fullPar = enframe(c(par_0, fixedPar)) %>% pivot_wider(),
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))
#
# fixedPar
# res_Abis$par
# # 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(res_Abis$par, fixedPar)) %>% 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 <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern = temperaturePattern %>% filter(metapop == "semelparous"),
# Nind = 100000,
# growPar = parFemale %>% mutate(sigmaDeltaLVonBert =.6),
# RNGseed = 1) %>%
# arrange(ind,age) %>%
# filter(season == "spring")
#
#
# 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"}
......@@ -721,14 +734,15 @@ 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
lFirstMaturityForFemale = 45,#XML - Alosa sapidissima
linfVonBert = 76, #XML - Alosa alosa
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
linfVonBert = 76, #XML - Alosa alosa
sigmaDeltaLVonBert = .2) #XML - Alosa alosa
lengthAtHatching = 2.8)# XML - Alosa sapidissima
#optimisation
# XML - Alosa sapidissia)
......
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