Commit f45992ea authored by Poulet Camille's avatar Poulet Camille
Browse files

Merge branch 'develop' into exploration_GR3D_process

parents 57624722 de111f56
......@@ -73,12 +73,12 @@ connections_trans <- connections %>%
##Select the temperature experienced by fish during the corresponding season
#in winter
wint<-connections_trans %>%
wint <- connections_trans %>%
filter(season %in% c("winter","spring"), offshore_basin_type == "wintering") %>%
select(inshore_basin_name,offshore_basin_name,year,season,temperature)
#in summer
sum<-connections_trans %>%
sum <- connections_trans %>%
filter(season %in% c("summer","fall"), offshore_basin_type == "summering") %>%
select(inshore_basin_name,offshore_basin_name,year,season,temperature)
......@@ -87,7 +87,7 @@ tempEffectInOffshore <- bind_rows(wint,sum) %>%
ungroup() %>%
arrange(inshore_basin_name, year,season) %>%
group_by(inshore_basin_name,offshore_basin_name,season) %>%
summarise(Tmean = mean(temperature))
summarise(temperature_O = mean(temperature), .groups = 'drop')
#mutate(temperatureEffect = temperatureEffect(Tmean, growPar$tempMinGrow, growPar$tempOptGrow, growPar$tempMaxGrow))
#delete intermediate files
......@@ -118,24 +118,24 @@ tempEfffectInRiver <- tempInriver %>%
#Optimzation
tempRiver <- tempEfffectInRiver %>%
filter(season %in% c("spring","summer") & zone == "river_temperature")
filter(season %in% c("spring","summer") & zone == "river_temperature")
tempInshore <- tempEfffectInRiver %>%
filter(season %in% c("fall","winter") & zone == "inshore_temperature")
filter(season %in% c("fall","winter") & zone == "inshore_temperature")
tempEffectInFreshWat <- bind_rows(tempRiver,tempInshore) %>%
arrange(basin_name,year,season) %>%
select(basin_name,year, season, temperature) %>%
group_by(basin_name,season) %>%
summarize(Tmean_RIB = mean(temperature))
summarize(temperature_RI = mean(temperature), .groups = 'drop')
# tempEfffectInRiver %>%
# {bind_rows(
# filter(season == "summer" & zone == "river_temperature"),
# filter(season == "fall" & zone == "inshore_temperature"))
# }
```
......@@ -168,7 +168,7 @@ To optimize the growth parameters included in GR3D for the Amercian shad, we use
vonBertalanffyLengthAtHatch = function(Linf, K, t0){
L0 = Linf * (1-exp(t0*K))
L0 = Linf * (1 - exp(t0*K))
return(L0)
}
......@@ -181,9 +181,7 @@ Stich2020_sel <- read_csv("../NEA_calibration_offline/Stich_Table 9.csv") %>%
dplyr::rename("catchment"="cachtment") %>%
select(parameter,catchment,mean) %>%
pivot_wider(names_from = parameter, values_from = mean) %>%
mutate(L0_theo = vonBertalanffyLengthAtHatch(Linf,K,t0)) %>%
mutate_at(vars(c(Linf,L0_theo)), funs(./10)) #from mm to cm
mutate(L0_theo = vonBertalanffyLengthAtHatch(Linf, K, t0))
#create a data frame and compute L at age with values from Stich
growthInOffshore <- expand_grid(basin_name = nea_riverBasinFeatures$basin_name , age = seq(0,9.75,.25)) %>%
......@@ -200,24 +198,66 @@ growthInOffshore <- expand_grid(basin_name = nea_riverBasinFeatures$basin_name ,
lat_outlet > 41.28793 ~ "northern iteroparous")) %>%
inner_join(Stich2020_sel, by = c("metapop" = "catchment")) %>%
group_by(basin_name,metapop) %>%
mutate(LStich = vonBertalanffyGrowth(age, L0_theo, Linf, K)) %>%
mutate(LStich = vonBertalanffyGrowth(age, L0_theo, Linf, K) / 10) %>% # in cm
ungroup()
growthInOffshore <-growthInOffshore %>%
inner_join(tempEffectInOffshore, by = c("basin_name" ="inshore_basin_name",
"season")) %>%
right_join(tempEffectInFreshWat, by = c("basin_name","season"))
growthInOffshore<-growthInOffshore %>%
mutate(temperature = ifelse (season == "spring" & age == 0, Tmean_RIB,
ifelse(season == "summer" & age == 0.25, Tmean_RIB,
ifelse(season == "fall" & age == 0.50, Tmean_RIB, Tmean)))) %>%
inner_join(tempEffectInOffshore,
by = c("basin_name" = "inshore_basin_name", "season")) %>%
right_join(tempEffectInFreshWat,
by = c("basin_name","season"))
growthInOffshore <- growthInOffshore %>%
mutate(temperature = ifelse(season == "spring" & age == 0, temperature_RI,
ifelse(season == "summer" & age == 0.25, temperature_RI,
ifelse(season == "fall" & age == 0.50, temperature_RI, temperature_O)))) %>%
arrange(basin_name,age,season) %>%
select(basin_name,offshore_basin_name,metapop,age, season,temperature,LStich)
select(basin_name, offshore_basin_name, metapop,age, season, temperature_O, LStich)
#group_by(metapop,age,season) %>%
#summarize(LStich = unique(LStich))
# NEW =====
# add presence and reorder level for metapop and season
growthInBasin <- growthInBasin %>%
inner_join(nea_presence %>% select(-basin_id), by ='basin_name') %>%
mutate(metapop = factor(metapop, levels = c('northern iteroparous', 'southern iteroparous', 'semelparous')),
season = factor(season, levels = c('winter', 'spring', 'summer', 'fall')))
```
There is 1 river ( North (MA)) in Stich et al paper that is not referenced in our database.
```{r join between Stich and GR3D database}
# NEW =====
Stich2020_sel <- Stich2020_sel %>% mutate(basin_name = catchment) %>%
mutate(basin_name = replace(basin_name, basin_name == "Tar-Pamlico", "Pamlico-Tar"),
basin_name = replace(basin_name, basin_name == "Albemarle", "Roanoke"),
basin_name = replace(basin_name, basin_name == "Upper Chesapeake Bay", "Susquehanna"))
Stich2020_sel %>% filter(!catchment %in% c('coastwide', 'northern iteroparous', 'semelparous', 'southern iteroparous' )) %>%
select(catchment, basin_name) %>%
left_join(growthInBasin %>% distinct(basin_name, metapop), by = c("basin_name") ) %>%
arrange(metapop) %>%
flextable() %>%
autofit()
```
```{r range of temperature}
# NEW =====
growthInBasin %>%
group_by(obs_1900_1950) %>%
summarise(Tmin = min(temperature_RIO), Tmax = max(temperature_RIO), .groups = 'drop')
growthInBasin %>%
inner_join(Stich2020_sel %>% select(basin_name), by ='basin_name') %>%
group_by(obs_1900_1950) %>%
summarise(Tmin = min(temperature_RIO), Tmax = max(temperature_RIO), .groups = 'drop')
growthInBasin %>%
filter(temperature_RIO < 1.6) %>%
distinct(basin_name, obs_1900_1950)
```
```{r Growth in GR3D with temperature effect, echo = FALSE, warning=FALSE}
......@@ -237,7 +277,19 @@ growParModifed <- growPar %>%
```
```{r growth optimisation, echo = FALSE, warning=FALSE}
```{r}
# new ====
save(growParModifed,growthInBasin, nea_presence, Stich2020_sel, file = 'SOS.rdata' )
rm(list = ls())
load('SOS.rdata')
source("../GR3D_Rdescription/GR3Dfunction.R")
```
```{r growth optimisation function Camille, echo = FALSE, warning=FALSE}
#vector on paramters to get optimized
#par[1] = Tmin
#par[2] = Topt
......@@ -246,12 +298,6 @@ growParModifed <- growPar %>%
#par[5] = Linf
#par[6] = k
vecPar = growParModifed %>%
select(tempMinGrow,tempOptGrow,tempMaxGrow,lengthAtHatching,Linf, kOpt) %>%
pivot_longer(tempMinGrow:kOpt, names_to = "parameter",values_to = "value")
vecPar = as.vector(vecPar$value)
#Return data
computeGrowOpt = function(par,data){
data <- data %>% arrange(age, basin_name) %>%
......@@ -294,8 +340,8 @@ computeGrowOpt = function(par,data){
#Function to be optimized
computeSSE = function(par,data){
data = computeGrowOpt(par,data)
data = computeGrowOpt(par,data)
SSE = data %>%
group_by(metapop) %>%
mutate(W = 1/n()) %>%
......@@ -313,34 +359,281 @@ computeSSE = function(par,data){
#computeSSE (vecPar, growthInOffshore)
```
```{r growth optimisation function Patrick, echo = FALSE, warning=FALSE}
computeGrowAllBasins = function(data, growPar){
data <- data %>% arrange(age, basin_name) %>%
mutate(L = if_else(age == 0, growPar$lengthAtHatching, 0),
temperatureEffect = temperatureEffect(temperature, growPar$tempMinGrow, growPar$tempOptGrow, growPar$tempMaxGrow))
ages <- data %>% distinct(age) %>% arrange(age) %>% unlist(use.names = FALSE)
for (i in 2:length(ages)) {
previousAge = ages[i - 1]
currentAge = ages[i]
tempEffect = data %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE)
previousL <- data %>% 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,
sigma = growPar$sigmaDeltaLVonBert,
tempEffect = tempEffect)
data = data %>% mutate(L =replace(L, age == ages[i], currentL) )
}
return(data)
}
par2parGrowth = function(par){
parGrowth = data.frame(tempMinGrow = unname(par['tempOptGrow'] - par['epsilonMinus']),
tempOptGrow = unname(par['tempOptGrow']),
tempMaxGrow = unname(par['tempOptGrow'] + par['epsilonPlus']),
lengthAtHatching = unname(par['lengthAtHatching']),
Linf = unname(par['Linf']),
kOpt = unname(par['kOpt']),
sigmaDeltaLVonBert = 0 )
return(parGrowth)
}
# avec Topt, Tmin = Topt - espsilonMinus, Tmax = Topt + epsilonPlus
objFn = function(par, data) {
growPar = par2parGrowth(par)
data = computeGrowAllBasins(data, growPar)
SSE <- data %>%
# fit only on one-season L to avoid smoothing of seasonal growth
# filter(season == 'spring') %>%
# mutate(age_year = floor(age)) %>%
# group_by(basin_name, metapop, age_year) %>%
# summarise( L = mean(L), LStich = mean(LStich), .groups = 'drop') %>%
# on pondère chaque basin versant pour garder des points équivalents entre les metapop
group_by(metapop) %>% mutate(W = 1/n()) %>% ungroup() %>%
transmute(squareError = W * ((L - LStich)^2)) %>% summarise(sum(squareError)) %>% unlist(use.names = FALSE)
return(SSE)
}
tic()
# équivalent de la fonction computeSSE (un peu plus rapide)
objFn_B = function(par, data) {
growPar = enframe(par) %>%
pivot_wider() %>% mutate(sigmaDeltaLVonBert = 0)
data = computeGrowAllBasins(data, growPar)
SSE <- data %>%
# fit only on one-season L to avoid smoothing of seasonal growth
# filter(season == 'spring') %>%
# mutate(age_year = floor(age)) %>%
# group_by(basin_name, metapop, age_year) %>%
# summarise( L = mean(L), LStich = mean(LStich), .groups = 'drop') %>%
# on pondère chaque basin versant pour garder des points équivalents entre les metapop
group_by(metapop) %>% mutate(W = 1/n()) %>% ungroup() %>%
transmute(squareError = W * ((L - LStich)^2)) %>% summarise(sum(squareError)) %>% unlist(use.names = FALSE)
return(SSE)
}
# équivalent de la fonction objFn_B mais avec la possibilité de fixer des parametres
objFn_C = function(par, data, fixedPar) {
growPar = enframe(c(par, fixedPar)) %>% pivot_wider() %>%
mutate(sigmaDeltaLVonBert = 0)
data = computeGrowAllBasins(data, growPar)
SSE <- data %>%
# fit only on one-season L to avoid smoothing of seasonal growth
# filter(season == 'spring') %>%
# mutate(age_year = floor(age)) %>%
# group_by(basin_name, metapop, age_year) %>%
# summarise( L = mean(L), LStich = mean(LStich), .groups = 'drop') %>%
# on pondère chaque basin versant pour garder des points équivalents entre les metapop
group_by(metapop) %>% mutate(W = 1/n()) %>% ungroup() %>%
transmute(squareError = W * ((L - LStich)^2)) %>% summarise(sum(squareError)) %>% unlist(use.names = FALSE)
return(SSE)
}
```
```{r starting parameters for optimisation}
# starting point from XML
vecPar = growParModifed %>%
select(tempMinGrow, tempOptGrow, tempMaxGrow, lengthAtHatching, Linf, kOpt) %>%
pivot_longer(tempMinGrow:kOpt, names_to = "parameter", values_to = "value")
vecPar = as.vector(vecPar$value)
par0 = c(tempOptGrow = vecPar[2],
epsilonMinus = vecPar[2] - vecPar[1],
epsilonPlus = vecPar[3] - vecPar[2],
lengthAtHatching = vecPar[4],
Linf = vecPar[5],
kOpt = vecPar[6])
par0_B = c(tempMinGrow = vecPar[1],
tempOptGrow = vecPar[2],
tempMaxGrow = vecPar[3],
lengthAtHatching = vecPar[4],
Linf = vecPar[5],
kOpt = vecPar[6])
par0_C = par0_B[c('tempOptGrow', 'lengthAtHatching', 'Linf', 'kOpt')]
#fixedPar_C = par0_B[c('tempMinGrow', 'tempMaxGrow')]
fixedPar_C = c(tempMinGrow = 1.6, tempMaxGrow = 27.9)
# # another starting point
# par0 = c(tempOptGrow = 7,
# epsilonMinus = 5,
# epsilonPlus = 10,
# lengthAtHatching = 5,
# Linf = 45,
# kOpt = .5)
#
# par0_B = par2parGrowth(par0) %>%
# select(-sigmaDeltaLVonBert) %>%
# unlist()
#
# par0_C = par0_B[c('tempOptGrow', 'lengthAtHatching', 'Linf', 'kOpt')]
# #fixedPar_C = par0_B[c('tempMinGrow', 'tempMaxGrow')]
# fixedPar_C = c(tempMinGrow = 1.6, tempMaxGrow = 27.9)
#
# vecPar = par2parGrowth(par0) %>%
# select(-sigmaDeltaLVonBert) %>%
# pivot_longer(tempMinGrow:kOpt, names_to = "parameter",values_to = "value")
# vecPar = as.vector(vecPar$value)
```
```{r data for calibration}
dataCalibration = growthInBasin %>%
inner_join(Stich2020_sel %>% select(basin_name), by ='basin_name') %>%
mutate(temperature = temperature_RIO)
```
Avec uniqument le jeu de Stich cela part en vrille avec objFn_B et computeSSE ( sans donner les mémes resultats argh) avec Topt < Tmin. On a des résulats corrects avec objFn ( avec Topt +- epsilon). Mais comme on va utiliser l'option en fixant Tmin et Tmax
```{r optimisation run}
tic('objFn')
# avec Topt, Tmin = Topt - espsilonMinus, Tmax = Topt + epsilonPlus
res <- optim(par = par0,
fn = objFn,
data = dataCalibration,
# method = "L-BFGS-B",
# lower = c(2, 2, 2, 2, 30, 0.2),
# upper = c(20, 20, 20, 20, 100, 2),
control = list(trace = 0,
#parscale = c(1,1,1,1,1,0.1),
maxit = 1000))
toc()
tic('objFn_B')
# équivalent de la fonction computeSSE (un peu plus rapide)
res_B <- optim(par = par0_B,
fn = objFn_B,
data = dataCalibration,
# method = "L-BFGS-B",
# lower = c(0,2,10,2,30,0.2),
# upper = c(5,10,30,20,100,2),
control = list(trace = 0,
maxit = 1000))
toc()
tic('objFn_C')
# équivalent de la fonction objFn_B mais avec la possibilité de fixer des parametres
res_C <- optim(par = par0_C,
fn = objFn_C,
data = dataCalibration,
fixedPar = fixedPar_C,
# method = "L-BFGS-B",
# lower = c(0,2,10,2,30,0.2),
# upper = c(5,10,30,20,100,2),
control = list(trace = 0,
maxit = 1000))
toc()
tic('computeSSE')
growth_optimal <- optim(c(Tmin = vecPar[1],
Topt = vecPar[2],
Tmax = vecPar[3],
L0 = vecPar[4],
Linf = vecPar[5],
K = vecPar[6]),
computeSSE,
data = growthInOffshore,
lower = c(0,0,0,0,0,0),
method = "L-BFGS-B",
control = list(trace = TRUE, parscale = c(1,1,1,0.1,1,0.01)))
fn = computeSSE,
data = dataCalibration,
# lower = c(0,0,0,0,0,0),
# method = "L-BFGS-B",
control = list(trace = 0,
#parscale = c(1,1,1,0.1,1,0.01),
maxit = 10000))
toc()
# results
par2parGrowth(res$par)
res_B$par
c(res_C$par, fixedPar_C)
growth_optimal$par
res$value
res_B$value
res_C$value
growth_optimal$value
#Optimized parameters
growParOptim <- data.frame(value = growth_optimal$par) %>%
rownames_to_column("parameter")
#Rdata
save(growth_optimal, file = "growth_optimal.RData")
save(growth_optimal, res, res_B, res_C, file = "growth_optimal.RData")
#csv
write_csv(growParOptim, "growParOptim.csv")
```
```{r graph to present results}
# NEW
# ===============================================================================
#parGrowthUsed = c(res_B$par) %>% enframe() %>% pivot_wider()
parGrowthUsed = c(res_C$par, fixedPar_C) %>% enframe() %>% pivot_wider()
# ================================================================================
dataVerif <- computeGrowAllBasins(data = dataCalibration,
growPar = parGrowthUsed %>% mutate(sigmaDeltaLVonBert = 0) )
# graph by meta pop
dataVerif %>% group_by(metapop, age) %>%
summarise(L = mean(L), LStich = mean(LStich), .groups = 'drop') %>%
mutate(metapop = factor(metapop, levels = c('northern iteroparous', 'southern iteroparous', 'semelparous'))) %>%
ggplot(aes(x = age, color = metapop)) +
geom_line(aes(y = LStich)) +
geom_line(aes(y = L))
# thermal response curve with average seasonal temperature experienced by fish
dataTemperatureMetapop <- dataVerif %>%
group_by(metapop, season) %>%
summarise(temperature = mean(temperature), .groups = 'drop') %>%
mutate(tempEffect = temperatureEffect(temperature, parGrowthUsed$tempMinGrow, parGrowthUsed$tempOptGrow, parGrowthUsed$tempMaxGrow ),
metapop = factor(metapop, levels = c('northern iteroparous', 'southern iteroparous', 'semelparous')),
season = factor(season, levels = c('winter', 'spring', 'summer', 'fall')))
data.frame(temperature = seq(0,40,0.1)) %>%
mutate(tempEffect = temperatureEffect(temperature, parGrowthUsed$tempMinGrow, parGrowthUsed$tempOptGrow, parGrowthUsed$tempMaxGrow )) %>%
ggplot() + ylim(0,1)+
geom_line(aes(x = temperature, y = tempEffect)) +
geom_point(data = dataTemperatureMetapop, aes(x = temperature, y = tempEffect, color = metapop, shape = season ), size=3)
```
```{r draw growth curves}
load("growth_optimal.RData")
......
parameter,value
<<<<<<< HEAD
Tmin,0
Topt,4.5125688473103995
Tmax,76.23141515590089
L0,6.5184840221819975
Linf,47.53851690994118
K,0.4964834345987972
=======
Tmin,-12.713186580960963
Topt,-35.999892199799476
Tmax,77.91401250769903
L0,6.437410947346945
Linf,47.55230062211272
K,1.2212313819420748
>>>>>>> develop
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