R/atlantisfmsy_core.R

Defines functions atlantisfmsy_fmaxcatch atlantisfmsy_inisimu atlantis_bachchange atlantis_fchange

Documented in atlantis_bachchange atlantis_fchange atlantisfmsy_fmaxcatch atlantisfmsy_inisimu

#' Change F in Atlantis.
#'
#' @description It modifies the harvest parameters file by applying a new value
#'   of F for the functional group \code{func_grp}. It looks for \code{mFC_XXX}
#'   parameter and and change the value of the vector using the new F value
#'   \code{f_test} that need to be simulated and the distribution of the fishing
#'   pressure accross fleets \code{f_prop} from the calibrated model.
#' @param func_grp The code of the Atlantis functional group for which Fmsy will
#'   be estimated.
#' @param path The directory of the batch or shell file.
#' @param harvest_filename The name of the harvest parameters file with its
#'   extension (ex: file.prm).
#' @param f_prop The proportion of the total fishing pressure applied per each
#'   fleet. The order of the fleets is the same as the one used in Atlantis.
#' @param f_test The new value of total fishing mortality to be tested in days.
#' @param fishing_para A dataframe containing data from the fishing fleet csv
#'   file plus two extra column, one to indicate if the fleet is active (1) or
#'   not (0) named \code{active_flt} and one named \code{effortmodel} containing
#'   the effortmodel option used for each fleet (output from the function
#'   \code{\link{atlantis_checkf}}).
#' @return Change value of F in the harvest parameters files.
#' @examples
#' atlantis_fchange("COD", "C:/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "AEEC_harvest.prm", c(0.2, 0.45, 0.05, 0.12, 0.08, 0.1), 0.0002739726,
#'   fishing_para)
#' atlantis_fchange("COD", "/home/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "AEEC_harvest.prm", c(0.2, 0.45, 0.05, 0.12, 0.08, 0.1), 0.0002739726,
#'   fishing_para)
#' @seealso \code{\link{atlantis_openfile}} to open a parameters file and select
#'   a parameter.
#' @export

# Function used:
# - atlantis_openfile (fileselect.R)

atlantis_fchange = function(func_grp, path, harvest_filename, f_prop, f_test, fishing_para) {
  #open harvest file.
  para <- atlantis_openfile(path, harvest_filename, paste("mFC_", func_grp, sep = ""))
  params <- para[[1]]
  idxline <- para[[2]] + 1

  #extract F parameter for species func_grp.
  mfc_sp <- as.numeric(unlist(strsplit(params[idxline], " ")))
  mfc_sp <- mfc_sp[!is.na(mfc_sp)]

  #distribute new F to fleets according to previous distribution of F in the calibrated model.
  mfc_sp <-  f_prop * f_test[1]

  params[idxline] <- paste(as.character(mfc_sp), collapse = " ")

  # writeLines(params, con = infilecon)
  write(params, file.path(path, harvest_filename))
  gc()
}

#' Change bach/shell Atlantis file.
#'
#' @description It modifies the name of the output files for the future run in
#'   the batch/shell file used to run Atlantis. The name is coded with the code
#'   of the functional group selected \code{func_grp} and the value of F that
#'   will be tested \code{f_test} (example: For \code{func_grp} = "COD" with a
#'   \code{f_test} = 0.0002739726 (0.1 y-1) the output file prefix will be
#'   MSY_COD_F010).
#' @section Warning:
#'   Only work if \code{fmax} < 10 y-1.
#' @param func_grp The code of the Atlantis functional group for which Fmsy will
#'   be estimated.
#' @param path The directory of the batch or shell file.
#' @param exe_name The name of the atlantis executable you used (ex:
#'   atlantismain, atlantisNew).
#' @param f_test The new value of total fishing mortality to be tested in days.
#' @param batch_file The name of the batch/shell file with extension you are using
#'   to run your model. If not provided, the function will search for the unique
#'   batch file in your \code{folder_path}. \strong{Default:} NULL.
#' @param msy Add msy to the outputs names or not. \strong{Default:} TRUE.
#' @param output Return \code{f_test} values without the one tested or not.
#' \strong{Default:} TRUE.
#' @return \code{f_test} The new value of total fishing mortality to be tested
#'   in days, and modify batch/shell Atlantis run file.
#' @examples
#' atlantis_bachchange("COD", "C:/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "atlantismain", 0.0002739726, "runAtlantis.bat")
#' atlantis_bachchange("COD", "/home/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "atlantisNew", 0.0002739726, "runAtlantis.sh")
#'
#' @seealso \code{\link{atlantis_openfile}} to open a parameters file and select
#'   a parameter.
#' @export

