Commit 3e87cc1a authored by patrick.lambert's avatar patrick.lambert
Browse files

with a loop on several simulation

parent 1a3d84fc
library(tidyverse)
library(rJava) # see J4R
#library(yaml)
library(tidyverse)
......@@ -5,40 +6,90 @@ library(jdx)
rm(list = ls())
# # parameters and prios=rs of parameters to be fitted
# inputData = yaml.load_file("ABC_seq_2param_3stats_NewAfterSA.yaml")
# # load prior
# prior = list()
# for (i in 1:length(inputData$parameters)) {
# p = inputData$parameters[[i]]
# prior[[i]] = unlist(p$prior)
# }
#
# parametersNames = list()
# for (i in 1:length(inputData$parameters)) {
# parametersNames = c(parametersNames, inputData$parameters[[i]]$javapath)
# }
# parametersNames = unlist(parametersNames)
# GR3D arguments -----------------------------------------
## to have the same working directory as GR3D
#setwd("~/Documents/workspace/GR3D/exploration/NEA_sensitivity_analysis")
setwd("../..")
# path to outputs
#outputDir = "simus/"
#outputFileName = paste0("SA/",seed)
# ====================================== LOCAL FUNCTIONS =====================================
convertFromSALtoR = function(java_out){
require(jdx)
# headers
headers = convertToR(java_out[[1]])
# data
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())
return(R_out)
}
runEasyRun = function(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV) {
# run the GR3D with updated parameter values in processes --------------------------------------------
.jcall("analysis.EasyRun", "V", "runSimulation", arguments,
.jarray(parametersNamesANG), .jarray(thetasANG),
.jarray(parametersNamesENV), .jarray(thetasENV))
}
exploreEasyRun = function(){
# call GR3D method to get model outputs ----------------------------------------------
# 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)
# 3. log-likelihood ----
# the proportion of reproductions leading to a recruitment > 50 over the 1920-1950 period was used
# as a proxy of predicted probability of presence in a basin.
# Then a log-likelihood of a binomial distribution is computed
logL = .jcall("analysis.EasyRun", "D", "getSingleValueFromAquanismGroupProcess", "processes.processesAtEnd.0.getLogLikelihood");
# 4. proportion of first-time spawners
# annual % of first-time spawners according to latitude (or basin_id)
# between 1920-1950
# 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)
out_spawners_features = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.9.exportToR")
spawners_features = convertFromSALtoR(out_spawners_features)
## ICI coder pente de la droite moyenne des ages
# ------------------------------------------------------------ #
# 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)
## ICI lancer metapopulationIdentification
return(list(northern_basin_id = range$northern_basin_id,
effective_centroid_latitude = range$effective_centroid_latitude,
logL = logL,
spawners_features = spawners_features,
exchanges = exchanges))
}
# GR3D arguments -----------------------------------------
jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
simDuration = 150
simBegin = 1
timeStepDuration = 1
seed = 1
arguments = c('-q','-simDuration', simDuration, '-simBegin',simBegin,
#'-q',
arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
'-timeStepDuration', timeStepDuration,
'-RNGStatusIndex', format(seed,scientific = FALSE),
'-groups',"data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml",
......@@ -47,111 +98,66 @@ arguments = c('-q','-simDuration', simDuration, '-simBegin',simBegin,
# -----------------------------------------------------
# initializes the Java Virtual Machine -------------------------------------
.jinit(classpath ="", force.init = TRUE)
.jinit(NULL, force.init = TRUE)
.jinit(classpath = jarfile, force.init = TRUE)
.jclassPath()
# modification of xlm parameters -----------------------
parametersNamesANG = c("processes.processesEachStep.11.Soffset")
thetasANG = c(0)
thetasANG = c(-1)
parametersNamesENV = c("simulationName")
thetasENV = c("noAllee")
runEasyRun = function(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV ){
# run the GR3D with updated parameter values in processes --------------------------------------------
.jclassPath()
.jcall("analysis.EasyRun", "V", "runSimulation", arguments,
.jarray(parametersNamesANG), .jarray(thetasANG),
.jarray(parametersNamesENV), .jarray(thetasENV))
# call GR3D method to get model outputs ----------------------------------------------
#1. northern limit and distribution centroid ----
out = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.10.exportToR")
range <- data.frame(matrix(convertToR(out[[2]]), nrow = 1, dimnames = list(NULL, convertToR(out[[1]])))) %>% type_convert()
return(list(range = range))
}
thetasENV = c("Allee")
# run the GR3D with updated parameter values in processes --------------------------------------------
runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV )
# run the GR3D with updated parameter values in processes --------------------------------------------
.jclassPath()
.jcall("analysis.EasyRun", "V", "runSimulation", arguments,
.jarray(parametersNamesANG), .jarray(thetasANG),
.jarray(parametersNamesENV), .jarray(thetasENV))
# .jcall("analysis.EasyRun", "V", "runSimulation", arguments,
# .jarray(parametersNamesANG), .jarray(thetasANG),
# .jarray(parametersNamesENV), .jarray(thetasENV))
# call GR3D method to get model outputs ----------------------------------------------
#TODO export [char] and then data.frame(toto) %>% type_convert()
# ----------------------------------------------------------------------------------#
# 1. northern limit and distribution centroid ----
# only basins with more than 50 recruits in average in the last 50 years (colonised basins)
# compute only in SPRING
# access through getRangeDistribution in DiadromousFishGroup ==== NOT CORRECT
# range[1] = latitude of the median distribution ( 50 % of effectives are northern)
# range[2] = southern limit latitude of the fish distribution
out = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.10.exportToR")
range <- data.frame(matrix(convertToR(out[[2]]), nrow = 1, dimnames = list(NULL, convertToR(out[[1]])))) %>% type_convert()
# ------------------------------------------------------------ #
# 2. log-likelihood ----
# the proportion of reproductions leading to a recruitment > 50 over the 1920-1950 period was used
# as a proxy of predicted probability of presence in a basin.
# Then a log-likelihood of a binomial distribution is computed
.jcall("analysis.EasyRun", "D", "getSingleValueFromAquanismGroupProcess", "processes.processesAtEnd.0.getLogLikelihood");
# ------------------------------------------------------------ #
# 3. proportion of first-time spawners
# annual % of first-time spawners according to latitude (or basin_id)
# between 1920-1950
# then in R: annual slope of the linear regression
# and average over the period
# 4. Mean age of the spawner for male and female ----
# annual mean age in each basin,
# then in R: average per year,
# and average over the period (1920 and 1950)
# mean (+sd) and slope after
out = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.8.exportToR")
convertFromSALtoR = function(java_out){
require(jdx)
out2 = matrix(unlist(lapply(java_out, convertToR)), nrow = length(java_out) , byrow = TRUE)
out3 <- out2[-1,]
colnames(out3) <- unlist(out2[1,])
R_out <- as.data.frame(out3) %>% type_convert()
return(R_out)
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),
'-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")
}
convertFromSALtoR(java_out = out)
# ------------------------------------------------------------ #
# ------------------------------------------------------------ #
# 5. number of metapopulation based on the exchange matrix ----
# annual exchange matrix between 1920 and 1950
# then in R: compute average matrix over the period,
# run metapopulationIdentification
# moyenne sur les bassins versants des moyennes sur la peride des moyennes des flux
# access through species.IdentifyPopulation in processesAtEnd
#TODO A FAIRE !!!
# =================================================================================
toto = listeSimulation[[i]]$parametersNamesANG
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