feat: Add RunModel.GRiwrmOutputsModel

WIP: Checks are still failing

Refs #173

NAMESPACE

@@ -21,8 +21,11 @@ S3method(CreateRunOptions,InputsModel)
 S3method(CreateRunOptions,character)
 S3method(RunModel,GR)
 S3method(RunModel,GRiwrmInputsModel)
+S3method(RunModel,GRiwrmOutputsModel)
 S3method(RunModel,InputsModel)
 S3method(RunModel,Supervisor)
+S3method(extractParam,GRiwrmOutputsCalib)
+S3method(extractParam,GRiwrmOutputsModel)
 S3method(isNodeDownstream,GRiwrm)
 S3method(isNodeDownstream,GRiwrmInputsModel)
 S3method(isNodeUpstream,GRiwrm)

R/CreateInputsModel.GRiwrm.R

@@ -116,57 +116,21 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
     warning("The usage of 'Qobs' is deprecated, use 'Qinf' instead")
     Qinf <- Qobs
   }
-  varNames <- c("Precip", "PotEvap", "TempMean", "Qinf", "Qmin",
-                "TempMin", "TempMax", "ZInputs", "HypsoData", "NLayers")
-  names(varNames) <- varNames
-  lapply(varNames, function(varName) {
-    v <- get(varName)
-    if (!is.null(v)) {
-      if (is.matrix(v) || is.data.frame(v)) {
-        if (is.null(colnames(v))) {
-          stop(sprintf(
-            "'%s' must have column names",
-            varName
-          ))
-        } else if (!all(colnames(v) %in% x$id)) {
-          stop(sprintf(
-            "'%s' column names must be included in 'id's of the GRiwrm object",
-            varName
-          ), "\n",
-          sprintf("These columns are not known: %s",
-                  paste(colnames(v)[!colnames(v) %in% x$id], collapse = ", ")))
-        } else if (any(duplicated(colnames(v)))) {
-          stop(sprintf(
-            "'%s' has duplicated column names: '%s'",
-            varName,
-            paste(colnames(v)[duplicated(colnames(v))], collapse = "', '")
-          ))
-        }
-        if (!varName %in% c("ZInputs", "NLayers", "HypsoData") && nrow(v) != length(DatesR)) {
-          stop(sprintf(
-            "'%s' number of rows and the length of 'DatesR' must be equal",
-             varName
-          ))
-        }
-        if (varName %in% c("Precip", "PotEvap", "Qmin")) {
-          if (any(is.na(v))) {
-            stop(sprintf(
-              "`NA` values detected in '%s'. Missing values are not allowed in InputsModel",
-              varName
-            ))
-          }
-          if (any(v < 0)) {
-            stop(sprintf(
-              "'%s' values must be positive or nul. Missing values are not allowed in InputsModel",
-              varName
-            ))
-          }
-        }
-      } else if (!varName %in% c("ZInputs", "NLayers")) {
-        stop(sprintf("'%s' must be a matrix or a data.frame", varName))
-      }
-    }
-  })
+
+  checkInputsModelArguments(
+    x,
+    DatesR,
+    Precip = Precip,
+    PotEvap = PotEvap,
+    TempMean = TempMean,
+    Qinf = Qinf,
+    Qmin = Qmin,
+    TempMin = TempMin,
+    TempMax = TempMax,
+    ZInputs = ZInputs,
+    HypsoData = HypsoData,
+    NLayers = NLayers
+  )
 
   if (is.null(Qinf)) Qinf <- matrix(0, ncol = 0, nrow = length(DatesR))
   if (is.null(Qrelease)) Qrelease <- matrix(0, ncol = 0, nrow = length(DatesR))

R/RunModel.GRiwrmOutputsModel.R

