Run SA

data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml

@@ -236,7 +236,7 @@
 				<species.Grow>
 					<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
 					<temperatureEffectGrow class="temperatureEffect.Rosso">
-						<Tmin>1.7</Tmin>
+						<Tmin>1.6</Tmin>
 						<Topt>13.37</Topt>
 						<Tmax>27.9</Tmax>
 					</temperatureEffectGrow>

exploration/NEA_sensitivity_analysis/prepareSensitivityAnalysis_CP.R

+require(tidyverse)
+
+# list of As factors with levels
+listFactorWithModalities = list()
+listFactorWithModalities[["Allee"]] =
+  list("sizeProportional" = list('processANG' = data.frame(process = 'species.ReproduceWithDiagnose', 
+                                                          par = 'Soffset',
+                                                          value = -1)),
+       "constant" = list('processANG' = data.frame(process = 'species.ReproduceWithDiagnose', 
+                                                  par = 'Soffset',
+                                                  value = 1000) ),
+       "noEffect" = list('processANG' = data.frame(process = 'species.ReproduceWithDiagnose', 
+                                                  par = 'Soffset',
+                                                  value = 0) ))
+
+listFactorWithModalities[["homing"]] =
+  list("high" = list('processANG' = data.frame(process = 'species.DisperseAndMigrateToRiverWithMultiNomDistriAndDeathBasin', 
+                                              par = 'pHomingAfterEquil',
+                                              value = .95)),
+       "low" = list('processANG' = data.frame(process = 'species.DisperseAndMigrateToRiverWithMultiNomDistriAndDeathBasin', 
+                                              par = 'pHomingAfterEquil',
+                                              value = 0.75)))
+
+
+listFactorWithModalities[["strayingDistance"]] =
+  list("high" = list('processANG' = data.frame(process = rep('species.DisperseAndMigrateToRiverWithMultiNomDistriAndDeathBasin',2), 
+                                               par = c('alpha0rep','alpha1Rep'),
+                                               value = c(-2.9, 19.7))),
+       "medium" = list('processANG' = data.frame(process = rep('species.DisperseAndMigrateToRiverWithMultiNomDistriAndDeathBasin',2), 
+                                                 par = c('alpha0rep','alpha1Rep'),
+                                                 value = c(-5.2, 21.2))),
+       "low" = list('processANG' = data.frame(process = rep('species.DisperseAndMigrateToRiverWithMultiNomDistriAndDeathBasin',2), 
+                                              par = c('alpha0Rep', 'alpha1Rep'),
+                                              value = c(-40.8, 35.8)))) 
+
+listFactorWithModalities[["survivalAfterReproduction"]] =
+  list("constant" = list('processANG' = data.frame(process = rep('species.SurviveAfterReproduction',2), 
+                                                   par = c('maximalSurvivalRate', 'temperatureEffectSurvivalAfterReproduction.alpha'),
+                                                   value = c(0.1*2, 0))),
+       "temperatureDependant" = list('processANG' = data.frame(process = rep('species.SurviveAfterReproduction',2), 
+                                                               par = c('maximalSurvivalRate', 'temperatureEffectSurvivalAfterReproduction.alpha'),
+                                                               value = c(0.18, -0.58))))
+
+# ---------------------------------------- sampling design ----
+#selectedListFactorWithModalities = listFactorWithModalities[c('Allee', 'homing')]
+ selectedListFactorWithModalities = listFactorWithModalities
+factorModalities = tibble(names(selectedListFactorWithModalities[[1]]))
+names(factorModalities) <- names(selectedListFactorWithModalities[1]) 
+for (i in 2:length(selectedListFactorWithModalities)) {
+  factor = tibble(names(selectedListFactorWithModalities[[i]]))
+  names(factor) <- names(selectedListFactorWithModalities[i]) 
+  factorModalities <-  factorModalities %>%  expand_grid(factor)
+}
+factorModalities <- factorModalities %>% expand_grid(replicate = 1)#nombre de rplicats
+
+# -------------------------------------  simulation list ----
+simulationList = list()
+for (simul in 1:nrow(factorModalities)) {
+  
+  # simulation names
+  nameSimul = tibble(name = paste(names(factorModalities), factorModalities[simul,], sep = ':')) %>% 
+    summarise(name = paste(name, collapse = '_')) %>% unlist(use.names = FALSE)
+  #cat(nameSimul, '\n')
+  
+  # parameters to be modified in class Environment 
+  parametersValueENV= c(par = 'simulationName', value = paste0('simul_', simul))  
+
+  
+  # parameters to be modified in class AquaNismGroup 
+  parametersValueANG = c() 
+  for (factor in names(factorModalities %>% select(-replicate))) {
+    modality = factorModalities[simul, factor] %>% unlist(use.names = FALSE)
+    #cat('\t',factor , ':', modality, '\n')
+    
+    parametersValueANG = parametersValueANG %>% 
+      bind_rows(listFactorWithModalities[[factor]][[modality]]$processANG %>% 
+                  inner_join(no_processesEachStep %>% select(process, no_java), by =c('process')) %>% 
+                  mutate(target = paste("processes.processesEachStep", no_java, par, sep ='.')) %>% 
+                  select(target, value, process))
+    
+    
+  }
+
+  simulationList[[simul]] = list(name = paste0('simul_', simul),
+                                 longName = nameSimul,
+                                 replicate = factorModalities[simul,'replicate'],
+                                 parametersValueANG = parametersValueANG,
+                                 parametersValueENV = parametersValueENV)
+}
+simulationList
+simulationList = unlist(map(simulationList, 'longName'))
+
+remove(factor, modality, nameSimul, simul, parametersValueANG, parametersValueENV, selectedListFactorWithModalities)

