Commit a2078fd2 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Merge branch '109-runmodel_lag-add-an-argument-qcontribdown-to-runmodel_lag' into 'dev'

Resolve "RunModel_Lag: Add an argument `QcontribDown` to `RunModel_Lag`"

Closes #109 and #108

See merge request !37
parents 22fbd89b ae70dd4b
Pipeline #22291 passed with stages
in 25 minutes and 25 seconds
......@@ -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) {
......
......@@ -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
}
......
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
}
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")) {
......
......@@ -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)
......
......@@ -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
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment