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

Merge branch 'SpawnerRunAnalysis' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D...

Merge branch 'SpawnerRunAnalysis' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D into exploration_GR3D_process
parents 1b6f6d7a d0ad927b
#=============================================================================================================
# 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))
}
library(tidyverse)
#TODO extract from xml
thermalRange = tibble(process = "Grow", Tmin = 1.7, Topt = 4.5, Tmax = 27.9, type = "Rosso") %>%
bind_rows(tibble(process = "SurvivalSpawnerInRiv", Tmin = 5.4, Topt = 16.7, Tmax = 27.5, type = "Rosso") ) %>%
bind_rows(tibble(process = "Reproduction", Tmin = 5.1, Topt = 13.5, Tmax = 24.5, type = "Rosso") ) %>%
bind_rows(tibble(process = "SurviveAfterReproduction", Tmin = -Inf, Topt = 19.58 - log(19)/0.58, Tmax = 19.58 + log(19)/0.58, type = "logit") )
thermalRange %>% mutate(numero = row_number(),
process = factor(process,
levels = c("Grow","SurvivalSpawnerInRiv", "SurviveAfterReproduction", "Reproduction" ))) %>%
ggplot() +
geom_segment(aes(x = Tmin, y = process, xend = Tmax, yend = process)) +
geom_point(aes(x = Topt, y = process))
library(tidyverse)
library(XML)
#source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
setwd("~/Documents/workspace/GR3D/exploration/NEA_sensitivity_analysis")
......
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) %>%
......
......@@ -3,15 +3,16 @@ library(rJava) # see J4R
#library(yaml)
library(tidyverse)
library(jdx)
library(tictoc)
rm(list = ls())
rm(list = ls())
## to have the same working directory as GR3D
#setwd("~/Documents/workspace/GR3D/exploration/NEA_sensitivity_analysis")
setwd("../..")
source("exploration/GR3D_Rdescription/GR3Dmetapopulation.R")
# ====================================== LOCAL FUNCTIONS =====================================
convertFromSALtoR = function(java_out){
require(jdx)
......@@ -22,8 +23,8 @@ convertFromSALtoR = function(java_out){
R_out = lapply(seq(2, length(java_out)), FUN = function(i){convertToR(java_out[[i]])}) %>%
unlist() %>%
matrix(nrow = length(java_out) - 1, byrow = TRUE, dimnames = list(NULL, headers) ) %>%
data.frame() %>%
suppressMessages(type_convert())
as_tibble() %>%
type_convert()
return(R_out)
}
......@@ -41,7 +42,7 @@ exploreEasyRun = function(){
# 1. northernmost populated basin_id
# 2. centroid latitude
out_range = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.10.exportToR")
range = convertFromSALtoR(out_range)
range = suppressMessages(convertFromSALtoR(out_range))
# 3. log-likelihood ----
# the proportion of reproductions leading to a recruitment > 50 over the 1920-1950 period was used
......@@ -55,66 +56,79 @@ exploreEasyRun = function(){
# then in R: annual slope of the linear regression
# and average over the period
# 5. mean age of female and age
# annual mean age in each basin,
# then in R: average per year,
# and average over the period (1920 and 1950)
# annual mean age over the period (1920 and 1950) in each basin,
# and average over basin
out_spawners_features = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.9.exportToR")
spawners_features = convertFromSALtoR(out_spawners_features)
spawners_features =suppressMessages (convertFromSALtoR(out_spawners_features))
## ICI coder pente de la droite moyenne des ages
slopePrimiparous = lm(pct_primiparous~latitude, data = spawners_features)$coeff[2]
meanAgeFemale = mean(spawners_features$mean_age_female, na.rm = TRUE)
meanAgeMale = mean(spawners_features$mean_age_male, na.rm = TRUE)
# ------------------------------------------------------------ #
# 5. number of metapopulation based on the exchange matrix ----
# mean annual exchange matrix between 1920 and 1950
out_exchanges = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.16.exportToR")
exchanges = convertFromSALtoR(out_exchanges)
exchanges = suppressMessages(convertFromSALtoR(out_exchanges))
# keep only basin with a totalRun (sum over migrationBasin) >0
selectedBasin <- exchanges %>% group_by(basin = migrationBasin) %>%
summarise(totalRun = sum(effective), .groups = 'drop') %>%
filter(totalRun > 0) %>% select(basin)
## ICI lancer metapopulationIdentification
metapopulations <- exchanges %>% inner_join(selectedBasin, by=c('migrationBasin' = 'basin')) %>%
inner_join(selectedBasin, by=c('natalBasin' = 'basin')) %>%
mutate(year =0) %>% rename(originBasin = natalBasin) %>%
identyMetapopulation(., threshold = .95, verbose = FALSE)
nbMetapopulation = metapopulations$metapopulation %>%
distinct(metapop) %>%
summarise(n = n()) %>%
unlist(use.names = FALSE)
return(list(northern_basin_id = range$northern_basin_id,
effective_centroid_latitude = range$effective_centroid_latitude,
logL = logL,
spawners_features = spawners_features,
exchanges = exchanges))
return(tibble(northern_basin_id = range$northern_basin_id,
effective_centroid_latitude = range$effective_centroid_latitude,
logL = logL,
slopePrimiparous = slopePrimiparous,
meanAgeFemale = meanAgeFemale,
meanAgeMale = meanAgeMale,
nbMetapopulation = nbMetapopulation))
}
# GR3D arguments -----------------------------------------
jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
simDuration = 150
simBegin = 1
timeStepDuration = 1
seed = 1
#'-q',
arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
'-timeStepDuration', timeStepDuration,
'-RNGStatusIndex', format(seed,scientific = FALSE),
'-groups',"data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml",
'-env',"data/input/northeastamerica/RIOBNneaBasins_Rjava.xml",
'-observers',"data/input/northeastamerica/RIO_obs_empty.xml")
#' jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
#' simDuration = 150
#' simBegin = 1
#' timeStepDuration = 1
#'
#' seed = 1
#' #'-q',
#' arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
#' '-timeStepDuration', timeStepDuration,
#' '-RNGStatusIndex', format(seed,scientific = FALSE),
#' '-groups',"data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml",
#' '-env',"data/input/northeastamerica/RIOBNneaBasins_Rjava.xml",
#' '-observers',"data/input/northeastamerica/RIO_obs_empty.xml")
# -----------------------------------------------------
# initializes the Java Virtual Machine -------------------------------------
jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
.jinit(NULL, force.init = TRUE)
.jinit(classpath = jarfile, force.init = TRUE)
.jclassPath()
# modification of xlm parameters -----------------------
parametersNamesANG = c("processes.processesEachStep.11.Soffset")
thetasANG = c(-1)
parametersNamesENV = c("simulationName")
thetasENV = c("Allee")
# parametersNamesANG = c("processes.processesEachStep.11.Soffset")
# thetasANG = c(-1)
#
# parametersNamesENV = c("simulationName")
# thetasENV = c("Allee")
# run the GR3D with updated parameter values in processes --------------------------------------------
runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV )
# runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV )
# .jcall("analysis.EasyRun", "V", "runSimulation", arguments,
# .jarray(parametersNamesANG), .jarray(thetasANG),
......@@ -122,42 +136,60 @@ runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasEN
# call GR3D method to get model outputs ----------------------------------------------
res1 = exploreEasyRun()
res1$spawners_features
#
listeSimulation <- list(list(replicat = 1,
parametersNamesANG = c("processes.processesEachStep.11.Soffset"),
thetasANG = c(-1),
parametersNamesENV = c("simulationName"),
thetasENV = c("Allee")),
list(replicat = 1,
parametersNamesANG = c("processes.processesEachStep.11.Soffset"),
thetasANG = c(0),
parametersNamesENV = c("simulationName"),
thetasENV = c("no Allee")))
for (i in 1:length(listeSimulation)) {
arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
'-timeStepDuration', timeStepDuration,
'-RNGStatusIndex', format(listeSimulation[[i]]$replicat,scientific = FALSE),
runSimulation = function(simulationParameter) {
arguments = c('-simDuration', 600, '-simBegin',2,
'-timeStepDuration', 1,
'-RNGStatusIndex', format(simulationParameter$replicat,scientific = FALSE),
'-groups',"data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml",
'-env',"data/input/northeastamerica/RIOBNneaBasins_Rjava.xml",
'-observers',"data/input/northeastamerica/RIO_obs_empty.xml")
parametersNamesANG = listeSimulation[[i]]$parametersNamesANG
thetasANG = listeSimulation[[i]]$thetasANG
parametersNamesENV = listeSimulation[[i]]$parametersNamesENV
thetasENV = listeSimulation[[i]]$thetasENV
runEasyRun(arguments,
parametersNamesANG,
thetasANG,
parametersNamesENV ,
thetasENV )
res1 = exploreEasyRun()
cat(res1$effective_centroid_latitude, "\n")
parametersNamesANG = simulationParameter$parametersNamesANG,
thetasANG = simulationParameter$thetasANG,
parametersNamesENV = simulationParameter$parametersNamesENV,
thetasENV = simulationParameter$thetasENV )
#tic()
result = exploreEasyRun()
#toc()
return(result)
}
toto = listeSimulation[[i]]$parametersNamesANG
#============================
listeSimulation <- list(list(replicat = 1,
parametersNamesANG = c("processes.processesEachStep.11.Soffset",
"processes.processesEachStep.6.pHomingAfterEquil"),
thetasANG = c(-1, 0.97),
parametersNamesENV = c("simulationName"),
thetasENV = c("Allee_highHoming")),
list(replicat = 1,
parametersNamesANG = c("processes.processesEachStep.11.Soffset",
"processes.processesEachStep.6.pHomingAfterEquil"),
thetasANG = c(-1, 0.75),
parametersNamesENV = c("simulationName"),
thetasENV = c("Allee_lowHoming")),
list(replicat = 1,
parametersNamesANG = c("processes.processesEachStep.11.Soffset",
"processes.processesEachStep.6.pHomingAfterEquil"),
thetasANG = c(0, 0.97),
parametersNamesENV = c("simulationName"),
thetasENV = c("noAllee_highHoming")),
list(replicat = 1,
parametersNamesANG = c("processes.processesEachStep.11.Soffset",
"processes.processesEachStep.6.pHomingAfterEquil"),
thetasANG = c(0, 0.75),
parametersNamesENV = c("simulationName"),
thetasENV = c("noAllee_lowHoming"))
)
tic()
resSimulation = lapply(listeSimulation, runSimulation)
toc()
save(listeSimulation, resSimulation, file = 'simulation.Rdata')
......@@ -254,40 +254,34 @@ public class AnalyseFishDistribution extends AquaNismsGroupProcess<DiadromousFis
// initialise the distribution range
double southernEdge = northernLimit;
double northernEdge = southernLimit;
int southernId = -1;
int northernId = -1;
String southernName = "";
String northernName = "";
// String southernName = "";
// String northernName = "";
RiverBasin northernBasin = null;
RiverBasin southernBasin = null;
// initialise for centroid
TreeMapForCentile latitudeEffective = new TreeMapForCentile();
TreeMapForCentile latitudePresence = new TreeMapForCentile();
String[][] export = new String[2][6];
for (RiverBasin riverBasin : riverBasins) {
double meanLastRecruitment = riverBasin.getLastRecruitments().getMean();
if (meanLastRecruitment >= minimumRecruitsForPopulatedBasin) {
// the river basin is considered populated
// NOTE : recruiit number is calulated after thermal tolerance impact
// CHECK : recruiit number is calulated after thermal tolerance impact
// (which is intreaged in the stock-recruitment relationship)
// southern edge
if (riverBasin.getLatitude() < southernEdge) {
// the basin is the new southern edge
southernEdge = riverBasin.getLatitude();
southernId = riverBasin.getBasin_id();
southernName = riverBasin.getName();
southernBasin = riverBasin;
}
// northern edge
if (riverBasin.getLatitude() > northernEdge) {
// the basin is the new northern edge
northernEdge = riverBasin.getLatitude();
northernId = riverBasin.getBasin_id();
northernName = riverBasin.getName();
northernBasin = riverBasin;
}
// for distribution centroide computation
......@@ -297,29 +291,41 @@ public class AnalyseFishDistribution extends AquaNismsGroupProcess<DiadromousFis
}
}
// if the universe is empty (none of basins is populated)
if (southernEdge == northernLimit & northernEdge == southernLimit) {
southernEdge = (northernLimit + southernLimit) / 2.;
northernEdge = southernEdge;
}
// distribution centroids computation
double basinCentroid = latitudePresence.calculateMedian();
double effectiveCentroid = latitudeEffective.calculateMedian();
// export
String[][] export = new String[2][8];
export[0][0] = "southern_basin_id";
export[0][1] = "southern_basin_name";
export[0][2] = "northern_basin_id";
export[0][3] = "northern_basin_name";
export[0][4] = "basin_centroid_latitude";
export[0][5] = "effective_centroid_latitude";
export[1][0] = String.valueOf(southernId);
export[1][1] = southernName;
export[1][2] = String.valueOf(northernId);
export[1][3] = northernName;
export[1][4] = String.valueOf(basinCentroid);
export[1][5] = String.valueOf(effectiveCentroid);
export[0][2] = "southern_latitude";
export[0][3] = "northern_basin_id";
export[0][4] = "northern_basin_name";
export[0][5] = "northern_latitude";
export[0][6] = "basin_centroid_latitude";
export[0][7] = "effective_centroid_latitude";
if (southernBasin != null & northernBasin != null) {
export[1][0] = String.valueOf(southernBasin.getBasin_id());
export[1][1] = southernBasin.getName();
export[1][2] = String.valueOf(southernBasin.getLatitude());
export[1][3] = String.valueOf(northernBasin.getBasin_id());
export[1][4] = northernBasin.getName();
export[1][5] = String.valueOf(northernBasin.getLatitude());