run GR3D from R

exploration/NEA_sensitivity_analysis/extractNumeroProcess.R

 library(tidyverse)
+library(XML)
 #source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
 
 setwd("~/Documents/workspace/GR3D/exploration/NEA_sensitivity_analysis")

exploration/NEA_sensitivity_analysis/tryFromJava.R

@@ -3,6 +3,7 @@ library(rJava) # see J4R
 #library(yaml)
 library(tidyverse)
 library(jdx)
+library(tictoc)
 
 rm(list = ls())
 
@@ -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,28 +56,40 @@ 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))
   
-  ## ICI  lancer metapopulationIdentification
+  # 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)
   
+  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)
   
-  return(list(northern_basin_id = range$northern_basin_id, 
-              effective_centroid_latitude = range$effective_centroid_latitude,
-              logL = logL,
-              spawners_features = spawners_features,
-              exchanges = exchanges))
+  nbMetapopulation = metapopulations$metapopulation %>% distinct(metapop) %>% summarise(n = n())
+  
+  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))
 }
 
 
@@ -97,9 +110,6 @@ arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
               '-observers',"data/input/northeastamerica/RIO_obs_empty.xml")
 
 
-
-
-
 # -----------------------------------------------------
 #  initializes the Java Virtual Machine -------------------------------------
 .jinit(NULL, force.init = TRUE)
@@ -122,11 +132,32 @@ runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasEN
 
 
 # call GR3D method to get model outputs ----------------------------------------------
-res1 = exploreEasyRun()
-res1$spawners_features
-
 
 #
+
+
+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"), 
                              thetasANG = c(-1),
@@ -138,26 +169,33 @@ listeSimulation <- list(list(replicat = 1,
                              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")
-}
+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.97),
+                             parametersNamesENV = c("simulationName"), 
+                             thetasENV = c("noAllee_lowHoming"))
+                        )
 
-toto = listeSimulation[[i]]$parametersNamesANG
+lapply(listeSimulation, runSimulation)
\ No newline at end of file