+#' RunModel for hot restart after a previous simulation period
+#'
+#' This function allows to restart a simulation at the end of a previous
+#' simulation period. Parameters `Qinf`, `Qrelease`, and `Qmin` can be
+#' redefined for this new simulation period.
+#'
+#' @details
+#' `IndPeriod_Run` or `DatesR` must must be continuous periods starting the
+#' time step after the last simulation time step of `OutputsModel`.
+#'
+#' `Qinf`, `Qmin`, and `Qrelease` are used for overwriting the corresponding
+#' arguments provided to [CreateInputsModel.GRiwrm] for the period to be simulated.
+#' Therefore, the number of rows of these arguments must correspond to
+#' `IndPeriod_Run` or `DatesR` lengths.
+#'
+#' @inheritParams getNextTimeSteps
+#' @inheritParams Calibration
+#' @inheritParams airGR::CreateRunOptions
+#' @param DatesR (optional) [POSIXt] vector of dates of period to be used for
+#' the model run. See details
+#' @param Qinf (optional) [matrix] or [data.frame] of [numeric] containing
+#'        observed flows. It must be provided only for nodes of type "Direct
+#'        injection" and "Diversion" \[m3 per time step\].
+#'        Column names correspond to node IDs. Negative flows are abstracted from
+#'        the model and positive flows are injected to the model. See details
+#' @param Qmin (optional) [matrix] or [data.frame] of [numeric] containing
+#'        minimum flows to let downstream of a node with a Diversion \[m3 per
+#'        time step\]. Default is zero. Column names correspond to node IDs.
+#'        See details
+#' @param Qrelease (optional) [matrix] or [data.frame] of [numeric] containing
+#'        release flows by nodes using the model `RunModel_Reservoir` \[m3 per
+#'        time step\]. See details
+#'
+#' @inherit RunModel.GRiwrmInputsModel return
+#' @export
+#'
+RunModel.GRiwrmOutputsModel <- function(OutputsModel,
+                                        InputsModel,
+                                        RunOptions,
+                                        IndPeriod_Run = which(InputsModel[[1]]$DatesR %in% DatesR),
+                                        DatesR = getNextTimeSteps(OutputsModel),
+                                        Qinf = NULL,
+                                        Qrelease = NULL,
+                                        Qmin = NULL) {
+  stopifnot(inherits(OutputsModel, "GRiwrmOutputsModel"),
+            inherits(InputsModel, "GRiwrmInputsModel"),
+            inherits(RunOptions, "GRiwrmRunOptions"))
+  # Check Run Period
+  next_time_step <- getNextTimeSteps(OutputsModel)
+  next_index <- which(InputsModel[[1]]$DatesR == next_time_step)
+  if (IndPeriod_Run[1] != next_index) {
+    stop("`IndPeriod_Run` should have its first element equal to ", next_index)
+  }
+
+  # State Initiation
+  for (id in names(RunOptions)) {
+    # Run model for the sub-basin and one time step
+    RunOptions[[id]]$IniResLevels <- NULL
+    RunOptions[[id]]$IniStates <- serializeIniStates(OutputsModel[[id]]$StateEnd)
+    RunOptions[[id]]$IndPeriod_WarmUp <- 0L
+    RunOptions[[id]]$IndPeriod_Run <- IndPeriod_Run
+  }
+
+  # Inputs change
+  checkInputsModelArguments(
+    attr(InputsModel, "GRiwrm"),
+    InputsModel[[1]]$DatesR[IndPeriod_Run],
+    Qinf = Qinf,
+    Qrelease = Qrelease,
+    Qmin = Qmin
+  )
+  inputs <- list(Qinf = Qinf, Qrelease = Qrelease, Qmin = Qmin)
+  inputs[sapply(inputs, is.null)] <- NULL
+  for (inputArg in names(inputs)) {
+    input <- inputs[[inputArg]]
+    if (length(IndPeriod_Run) != nrow(input)) {
+      stop("The Argument ", inputArg,
+           " must have a number of rows identical to the lenght of `IndPeriod_Run`")
+    }
+    for (id in colnames(input)) {
+      v <- input[, id, drop = TRUE]
+      if (inputArg %in% c("Qrelease", "Qmin")) {
+        if (inputArg == "Qrelease" && !InputsModel[[id]]$isReservoir) {
+            stop("The column ", id, " of the argument `Qrelease` does not refer to a Reservoir node")
+        }
+        if (inputArg == "Qmin" && !InputsModel[[id]]$isDiversion) {
+          stop("The column ", id, " of the argument `Qmin` does not refer to a Diversion node")
+        }
+        if (is.null(InputsModel[[id]][[inputArg]])) {
+            stop("InputsModel[['", id, "']] should contain a `", inputArg, "` item")
+        }
+        InputsModel[[id]][[inputArg]][IndPeriod_Run] <- v
+      }
+      if (inputArg == "Qinf") {
+        if (is.null(InputsModel[[id]])) {
+          # Direct Injection
+          g <- attr(InputsModel, "GRiwrm")
+          id_down <- g$down[g$id == id]
+          InputsModel[[id_down]]$Qupstream[IndPeriod_Run, id] <- v
+        } else {
+          if (!InputsModel[[id]]$isDiversion) {
+            stop("The column ", id, " of the argument `Qinf` does not refer to a DirectInjection or a Diversion node")
+          }
+          InputsModel[[id]]$Qdiv[IndPeriod_Run] <- v
+        }
+      }
+    }
+  }
+
+  # Run the model
+  return(suppressMessages(
+    RunModel(
+      InputsModel,
+      RunOptions = RunOptions,
+      Param = extractParam(OutputsModel)
+    )
+  ))
+}