exploration/NEA_sensitivity_analysis/tryFromJava_CP.R

+library(tidyverse)
+library(rJava) # see J4R
+#library(yaml)
+library(jdx)
+library(tictoc)
+library(XML)
+library(doParallel)
+
+
+## to have the same working directory as GR3D
+#setwd("~/Documents/workspace/GR3D/exploration/NEA_sensitivity_analysis")#A adpater sous Windows 
+setwd("C:/Users/camille.poulet/workspace/GR3D/exploration/NEA_sensitivity_analysis") 
+setwd("../..")
+
+rm(list = ls())
+
+
+# ======================================================== 
+# load the function to calculate the metapopulations
+#source("exploration/GR3D_Rdescription/GR3Dmetapopulation.R")
+
+source("exploration/GR3D_Rdescription/GR3Dmetapopulation.R")
+
+# load the xml
+#fishXML <- xmlToList(xmlParse("data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml"))
+fishXML <- xmlToList(xmlParse("data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml"))
+
+
+#Lorsqu'un process est commente dans JAVA, cration d'une ligne R 
+no_processesEachStep = tibble( process = names(fishXML[["species.DiadromousFishGroup"]][["processes"]][["processesEachStep"]][])) %>% 
+  mutate(no_R = row_number()) %>% 
+  filter(process != "comment") %>% 
+  mutate(no_java = row_number()-1) %>% 
+  select(no_R, no_java, process)
+
+no_processesAtEnd = tibble( process = names(fishXML[["species.DiadromousFishGroup"]][["processes"]][["processesAtEnd"]][])) %>% 
+  mutate(no_R = row_number()) %>% 
+  filter(process != "comment") %>% 
+  mutate(no_java = row_number()-1) %>% 
+  select(no_R, no_java, process)
+
+
+
+# ====================================== LOCAL FUNCTIONS =====================================
+#renvoie un R object en convertissant un fichier Java en df R 
+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) ) %>%
+    as_tibble() %>% 
+    type_convert()#convertit directement au bon format, et non en String 
+  
+  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))
+# }
+
+#fonction locale qui permet d'aller explorer les processus qui nous intresses 
+exploreEasyRun = function(simulationName){
+  # call GR3D method to get model outputs ----------------------------------------------
+  
+  # 00. simulation name 
+  # NE MARCHE PAS
+  # simulationName = .jcall("analysis.EasyRun","[Ljava/lang/String;", "getStringsFromEnvironment",
+  #                         "simulationName")
+  
+  # 0. abundance (caution spawner abundance of the last reproduction (no averaging) )
+  lastSpawnerAdundance = .jcall("analysis.EasyRun","D", "getSingleValueFromAquanismGroupProcess", 
+                                "processes.processesEachStep.10.getSpawnerEffective")
+  
+  # 1. northernmost populated basin_id   
+  # 2. centroid latitude
+  #r[[Ljava/lang/String; renvoie un tableau de String 
+  out_range = .jcall("analysis.EasyRun","[[Ljava/lang/String;", "getValuesFromAquanismGroupProcess", 
+                     "processes.processesEachStep.10.exportToR")
+  range = suppressMessages(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 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 = suppressMessages (convertFromSALtoR(out_spawners_features))
+  
+  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)
+  
+  # ------------------------------------------------------------ #
+  # 6. 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 = 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)
+  if (nrow(selectedBasin > 0)) {
+    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)
+  } else
+    nbMetapopulation =NA
+  return(tibble(simulationName = simulationName,
+                lastSpawnerAdundance=lastSpawnerAdundance,
+                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")
+
+
+# -----------------------------------------------------
+#  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")
+
+# run the GR3D with updated parameter values in processes --------------------------------------------
+# runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV )
+
+# .jcall("analysis.EasyRun", "V", "runSimulation", arguments, 
+#        .jarray(parametersNamesANG), .jarray(thetasANG),
+#        .jarray(parametersNamesENV), .jarray(thetasENV))
+
+
+# call GR3D method to get model outputs ----------------------------------------------
+
+# #==================== run simulation 1 =====
+# 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")
+#   
+#   runEasyRun(arguments, 
+#              parametersNamesANG = simulationParameter$parametersNamesANG,
+#              thetasANG = simulationParameter$thetasANG, 
+#              parametersNamesENV = simulationParameter$parametersNamesENV, 
+#              thetasENV = simulationParameter$thetasENV )
+#   
+#   #tic()
+#   result = exploreEasyRun()
+#   #toc()
+#   return(result)
+# }
+# 
+# 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()
+
+# simulationParameter = simulationList[[35]]
+# toto = simulationParameter$parametersValueANG[,'value'] 
+# dim(toto)
+# titi= c(0, 0.75)
+# dim(titi)
+#==================== run simulation 2 =====
+#lance la machine virtuelle qui permet de connecter le point Jar en R 
+#mis les pramtres correctement et executer la simulation pour explorer les rsultats du modle 
+runSimulation2 = function(simulationParameter, simSDuration = 600) {
+  cat(simulationParameter$name, ' = ', simulationParameter$longName,'\n')
+  # ---------------------------------------------------------- start the VM
+  jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
+  .jinit(classpath = jarfile, force.init = FALSE)
+  .jclassPath()
+  
+  #---------------------------------------------- parameters for simulation 
+  # replicate and default values
+  arguments = c('-simDuration', simSDuration, '-simBegin',2,
+                '-timeStepDuration', 1, 
+                '-RNGStatusIndex', format(simulationParameter$replicat %>% unlist(use.names = FALSE), 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")
+  arguments
+  # new values in the XLMs
+  parametersNamesANG = simulationParameter$parametersValueANG[,'target']
+  thetasANG = simulationParameter$parametersValueANG[,'value']
+  parametersNamesENV = simulationParameter$parametersValueENV['par']
+  thetasENV = simulationParameter$parametersValueENV['value']
+  # cat(parametersNamesANG, '\n')
+  # cat(thetasANG, '\n')
+  # cat(parametersNamesENV,'\n')
+  # cat(thetasENV, '\n')
+  # --------------------------------------- run the simulation with new parameters
+  #tic("simulationRun")
+  .jcall("analysis.EasyRun", "V", "runSimulationANG_ENV", arguments,
+         .jarray(parametersNamesANG), .jarray(thetasANG),#la fonction Jarray traduit le format R en Java
+         .jarray(parametersNamesENV), .jarray(thetasENV))
+  #toc()
+  
+  #tic("simulationExploration")
+  result = exploreEasyRun(simulationParameter$name)
+  #toc()
+  
+  return(result)
+}
+
+
+runSimulation2(simulationList[[35]])
+
+
+#en squentielle excute toutes les simulations de la liste
+#resSimulation = lapply(simulationList, runSimulation2)
+
+
+# registerDoParallel(NULL)
+# cl <- detectCores() %>% -1 %>% makeCluster
+#registerDoParallel(cores = 1)
+
+stopImplicitCluster()#stopper le coeur des processuers 
+#cl <- makeCluster(3)
+cl <- detectCores() %>% -1 %>% makeCluster #je fais un cluster avec mes 7 coeur de processeur, pour en laisser un disponible pour faire les autres taches pour que je puisse travailler
+registerDoParallel(cl) #assigner des taches  chaque processeru/retablir les taches 
+rm(resSimulations)
+
+#unlist(map(simulationList, 'longName'))[c(1,7,13,19 )]
+tic('parralel')
+#foreach
+#%dopar : exectute en parallle les tache du foreach 
+#quivalent ici du lapply en plus puissant car permet de mobiliser plusieurs coeur en mme temps.
+#dans simul = c(...), on prcise les lignes de la lignes de simulations qu'on souhaite oprer. 
+#Ici, on veut les 36l du tableau donc on mettra 1:35
+resSimulations_2 = foreach(simul = seq(), 
+                         .combine = rbind, .verbose = FALSE, .multicombine = TRUE,
+                         .packages = c('tidyverse', 'rJava', 'jdx'),
+                         .inorder = FALSE) %dopar% {
+                           res = runSimulation2(simulationList[[simul]],
+                                                )
+                           #res = runSimulation3(simul)
+                         }
+toc()
+resSimulationsSequentielle
+
+simulationList[[1]]$parametersValueANG$value
+simulationList[[17]]$parametersValueANG$value
+
+save(resSimulations, file = 'exploration/NEA_sensitivity_analysis/simulationSequentielle.Rdata')
+
+
+