Merge branch '48-handling-correctly-initial-conditions' into 'dev'

Resolve "Handling correctly initial conditions"

Closes #48

See merge request !21

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

R/Calibration.GRiwrmInputsModel.R

@@ -19,7 +19,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
 
     if(useUpstreamQsim && any(IM$UpstreamIsRunoff)) {
       # Update InputsModel$Qupstream with simulated upstream flows
-      IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]]$IndPeriod_Run, OutputsModel)
+      IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]], OutputsModel)
     }
 
     OutputsCalib[[IM$id]] <- Calibration(

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.GRiwrmInputsModel.R

@@ -80,9 +80,8 @@ RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
 
     # Update x[[id]]$Qupstream with simulated upstream flows
     if(any(x[[id]]$UpstreamIsRunoff)) {
-      x[[id]] <- UpdateQsimUpstream(x[[id]], RunOptions[[id]]$IndPeriod_Run, OutputsModel)
+      x[[id]] <- UpdateQsimUpstream(x[[id]], RunOptions[[id]], OutputsModel)
     }
-
     # Run the model for the sub-basin
     OutputsModel[[id]] <- RunModel.InputsModel(
       x[[id]],

R/RunModel.InputsModel.R

@@ -14,5 +14,8 @@ RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) {
     # Add Qsim_m3 in m3/timestep
     OutputsModel$Qsim_m3 <- OutputsModel$Qsim * sum(x$BasinAreas) * 1e3
   }
+  if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$WarmUpQsim_m3 <- OutputsModel$WarmUpQsim * sum(x$BasinAreas) * 1e3
+  }
   return(OutputsModel)
 }

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/UpdateQsimUpstream.R

 #' Update of InputsModel$Qupstream with simulated upstream flows provided by a GRiwrmOutputsModels object
 #'
-#' @param InputsModel \[\emph{InputsModel} object\] see [airGR::CreateInputsModel] for details
-#' @param IndPeriod_Run [numeric] index of period to be used for the model run (-)
-#' @param OutputsModel \[\emph{GRiwrmOutputsModel} object\] output from [RunModel.GRiwrmInputsModel]
+#' @param InputsModel \emph{InputsModel} object. See [airGR::CreateInputsModel]
+#' @param RunOptions \emph{RunOptions} object. See [airGR::CreateRunOptions]
+#' @param OutputsModel \emph{GRiwrmOutputsModel} object provided by [RunModel.GRiwrmInputsModel].
 #'
-#' @description This function is used by [RunModel.GRiwrmInputsModel] and [Calibration.GRiwrmInputsModel] in order to provide upstream simulated flows to a node
+#' @description This function is used by [RunModel.GRiwrmInputsModel] and [Calibration.GRiwrmInputsModel]
+#' in order to provide upstream simulated flows to a node.
 #'
 #' @return `InputsModel` object with updated `Qupstream`
 #' @noRd
-UpdateQsimUpstream <- function(InputsModel, IndPeriod_Run, OutputsModel) {
+#'
+UpdateQsimUpstream <- function(InputsModel, Runoptions, OutputsModel) {
   iQ <- which(InputsModel$UpstreamIsRunoff)
   for(i in iQ) {
-      InputsModel$Qupstream[IndPeriod_Run, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]]$Qsim_m3
+      InputsModel$Qupstream[Runoptions$IndPeriod_Run, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]]$Qsim_m3
+      if (!is.null(OutputsModel[[InputsModel$UpstreamNodes[i]]]$WarmUpQsim_m3)) {
+        InputsModel$Qupstream[Runoptions$IndPeriod_WarmUp, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]]$WarmUpQsim_m3
+      }
   }
   return(InputsModel)
 }

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

-context("RunModel.Supervisor")
+# data set up
 
-# Load data
 data(Severn)
 
-# Network configuration
-nodes <- Severn$BasinsInfo[c(1,2,5), c("gauge_id", "downstream_id", "distance_downstream", "area")]
-nodes$distance_downstream <- nodes$distance_downstream # Conversion km -> m
-nodes$model <- NA
-nodes$model[1] <- "RunModel_GR4J"
-griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
-
-# InputsModel
-DatesR <- Severn$BasinsObs[[1]]$DatesR
-BasinsObs <- Severn$BasinsObs[griwrm$id]
+# Format observation
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
 PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
 PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