R/RunModel_Reservoir.R

@@ -146,6 +146,7 @@ RunModel_Reservoir <- function(InputsModel, RunOptions, Param) {
       OutputsModel$RunOptions$WarmUpQdiv_m3 <- Qdiv_m3[iWarmUp]
     }
   }
+  OutputsModel$RunOptions$Param <- Param
   iRun <- length(IndPerWarmUp) + seq(length(RunOptions$IndPeriod_Run))
   OutputsModel$Qsim_m3 <- Qsim_m3[iRun]
   OutputsModel$Vsim <- Vsim[iRun]

R/utils.Calibration.R

@@ -278,6 +278,18 @@ transferGRparams <- function(InputsModel, Param, donor, receiver, default_param
 #' @export
 #'
 extractParam <- function(x) {
-  stopifnot(inherits(x, "GRiwrmOutputsCalib"))
+  UseMethod("extractParam")
+}
+
+#' @export
+#' @rdname extractParam
+extractParam.GRiwrmOutputsCalib <- function(x) {
   lapply(x, "[[", "ParamFinalR")
 }
+
+#' @export
+#' @rdname extractParam
+extractParam.GRiwrmOutputsModel <- function(x) {
+  lapply(x, function(o) o$RunOptions$Param)
+}
+

R/utils.CreateInputsModel.R

@@ -57,3 +57,65 @@ checkQinfQrelease <- function(g, varname, Q) {
   }
   return(Q)
 }
+
+#' Check the parameters provided to CreateInputsModel.GRiwrm
+#'
+#' @param x GRiwrm
+#' @param DatesR DatesR
+#' @param ... Parameters to check
+#'
+#' @return Nothing
+#' @noRd
+#'
+checkInputsModelArguments <- function(x, DatesR, ...) {
+  dots <- list(...)
+
+  lapply(names(dots), function(varName) {
+    v <- dots[[varName]]
+    if (!is.null(v)) {
+      if (is.matrix(v) || is.data.frame(v)) {
+        if (is.null(colnames(v))) {
+          stop(sprintf(
+            "'%s' must have column names",
+            varName
+          ))
+        } else if (!all(colnames(v) %in% x$id)) {
+          stop(sprintf(
+            "'%s' column names must be included in 'id's of the GRiwrm object",
+            varName
+          ), "\n",
+          sprintf("These columns are not known: %s",
+                  paste(colnames(v)[!colnames(v) %in% x$id], collapse = ", ")))
+        } else if (any(duplicated(colnames(v)))) {
+          stop(sprintf(
+            "'%s' has duplicated column names: '%s'",
+            varName,
+            paste(colnames(v)[duplicated(colnames(v))], collapse = "', '")
+          ))
+        }
+        if (!varName %in% c("ZInputs", "NLayers", "HypsoData") && nrow(v) != length(DatesR)) {
+          stop(sprintf(
+            "'%s' number of rows and the length of 'DatesR' must be equal",
+            varName
+          ))
+        }
+        if (varName %in% c("Precip", "PotEvap", "Qmin")) {
+          if (any(is.na(v))) {
+            stop(sprintf(
+              "`NA` values detected in '%s'. Missing values are not allowed in InputsModel",
+              varName
+            ))
+          }
+          if (any(v < 0)) {
+            stop(sprintf(
+              "'%s' values must be positive or nul. Missing values are not allowed in InputsModel",
+              varName
+            ))
+          }
+        }
+      } else if (!varName %in% c("ZInputs", "NLayers")) {
+        stop(sprintf("'%s' must be a matrix or a data.frame", varName))
+      }
+    }
+  })
+}

