diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 69a3da22b2b312e01e80b37a0103d8fe1663dba8..0dd1eec50da9aca74f4c9353c0b8e6c878d9207d 100644
--- a/R/Calibration_Michel.R
+++ b/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) {
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 78302745a0d96af3512c53068cacf531d964a566..1d13f2ae29cf394512b220986941176b575b2bf5 100644
--- a/R/CreateCalibOptions.R
+++ b/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
}
diff --git a/R/RunModel.R b/R/RunModel.R
index 1a6104d15fcceaa8195dc3944923d7110ace56dc..c21b803cfb83b02de22136b46a0ad95f38d94424 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -1,21 +1,20 @@
-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
+}
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index f1ed6a2cd40bdd147d852b25b62a3a52f53eb01a..b7f4398897b991ff7cb54980f12aad0f930ef98f 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -1,4 +1,4 @@
-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")) {
diff --git a/man/RunModel_Lag.Rd b/man/RunModel_Lag.Rd
index b37487679c2734efdddda7f066a4ad9133e07cff..5819f083a7c52e5ba3a6587db83c63daf64252a9 100644
--- a/man/RunModel_Lag.Rd
+++ b/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)
diff --git a/tests/testthat/test-RunModel_Lag.R b/tests/testthat/test-RunModel_Lag.R
index 0432fa091351418eaa4bfeee93a906dd805183f6..8b83924c280147c3a826e940a66afaffa20de40b 100644
--- a/tests/testthat/test-RunModel_Lag.R
+++ b/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