# Function used:
# - atlantis_openfile (fileselect.R)

atlantis_bachchange = function(func_grp, path, exe_name, f_test, batch_file = NULL, msy=TRUE, output = TRUE) {
  #check os used
  os <- Sys.info()['sysname']

  #find the Atlantis run bach/shell file.
  if(is.null(batch_file)){
    if (os == "Windows") {
      infilename <- list.files(path)[regexpr(".bat", list.files(path), fixed = T) != -1]
    } else {
      infilename <- list.files(path)[regexpr(".sh", list.files(path), fixed = T) != -1]
    }
  } else {
    infilename <- batch_file
  }

  para <- atlantis_openfile(path, infilename, exe_name)
  params <- para[[1]]
  idxline <- para[[2]]

  line_para <- unlist(strsplit(params[idxline], " "))

  output_name <- line_para[which(line_para == "-o") + 1]
  output_name <- unlist(strsplit(output_name, "/"))

  morta <- as.character(round(f_test[1] * 365, digits = 2))
  morta <- gsub("[.]", "", morta)

  while(nchar(morta) != 3) {
    morta <- paste(morta, "0", sep = "")
  }
  if(msy){
    output_name[length(output_name)] <- paste("MSY_", func_grp, "_F", morta, ".nc", sep = "")
  } else {
    output_name[length(output_name)] <- paste(func_grp, "_F", morta, ".nc", sep = "")
  }

  line_para[which(line_para == "-o") + 1] <- paste(output_name, collapse = "/")
  params <- paste(line_para, collapse = " ")


  write(params, file.path(path, infilename))
  f_test <- f_test[-1]

  gc()
  if(output) return(f_test)
}