R/utils.RunModel.R

@@ -61,7 +61,9 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
 #' @noRd
 #'
 serializeIniStates <- function(IniStates) {
-  unlist(IniStates)
+  IniStates <- unlist(IniStates)
+  IniStates[is.na(IniStates)] <- 0
+  return(IniStates)
 }
 
 

man/RunModel.GRiwrmOutputsModel.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RunModel.GRiwrmOutputsModel.R
+\name{RunModel.GRiwrmOutputsModel}
+\alias{RunModel.GRiwrmOutputsModel}
+\title{RunModel for hot restart after a previous simulation period}
+\usage{
+\method{RunModel}{GRiwrmOutputsModel}(
+  OutputsModel,
+  InputsModel,
+  RunOptions,
+  IndPeriod_Run = which(InputsModel[[1]]$DatesR \%in\% DatesR),
+  DatesR = getNextTimeSteps(OutputsModel),
+  Qinf = NULL,
+  Qrelease = NULL,
+  Qmin = NULL
+)
+}
+\arguments{
+\item{OutputsModel}{Object returned by \link{RunModel.GRiwrmInputsModel},
+\link{RunModel.Supervisor}, or \link{RunModel.GRiwrmOutputsModel}}
+
+\item{InputsModel}{[object of class \emph{InputsModel} or \emph{GRiwrmInputsModel}] see \link{CreateInputsModel}}
+
+\item{RunOptions}{[object of class \emph{RunOptions} or \emph{GRiwrmRunOptions}] see \link{CreateRunOptions}}
+
+\item{IndPeriod_Run}{[numeric] index of period to be used for the model run [-]. See details}
+
+\item{DatesR}{(optional) \link{POSIXt} vector of dates of period to be used for
+the model run. See details}
+
+\item{Qinf}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+observed flows. It must be provided only for nodes of type "Direct
+injection" and "Diversion" [m3 per time step].
+Column names correspond to node IDs. Negative flows are abstracted from
+the model and positive flows are injected to the model. See details}
+
+\item{Qrelease}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+release flows by nodes using the model \code{RunModel_Reservoir} [m3 per
+time step]. See details}
+
+\item{Qmin}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+minimum flows to let downstream of a node with a Diversion [m3 per
+time step]. Default is zero. Column names correspond to node IDs.
+See details}
+}
+\value{
+An object of class \emph{GRiwrmOutputsModel}.
+This object is a \link{list} of \emph{OutputsModel} objects produced by \link{RunModel.InputsModel}
+for each node of the semi-distributed model.
+
+It also contains the following attributes (see \link{attr}):
+\itemize{
+\item "Qm3s": a \link{data.frame} containing the dates of simulation and one column by node
+with the simulated flows in cubic meters per seconds (See \link{plot.Qm3s})
+\item "GRiwrm":  a copy of the \emph{GRiwrm} object produced by \link{CreateGRiwrm} and used for the simulation
+\item "TimeStep":  time step of the simulation in seconds
+}
+}
+\description{
+This function allows to restart a simulation at the end of a previous
+simulation period. Parameters \code{Qinf}, \code{Qrelease}, and \code{Qmin} can be
+redefined for this new simulation period.
+}
+\details{
+\code{IndPeriod_Run} or \code{DatesR} must must be continuous periods starting the
+time step after the last simulation time step of \code{OutputsModel}.
+
+\code{Qinf}, \code{Qmin}, and \code{Qrelease} are used for overwriting the corresponding
+arguments provided to \link{CreateInputsModel.GRiwrm} for the period to be simulated.
+Therefore, the number of rows of these arguments must correspond to
+\code{IndPeriod_Run} or \code{DatesR} lengths.
+}

