diff --git a/exploration/GR3D_Rdescription/GR3Dmetapopulation.R b/exploration/GR3D_Rdescription/GR3Dmetapopulation.R
new file mode 100644
index 0000000000000000000000000000000000000000..cd1dbe14251d668fe6a760618649a1ac750ab9d8
--- /dev/null
+++ b/exploration/GR3D_Rdescription/GR3Dmetapopulation.R
@@ -0,0 +1,101 @@
+
+#=============================================================================================================
+# metapopulation identification
+# =============================================================================================================
+computeAutochtonousRate = function(data){
+ data %>% filter(migrationBasin == originBasin) %>%
+ select(-originBasin) %>%
+ inner_join(data %>% group_by(year, migrationBasin) %>%
+ summarise(totalRun = sum(effective), .groups = 'drop'), by = c('migrationBasin', 'year')) %>%
+ mutate(autochtonousRate = effective / totalRun) %>%
+ select(year, migrationBasin, autochtonousRate)
+}
+
+computeHomingRate = function(data){
+ data %>% filter(migrationBasin == originBasin) %>%
+ select(-migrationBasin) %>%
+ inner_join(data %>% group_by(year, originBasin) %>%
+ summarise(production = sum(effective), .groups = 'drop'), by = c('originBasin', 'year')) %>%
+ mutate(homingRate = effective / production) %>%
+ select(year, originBasin, homingRate)
+}
+
+# comparison between autochtonous rate and homing rate
+# computeAutochtonousRate(exchangesData) %>%
+# inner_join(computeHomingRate(data = exchangesData), by = c('migrationBasin' = 'originBasin', 'year'))
+
+
+
+identyMetapopulation = function(exchangesData, threshold = .95, verbose = FALSE) {
+ # exchangesData {year, migrationBasin, originBasin)}
+
+ exchangesDataUpdated = exchangesData
+
+ # initialise the connection table between basin and metapopulation
+ metapopulation = exchangesDataUpdated %>%
+ distinct(year, basin = migrationBasin) %>%
+ mutate(metapop = basin)
+
+ iteration = 1
+
+ # find the basin with the minimum autochtonous rate for the first time,
+ basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
+ group_by(year) %>%
+ slice(which.min(autochtonousRate))
+
+ # loop while autochtonousRate is still < 0.95
+
+ while (basinWithMinAutochtonousRate$autochtonousRate < threshold ) {
+ # while (iteration <= 4 ) {
+ if (verbose) cat(iteration, ": ", basinWithMinAutochtonousRate$autochtonousRate, '\n' )
+
+ # basinWithMinAutochtonousRate will be merged with origin basin sending the maximum number of fish in this basin
+ metapopsToBeMerged <- exchangesDataUpdated %>%
+ inner_join(basinWithMinAutochtonousRate %>% select(-autochtonousRate),
+ by = c('year', 'migrationBasin' )) %>%
+ filter(migrationBasin != originBasin) %>% # to avoid self merging
+ group_by(year) %>%
+ slice(which.max(effective)) %>%
+ ungroup() %>%
+ select(year, migrationBasin, originBasin)
+
+ # update the connection table between basin and metapopulation by merging metapopulation
+ for (i in 1:nrow(metapopsToBeMerged)) {
+ mergedName = paste0("M", iteration, '_',i)
+ metapopsToBeMerged_i <- metapopsToBeMerged %>% slice(i)
+ if (verbose) cat("\t", metapopsToBeMerged_i$migrationBasin, ' + ', metapopsToBeMerged_i$originBasin, ' = ', mergedName, "\n")
+
+ metapopulation <- metapopulation %>%
+ mutate(metapop = if_else(year == metapopsToBeMerged_i$year &
+ metapop %in% (metapopsToBeMerged_i %>% select(ends_with("Basin")) %>% unlist(use.names = FALSE)),
+ mergedName,
+ metapop))
+ }
+ # sum by migration and origin basins according to updated metapopulations
+ # computed with initial exchangesData
+ exchangesDataUpdated <- exchangesData %>%
+ #sum up on migration basins
+ inner_join(metapopulation,
+ by = c('year', 'migrationBasin' = 'basin')) %>%
+ group_by(year, metapop, originBasin) %>%
+ summarise(effective = sum(effective), .groups = 'drop') %>%
+ rename(migrationBasin = metapop) %>%
+ # sum on origine basins
+ inner_join(metapopulation,
+ by = c('year', 'originBasin' = 'basin')) %>%
+ group_by(year, metapop, migrationBasin) %>%
+ summarise(effective = sum(effective), .groups = 'drop') %>%
+ rename(originBasin = metapop)
+
+ # minimum of autochtonous rate
+ basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
+ group_by(year) %>%
+ slice(which.min(autochtonousRate))
+
+ iteration = iteration + 1
+ }
+ return(list(metapopulation = metapopulation, exchangesData = exchangesDataUpdated))
+}
+
+
+
diff --git a/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R b/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R
index 61ddeb400ff91bf41b8badc9f4859fe62288aaee..03b1c22efb62721f7255a9082822393dade46e3e 100644
--- a/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R
+++ b/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R
@@ -1,102 +1,11 @@
-library(tidyverse)
-
-computeAutochtonousRate = function(data){
- data %>% filter(migrationBasin == originBasin) %>%
- select(-originBasin) %>%
- inner_join(data %>% group_by(year, migrationBasin) %>%
- summarise(totalRun = sum(effective), .groups = 'drop'), by = c('migrationBasin', 'year')) %>%
- mutate(autochtonousRate = effective / totalRun) %>%
- select(year, migrationBasin, autochtonousRate)
-}
-
-computeHomingRate = function(data){
- data %>% filter(migrationBasin == originBasin) %>%
- select(-migrationBasin) %>%
- inner_join(data %>% group_by(year, originBasin) %>%
- summarise(production = sum(effective), .groups = 'drop'), by = c('originBasin', 'year')) %>%
- mutate(homingRate = effective / production) %>%
- select(year, originBasin, homingRate)
-}
-
-# comparison between autochtonous rate and homing rate
-# computeAutochtonousRate(exchangesData) %>%
-# inner_join(computeHomingRate(data = exchangesData), by = c('migrationBasin' = 'originBasin', 'year'))
-
-
-
-identyMetapopulation = function(exchangesData, threshold = .95, verbose = FALSE) {
- exchangesDataUpdated = exchangesData
-
- # initialise the connection table between basin and metapopulation
- metapopulation = exchangesDataUpdated %>%
- distinct(year, basin = migrationBasin) %>%
- mutate(metapop = basin)
-
- iteration = 1
-
- # find the basin with the minimum autochtonous rate for the first time,
- basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
- group_by(year) %>%
- slice(which.min(autochtonousRate))
-
- # loop while autochtonousRate is still < 0.95
-
- while (basinWithMinAutochtonousRate$autochtonousRate < threshold ) {
- # while (iteration <= 4 ) {
- if (verbose) cat(iteration, ": ", basinWithMinAutochtonousRate$autochtonousRate, '\n' )
-
- # basinWithMinAutochtonousRate will be merged with origin basin sending the maximum number of fish in this basin
- metapopsToBeMerged <- exchangesDataUpdated %>%
- inner_join(basinWithMinAutochtonousRate %>% select(-autochtonousRate),
- by = c('year', 'migrationBasin' )) %>%
- filter(migrationBasin != originBasin) %>% # to avoid self merging
- group_by(year) %>%
- slice(which.max(effective)) %>%
- ungroup() %>%
- select(year, migrationBasin, originBasin)
-
- # update the connection table between basin and metapopulation by merging metapopulation
- for (i in 1:nrow(metapopsToBeMerged)) {
- mergedName = paste0("M", iteration, '_',i)
- metapopsToBeMerged_i <- metapopsToBeMerged %>% slice(i)
- if (verbose) cat("\t", metapopsToBeMerged_i$migrationBasin, ' + ', metapopsToBeMerged_i$originBasin, ' = ', mergedName, "\n")
-
- metapopulation <- metapopulation %>%
- mutate(metapop = if_else(year == metapopsToBeMerged_i$year &
- metapop %in% (metapopsToBeMerged_i %>% select(ends_with("Basin")) %>% unlist(use.names = FALSE)),
- mergedName,
- metapop))
- }
- # sum by migration and origin basins according to updated metapopulations
- # computed with initial exchangesData
- exchangesDataUpdated <- exchangesData %>%
- #sum up on migration basins
- inner_join(metapopulation,
- by = c('year', 'migrationBasin' = 'basin')) %>%
- group_by(year, metapop, originBasin) %>%
- summarise(effective = sum(effective), .groups = 'drop') %>%
- rename(migrationBasin = metapop) %>%
- # sum on origine basins
- inner_join(metapopulation,
- by = c('year', 'originBasin' = 'basin')) %>%
- group_by(year, metapop, migrationBasin) %>%
- summarise(effective = sum(effective), .groups = 'drop') %>%
- rename(originBasin = metapop)
-
- # minimum of autochtonous rate
- basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
- group_by(year) %>%
- slice(which.min(autochtonousRate))
-
- iteration = iteration + 1
- }
- return(list(metapopulation = metapopulation, exchangesData = exchangesDataUpdated))
-}
+require(tidyverse)
+source(GR3Dmetapopulation)
# upload data =====
# longer table of effective of different origin basins in each migration basin
exchangesData <- read_csv("../../data/output/northeastamerica/effectiveFluxes_1-observed.csv")
+
exchangesData %>% group_by(year, basin = originBasin) %>%
summarise(production = sum(effective), .groups = 'drop') %>%
inner_join(exchangesData %>% group_by(year, basin = migrationBasin) %>%