diff --git a/exploration/NEA_calibration_offline/maturation.Rmd b/exploration/NEA_calibration_offline/maturation.Rmd
index e8cdb12cbd838bd173a70b81a6a793d499ee2b9c..87ca4b7389e5e717765b558fa74af83b3a536235 100644
--- a/exploration/NEA_calibration_offline/maturation.Rmd
+++ b/exploration/NEA_calibration_offline/maturation.Rmd
@@ -28,10 +28,19 @@ 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(.)))
-growPar <- Stich2020_sel %>% filter(catchment == "semelparous") %>%
+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,
@@ -49,6 +58,10 @@ temperaturePattern <- growthInBasin %>%
# )) %>%
select(age, season,temperature)
+```
+
+```{r local function}
+
# generate a cohort of length trajectories
computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar){
ages = temperaturePattern$age
@@ -80,19 +93,64 @@ computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, grow
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)
+}
-tic()
-computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growPar) %>%
- filter(season == 'spring')
-toc()
-
+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(yearlyLengthTrajectories, lengthAtMaturity) {
- ogive <- yearlyLengthTrajectories %>%
+computeOgive = function(lengthTrajectories, lengthAtMaturity) {
+ ogive <- lengthTrajectories %>%
distinct(age) %>% arrange(age) %>%
- left_join(yearlyLengthTrajectories %>% filter(L > lengthAtMaturity) %>%
+ left_join(lengthTrajectories %>% filter(L > lengthAtMaturity) %>%
group_by(ind) %>%
summarise(age = min(age), .groups = 'drop') %>%
group_by(age) %>%
@@ -100,14 +158,26 @@ computeOgive = function(yearlyLengthTrajectories, lengthAtMaturity) {
by = 'age') %>%
replace_na(list(mature = 0)) %>%
mutate(immature = sum(mature) - cumsum(mature),
- p = if_else(mature + immature == 0, 1, mature/(mature + immature))) %>%
+ 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) %>%
@@ -116,21 +186,18 @@ computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
```
-```{r ogive varaibility}
-
-
-
+```{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 envelop}
-# envelop of growth curves
-computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growPar) %>%
+```{r growth curve envelope}
+# envelope of growth curves
+computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growParStich) %>%
group_by(age) %>%
summarise(Lmin = min(L),
L025 = quantile(L, .25),
@@ -149,11 +216,11 @@ computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growPar) %>%
# verify that all L < Linf
-computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growPar) %>%
- filter(L>growPar$Linf)
+computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growParStich) %>%
+ filter(L>growParStich$Linf)
```
-```{r}
+```{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) %>%
@@ -161,18 +228,26 @@ lengthTrajectoriesInSpring = computeMultipleLengthTrajectories(temperaturePatter
# calibrate only lengthAtMaturity
objFn_A = function(par, ogiveObs, yearlyLengthTrajectories){
- SSE = ogiveObs %>%
+ 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)) %>%
+ 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 = 35),
+objFn_A(par = c(lengthAtMaturity = 40),
ogiveObs,
yearlyLengthTrajectories = lengthTrajectoriesInSpring)
@@ -191,12 +266,12 @@ computeOgive(yearlyLengthTrajectories = lengthTrajectoriesInSpring,
# ----------------------------------------------------------------------------------------------------
# calibrate lengthAtMaturity and Kopt
-objFn_B = function(par, growPar, ogiveObs, temperaturePattern, Nind = 100000, season = 'spring', sigmaDeltaLVonBert = .2){
+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 == season) %>%
+ filter(season == seasonSelect) %>%
computeOgive(lengthAtMaturity = par[1])
SSE = ogiveObs %>%
@@ -252,4 +327,129 @@ computeMultipleLengthTrajectories(temperaturePattern,
geom_line(aes(y = Lmed)) +
labs(x = 'Age (y)', y = 'Length (cm)') +
geom_abline(slope =0, intercept = 44)
-```
\ No newline at end of file
+```
+
+```{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))
+
+```