R/Calibration_Michel.R

@@ -5,7 +5,8 @@ Calibration_Michel <- function(InputsModel,
                                FUN_MOD,
                                FUN_CRIT,           # deprecated
                                FUN_TRANSFO = NULL,
-                               verbose = TRUE) {
+                               verbose = TRUE,
+                               ...) {
 
 
   FUN_MOD  <- match.fun(FUN_MOD)
@@ -145,7 +146,7 @@ Calibration_Michel <- function(InputsModel,
     }
     ##Model_run
     Param <- CandidatesParamR[iNew, ]
-    OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD)
+    OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
 
     ##Calibration_criterion_computation
     OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
@@ -294,7 +295,7 @@ Calibration_Michel <- function(InputsModel,
     for (iNew in 1:nrow(CandidatesParamR)) {
       ##Model_run
       Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD)
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
       ##Calibration_criterion_computation
       OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
       if (!is.na(OutputsCrit$CritValue)) {
@@ -355,7 +356,7 @@ Calibration_Michel <- function(InputsModel,
       CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
       ##Model_run
       Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD)
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
       ##Calibration_criterion_computation
       OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
       if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {

R/CreateCalibOptions.R

@@ -76,6 +76,13 @@ CreateCalibOptions <- function(FUN_MOD,
     ObjectClass <- c(ObjectClass, "CemaNeigeGR6J")
     BOOL <- TRUE
   }
+  if (identical(FUN_MOD, RunModel_Lag)) {
+    ObjectClass <- c(ObjectClass, "Lag")
+    if (IsSD) {
+      stop("RunModel_Lag should not be used with 'isSD=TRUE'")
+    }
+    BOOL <- TRUE
+  }
   if (IsHyst) {
     ObjectClass <- c(ObjectClass, "hysteresis")
   }
@@ -136,6 +143,9 @@ CreateCalibOptions <- function(FUN_MOD,
         FUN_GR <- TransfoParam_CemaNeige
       }
     }
+    if (identical(FUN_MOD, RunModel_Lag)) {
+      FUN_GR <- TransfoParam_Lag
+    }
     if (is.null(FUN_GR)) {
       stop("'FUN_GR' was not found")
       return(NULL)
@@ -151,7 +161,7 @@ CreateCalibOptions <- function(FUN_MOD,
       FUN_LAG <- TransfoParam_Lag
     }
     ## set FUN_TRANSFO
-    if (sum(ObjectClass %in% c("GR4H", "GR5H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige")) > 0) {
+    if (sum(ObjectClass %in% c("GR4H", "GR5H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige", "Lag")) > 0) {
       if (!IsSD) {
         FUN_TRANSFO <- FUN_GR
       } else {
@@ -292,6 +302,10 @@ CreateCalibOptions <- function(FUN_MOD,
   if ("CemaNeigeGR6J" %in% ObjectClass) {
     NParam <- 8
   }
+  if ("Lag" %in% ObjectClass) {
+    NParam <- 1
+  }
+
   if (IsHyst) {
     NParam <- NParam + 2
   }

R/RunModel.R

-RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD) {
-  
+RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, ...) {
+
   FUN_MOD <- match.fun(FUN_MOD)
-  
-  if (inherits(InputsModel, "SD")) {
+
+  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
     # Lag model take one parameter at the beginning of the vector
     iFirstParamRunOffModel <- 2
   } else {
     # All parameters
     iFirstParamRunOffModel <- 1
   }
-  
+
   OutputsModel <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
-                          Param = Param[iFirstParamRunOffModel:length(Param)])
-  
-  if (inherits(InputsModel, "SD")) {
-    InputsModel$OutputsModel <- OutputsModel
-    OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1])
+                          Param = Param[iFirstParamRunOffModel:length(Param)], ...)
+
+  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
+    OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel)
   }
   return(OutputsModel)
-}
\ No newline at end of file
+}

R/RunModel_Lag.R

