Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils Loading required package: utils > #!/usr/bin/env Rscript > > ############################################################################### > # script pour lancer recherche LHS sur Datarmor avec PBS qsub > # toutes les variables d’environnement sont définies dans le script PBS > ############################################################################### > > # normalement ce script doit être exécuté dans $PBS_O_WORKDIR > > # normalement 2 packages chargés: {Rmpi} et {snow} qui mettent à dispo: > # `makeCluster`, `clusterEvalQ`, `clusterExport`, `stopCluster` > sessionInfo() R version 3.6.3 (2020-02-29) Platform: x86_64-pc-linux-gnu (64-bit) Running under: SUSE Linux Enterprise Server 12 SP1 Matrix products: default BLAS/LAPACK: /appli/intel/parallel_studio/parallel_studio_xe_2017_update2/compilers_and_libraries_2017.2.174/linux/mkl/lib/intel64_lin/libmkl_rt.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices datasets utils methods base other attached packages: [1] snow_0.4-3 Rmpi_0.6-9 loaded via a namespace (and not attached): [1] compiler_3.6.3 parallel_3.6.3 > > library(doParallel, include.only = c("registerDoParallel")) # pour éviter d’écraser `makeCluster` > library(foreach) > > # utilisés par la partie non parallélisée (hors {foreach}) > library(dplyr, warn.conflicts = FALSE) > library(data.table, include.only = c("fread")) > > cl <- makeCluster() > registerDoParallel(cl) > > # utilisés par la partie parallélisée (avec {foreach}) > clusterEvalQ(cl, { + library(dplyr, warn.conflicts = FALSE) + library(data.table, include.only = c("fread")) + library(tidyr, include.only = c("separate", "pivot_longer")) + library(stringi, include.only = c("stri_extract_last_regex", "stri_join")) + options(dplyr.summarise.inform = FALSE) + }) %>% # cacher des messages encombrants + suppressMessages() %>% + suppressWarnings() %>% + capture.output() %>% + invisible() > > RERUN_SIMUL <- Sys.getenv("RERUN_SIMUL") %>% as.logical() > > ############################################################################### > > lesq <- c(paste0("q", 0:6), paste0("tar", 0:1)) > parameter_names <- lesq > > # capture en nombre par âge × année (France) > obs_catch_age_fr <- readRDS("Inputs/sardine_catch_age.RDS") %>% + rename(obs = data) %>% + select(names, obs) > > # capture en nombre par âge × année (Espagne) > obs_catch_age_esp <- readRDS("Inputs/sardine_catch_age_esp.RDS") %>% + rename(obs_esp = data, + names = name) %>% + select(names, obs_esp) > > # capture en nombre par âge × année (France + Espagne) > obs_catch_age <- obs_catch_age_fr %>% + left_join(obs_catch_age_esp) %>% + mutate(obs = obs + obs_esp) %>% + select(-obs_esp) Joining, by = "names" > > # capture en poids par flottille > obs_catch_flt <- readRDS("Inputs/sardine_catch_fleet.RDS") %>% + rename(obs = data) %>% + select(names, obs) > > # capture en poids par pays (France) > obs_catch_ctr_fr <- readRDS("Inputs/sardine_catch_country.RDS") %>% + mutate(year = as.character(year)) %>% + rename(obs = france, + names = year) %>% + select(names, obs) > > # capture en poids par pays (Espagne) > obs_catch_ctr_esp <- readRDS("Inputs/sardine_catch_country.RDS") %>% + mutate(year = as.character(year)) %>% + rename(obs = spain, + names = year) %>% + select(names, obs) > > # capture en poids par mois (France) > obs_catch_month_fr <- readRDS("Inputs/sardine_catch_month.RDS") %>% + rename(obs = data) %>% + select(names, obs) > > # capture en poids par mois (Espagne) > obs_catch_month_esp <- readRDS("Inputs/sardine_catch_month_esp.RDS") %>% + rename(obs = data) %>% + select(names, obs) > > # capture en poids (extra) > obs_catch_month_extra <- readRDS("Inputs/sardine_catch_extra.RDS") %>% + mutate(year = as.character(year)) %>% + group_by(year) %>% + summarise(catches_extra = sum(catches_extra)) %>% + ungroup() %>% + rename(obs = catches_extra, + names = year) %>% + select(names, obs) `summarise()` ungrouping output (override with `.groups` argument) > > planExperiences <- "LHS_q_et_al/LHS_q_tarEsp_pil/LHS_q_tarEsp.csv" %>% + fread(data.table = FALSE, col.names = parameter_names) > # le fichier a pas de nom colonnes pour au cas où utiliser avec script java plan simul > > # nombre de simulations > N_SIM <- Sys.getenv("PARAMS") %>% + paste0("LHS_q_et_al/LHS_", .) %>% + list.files(path = ., pattern = "row") %>% + length() > stopifnot(N_SIM == nrow(planExperiences)) > > # Only for tests : > # N_SIM = 4 > > # créer les dossiers ayant la même structure que quand on calibre > q_dirs <- sprintf("RUN_%s/i%04d", Sys.getenv("SIMUL_BASENAME"), 1:N_SIM) > for (i in q_dirs) dir.create(i, recursive = TRUE, showWarnings = FALSE) > > if (RERUN_SIMUL) invisible(file.copy( + from = sprintf("LHS_q_et_al/LHS_q_tarEsp_pil/row_%04d.csv", 1:N_SIM), + to = paste0(q_dirs, "/", Sys.getenv("PARAMS"), ".csv") + )) > > #' @title residual sum of squares function > #' @details traiter aussi le cas où obs = 0, > #' on s’en fout le cas NA car obs inexistantes sont enlevées avec merge(sim_data, obs_data) > err_calc <- function(obs, sim) { + res <- if_else(obs == 0, sim, (obs - sim) / obs) + return(res ^ 2) + } > > # export “biomasses fécondées” manque nom des colonnes > biomass_cols <- c("pop", "fish_age", "zonePop", "step", "value") > > start_year <- 2010 > sep <- "_" > > ############################################################################### > # partie calculs parallélisés > # les messages d’erreur peuvent être très cryptiques, donc vérifier bien les codes > # remarque: éviter d’avoir plusieurs objets portant même nom, des fois MPI aime pas ça > > my_progress <- function(...) cat(..., "\n", format(Sys.time(), format = "%F %Hh%M", tz = "CET", usetz = TRUE), "\n") > > my_progress("Step 0:") Step 0: 2024-02-19 22h50 CET > > # lancer simul, récupérer valeurs simulées ========================== > tmp <- foreach(i = 1:N_SIM) %dopar% { + + # run ISIS -------------------------------------------- + # on est pas dans les dossiers RUN_* comme {calibrar} + # mais toujours dans $PBS_O_WORKDIR + + # lancer simul + nomSimul_i <- sprintf("%s_i%04d", Sys.getenv("SIMUL_BASENAME"), i) + debug_file <- file.path(Sys.getenv("PBS_O_WORKDIR"), q_dirs[i], "debugISIS.log") + + if (RERUN_SIMUL) system2("myIsis", args = nomSimul_i, stdout = debug_file, stderr = debug_file) + + dossier_simul <- file.path(Sys.getenv("ISIS_BD"), "isis-database/simulations", nomSimul_i, "resultExports") + + # read isis csv --------------------------------------- + + biomasses_fecondees <- file.path(dossier_simul, "BiomasseFeconde.csv.gz") %>% + fread(data.table = FALSE, col.names = biomass_cols) %>% + filter(step %% 12 == 0) %>% # janvier + mutate( + year = as.integer(step / 12 + start_year), + ) %>% + group_by(year, fish_age) %>% + summarise(biomass = sum(value)) %>% + ungroup() + + # capture en poids par âge × année [anchois et sardine] + sim_catch_pds_age <- file.path(dossier_simul, "CapturesWeight_AgeStructure_3pop.csv.gz") %>% + fread(data.table = FALSE) %>% # regarder nom colonnes dans fichier script java export + #mutate(comp_names = stri_join("gp", group, sep, year + start_year)) %>% + mutate(comp_names = stri_join(year + start_year, sep, group)) %>% + rename(sim = value) %>% + select(comp_names, sim) + + # capture en poids par flottille × année + sim_catch_pds_flt <- file.path(dossier_simul, "CatchWeightSpFltYear.csv.gz") %>% + fread(data.table = FALSE) %>% # regarder nom colonnes dans fichier script java export + mutate(year = floor(yearQuarter/4), + comp_names = stri_join(year + start_year, sep, FLEET) + ) %>% + group_by(comp_names) %>% + summarise(sim = sum(valSim)) %>% + ungroup() + + sim_catch_pds_fr <- file.path(dossier_simul, "CatchWeightSpFltYear.csv.gz") %>% + fread(data.table = FALSE) %>% # regarder nom colonnes dans fichier script java export + filter(!FLEET %in% c("Espagnols","extra")) %>% + mutate(year = floor(yearQuarter/4),comp_names = stri_join(year + start_year)) %>% + group_by(comp_names) %>% + summarise(sim = sum(valSim)) %>% + ungroup() + + sim_catch_pds_esp <- file.path(dossier_simul, "CatchWeightSpFltYear.csv.gz") %>% + fread(data.table = FALSE) %>% # regarder nom colonnes dans fichier script java export + filter(FLEET == "Espagnols") %>% + mutate(year = floor(yearQuarter/4),comp_names = stri_join(year + start_year)) %>% + group_by(comp_names) %>% + summarise(sim = sum(valSim)) %>% + ungroup() + + sim_catch_pds_extra <- file.path(dossier_simul, "CatchWeightSpFltYear.csv.gz") %>% + fread(data.table = FALSE) %>% # regarder nom colonnes dans fichier script java export + filter(FLEET == "extra") %>% + mutate(year = floor(yearQuarter/4),comp_names = stri_join(year + start_year)) %>% + group_by(comp_names) %>% + summarise(sim = sum(valSim)) %>% + ungroup() + + # capture en poids par mois (France) + #sim_catch_pds_month <- file.path(dossier_simul, "CatchWeightSpMonth.csv") %>% + sim_catch_pds_month_fr <- file.path(dossier_simul, "CatchWeightSpFltMonth.csv.gz") %>% + fread(data.table = FALSE) %>% # regarder nom colonnes dans fichier script java export + filter(fleet != "Espagnols") %>% + group_by(month) %>% + summarise(valSim = sum(valSim)) %>% + ungroup() %>% + mutate( + year = month %/% 12, + month = month %% 12 + 1, + comp_names = stri_join(year + start_year, sep, month) + ) %>% + group_by(comp_names) %>% + summarise(sim = sum(valSim)) %>% + ungroup() + + # capture en poids par mois (Espagne) + #sim_catch_pds_month <- file.path(dossier_simul, "CatchWeightSpMonth.csv") %>% + sim_catch_pds_month_esp <- file.path(dossier_simul, "CatchWeightSpFltMonth.csv.gz") %>% + fread(data.table = FALSE) %>% # regarder nom colonnes dans fichier script java export + filter(fleet == "Espagnols") %>% + group_by(month) %>% + summarise(valSim = sum(valSim)) %>% + ungroup() %>% + mutate( + year = month %/% 12, + month = month %% 12 + 1, + comp_names = stri_join(year + start_year, sep, month) + ) %>% + group_by(comp_names) %>% + summarise(sim = sum(valSim)) %>% + ungroup() + + # finish ---------------------------------------------- + + return(list( + params = unlist(planExperiences[i,]), # convert from row to named vector + biomasses_fecondees = biomasses_fecondees, + sim_catch_pds_age = sim_catch_pds_age, + sim_catch_pds_flt = sim_catch_pds_flt, + sim_catch_pds_fr = sim_catch_pds_fr, + sim_catch_pds_esp = sim_catch_pds_esp, + sim_catch_pds_extra = sim_catch_pds_extra, + sim_catch_pds_month_fr = sim_catch_pds_month_fr, + sim_catch_pds_month_esp = sim_catch_pds_month_esp + )) + + }