+
+# Set network
+nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes$distance_downstream <- nodes$distance_downstream
+nodes$model <- "RunModel_GR4J"
+griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+
+# Convert meteo data to SD (remove upstream areas)
 Precip <- ConvertMeteoSD(griwrm, PrecipTot)
 PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
-Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
+
+# Calibration parameters
+ParamMichel <- list(
+  `54057` = c(0.777323612737634, 146.867455906443, -0.100520484469163, 0.0891311328631119, 8.61545521845828),
+  `54032` = c(1.10493336411855, 1925.34891922787, -0.144768939617915, 6.3210752928989, 1.81732963586905),
+  `54001` = c(2.3483563315249, 4402.81769423169, -10903.649376187, 35.1631971451066, 17.4996745805422),
+  `54095` = c(252.52391159181, 0.0314385375487947, 55.0975274220997, 3.2928681361255),
+  `54002` = c(223.631587680546, -0.0200013333600003, 18.915846312255, 2.1981981981982),
+  `54029` = c(220.600635340659, -0.0843239834153641, 37.743931332494, 2.11619461485516)
+)
+
+# set up inputs
 InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
 
 # RunOptions
@@ -26,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,
@@ -34,21 +46,54 @@ RunOptions <- CreateRunOptions(
 )
 
 # RunModel.GRiwrmInputsModel
-Param <- list("54057" = c(0.727,  175.493,   -0.082,    0.029,    4.654))
 OM_GriwrmInputs <- RunModel(
   InputsModel,
   RunOptions = RunOptions,
-  Param = Param
+  Param = ParamMichel
 )
 
-test_that("RunModelSupervisor with no regulation should returns same results as RunModel.GRiwrmInputsModel", {
+context("RunModel.GRiwrmInputsModel")
+
+test_that("RunModel.GRiwrmInputsModel should return same result with separated warm-up", {
+  RO_WarmUp <- CreateRunOptions(
+    InputsModel = InputsModel,
+    IndPeriod_WarmUp = 0L,
+    IndPeriod_Run = IndPeriod_WarmUp
+  )
+  OM_WarmUp <- RunModel(
+    InputsModel,
+    RunOptions = RO_WarmUp,
+    Param = ParamMichel
+  )
+  RO_Run <- CreateRunOptions(
+    InputsModel = InputsModel,
+    IndPeriod_WarmUp = 0L,
+    IndPeriod_Run = IndPeriod_Run,
+    IniStates = lapply(OM_WarmUp, "[[", "StateEnd")
+  )
+  OM_Run <- RunModel(
+    InputsModel,
+    RunOptions = RO_Run,
+    Param = ParamMichel
+  )
+  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")
+
+test_that("RunModel.Supervisor with no regulation should returns same results as RunModel.GRiwrmInputsModel", {
   sv <- CreateSupervisor(InputsModel)
   OM_Supervisor <- RunModel(
     sv,
     RunOptions = RunOptions,
-    Param = Param
+    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
@@ -65,7 +110,7 @@ Qobs2 <- cbind(Qobs, matrix(data = rep(0, 2*nrow(Qobs)), ncol = 2))
 colnames(Qobs2) <- c(colnames(Qobs2)[1:6], "R1", "R2")
 InputsModel <- CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qobs2)
 
-test_that("RunModelSupervisor with two regulations that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
+test_that("RunModel.Supervisor with two regulations that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
   # Create Supervisor
   sv <- CreateSupervisor(InputsModel)
   # Function to withdraw half of the measured flow
@@ -80,12 +125,12 @@ test_that("RunModelSupervisor with two regulations that cancel each other out sh
   OM_Supervisor <- RunModel(
     sv,
     RunOptions = RunOptions,
-    Param = Param
+    Param = ParamMichel
   )
   expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
 })
 
-test_that("RunModelSupervisor with multi time steps controller, two regulations in 1 centralised controller that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
+test_that("RunModel.Supervisor with multi time steps controller, two regulations in 1 centralised controller that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
   sv <- CreateSupervisor(InputsModel, TimeStep = 10L)
   fEverything <- function(y) {
     matrix(c(y[,1]/2, -y[,1]/2), ncol = 2)
@@ -94,7 +139,7 @@ test_that("RunModelSupervisor with multi time steps controller, two regulations
   OM_Supervisor <- RunModel(
     sv,
     RunOptions = RunOptions,
-    Param = Param
+    Param = ParamMichel
   )
   expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
 })

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,])
 ```