watervalues: Generating water values for Antares

Documented in Grid_Matrix

#' Calculate grid layer matrix of Bellman values and water values
#'
#' @param area An 'antares' area.
#' @param simulation_names Names of simulations to retrieve.
#' @param simulation_values Values for the simulation.
#' @param reward_db a table contains the reward values generated by the function \code{get_Reward}
#' If it's NULL it auto calculated.
#' @param inflow a table contains the hydro storage generated by the function \code{readAntares}
#' with th option \code{hydrostorage = TRUE}. If it's NULL it auto calculated.
#' @param nb_cycle Number of times to run the algorithm.
#' @param district_name Name of the district used to store output.
#' @param mcyears MC years to consider, by default all of them.
#' @param week_53 Water values for week 53, by default 0.
#' @param method Perform mean of grids algorithm or grid of means algorithm or
#'  grid of quantile algorithm.
#' @param states_step_ratio Discretization ratio to generate steps levels
#' between the reservoir capacity and zero . Defaults to 0.05
#' @param q_ratio from 0 to 1. the probability used in quantile method
#' to determine a bellman value which q_ratio all bellman values are equal or
#' less to it. (quantile(q_ratio))
#' @param monotonic_bellman force increasing bellman values with the stock
#' level in the calculation.
#' @param test_week the week number u want to see it's calculation information
#' in the console
#' @param reservoir_capacity Reservoir capacity for the given area in MWh,
#'  if \code{NULL} (the default), value in Antares is used if available else
#'  a prompt ask the user the value to be used.
#' @param na_rm Remove NAs
#' @param correct_outliers If TRUE, outliers in Bellman values are replaced
#' by spline interpolations. Defaults to FALSE.
#' @param only_input if TRUE skip bellman values calculation and return the input
#' @param inaccessible_states Numeric in [0,1]. Tolerance of inaccessible states.
#' For example if equal to 0.9 we delete the state if this states is inaccessible by 90\% of scenarios.
#' @param until_convergence Boolean. TRUE to repeat cycle until convergence or
#'  attending the limit.
#' @param convergence_rate from 0 to 1. Define the convergence level from which
#'  we suppose that no need to continue another cycle..
#' @param convergence_criteria the value define convergence. if the difference
#' between two water values is less then this value those values are converged.
#' @param cycle_limit Define cycles limit when you are in the until_convergence mod.
#' @param pumping Boolean. True to take into account the pumping.
#' @param efficiency in [0,1]. efficient ratio of pumping.
#' @param opts
#'   List of simulation parameters returned by the function
#'   \code{antaresRead::setSimulationPath}
#' @param shiny Boolean. True to run the script in shiny mode.
#' @param stop_rate the percent from which the calculation stop. for example
#' \code{stop_rate=5} means the calculation stop if there is a week with less then
#' 5\% accessibles states.
#' @param debug_week the number of the week to open the process in debug mode
#' @param ... further arguments passed to or from other methods.
#' @return a \code{data.table}
#' @export
#'
#' @importFrom assertthat assert_that
#' @importFrom antaresRead readAntares setSimulationPath
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom dplyr left_join
#' @import data.table
#' @importFrom shinybusy show_modal_spinner remove_modal_spinner
#'


  Grid_Matrix <- function(area, simulation_names,reward_db=NULL,inflow=NULL,
                             simulation_values = NULL, nb_cycle = 1L,
                             district_name = "water values district", mcyears = NULL,
                             week_53 = 0,
                             states_step_ratio = 0.01,
                             reservoir_capacity = NULL,
                             na_rm = FALSE,
                             correct_outliers = FALSE,
                             method ,
                             only_input=FALSE,
                             q_ratio=0.5,
                             monotonic_bellman=FALSE,
                             test_week=NULL,
                             opts = antaresRead::simOptions(),
                             shiny=F,
                             inaccessible_states=1,
                             until_convergence=F,
                             convergence_rate=0.9,
                             convergence_criteria=1,
                             cycle_limit=10,
                             pumping=F,efficiency=1,stop_rate=0,
                             debug_week=54,
                        ...) {



  #----- shiny Loader


  methods <- c("mean-grid","grid-mean","quantile")
  if (!method %in% methods){
    stop("Unknown method. available methods: ('mean-grid','grid-mean','quantile') ")
  }

  assertthat::assert_that(class(opts) == "simOptions")

  # Number of weeks
  n_week <- 52

  # max hydro
  max_hydro <- get_max_hydro(area)


  if(is.null(reward_db)){
    if(!is.null(simulation_names)){
      if (is.null(simulation_values)) {
        simulation_values <- constraint_generator(area,length(simulation_names),pumping,opts)
        message(paste0("Using simulation_values: ", paste(simulation_values, collapse = ", ")))
      }}
    }else
      {simulation_names <- reward_db$simulation_names
      simulation_values <- reward_db$simulation_values
      reward_db <- as.data.table(reward_db$reward)
      }

  if (is.null(mcyears)) {
    mcyears <- opts$parameters$general$nbyears
    mcyears <- seq(0,mcyears)
  }



  # Niveau max
  {if (is.null(reservoir_capacity)) {
    niveau_max <- get_reservoir_capacity(area = area)
    if (length(niveau_max) == 0) {
      ask_niveau_max <- "Failed to retrieve reservoir capacity from Antares, please specify value (in GWh, e.g. 10000 for France):\n"
      niveau_max <- readline(prompt = ask_niveau_max)
      niveau_max <- as.numeric(niveau_max)*1000
    }
    if (length(niveau_max) == 0)
      stop("Failed to retrieve reservoir capacity, please specify it explicitly with 'reservoir_capacity'.")
  } else {
    niveau_max <- reservoir_capacity
  }}


  # synchronizing between the simulations and the states discretisation
  if ((!is.numeric(states_step_ratio))|(states_step_ratio<0)|(states_step_ratio<0))
    stop("Failed Not valid States_step_ratio, please change it between 0 and 1.")

  states_steps <- niveau_max*states_step_ratio


  # States matrix
  states <- matrix( rep(seq(from = niveau_max, to = 0, by = -states_steps), n_week + 1), byrow = FALSE, ncol = n_week + 1)


  if(is.null(efficiency))
  {
    efficiency <- getPumpEfficiency(area = area,opts = opts)
  }
  decision_space <- simulation_values
  # decision_space <- unlist(lapply(decision_space,FUN = function(x) efficiency_effect(x,efficiency)))
  decision_space <- round(decision_space)

  decimals <- 6
  {
    if(is.null(inflow)){
      tmp_name <- getSimulationNames(pattern = simulation_names[1], opts = opts)[1]
      tmp_opt <- antaresRead::setSimulationPath(path = opts$studyPath, simulation = tmp_name)
      inflow <- antaresRead::readAntares(areas = area, hydroStorage = TRUE, timeStep = "weekly", mcYears = mcyears, opts = tmp_opt)
    }
    inflow[with(inflow, order(mcYear, timeId)),]
    inflow <- inflow[, list(area, tsId = mcYear, timeId, time, hydroStorage)]
    # inflow[, timeId := gsub(pattern = "\\d{4}-w", replacement = "", x = time)]
    inflow[, timeId := as.numeric(timeId)]
    inflow <- inflow[, list(hydroStorage = sum(hydroStorage, na.rm = TRUE)), by = list(area, timeId, tsId)] # sum

  } # get the table (area,time,tsid,hydroStorage)

  options("antares" = opts)

  # Reward
  {
    if(is.null(reward_db))
    {

    reward_db <- get_Reward(simulation_names = simulation_names, district_name = district_name, opts = opts)$reward
    }

    reward_db <- reward_db[timeId %in% seq_len(n_week)]}

  # Reservoir (calque)
  {
    reservoir <- readReservoirLevels(area, timeStep = "weekly", byReservoirCapacity = FALSE, opts = opts)
    vars <- c("level_low", "level_avg", "level_high")
    reservoir[,
              (vars) := lapply(.SD, function(x) {round(x * max(states))}),
              .SDcols = vars
    ]
  }

  # preparation DATA (generate a table of weeks and years)
  watervalues <- data.table(expand.grid(weeks = seq_len(n_week+1), years = mcyears))

  # add states
  {
    statesdt <- as.data.table(states)  #convert states matrix to data.table
    statesdt <- melt(data = statesdt, measure.vars = seq_len(ncol(states)), variable.name = "weeks", value.name = "states")
    statesdt[, weeks := as.numeric(gsub("V", "", weeks))] #turn weeks to numbers V1==> 1
    statesdt[, statesid := seq_along(states), by = weeks] # add id to refer to the state
    statesdt[, states := round(states)]
  }

  # add states plus 1
  {
    statesplus1 <- copy(statesdt)
    statesplus1[, weeks := weeks - 1]
    statesplus1 <- statesplus1[, list(states_next = list(unlist(states))), by = weeks]
    statesplus1 <- dplyr::left_join(x = statesdt, y = statesplus1, by = c("weeks"), all.x = TRUE)
    watervalues <- dplyr::left_join(x = watervalues, y = statesplus1, by = "weeks", all.x = TRUE, all.y = FALSE, allow.cartesian = TRUE)
  }

  # add inflow
  watervalues <- dplyr::left_join(x = watervalues, y = inflow[, list(weeks = timeId, years = tsId, hydroStorage)], by = c("weeks", "years"))
  #at this point water values is the table containing (weeks,year,states,statesid;states_next,hydroStorage)

  #add reward
  col_names <- make.names(colnames(reward_db))
  setnames(reward_db,col_names)
  reward_l <- reward_db[, list(reward_db = list(unlist(.SD))), .SDcols =col_names[c(-1,-2)], by = list(weeks = timeId, years = mcYear)]


  watervalues <- dplyr::left_join(x = watervalues, y = reward_l, by = c("weeks","years"))


  #at this point we added the rewards for each weekly_amount

  # add reservoir
  watervalues <- dplyr::left_join(x = watervalues, y = reservoir[, list(weeks = timeId, level_low, level_high)], by = "weeks", all = TRUE)
  #here we added the lvl_high and low of the reservoir

  # add empty columns ---------------------
  watervalues$value_node <- NA_real_
  watervalues$transition <- NA_real_
  watervalues$transition_reward <- NA_real_
  watervalues$next_bellman_value <- NA_real_





  # prepare next function inputs
  {
  if (length(week_53) == 1) week_53 <- rep_len(week_53, length(states))
  next_week_values <- (week_53 * niveau_max)/2   # approximation to get initial bellman values from initial water values
  niveau_max = niveau_max
  E_max <-max_hydro$turb
  P_max <- max_hydro$pump
  max_mcyear <- length(mcyears)
  counter <- 0
  if(!pumping)P_max <- 0

  }

  ####

  if (only_input) return(watervalues)

  if(!until_convergence){
    next_week_values <- rep_len(next_week_values, nrow(watervalues[weeks==52]))

    for (n_cycl in seq_len(nb_cycle)) {

      cat("Calculating value nodes, cycle number:", n_cycl, "\n")

      pb <- txtProgressBar(min = 0, max = 51, style = 3)

      for (i in rev(seq_len(52))) { # rep(52:1, times = nb_cycle)


        temp <- watervalues[weeks==i]
        if(i==52){
        next_level_high <- watervalues$level_high[1]
        next_level_low <- watervalues$level_low[1]
        }else{
          next_level_high <- watervalues$level_high[i+1]
          next_level_low <- watervalues$level_low[i+1]
        }

        if(debug_week==i)browser()

        temp <- Bellman(Data_week=temp,
                        next_week_values_l = next_week_values,
                        decision_space=decision_space,
                        E_max=E_max,
                        P_max=P_max,
                        niveau_max=niveau_max,
                        method=method,
                        max_mcyear = max_mcyear,
                        q_ratio= q_ratio,
                        correct_outliers = correct_outliers,
                        test_week = test_week,
                        counter = i,
                        inaccessible_states=inaccessible_states,
                        next_level_high=next_level_high,
                        next_level_low=next_level_low,
                        stop_rate=stop_rate)



        if(shiny&n_cycl==1&i==52){
          shinybusy::show_modal_spinner(spin = "atom",color = "#0039f5")
        }







        # monotonic Bellman
        {

          if(monotonic_bellman){
            for (k in 1:max_mcyear){
              temp1 <- temp[weeks==i&years==k]
              m <- 0
              M <- 0

              for (j in 1:nrow(temp1)){

                if (is.na(temp1$value_node[j])|!is.finite(temp1$value_node[j])) next

                if(m==0)m <- j

                M <- j

              }


              temp1$value_node[m:M]<- temp1$value_node[m:M][order(temp1$value_node[m:M],decreasing = FALSE)]
              temp[(weeks==i&years==k),value_node :=temp1$value_node]
            }}



        }


        watervalues[weeks==i,value_node :=temp$value_node]
        watervalues[weeks==i,transition :=temp$transition]
        watervalues[weeks==i,transition_reward :=temp$transition_reward]
        watervalues[weeks==i,next_bellman_value :=temp$next_bellman_value]
        watervalues[weeks==i,accessibility :=temp$accessibility]
        watervalues[weeks==i,max_acc :=temp$max_acc]


        if (correct_outliers) {
          watervalues[weeks == i, value_node := correct_outliers(value_node), by = years]
          watervalues[weeks==i&value_node<0&is.finite(value_node),value_node:=NaN]
        }



        # next_week_values <- correct_outliers(temp$value_node)
        next_week_values <- temp$value_node
        setTxtProgressBar(pb = pb, value = 52 - i)

      }
      close(pb)
      next_week_values <- temp[weeks==1]$value_node
      watervalues[!is.finite(value_node),value_node:=NaN]
      value_nodes_dt <- build_data_watervalues(watervalues,inaccessible_states,statesdt,reservoir)

    }

  }else{

    for (n_cycl in seq_len(cycle_limit)) {

      cat("Calculating value nodes, cycle number:", n_cycl, "\n")

      pb <- txtProgressBar(min = 0, max = 51, style = 3)

      for (i in rev(seq_len(52))) { # rep(52:1, times = nb_cycle)


        temp <- watervalues[weeks==i]

          temp <- Bellman(Data_week=temp,
                          next_week_values_l = next_week_values,
                          decision_space=decision_space,
                          E_max=E_max,
                          P_max=P_max,
                          niveau_max=niveau_max,
                          method=method,
                          max_mcyear = max_mcyear,
                          q_ratio= q_ratio,
                          correct_outliers = correct_outliers,
                          test_week = test_week,
                          counter = i,
                          inaccessible_states=inaccessible_states)



        if(shiny&n_cycl==1&i==52){
          shinybusy::show_modal_spinner(spin = "atom",color = "#0039f5")
        }







        # monotonic Bellman
        {

          if(monotonic_bellman){
            for (k in 1:max_mcyear){
              temp1 <- temp[weeks==i&years==k]
              m <- 0
              M <- 0

              for (j in 1:nrow(temp1)){

                if (is.na(temp1$value_node[j])|!is.finite(temp1$value_node[j])) next

                if(m==0)m <- j

                M <- j

              }


              temp1$value_node[m:M]<- temp1$value_node[m:M][order(temp1$value_node[m:M],decreasing = FALSE)]
              temp[(weeks==i&years==k),value_node :=temp1$value_node]
            }}



        }


        watervalues[weeks==i,value_node :=temp$value_node]
        watervalues[weeks==i,transition :=temp$transition]
        watervalues[weeks==i,accessibility :=temp$accessibility]
        watervalues[weeks==i,max_acc :=temp$max_acc]



        if (correct_outliers) {
          watervalues[weeks == i, value_node := correct_outliers(value_node), by = years]
          watervalues[weeks==i&value_node<0&is.finite(value_node),value_node:=NaN]
        }



        # next_week_values <- correct_outliers(temp$value_node)
        next_week_values <- temp$value_node
        setTxtProgressBar(pb = pb, value = 52 - i)

      }
      close(pb)
      next_week_values <- temp[weeks==1]$value_node
      watervalues[!is.finite(value_node),value_node:=NaN]
      value_nodes_dt <- build_data_watervalues(watervalues,inaccessible_states,statesdt,reservoir)


      if(n_cycl>1){
        diff_vect <- last_wv -value_node_gen(watervalues,inaccessible_states,statesdt,reservoir)$vu
        convergence_value <- converged(diff_vect,conv=convergence_criteria)
        convergence_percent <- sprintf((convergence_value*100), fmt = '%#.2f')

        if(n_cycl>2){
        if(convergence_value>convergence_rate){
            cat(paste0("\033[0;42m", "        Cycle number:", n_cycl, ",    ==>",convergence_percent,"% Converged <==", "\033[0m \n"))
          break}


          if (identical(diff_vect,last_diff_vect)){
          cat(paste0("\033[0;43m", "        Cycle number:", n_cycl, ",    ==>", "Converged with:",convergence_percent,"% <==", "\033[0m \n"))
          break}


        if(convergence_value==last_conv){
          counter <- counter+1
          if(counter>1){
            cat(paste0("\033[0;43m", "        Cycle number:", n_cycl, ",    ==>", "Converged with:",convergence_percent,"% <==", "\033[0m \n"))
            break

          }
        }
        }


        cat(paste0("\033[0;41m", "        Cycle number:", n_cycl, ",    ==>",convergence_percent,"% Converged <==", "\033[0m \n"))
        last_diff_vect <- diff_vect
        last_conv <- convergence_value

      }

      last_wv <- value_node_gen(watervalues,inaccessible_states,statesdt,reservoir)$vu

      }
    }# end else

  #scenario elimination criteria Info
  elim_ratio <- 1-(min(watervalues$max_acc,na.rm = T)/max_mcyear)
  if(elim_ratio>0){
    message <- sprintf("Please select an elimnation criterea greater then %0.f%% to avoid infesable week.",elim_ratio*100)
    message(message)
  }


  if(min(value_nodes_dt,na.rm = T)<0.5)
  {
    value_nodes_dt[is.finite(vu)&vu<0.5,vu:=0.5]
    message("Minimal water value is 0.5 lesser valuers are changed automatically to 0.5.
            This change is done to assure the well functioning with Antares hydro model")
  }


  if(shiny){

    shinybusy::remove_modal_spinner()

  }


  result <- list()
  result$watervalues <- watervalues
  result$aggregated_results <- value_nodes_dt
  class(result) <- "detailled and aggregated results of watervalues calculation"
  return(result)
}