#!/usr/bin/env Rscript

###############################################################################
# script pour lancer calibration sur Datarmor avec PBS qsub
# toutes les variables d’environnement sont définies dans le script PBS
###############################################################################

# normalement 2 packages chargés: {Rmpi} et {snow} qui mettent à dispo:
# `makeCluster`, `clusterEvalQ`, `clusterExport`, `stopCluster`
sessionInfo()

library(doParallel, include.only = c("registerDoParallel")) # pour éviter d’écraser `makeCluster`

# la façon dont {calibrar} utilise {foreach} perturbe {doParallel}
# il faut toujours charger tous les pkgs nécessaires ici et dans `clusterEvalQ` dessous
# si vous le faites pas en double, erreur “function not found” directe!
library(calibrar, lib.loc = Sys.getenv("NICO_REPO"))
library(dplyr, warn.conflicts = FALSE)
library(data.table, include.only = c("fread"))
library(stringi, include.only = c("stri_extract_last_regex", "stri_join"))
options(dplyr.summarise.inform = FALSE)

cl <- makeCluster()
registerDoParallel(cl)

clusterEvalQ(cl, {
  library(calibrar, lib.loc = Sys.getenv("NICO_REPO"))

	library(dplyr, warn.conflicts = FALSE)
	library(data.table, include.only = c("fread"))
	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()

###############################################################################

runModel <- function(par) { # "par" provided by the calibrate function

	# run ISIS ======================================================
	# ATTENTION we are in RUN_*/i* folder

	# if necessary, manually save q of each generation
	current_time <- format(Sys.time(), format = "%F %Hh%M", tz = "CET", usetz = TRUE)
	current_gen <- 1L + length(list.files(pattern = "^params .+\\.csv$"))
	info_suffix <- sprintf(" gen %03d start %s", current_gen, current_time)

	# get back the capturability parameters (par) and write them in the format the model is able to read (csv, semantic)
	textLines <- readLines(Sys.getenv("PARAMS_FILE"))
	for (i in seq_along(par)) {
		id <- 2 + i
		tmp <- strsplit(textLines[id], ";")[[1]]
		textLines[id] <- paste0(tmp[1], ";", par[i])
	}
	write(textLines, file = Sys.getenv("PARAMS_FILE"))
	file.copy(from = Sys.getenv("PARAMS_FILE"), to = paste0("params", info_suffix, ".csv"))

	# launch ISIS from R, cf le script PBS
	current_ind <- basename(getwd()) # current individual of the population
	nomSimul_i <- Sys.getenv("SIMUL_BASENAME") %>% paste(., current_ind, sep = "_")
	# debug_file <- paste0("debugISIS", info_suffix, ".log")
	system2("myIsis", args = nomSimul_i, stdout = "debugISIS_octopus2017.log", stderr = "debugISIS_octopus2017.log")

	dossier_simul <- file.path(Sys.getenv("ISIS_BD"), "isis-database/simulations", nomSimul_i, "resultExports")

	# read the model outputs and tranform them into a list ==========

	# catch at age
	export1 <- file.path(dossier_simul, "CapturesPoids.csv") %>%
		fread(data.table = FALSE) %>%
		mutate(names = paste0("gp",group,"_Octopus2017")) %>%
		select(names, value)

	
	# combine (if several inputs)=======================================================
	# export <- rbind(export1, export2)
	export <- export1
	output <- as.list(export$value)
	names(output) <- export$names
	return(output)
	
	
	
	# path_export <- file.path(ifelse(Mon_disk == "~", "~", "c:/Users/Dedah"), paste0("isis-fish-4/isis-database/simulations/", nomSimu, "/resultExports/"))
	export1 <- read.table(paste(dossier_simul, "Captures_AgeStructure_pop.csv", sep = '/'), header = TRUE, sep = ';')
	export2 <- read.table(paste(dossier_simul, "Captures_MoisStructure_pop.csv", sep = '/'), header = TRUE, sep = ';')
	
	
	export2 = export2[-c(1:5, 11:12),]
	df.output <- data.frame('step'=c(paste("gp", c(0, 1:2) ,"_" ,"2018" , sep =""), paste("Mois", c(6:10) ,"_" ,"2018" , sep =""))
	                        , 'value'=NA)
	
	for (i in 1:3) {
	  df.output[i,2] <- export1[i, 4]
	}
	
	
	for (i in 4:8) {
	  df.output[i,2] <- export2[i-3, 3]
	}
	
	output <- list()
	for (i in 1:dim(df.output)[1]) {
	  output[[i]] <- df.output[i,2]
	}
	
	names(output) <- c(df.output[,1])
	
	return(output) ## List of nb year X nb pop group  elements : captures per age, per year
	

}

###############################################################################

#' @title residual sum of squares function
#' @description name of this fx must match what is declared in calibrationInfo.csv
#' @param obs,sim numeric values of LENGTH 1
userFunction <- function(obs, sim) {
	# obs peut être NA, c-a-d pas de valeurs observées mais valeurs simulées existent dans les exports
	# dû à des combinaisons irréelles entre flottille avec GSA zone pop
	res <-
		if (is.na(obs))    0
		else if (obs == 0) sim
		else               (obs - sim) / obs
	return(res ^ 2)
}

# calibrationInfo
calib_set <- calibration_setup(
	file = "data_calib/data_calib_Octopus2017/calibrationInfo.csv",
	control = list(col_skip = 1) # pour compatibilité avec ancienne version
)

# Observed data
calib_obs <- calibration_data(setup = calib_set, path = "data_calib/data_calib_Octopus2017")

# creation of the objective function
obj_fn <- calibration_objFn(model = runModel, setup = calib_set, observed = calib_obs)
# aggFn = calibrar:::.weighted.sum # by default

clusterExport(cl, c("userFunction", "obj_fn", "calib_set", "calib_obs"))

###############################################################################

ncpus <- as.integer(Sys.getenv("mpiproc"))

# control argument of the calibrate function
control <- list(
	master = "master/master_octopus2017", # a folder
	run = paste0("RUN_", Sys.getenv("SIMUL_BASENAME")), # a folder
	restart.file = paste0("ckpt_", Sys.getenv("SIMUL_BASENAME")),
	convergence = 0.001, maxgen = 500,
	popsize = if (ncpus < 50) 2 * ncpus else ncpus,
	REPORT = 1, trace = 5,
	verbose = TRUE, parallel = TRUE, nCores = ncpus
)

# N paramètres
calib <- calibrate(
	fn = obj_fn, method = "default", control = control,
	par = 9.87e-3,
	lower = 1e-5,
	upper = 1e-1
)

stopCluster(cl)