#' F simulations up to the collapse of the stock.
#'
#' @description In the directory \code{folder_path}, it creates the folder
#'   \code{AtlantisMSY} and the sub-folder named after the functional group
#'   \code{func_grp} where all simulations will take place. The calibrated model
#'   files from the directory \code{model_path} is copied into this newly
#'   created folder. Then a series of values of F between F = 0 to F =
#'   \code{fmax} with an incrementation of 0.4 y-1 are tested. When the stock
#'   collapse, the function stop and return the directory of the simulation
#'   folder.
#' @section Warning:
#'   Only work if fmax < 10 y-1, see \code{\link{atlantis_bachchange}}).
#' @param func_grp The code of the Atlantis functional group for which Fmsy will
#'   be estimated.
#' @param folder_path The working directory: where the Fmsy simulations will be
#'   run and stored.
#' @param model_path The directory of the calibrated model (containing all the
#'   parameters files and one bach file. Forcing files can be stored in a direct
#'   parent directory of model_path). \strong{WARNING:} Only working if the forcing
#'   folder is in the main model directory \code{model_path} or if it is in the
#'   direct parent directory. If not please either modify this package or modify
#'   the path structure of your Atlantis input forcing parameters file.
#' @param exe_name The name of the atlantis executable you used (ex:
#'   atlantismain, atlantisNew).
#' @param harvest_filename The name of the harvest parameters file with its
#'   extension (ex: file.prm).
#' @param run_time The total duration of the simulation in days
#'   (\code{burnin_time} + 10950).
#' @param f_prop The proportion of the total fishing pressure applied per each
#'   fleet. The order of the fleets is the same as the one used in Atlantis.
#' @param fmax The maximal value of F the simulation is going to test in y-1.
#'   \strong{WARNING:} only work if \code{fmax} < 10 y-1 (see
#'   \code{\link{atlantis_bachchange}}).
#' @param fmin The minimum value of F to test in y-1 (only used when the
#'   simulation are restarted after an interruption). The default value is 0.
#' @param fishing_para A dataframe containing data from the fishing fleet csv
#'   file plus two extra columns, one to indicate if the fleet is active (1) or
#'   not (0) named \code{active_flt} and one named \code{effortmodel} containing
#'   the effortmodel option used for each fleet (output from the function
#'   \code{\link{atlantis_checkf}}).
#' @param save_time The number of days to save from the end of the simulation.
#'   \strong{Default:} NULL.
#'
#' @return \code{simu_path} The directory of the simulation folder, and run
#'   first set of simulations that are stored in the output folder.
#' @examples
#' atlantisfmsy_inisimu("COD", "C:/Atlantis/AtlantisEEC/",
#'  "C:/Atlantis/AtlantisEEC/AtlantisEECF_v3", "atlantismain",
#'   "AEEC_harvest.prm", 18250, c(0.2, 0.45, 0.05, 0.12, 0.08, 0.1), 4, 0,
#'    fishing_para, "runAtlantis.bat")
#' atlantisfmsy_inisimu("COD", "/home/Atlantis/AtlantisEEC/",
#'  "/home/Atlantis/AtlantisEEC/AtlantisEECF_v3", "atlantisNew",
#'   "AEEC_harvest.prm", 18250, c(0.2, 0.45, 0.05, 0.12, 0.08, 0.1), 4, 0,
#'    fishing_para, "runAtlantis.sh")
#'
#' \dontrun{atlantisfmsy_inisimu("COD", "C:/Atlantis/AtlantisEEC/",
#'  "C:/Atlantis/AtlantisEEC/AtlantisEECF_v3", "atlantismain",
#'   "AEEC_harvest.prm", 18250, c(0.2, 0.45, 0.05, 0.12, 0.08, 0.1), 10, 0,
#'    fishing_para)}
#' @seealso \code{\link{atlantisfmsy_modelcopy}} to create the simulation
#'   directory and copy the calibrated model inside it,
#'   \code{\link{atlantis_runtime}} to modify the runtime in the run parameters
#'   file, \code{\link{atlantis_wsummary}} to modify periodicity of summary
#'   outputs, \code{\link{atlantis_fchange}} to change the value of F in the
#'   harvest pramaters file, \code{\link{atlantis_bachchange}} to change output
#'   names in the batch/shell file, and \code{\link{atlantis_avbiomsp}} to
#'   estimate the functional group biomass.
#' @export

# Function used:
# - atlantisfmsy_modelcopy (simudir.R)
# - atlantis_runtime (ini.R)
# - atlantis_fchange (core.R)
# - atlantis_bachchange (core.R)
# - atlantis_avbiomsp (check.R)