man/extractParam.Rd

@@ -2,9 +2,15 @@
 % Please edit documentation in R/utils.Calibration.R
 \name{extractParam}
 \alias{extractParam}
+\alias{extractParam.GRiwrmOutputsCalib}
+\alias{extractParam.GRiwrmOutputsModel}
 \title{Extract calibrated parameters}
 \usage{
 extractParam(x)
+
+\method{extractParam}{GRiwrmOutputsCalib}(x)
+
+\method{extractParam}{GRiwrmOutputsModel}(x)
 }
 \arguments{
 \item{x}{A \emph{GRiwrmOutputsModel} object returned by \link{Calibration.GRiwrmInputsModel}}

tests/testthat/helper_1_RunModel.R

@@ -17,6 +17,7 @@ setupRunModel <-
            griwrm = NULL,
            Qinf = NULL,
            Qrelease = NULL,
+           Qmin = NULL,
            IsHyst = FALSE) {
 
     data(Severn)
@@ -100,8 +101,9 @@ setupRunModel <-
     InputsModel <-
       suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap,
                                          TempMean = TempMean,
-                                         Qobs = Qinf,
+                                         Qinf = Qinf,
                                          Qrelease = Qrelease,
+                                         Qmin = Qmin,
                                          IsHyst = IsHyst))
 
     # RunOptions

tests/testthat/test-RunModel.GRiwrmOutputsModel.R

+skip_on_cran()
+test_that("RunModel.GRiwrmOutputsModel works", {
+  # Setup model
+  griwrm <- CreateGRiwrm(rbind(
+    n_derived_rsrvr,
+    data.frame(
+      id = "WD",
+      down = "Dam",
+      length = 0,
+      area = NA,
+      model = NA
+    )
+  ))
+  data(Severn)
+  DatesR <-  Severn$BasinsObs[[1]]$DatesR
+  Qinf <- data.frame(
+    `54095` = rep(-1500000, length(DatesR)),
+    WD = rep(-250000, length(DatesR))
+  )
+  names(Qinf)[1] <- "54095"
+  Qrelease <- data.frame(Dam = rep(100000, length(DatesR)))
+  Qmin <- data.frame("54095" = rep(1000000, length(DatesR)))
+  names(Qmin) <- "54095"
+  e <- setupRunModel(
+    griwrm = griwrm,
+    Qinf = Qinf,
+    Qrelease = Qrelease,
+    Qmin = Qmin,
+    runRunOptions = FALSE
+  )
+  for (x in ls(e)) assign(x, get(x, e))
+
+  # Set up initial conditions
+  RunOptions <- CreateRunOptions(InputsModel, IndPeriod_WarmUp = 1:364, IndPeriod_Run = 365L)
+  Param <- c(ParamMichel[names(ParamMichel) %in% griwrm$id], list(Dam = c(100E6, 1)))
+  OM <- RunModel(InputsModel, RunOptions, Param)
+
+  # Loop over periods months periods
+  dfTS <- data.frame(
+    DatesR = DatesR,
+    yearmonth = format(DatesR, "%Y-%m")
+  )
+  Qm3s <- attr(OM, "Qm3s")
+  for(ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) {
+    OM <- RunModel(OM,
+                   InputsModel = InputsModel,
+                   RunOptions = RunOptions,
+                   IndPeriod_Run = which(dfTS$yearmonth == ym))
+    Qm3s <- rbind(Qm3s, attr(OM, "Qm3s"))
+  }
+  expect_equal(nrow(Qm3s), length(DatesR) - 364)
+})