LakeEnsemblR: Run Ensemble of Lake Models

#' douplicate section in gotm.yaml
#'@description
#' doublicate a section in the gotm.yaml file. typically used for in- or outflow
#'
#' @name doubl_yaml_sec
#' @param yaml_file filepath; to gotm.yaml
#' @param sec_name name of the section to doublicate
#' @param ap appendix for the name of the section
#' @noRd


doubl_yaml_sec <- function(yaml_file, sec_name, ap) {

  # read yaml file
  yml <- readLines(yaml_file)

  # Prevent from finding labels/keys in comments
  yml_no_comments <- unname(sapply(yml, function(x) strsplit(x, "#")[[1]][1]))

  #Find index of section
  sec_id <- paste0(sec_name, ":")
  ind_sec <- grep(paste0("\\b", sec_id), yml_no_comments)

  if(length(ind_sec) == 0){
    stop(sec_name, " not found in ", yaml_file)
  }

  # find index of next section
  find <- TRUE
  # index
  i <- 1
  # number of whitespace in section
  lws <- nchar(strsplit(yml_no_comments[ind_sec], sec_id)[[1]][1])
  while(find) {

    ws <- nchar(strsplit(yml_no_comments[ind_sec + i], "\\w+\\:")[[1]][1])
    if(ws > lws) {
      i <- i + 1
    } else {
      find <- FALSE
    }

  }

  # whole sdction
  sec_yml <- yml[ind_sec:(ind_sec + i - 1)]
  # change name by adding a number
  sec_yml[1] <- gsub(pattern = sec_name, replacement = paste0(sec_name, ap), sec_yml[1])

  # insert section
  yml_out <- c(yml[1:(ind_sec + i - 1)], sec_yml, yml[(ind_sec + i):length(yml)])

  #Write to file
  writeLines(yml_out, yaml_file)
}


#' remove section in gotm.yaml
#'@description
#' remove a section  in the gotm.yaml file. typically used for in- or outflow
#'
#' @name doubl_yaml_sec
#' @param yaml_file filepath; to gotm.yaml
#' @param sec_name name of the section to be removed
#' @noRd


rm_yaml_sec <- function(yaml_file, sec_name) {

  # read yaml file
  yml <- readLines(yaml_file)

  # Prevent from finding labels/keys in comments
  yml_no_comments <- unname(sapply(yml, function(x) strsplit(x, "#")[[1]][1]))

  #Find index of section
  sec_id <- paste0(sec_name, ":")
  ind_sec <- grep(paste0("\\b", sec_id), yml_no_comments)

  if(length(ind_sec) == 0){
    stop(sec_name, " not found in ", yaml_file)
  }

  # find index of next section
  find <- TRUE
  # index
  i <- 1
  # number of whitespace in section
  lws <- nchar(strsplit(yml_no_comments[ind_sec], sec_id)[[1]][1])
  while(find) {

    ws <- nchar(strsplit(yml_no_comments[ind_sec + i], "\\w+\\:")[[1]][1])
    if(ws > lws) {
      i <- i + 1
    } else {
      find <- FALSE
    }

  }

  # remove section
  yml_out <- yml[c(1:(ind_sec - 1), (ind_sec + i):length(yml))]

  #Write to file
  writeLines(yml_out, yaml_file)
}

#' check naming convention for in/outflows
#'@description
#'check if the headder in in/outflow files follow the naming convention
#'
#' @name chk_names_flow
#' @param flow data.frame of the read in flow data
#' @param num_flows number of in/outflows
#' @param file_n file name
#' @noRd
chk_names_flow <- function(flow, num_flows, file_n) {

  # colnames of the data
  cln <- colnames(flow)
  # remove numbers if multiple in/outflows are there
  if(num_flows > 1) {
    cln <- gsub("(\\w+)\\_\\d+\\>", "\\1", cln)
  }
  # test if names are right
  chck_flow <- sapply(list(cln), function(x) x %in% lake_var_dic$standard_name)
  if(any(!chck_flow)){
    chck_flow[which(chck_flow == FALSE)] <- sapply(list(cln[which(
      chck_flow == FALSE)]), function(x) x %in% met_var_dic$standard_name)

    if(any(!chck_flow)){
      stop(paste0("Colnames of",  file_n, " file are not in standard notation! ",
           "They should be one of: \ndatetime\nFlow_metersCubedPerSecond\n",
           "Water_Temperature_celsius\nSalinity_practicalSalinityUnits"))
    }
  }
}

#' scale in/outflows
#'@description
#' scale the flow of in/outflow data.frames
#'
#' @name scale_flow
#' @param flow data.frame of the read in flow data
#' @param num_flows number of in/outflows
#' @param scale_param scaling parameters
#' @noRd
scale_flow <- function(flow, num_flows, scale_param) {

  if(num_flows == 1) {
    flow[["Flow_metersCubedPerSecond"]] <- flow[["Flow_metersCubedPerSecond"]] *
      scale_param
  } else if(num_flows > 1) {
    for (i in 1:num_flows) {
      flow[[paste0("Flow_metersCubedPerSecond_", i)]] <-
        flow[[paste0("Flow_metersCubedPerSecond_", i)]] * scale_param[i]
    }
  }
  return(flow)
}

#' function to generate a vector of strings for the Simstrat flow files
#'@description
#' returns a vector of strings that can be written using writeLines
#'
#' @name format_flow_simstrat
#' @param flow_file data.frame of the flow data, generated by format_inflow
#'  or format_outflow
#' @param levels vector of depths of the flows
#' @param surf_flow boolean; vector what flows are surface flows
#' @param in_or_out character; 'inflow' or 'outflow'
#' @param type character; 'discharge', 'temp', or 'salt'
#' @param sim_par filepath; path to the simstrat .par file
#' @noRd

# Based on tests with Simstrat v3.0.1 and the manual:
## Surface flows (i.e. relative to current wlvl) at the same depth will need to
## be summed. Deep flows (i.e. relative to initial wlvl) work differently, and
## do not need to be summed. Temp and salt need to match in terms of columns.
## InflowMode 2 (density-driven inflows) makes that the unit of inflows (m3) is
## different from that of the outflows (m2). We could not create a constant
## water balance using this option, and therefore we do not support this option
## yet.