-RunModel_Lag <- function(InputsModel, RunOptions, Param) {
+RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
   NParam <- 1
 
   ##Arguments_check
@@ -17,19 +17,23 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
   if (sum(!is.na(Param)) != NParam) {
     stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
   }
-  if (is.null(InputsModel$OutputsModel)) {
-    stop("'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment")
-  }
-  if (is.null(InputsModel$OutputsModel$Qsim)) {
-    stop("'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
+  if (inherits(QcontribDown, "OutputsModel")) {
+    if (is.null(QcontribDown$Qsim)) {
+      stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
+    }
+    OutputsModel <- QcontribDown
+    OutputsModel$QsimDown <- OutputsModel$Qsim
+  } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
+    OutputsModel <- list()
+    class(OutputsModel) <- c("OutputsModel", class(OutputsModel))
+    OutputsModel$QsimDown <- QcontribDown
+  } else {
+    stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
   }
-  if (sum(!is.na(InputsModel$OutputsModel$Qsim)) != length(RunOptions$IndPeriod_Run)) {
-    stop("'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA")
+  if (length(OutputsModel$QsimDown) != length(RunOptions$IndPeriod_Run)) {
+    stop("Time series in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
   }
 
-  OutputsModel <- InputsModel$OutputsModel
-  OutputsModel$QsimDown <- OutputsModel$Qsim
-
   if (inherits(InputsModel, "hourly")) {
     TimeStep <- 60 * 60
   } else if (inherits(InputsModel, "daily")) {

man/RunModel_Lag.Rd

@@ -14,26 +14,29 @@ Function which performs a single run for the Lag model over the test period.
 
 
 \usage{
-RunModel_Lag(InputsModel, RunOptions, Param)
+RunModel_Lag(InputsModel, RunOptions, Param, QcontribDown)
 }
 
 
 \arguments{
-\item{InputsModel}{[object of class \emph{InputsModel}] created with SD model inputs, see \code{\link{CreateInputsModel}} for details. The object should also contain a key \emph{OutputsModel} of class \code{\link{CreateInputsModel}} coming from the simulation of the downstream subcatchment runoff.}
+  \item{InputsModel}{[object of class \emph{InputsModel}] created with SD model inputs, see \code{\link{CreateInputsModel}} for details. The object should also contain a key \emph{OutputsModel} of class \code{\link{CreateInputsModel}} coming from the simulation of the downstream subcatchment runoff.}
 
-\item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
+  \item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
+
+  \item{Param}{[numeric] vector of 1 parameter
+    \tabular{ll}{
+      Velocity \tab Mean flow velocity [m/s]
+    }
+  }
+  \item{QcontribDown}{[numeric] vector or [OutputsModel] containing the time series of the runoff contribution of the downstream sub-basin}
 
-\item{Param}{[numeric] vector of 1 parameter
-  \tabular{ll}{
-    Velocity \tab Mean flow velocity [m/s]
-  }}
 }
 
 
 \value{
 [list] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details.
 
-The list value contains an extra item named \code{QsimDown} which is a copy of \code{InputsModel$OutputsModel$Qsim}, a numeric series of simulated discharge [mm/time step] related to the runoff contribution of the downstream sub-catchment.
+The list value contains an extra item named \code{QsimDown} which is a copy of the runoff contribution of the downstream sub-basin contained in argument \code{QcontribDown} in [mm/time step].
 }
 
 
@@ -88,12 +91,11 @@ OutputsModelDown <- RunModel_GR4J(InputsModel = InputsModel,
 ## with a delay of 2 days for 150 km, the flow velocity is 75 km per day
 Velocity <- (LengthHydro * 1e3 / 2) / (24 * 60 * 60) ## Conversion km/day -> m/s
 
-## add the output of the precipitation-runoff model in the InputsModel object
-InputsModel$OutputsModel <- OutputsModelDown
-
 ## run the lag model which routes precipitation-runoff model and upstream flows
 OutputsModel <- RunModel_Lag(InputsModel = InputsModel,
-                             RunOptions = RunOptions, Param = Velocity)
+                             RunOptions = RunOptions,
+                             Param = Velocity,
+                             QcontribDown = OutputsModelDown)
 
 ## results preview of comparison between naturalised (observed) and influenced flow (simulated)
 plot(OutputsModel, Qobs = OutputsModel$QsimDown)

tests/testthat/test-RunModel_Lag.R

@@ -39,10 +39,19 @@ RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                                 InputsModel = InputsModel,
                                                 IndPeriod_Run = Ind_Run))
 
-test_that("InputsModel parameter should contain an OutputsModel key", {
+test_that("QcontribDown parameter should be a numeric vector or an OutputModel object", {
+  regexp = "'QcontribDown' must be a numeric vector or a 'OutputsModel' object"
   expect_error(
-    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
-    regexp = "'InputsModel' should contain an 'OutputsModel' key"
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1, QcontribDown = "A"),
+    regexp = regexp
+  )
+  expect_error(
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1, QcontribDown = NULL),
+    regexp = regexp
+  )
+  expect_error(
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1, QcontribDown = matrix(1, ncol = 1)),
+    regexp = regexp
   )
 })
 
@@ -52,31 +61,36 @@ OutputsGR4JOnly <- RunModel_GR4J(InputsModel = InputsModel,
                                  RunOptions = RunOptions,
                                  Param = Param)
 
-test_that("InputsModel$OutputsModel should contain a Qsim key", {
-  InputsModel$OutputsModel <- OutputsGR4JOnly
-  InputsModel$OutputsModel$Qsim <- NULL
+test_that("QcontribDown should contain a Qsim key", {
+  QcontribDown <- OutputsGR4JOnly
+  QcontribDown$Qsim <- NULL
   expect_error(
-    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1, QcontribDown = QcontribDown),
     regexp = "should contain a key 'Qsim'"
   )
 })
 
-test_that("'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run'", {
-  InputsModel$OutputsModel <- OutputsGR4JOnly
-  InputsModel$OutputsModel$Qsim <- c(InputsModel$OutputsModel$Qsim, 0)
+test_that("'QcontribDown$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run'", {
+  QcontribDown <- OutputsGR4JOnly
+  QcontribDown$Qsim <- c(QcontribDown$Qsim, 0)
   expect_error(
-    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1, QcontribDown = QcontribDown),
     regexp = "should have the same lenght as"
   )
 })
 
-test_that("'InputsModel$OutputsModel$Qsim' should contain no NA'", {
-  InputsModel$OutputsModel <- OutputsGR4JOnly
-  InputsModel$OutputsModel$Qsim[10L] <- NA
-  expect_error(
-    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
-    regexp = "contain no NA"
-  )
+test_that("RunModel(FUN=RunModel_Lag) should give same result as RunModel_Lag", {
+  QcontribDown <- OutputsGR4JOnly
+  Output_RunModel_Lag <- RunModel_Lag(InputsModel = InputsModel,
+                                      RunOptions = RunOptions,
+                                      Param = 1,
+                                      QcontribDown = QcontribDown)
+  Output_RunModel <- RunModel(InputsModel = InputsModel,
+                              RunOptions = RunOptions,
+                              Param = 1,
+                              FUN_MOD = RunModel_Lag,
+                              QcontribDown = QcontribDown)
+  expect_equal(Output_RunModel, Output_RunModel_Lag)
 })
 
 test_that("'Qupstream' contain NA values", {
@@ -96,16 +110,16 @@ test_that("'Qupstream' contain NA values", {
   RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                                   InputsModel = InputsModel,
                                                   IndPeriod_Run = Ind_Run))
-  InputsModel$OutputsModel <- OutputsGR4JOnly
-  # Warning with RunModel
+  QcontribDown <- OutputsGR4JOnly
+  # Warning with RunModel_Lag
   expect_warning(
-    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1, QcontribDown = QcontribDown),
     regexp = "time steps with NA values"
   )
   # No warning during calibration
   RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal
   expect_warning(
-    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1, QcontribDown = QcontribDown),
     regexp = NA
   )
 })
@@ -150,15 +164,16 @@ test_that("1 input with lag of 0.5 time step delay out gives an output delayed o
   expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
 })
 
+InputsCrit <- CreateInputsCrit(
+  FUN_CRIT = ErrorCrit_NSE,
+  InputsModel = InputsModel,
+  RunOptions = RunOptions,
+  VarObs = "Q",
+  Obs = (c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
+           BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas)
+)
+
 test_that("Params from calibration with simulated data should be similar to initial params", {
-  InputsCrit <- CreateInputsCrit(
-    FUN_CRIT = ErrorCrit_NSE,
-    InputsModel = InputsModel,
-    RunOptions = RunOptions,
-    VarObs = "Q",
-    Obs = (c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
-             BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas)
-  )
   CalibOptions <- CreateCalibOptions(
     FUN_MOD = RunModel_GR4J,
     FUN_CALIB = Calibration_Michel,
@@ -175,6 +190,23 @@ test_that("Params from calibration with simulated data should be similar to init
   expect_equal(OutputsCalib$ParamFinalR[1L], ParamSD[1L], tolerance = 2e-3)
 })
 
+test_that("Params from calibration with simulated data should be similar to initial params", {
+  CalibOptions <- CreateCalibOptions(
+    FUN_MOD = RunModel_Lag,
+    FUN_CALIB = Calibration_Michel,
+    IsSD = FALSE
+  )
+  OutputsCalib <- Calibration_Michel(
+    InputsModel = InputsModel,
+    RunOptions = RunOptions,
+    InputsCrit = InputsCrit,
+    CalibOptions = CalibOptions,
+    FUN_MOD = RunModel_Lag,
+    QcontribDown = OutputsGR4JOnly
+  )
+  expect_equal(OutputsCalib$ParamFinalR[1L], ParamSD[1L], tolerance = 2e-3)
+})
+
 test_that("1 no area input with lag of 1 time step delay out gives an output delayed of one time step converted to mm", {
   Qm3GR4Only <- OutputsGR4JOnly$Qsim * BasinAreas[2L] * 1e3
   # Specify that upstream flow is not related to an area