atlantisfmsy_inisimu = function(func_grp, folder_path, model_path, exe_name, harvest_filename, run_time, f_prop, fmax, fmin, fishing_para, save_time = NULL, batch_file = NULL) {
  # check os used
  os <- Sys.info()['sysname']

  gwd_ini <- getwd()
  if (fmax >= 10) stop("You need to modify the code to use F higher than or equal to 10: problem with automatization of Atlantis output names for each simulation")

  #create simulation folder for functional group func_grp and copy model parameters files into the new directory.
  if(is.null(fmin)){
    simu_path <- atlantisfmsy_modelcopy(func_grp, folder_path, model_path, exe_name, batch_file)
    fmin <- 0
  } else {
    simu_path <- file.path(folder_path, "AtlantisMSY", func_grp) #if restart the simulation after computer shutdown.
  }

  f_test <- seq(fmin, fmax, by = 0.4) / 365 #0.4 previously.

  #change running time in Atlantis.
  atlantis_runtime(simu_path, exe_name, run_time, batch_file)

  #change stock state summary periodicity.
  atlantis_wsummary(simu_path, exe_name, batch_file)

  #change the number of days to save
  if(!is.null(save_time)) atlantis_savestart(simu_path, exe_name, run_time, save_time, batch_file)

  biom_sp <- 1
  #set F for each simulation.
  setwd(simu_path)
  while (length(f_test) > 0 & biom_sp != 0) {
    atlantis_fchange(func_grp, simu_path, harvest_filename, f_prop, f_test, fishing_para) #change F in harvest file.
    f_test <- atlantis_bachchange(func_grp, simu_path, exe_name, f_test, batch_file) #renamed the output files.

    if (os == "Windows") {
      shell(list.files()[regexpr(".bat", list.files(), fixed = T) != -1]) #run Atlantis on Windows.
    } else {
      system(paste("sh",list.files()[regexpr(".sh", list.files(), fixed = T) != -1], sep= " ")) #run Atlantis on Linux.
    }
    biom_sp <- atlantis_avbiomsp(func_grp, simu_path, exe_name, run_time, batch_file) #Check if functional group collapsed.
    gc()
  }
  setwd(gwd_ini)
  return(simu_path)
}

#' Find F with the maximum catch and estimate the next F values to simulate.
#'
#' @description It looks into the output files \code{...TOTCATCH.nc} and
#'   determine the value of F (which is coded in the output files name, see
#'   \code{\link{atlantis_bachchange}}) that result on the maximum catch for the
#'   functional group considered \code{func_grp}. Then it determines the two
#'   next values of F that need to be tested.
#' @param func_grp The code of the Atlantis functional group for which Fmsy will
#'   be estimated.
#' @param path The directory of the batch or shell file.
#' @param exe_name The name of the atlantis executable you used (ex:
#'   atlantismain, atlantisNew).
#' @param run_time The total duration of the simulation in days
#'   (\code{burnin_time} + 10950).
#' @param fmax The maximal value of F the simulation is going to test in y-1.
#'   \strong{WARNING:} only work if \code{fmax} < 10 y-1 (see
#'   \code{\link{atlantis_bachchange}}).
#' @param batch_file The name of the batch/shell file with extension you are using
#'   to run your model. If not provided, the function will search for the unique
#'   batch file in your \code{folder_path}. \strong{Default:} NULL.
#' @return It writes on the drive, in the output folder, the file Fnext_simu.txt
#'   containing the functional group code \code{func_grp} (sp), the current F
#'   value with a maximum catch (F), and the two next F value to run (Ft1 and
#'   Ft2). It writes as well the file Fmsy_XXX.txt containing the functional
#'   group code \code{func_grp} (sp), the current F value with a maximum catch
#'   (F), and the yearly average catch over the last 5 years of simulation
#'   (yield).
#' @examples
#' atlantisfmsy_fmaxcatch("COD", "C:/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "atlantismain", 18250, 4, "runAtlantis.bat")
#' atlantisfmsy_fmaxcatch("COD", "/home/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "atlantisNew", 18250, 4, "runAtlantis.sh")
#'
#' \dontrun{atlantisfmsy_fmaxcatch("COD", "C:/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "atlantismain", 18250, 10, "runAtlantis.bat")}
#' \dontrun{atlantisfmsy_fmaxcatch("COD", "/home/Atlantis/AtlantisEEC/AtlantisMSY/COD",
#'  "atlantisNew", 18250, 10, "runAtlantis.sh")}
#' @seealso \code{\link{atlantis_paraselect}} for parameters file selection,
#'   \code{\link[ncdf4]{nc_open}} to open netCDF4 files,
#'   \code{\link[ncdf4]{ncvar_get}} to extract variables from netCDF4 files,
#'   \code{\link[ncdf4]{nc_close}} to close netCDF4 files,
#'   \code{\link[stats]{aggregate}} to aggregate data and
#'   \code{\link[utils]{write.table}} to write table on the drive.
#' @export