format_flow_simstrat <- function(flow_file, levels, surf_flow, in_or_out, type,
                                 sim_par = "simstrat.par"){
  inflow_mode <- get_json_value(sim_par, "ModelConfig", "InflowMode")
  
  # In case you want to generate a 0 outflow file despite using inflows
  if(is.null(flow_file)){
    inflow_mode <- 0
  }
  
  ## Error checking
  # InflowMode != 1
  if(inflow_mode == 2L){
    stop("InflowMode 2 not yet supported by format_flow_simstrat...")
  }else if(!(inflow_mode %in% c(0L, 1L))){
    stop("Invalid inflow_mode argument for format_flow_simstrat...")
  }
  
  # flow depths
  if(length(levels) != length(surf_flow)){
    stop("Not the same length of levels and surf_flow!")
  }
  
  # in_or_out
  if(!(in_or_out %in% c("inflow", "outflow"))){
    stop("Invalid argument in_or_out: ", in_or_out)
  }
  
  # type
  if(!(type %in% c("discharge", "temp", "salt"))){
    stop("Invalid argument type: ", type)
  }
  if(in_or_out == "outflow" & type %in% c("temp", "salt")){
    stop("Simstrat outflow does not need temperature or salt!")
  }
  
  ## Line 1
  if(in_or_out == "inflow"){
    if(type == "discharge"){
      flow_line_1 <- "Time [d]\tQ_in [m2/s]"
    }else if(type == "temp"){
      flow_line_1 <- "Time [d]\tT_in [degC]"
    }else if(type == "salt"){
      flow_line_1 <- "Time [d]\tS_in [perMille]"
    }
  }else if(in_or_out == "outflow" & type == "discharge"){
    flow_line_1 <- "Time [d]\tQ_out [m3/s]"
  }
  
  if(inflow_mode == 0L){
    flow_line_2 <- "1"
    flow_line_3 <- "-1\t0.00"
    start_sim <- get_json_value(sim_par, "Simulation", "Start d")
    end_sim <- get_json_value(sim_par, "Simulation", "End d")
    flow_line_4 <- paste0(start_sim, "\t", 0.000)
    flow_line_5 <- paste0(end_sim, "\t", 0.000)
    
    return(c(flow_line_1, flow_line_2, flow_line_3, flow_line_4, flow_line_5))
  }
  
  ## Line 2
  num_deep_flows <- length(levels) - sum(surf_flow)
  num_surf_flows <- sum(surf_flow)
  flow_line_2 <- paste0(num_deep_flows * 3, "\t",
                        ifelse(num_surf_flows == 0L, 0, 2))
  
  ## Line 3
  flow_line_3 <- "-1"
  for(i in levels[!surf_flow]){
    flow_line_3 <- paste(flow_line_3, i - 1, i, i + 1, sep = "\t")
  }
  if(any(surf_flow)){
    flow_line_3 <- paste(flow_line_3, "-1.00\t0.00", sep = "\t")
  }
  
  ## Line 4
  flow_line_4 <- seq_len(length(flow_file$datetime))
  if(type == "discharge"){
    col_header <- "Flow_metersCubedPerSecond_"
  }else if(type == "temp"){
    col_header <- "Water_Temperature_celsius_"
  }else if(type == "salt"){
    col_header <- "Salinity_practicalSalinityUnits_"
  }
  
  # first the deep flows
  # deep flow discharges are bordered by 0s, temp and salt not
  for(i in seq_len(length(levels))){
    if(surf_flow[i]) next
    if(type == "discharge"){
      flow_line_4 <- paste(flow_line_4,
                           rep(0, length(flow_file$datetime)),
                           flow_file[, paste0(col_header, i)],
                           rep(0, length(flow_file$datetime)),
                           sep = "\t")
    }else{
      flow_line_4 <- paste(flow_line_4,
                           flow_file[, paste0(col_header, i)],
                           flow_file[, paste0(col_header, i)],
                           flow_file[, paste0(col_header, i)],
                           sep = "\t")
    }
  }
  
  # then the surface flows
  # all surface flows need to be summed (discharge) or weight-averaged (temp/salt)!
  if(any(surf_flow)){
    inds <- which(surf_flow)
    
    if(type == "discharge"){
      vals_to_enter <- rowSums(data.frame(flow_file[, paste0(col_header,
                                                             inds)]))
    }else{
      # Weighted average: Sum(Flow_i * Temp/Salt_i) / Total_Flow
      total_discharge <- rowSums(data.frame(flow_file[,
                                                      paste0("Flow_metersCubedPerSecond_",
                                                             inds)]))
      value_times_discharge <- data.frame(flow_file[, paste0(col_header, inds)]) *
        data.frame(flow_file[, paste0("Flow_metersCubedPerSecond_", inds)])
      
      vals_to_enter <- rowSums(value_times_discharge) / total_discharge
    }
    # if there were 0s remove NAs
    vals_to_enter[is.na(vals_to_enter)] <- 0
    # Enter values twice; for -1.00 and 0.00 m below surface, as the unit is
    #   m2/s; degC*m2/s; or perMille*m2/s. 
    flow_line_4 <- paste(flow_line_4,
                         vals_to_enter,
                         vals_to_enter,
                         sep = "\t")
  }
  
  return(c(flow_line_1, flow_line_2, flow_line_3, flow_line_4))
}