Resolve "Handling correctly initial conditions"

Closes #48 (closed)

R/CreateRunOptions.GRiwrmInputsModel.R

+#' @param IniStates (optional) [numeric] object or [list] of [numeric] object of class \emph{IniStates}, see [airGR::CreateIniStates] for details
 #' @rdname CreateRunOptions
 #' @export
-CreateRunOptions.GRiwrmInputsModel <- function(InputsModel, ...) {
+CreateRunOptions.GRiwrmInputsModel <- function(InputsModel, IniStates = NULL, ...) {
 
   RunOptions <- list()
   class(RunOptions) <- append(class(RunOptions), "GRiwrmRunOptions")
 
-  for(InputsModelBasin in InputsModel) {
-    RunOptions[[InputsModelBasin$id]] <- CreateRunOptions(InputsModel = InputsModelBasin, ...)
+  for(id in names(InputsModel)) {
+    RunOptions[[id]] <- CreateRunOptions(InputsModel = InputsModel[[id]], IniStates = IniStates[[id]], ...)
   }
   return(RunOptions)
 }

R/CreateRunOptions.R

@@ -7,9 +7,11 @@
 #'
 #' @details See [airGR::CreateRunOptions] documentation for a complete list of arguments.
 #'
+#' If `InputsModel` argument is a \emph{GRiwrmInputsModel} object, `IniStates` must be a list of [numeric] object of class \emph{IniStates} with one item per modelled sub-catchment.
+#'
 #' With a \emph{GRiwrmInputsModel} object, all arguments are applied on each sub-catchments of the network.
 #'
-#' @return Depending on the class of `InputsModel` argument (respectively `InputsModel` and `GRiwrmInputsModel` object), the returned value is respectively:
+#' @return Depending on the class of `InputsModel` argument (respectively \emph{InputsModel} and \emph{GRiwrmInputsModel} object), the returned value is respectively:
 #' - a `RunOptions` object (See [airGR::CreateRunOptions])
 #' - a `GRiwrmRunOptions` object which is a [list] of `RunOptions` object with one item per modelled sub-catchment
 #'

R/RunModel.GR.R

@@ -11,6 +11,7 @@ RunModel.GR <- function(x, RunOptions, Param, ...) {
   if (inherits(x, "SD")) {
     # Lag model take one parameter at the beginning of the vector
     iFirstParamRunOffModel <- 2
+    RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - 1
   } else {
     # All parameters
     iFirstParamRunOffModel <- 1

R/RunModel.Supervisor.R

@@ -43,16 +43,12 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   # Save Qsim for step by step simulation
   QcontribDown <- do.call(
     cbind,
-    lapply(x$OutputsModel, function(OM) {
-      OM$Qsim
-    })
+    lapply(x$OutputsModel, "[[", "Qsim")
   )
 
   Qsim_m3 <- do.call(
     cbind,
-    lapply(x$OutputsModel, function(OM) {
-      OM$Qsim_m3
-    })
+    lapply(x$OutputsModel, "[[", "Qsim_m3")
   )
 
   # Initialisation of model states by running the model with no supervision on warm-up period
@@ -62,16 +58,17 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
     RunOptionsWarmUp[[id]]$IndPeriod_WarmUp <- 0L
     RunOptionsWarmUp[[id]]$Outputs_Sim <- c("StateEnd", "Qsim")
   }
-  x$OutputsModel <- suppressMessages(
+  OM_WarmUp <- suppressMessages(
     RunModel.GRiwrmInputsModel(x$InputsModel,
                                RunOptions = RunOptionsWarmUp,
                                Param = Param)
   )
 
-  # Adapt RunOptions to step by step simulation
+  # Adapt RunOptions to step by step simulation and copy states
   for(id in getSD_Ids(x$InputsModel)) {
     RunOptions[[id]]$IndPeriod_WarmUp <- 0L
-    RunOptions[[id]]$Outputs_Sim <- "StateEnd"
+    RunOptions[[id]]$Outputs_Sim <- c("Qsim_m3", "StateEnd")
+    x$OutputsModel[[id]]$StateEnd <- serializeIniStates(OM_WarmUp[[id]]$StateEnd)
   }
 
   # Loop over time steps with a step equal to the supervision time step
@@ -88,7 +85,7 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
     for(id in getSD_Ids(x$InputsModel)) {
       # Run the SD model for the sub-basin and one time step
       RunOptions[[id]]$IndPeriod_Run <- iTS
-      RunOptions[[id]]$IniStates <- unlist(x$OutputsModel[[id]]$StateEnd)
+      RunOptions[[id]]$IniStates <- serializeIniStates(x$OutputsModel[[id]]$StateEnd)
       x$OutputsModel[[id]] <- RunModel.SD(
         x$InputsModel[[id]],
         RunOptions = RunOptions[[id]],

R/utils.R

@@ -151,3 +151,16 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
   class(dfQsim) <- c("Qm3s", class(dfQsim)) # For S3 methods
   return(dfQsim)
 }
+
+#' Convert IniStates list into a vector
+#'
+#' @param IniStates see [CreateIniStates]
+#'
+#' @return A vector as in `RunOptions$IniStates`
+#' @noRd
+#'
+serializeIniStates <- function(IniStates) {
+  IniStates <- unlist(IniStates)
+  IniStates[is.na(IniStates)] <- 0
+  return(IniStates)
+}

man/CreateRunOptions.Rd

@@ -7,7 +7,7 @@
 \alias{CreateRunOptions}
 \title{Creation of the CalibOptions object}
 \usage{
-\method{CreateRunOptions}{GRiwrmInputsModel}(InputsModel, ...)
+\method{CreateRunOptions}{GRiwrmInputsModel}(InputsModel, IniStates = NULL, ...)
 
 \method{CreateRunOptions}{InputsModel}(InputsModel, ...)
 
@@ -16,10 +16,12 @@ CreateRunOptions(InputsModel, ...)
 \arguments{
 \item{InputsModel}{object of class \emph{InputsModel} or \emph{GRiwrmInputsModel}. See \link{CreateInputsModel} for details}
 
+\item{IniStates}{(optional) \link{numeric} object or \link{list} of \link{numeric} object of class \emph{IniStates}, see \link[airGR:CreateIniStates]{airGR::CreateIniStates} for details}
+
 \item{...}{arguments passed to \link[airGR:CreateRunOptions]{airGR::CreateRunOptions}, see details}
 }
 \value{
-Depending on the class of \code{InputsModel} argument (respectively \code{InputsModel} and \code{GRiwrmInputsModel} object), the returned value is respectively:
+Depending on the class of \code{InputsModel} argument (respectively \emph{InputsModel} and \emph{GRiwrmInputsModel} object), the returned value is respectively:
 \itemize{
 \item a \code{RunOptions} object (See \link[airGR:CreateRunOptions]{airGR::CreateRunOptions})
 \item a \code{GRiwrmRunOptions} object which is a \link{list} of \code{RunOptions} object with one item per modelled sub-catchment
@@ -31,6 +33,8 @@ This function can be used either for a catchment (with an \emph{InputsModel} obj
 \details{
 See \link[airGR:CreateRunOptions]{airGR::CreateRunOptions} documentation for a complete list of arguments.
 
+If \code{InputsModel} argument is a \emph{GRiwrmInputsModel} object, \code{IniStates} must be a list of \link{numeric} object of class \emph{IniStates} with one item per modelled sub-catchment.
+
 With a \emph{GRiwrmInputsModel} object, all arguments are applied on each sub-catchments of the network.
 }
 \examples{

tests/testthat/test-RunModel.R

@@ -38,7 +38,7 @@ IndPeriod_Run <- seq(
   length(InputsModel[[1]]$DatesR) - nTS + 1,
   length(InputsModel[[1]]$DatesR)
 )
-IndPeriod_WarmUp = seq(IndPeriod_Run[1]-366,IndPeriod_Run[1]-1)
+IndPeriod_WarmUp = seq(IndPeriod_Run[1]-365,IndPeriod_Run[1]-1)
 RunOptions <- CreateRunOptions(
   InputsModel = InputsModel,
   IndPeriod_WarmUp = IndPeriod_WarmUp,
@@ -68,18 +68,18 @@ test_that("RunModel.GRiwrmInputsModel should return same result with separated w
   RO_Run <- CreateRunOptions(
     InputsModel = InputsModel,
     IndPeriod_WarmUp = 0L,
-    IndPeriod_Run = IndPeriod_Run
+    IndPeriod_Run = IndPeriod_Run,
+    IniStates = lapply(OM_WarmUp, "[[", "StateEnd")
   )
-  for(id in names(RO_Run)) {
-    RO_Run[[id]]$IniResLevels <- NULL
-    RO_Run[[id]]$IniStates <- airGRiwrm:::serializeIniStates(OM_WarmUp[[id]]$StateEnd)
-  }
   OM_Run <- RunModel(
     InputsModel,
     RunOptions = RO_Run,
     Param = ParamMichel
   )
-  expect_equal(OM_GriwrmInputs[["54057"]]$Qsim, OM_Run[["54057"]]$Qsim)
+  lapply(griwrm$id, function(id) {
+    # The 2 exclamation marks are for seeing the id in the test result (See ?quasi_label)
+    expect_equal(OM_GriwrmInputs[[!!id]]$Qsim, OM_Run[[!!id]]$Qsim)
+  })
 })
 
 context("RunModel.Supervisor")
@@ -91,7 +91,9 @@ test_that("RunModel.Supervisor with no regulation should returns same results as
     RunOptions = RunOptions,
     Param = ParamMichel
   )
-  expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+  lapply(griwrm$id, function(id) {
+    expect_equal(OM_Supervisor[[!!id]]$Qsim, OM_GriwrmInputs[[!!id]]$Qsim)
+  })
 })
 
 # Add 2 nodes to the network

vignettes/V02_Calibration_SD_model.Rmd

@@ -124,7 +124,7 @@ The **airGR** calibration process is applied on each node of the `GRiwrm` networ
 ```{r Calibration}
 OutputsCalib <- suppressWarnings(
   Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions))
-ParamMichel <- sapply(griwrm$id, function(x) {OutputsCalib[[x]]$Param})
+ParamMichel <- sapply(OutputsCalib, "[[", "ParamFinalR")
 ```
 
 ## Run the model with the optimised model parameters
@@ -147,5 +147,5 @@ The resulting flows of each node in m<sup>3</sup>/s is directly available and ca
 
 ```{r}
 Qm3s <- attr(OutputsModels, "Qm3s")
-plot(Qm3s[1:150,]) #on ne voit pas la courbe jaune je crois, est-ce normal ?
+plot(Qm3s[1:150,])
 ```

DESCRIPTION

@@ -31,4 +31,4 @@ BugReports: https://gitlab.irstea.fr/in-wop/airGRiwrm/-/issues/
 Depends:
     R (>= 2.10)
 Remotes:
-    url::https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/archive/dev/airgr-dev.zip
+    url::https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/archive/132-runmodel_lag-add-warmupqsim-to-outputsmodel/airgr-132-runmodel_lag-add-warmupqsim-to-outputsmodel.zip