# Function used:
# - atlantis_paraselect (fileselect.R)

atlantisfmsy_fmaxcatch = function(func_grp, path, exe_name, run_time, fmax, batch_file = NULL) {
  if (fmax >= 10) stop("You need to modify the code to use F higher than or equal to 10: problem with Atlantis output names")

  fyield <- data.frame(sp = NA, f = NA, yield = NA)

  output_path <- unlist(strsplit(atlantis_paraselect(path, exe_name, "-o", batch_file), "/"))  #search for the output directory in bach file.
  output_path <- file.path(path, output_path[-length(output_path)])
  files_name <- list.files(output_path, full.names = T)[regexpr("TOTCATCH.nc", list.files(output_path, full.names = T), fixed = T) != -1 & regexpr("MSY", list.files(output_path), fixed = T) != -1] #list of catch output files from Atlantis previous simulations.

  # setwd(paste(path, output_path, sep = "/"))
  i <- 1
  for (f in 1:length(files_name)) {
    catch.nc <- ncdf4::nc_open(files_name[f])
    time <- round(catch.nc$dim$t$vals / (3600 * 24), 0) #in days.

    catch <- ncdf4::ncvar_get(catch.nc, paste("Tot_", func_grp, "_Catch", sep = ""))  #extract total catch for the functional group considered.

    if (max(time) == run_time) {
      catch_tot <- mean(apply(catch[, time > (run_time - (365 * 5))], 2, sum)) #sum of total catch in each polygon + average catch of the 5 last years.

      if(i == 1) {
        fyield$sp <- func_grp
        fyield$f <- as.numeric(substr(tail(strsplit(files_name[f], "/")[[1]], 1), 10, 12)) / 100 #F in y-1 (extract from the file name).
        fyield$yield <- catch_tot  #total catch output.
      } else {
        fyield <- rbind(fyield, data.frame(sp = func_grp, f = as.numeric(substr(tail(strsplit(files_name[f], "/")[[1]], 1), 10, 12)) / 100, yield = catch_tot))
      }
      i <- i + 1
      ncdf4::nc_close(catch.nc)
    }
    gc()
  }

  utils::write.table(fyield, file.path(output_path, paste("Fmsy_", func_grp, ".txt", sep = "")), sep = ",", dec = ".", row.names = F) # yield for each F scenario already computed.
  fyield_final <- stats::aggregate(fyield$yield, list(sp = fyield$sp), max) #look for the maximum yield in the set of scenarios.
  names(fyield_final) <- gsub("x", "yield_max", names(fyield_final))
  fyield_ini <- merge(fyield_final, fyield)
  names(fyield_final) <- gsub("yield_max", "yield", names(fyield_final))
  fyield_final <- merge(fyield_final, fyield, x.all = T, y.all = F)

  #Dertermine the next simulation to run.
  fyield_final$pos <- which(fyield_ini$yield == fyield_ini$yield_max)
  fyield_final$f1 <- fyield$f[fyield_final$pos - 1]
  fyield_final$f2 <- fyield$f[fyield_final$pos + 1]
  fyield_final$ft1 <- round(fyield_final$f - abs(fyield_final$f1 - fyield_final$f) / 2, 2)
  fyield_final$ft2 <- round(fyield_final$f + abs(fyield_final$f2 - fyield_final$f) / 2, 2)

  fyield_final$ft2[fyield_final$ft2 >= fmax] <- NA
  fyield_ini <- fyield_final[, c("sp", "f", "ft1", "ft2")]
  fyield_ini <- fyield_ini[,!is.na(fyield_ini)]
  utils::write.table(fyield_ini, file.path(output_path, "Fnext_simu.txt"), sep = ",", dec = ".", row.names = F)
  gc()
}
rgirardi/atlantisfmsy documentation built on Aug. 30, 2022, 9:11 a.m.