Commit a501fad5 authored by Dorchies David's avatar Dorchies David
Browse files

fix: Cost calculation

parent 3ba8e5c1
Package: irmara
Title: Interactive Reservoir MAnagement Risk Assessment
Version: 0.0.1
Version: 0.0.2
Authors@R: person('David', 'Dorchies', email = 'david.dorchies@inrae.fr', role = c('cre', 'aut'))
Description: IRMaRA is a R Shiny interface providing probability of failure of flood and drought objective at key locations downstream of the 4 lakes regulating the Seine River.
License: AGPL-3
......@@ -24,5 +24,7 @@ BugReports: https://gitlab.irstea.fr/in-wop/irmara/issues
Suggests:
testthat,
knitr,
rmarkdown
rmarkdown,
zoo,
dplyr
VignetteBuilder: knitr
#' Calculate the cost for one lake
#' Calculate the cost function for one simulation (One objective at one station)
#'
#' It's the mean daily storage for a low-flow support objective and the mean daily available storage capacity for a high flow mitigation objective.
#' Uses ouput of PaChrono.txt or VObj\[1-4\].dat for the calculation.
#'
#' See [Vgest_cost.Objectives], [vgest_cost.ListVobj] and [vgest_obj.Chrono] for details.
#'
#' @param data Data source, see details
#' @param ... Other arguments passed to the methods
#'
#' @return
#' @export
#'
#' @examples
vgest_cost <- function(data, ...) {
UseMethod("vgest_cost", data)
}
#' Calculate the total cost for all the lakes
#'
#' @param dfVobj one item of the list given by [vgest_read_all()]
#' @param floodMaxCapacity the net storage capacity of the lake if the objective is high flow mitigation, set to `NULL` if the objective is low-flow support
#' For each lake, it's the mean daily storage for a low-flow support objective and the mean daily available storage capacity for a high flow mitigation objective.
#'
#' @return cost of the lake in m3/day
#' @param Vobj A [matrix] or a [data.frame] with one column by lake, in the same order as `objective$lakes`
#' @param objective one row of the dataframe given by [get_objectives()]
#'
#' @return total cost of all the lakes in m3/day
#'
vgest_cost_lakes <- function(Vobj, objective) {
lakes <- objective$lakes[[1]]
sum(
sapply(
1:nrow(lakes), function(i) {
V <- Vobj[,i]
if(objective$flood) {
V <- (lakes$max[i] - lakes$min[i]) * 1E6 - V
}
mean(V, na.rm = TRUE)
}
)
)
}
#' Calculate the total cost for one objective at one station
#'
#'
#' @param data the list given by [vgest_read_all()]
#' @param objective one row of the dataframe given by [get_objectives()]
#'
#' @return the total cost for one objective at one station in m3/day
#' @export
#'
#' @examples
......@@ -13,25 +55,20 @@
#' # This should be done after the execution of vgest for the concerned objective
#' objective <- get_objectives()[1,]
#' lResultVobj <- vgest_read_all(objective)
#' vgest_cost_lake(lResultVobj[[1]])
#' vgest_cost(lResultVobj, objective)
#' }
vgest_cost_lake <- function(dfVobj, floodMaxCapacity = NULL) {
# Calculate cost for each time step
cost <- unlist(dfVobj[,grep(pattern = "Y[0-9]{4}", names(dfVobj))])
cost <- cost[!is.na(cost)]
# For flood the cost is calculated from the free storage capacity
if(!is.null(floodMaxCapacity)) {
cost <- floodMaxCapacity - cost
}
mean(cost)
vgest_cost.ListVobj <- function(data, objective) {
lVobj <- lapply(data, function(x) {unlist(x[,grep(pattern = "Y[0-9]{4}", names(x))])})
Vobj <- matrix(unlist(lVobj), ncol = length(lVobj))
vgest_cost_lakes(Vobj, objective)
}
#' Calculate the total cost for one objective at one station
#'
#'
#' @param objective a row of the dataframe given by [get_objectives()]
#' @param lResultVobj the list given by [vgest_read_all()]
#'
#' @param data the [data.frame] given by one item of the [list] provided by [vgest_read_chrono()]
#' @param objective one row of the [data.frame] given by [get_objectives()]
#'
#' @return the total cost for one objective at one station in m3/day
#' @export
......@@ -40,19 +77,32 @@ vgest_cost_lake <- function(dfVobj, floodMaxCapacity = NULL) {
#' \donttest{
#' # This should be done after the execution of vgest for the concerned objective
#' objective <- get_objectives()[1,]
#' lResultVobj <- vgest_read_all(objective)
#' vgest_cost(objective, lResultVobj)
#' lChrono <- vgest_read_chrono(objective)
#' vgest_cost(lChrono[[1]], objective)
#' }
vgest_cost <- function(objective, lResultVobj) {
lakes <- objective$lakes[[1]]
lakes$floodMaxCapacity <- NULL
if(objective$flood) {
lakes$floodMaxCapacity <- lakes$max - lakes$min
}
lakes$Vobj <- sapply(lResultVobj, function(x) {list(x)})
sum(
apply(
lakes, 1, function(x) {vgest_cost_lake(x$Vobj, x$floodMaxCapacity)}
)
)
vgest_cost.Chrono <- function(data, objective) {
vgest_cost_lakes(as.data.frame(data[,grep("Vobj\\.[1-9]\\.", names(data))]), objective)
}
#' Calculate the total cost for one objective at one station
#'
#'
#' @param data the [list] given by [vgest_read_chrono()]
#' @param objective [data.frame] given by [get_objectives()]
#'
#' @return A list with items named [station]_[high/low]_[threshold] containing the total cost for a list of objectives and stations in m3/day
#' @export
#'
#' @examples
#' \donttest{
#' # This should be done after the execution of vgest for the concerned objective
#' objective <- get_objectives()[1,]
#' lChronos <- vgest_read_chrono(objective, distributionType = 2)
#' vgest_cost(lChronos, objective)
#' }
vgest_cost.ListChrono <- function(data, objective) {
lCosts <- lapply(1:length(data), function(i) {vgest_cost.Chrono (data[[i]], objective[i,])})
names(lCosts) <- names(data)
return(lCosts)
}
\ No newline at end of file
#' Read Chrono.txt or PaChrono.txt
#'
#'
#' See [vgest_read_chrono.default] for details.
#'
#' @param x
#' @param ...
#' @param x
#' @param ...
#'
#' @return
#' @export
......@@ -15,30 +15,30 @@ vgest_read_chrono <- function(x, ...) {
#' Read Chrono.txt or PaChrono.txt
#'
#'
#' Read the time series provided by VGEST for the forward calculation (Chrono.txt) and backward calculation (PaChrono.txt)
#'
#' The format of the file is has follow. Headers are in line 53 followed by the complete time series which is backward in time for PaChrono.txt. The columns and respected widths are:
#' - Downstream date (10)
#' - QXsous (22)
#'
#'
#' For each lake, we next have:
#' - Upstream date (14)
#' - 25 columns of flow and storage data (22)
#' - if `DistributionType %in% c(4,5)`, 1 extra column (22)
#' - 5 columns of codes (6)
#'
#'
#' And then at the end:
#' - 7 columns of flow (22)
#' - 3 columns of codes (8)
#' - 1 extra column of code for `DistributionType==5` (8)
#'
#'
#' The file is saved in RDS format for quicker reading the next time.
#'
#' @param x file to read with txt extension
#' @param nLakes number of lakes in the file
#' @param distributionType Distribution type. See [vgest_write_batch] details
#'
#'
#'
#' @return
#' @export
......@@ -64,8 +64,8 @@ vgest_read_chrono.default <- function(x, nLakes, distributionType) {
} else {
headLine <- 53
}
s <- readLines(x, n = headLine)[headLine]
headers = trimws(unlist(read.fwf(textConnection(s), widths = widths, as.is = TRUE), use.names = FALSE))
# Read file in fixed format
......@@ -86,16 +86,24 @@ vgest_read_chrono.default <- function(x, nLakes, distributionType) {
#' Read Chrono.txt or PaChrono.txt
#'
#' @param x
#' @param x
#' @param result.dir path for storing the result of vgest run. The result is stored in a subfolder named high or low (depending on \code{bFlood}) followed by the threshold
#' @param distributionType Distribution type. See [vgest_write_batch] details
#' @param backward boolean `TRUE` for reading "PaChrono.txt", `FALSE` for reading "Chrono.txt"
#'
#' @return
#' @return A list with items named [station]_[high/low]_[threshold] containing [data.frame] with the content of each file
#' @export
#'
#' @examples
vgest_read_chrono.Objectives <- function(x, result.dir, distributionType, backward = TRUE) {
#' \donttest{
#' objective <- get_objectives()[1,]
#' distributionType <- 2
#' vgest_run_store(objective,
#' 1, 1, "Q_NAT_1900-2009.txt",
#' "01/01/1900", "31/12/2009", distributionType)
#' lChrono <- vgest_read_chrono(objective, distributionType)
#' }
vgest_read_chrono.Objectives <- function(x, distributionType, result.dir = "database", backward = TRUE) {
if(backward) {
fileName <- "PaChrono.txt"
} else {
......@@ -112,6 +120,6 @@ vgest_read_chrono.Objectives <- function(x, result.dir, distributionType, backwa
)
})
names(l) <- paste(x$station, sLowHigh[x$flood + 1], x$threshold, sep = "_")
class(l) <- c("ListChrono", class(l))
return(l)
}
\ No newline at end of file
Version: 1.0
RestoreWorkspace: Default
SaveWorkspace: Default
RestoreWorkspace: No
SaveWorkspace: No
AlwaysSaveHistory: Default
EnableCodeIndexing: Yes
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment