Commit 3be3e839 authored by Dorchies David's avatar Dorchies David
Browse files

feat(RunModel_Lag): add error if IniStates has a wrong lenght

- Add error IniStates length in RunModel_Lag
- Add corresponding expectation in testhat

Refs #78
parent 87ae9ee4
Pipeline #18808 passed with stages
in 12 minutes and 37 seconds
RunModel_Lag <- function(InputsModel, RunOptions, Param) {
NParam <- 1
##Arguments_check
......@@ -19,23 +18,30 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
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")
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")
stop(
"'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment"
)
}
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")
stop(
"'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA"
)
}
OutputsModel <- InputsModel$OutputsModel
OutputsModel$QsimDown <- OutputsModel$Qsim
if (inherits(InputsModel, "daily")) {
TimeStep <- 60 * 60 * 24
}
if (inherits(InputsModel, "hourly")) {
TimeStep <- 60 * 60
} else if (inherits(InputsModel, "daily")) {
TimeStep <- 60 * 60 * 24
} else {
stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
}
# propagation time from upstream meshes to outlet
......@@ -44,27 +50,40 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
NbUpBasins <- length(InputsModel$LengthHydro)
LengthTs <- length(OutputsModel$QsimDown)
OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
OutputsModel$Qsim <-
OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
if(length(IniSD) > 0) {
if (length(IniSD) > 0) {
if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
stop(
sprintf(
"SD initial states has a length of %i and a length of %i is required",
length(IniSD),
sum(floor(PT)) + NbUpBasins
)
)
}
IniStates <- lapply(seq_len(NbUpBasins), function(x) {
iStart <- 1
if(x > 1) {
iStart <- iStart + sum(floor(PT[1:x-1]) + 1)
if (x > 1) {
iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
}
IniSD[iStart:(iStart + PT[x])]
})
} else {
IniStates <- lapply(seq_len(NbUpBasins), function(x) {
rep(0, floor(PT[x] + 1))})
rep(0, floor(PT[x] + 1))
})
}
for (upstream_basin in seq_len(NbUpBasins)) {
Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
Qupstream <-
InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
# Upstream flow with area needs to be converted to m3 by time step
Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
Qupstream <-
Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
}
OutputsModel$Qsim <- OutputsModel$Qsim +
c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])],
......@@ -89,4 +108,4 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
}
return(OutputsModel)
}
\ No newline at end of file
}
......@@ -56,7 +56,7 @@ RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
test_that("InputsModel parameter should contain an OutputsModel key", {
expect_error(
RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
regexp = "'InputsModel' should contain an 'OutputsModel' key"
)
})
......@@ -134,7 +134,6 @@ test_that("1 input with lag of 0.5 time step delay out gives an output delayed o
expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
})
test_that("Params from calibration with simulated data should be similar to initial params", {
InputsCrit <- CreateInputsCrit(
FUN_CRIT = ErrorCrit_NSE,
......@@ -168,45 +167,60 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
InputsModel$BasinAreas <- c(NA, BasinAreas[2L])
# Convert upstream flow to m3/day
InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinAreas[1L] * 1e3
OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
expect_false(any(is.na(OutputsSD$Qsim)))
Qm3SdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1e3
Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
# *** IniStates tests ***
IM <- InputsModel
IM$BasinAreas <- rep(BasinInfo$BasinArea, 3)
IM$Qupstream <- cbind(IM$Qupstream, IM$Qupstream)
IM$LengthHydro <- c(1000, 1500)
PSDini <- ParamSD
PSDini[1] <- PSDini[1] / 2 # 2 time step delay
Ind_Run1 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-12-31"))
Ind_Run2 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-12-31"))
# 1990
RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = IM, IndPeriod_Run = Ind_Run1))
OutputsModel1 <- RunModel(InputsModel = IM,
RunOptions = RunOptions1, Param = PSDini, FUN_MOD = RunModel_GR4J)
# 1990-1991
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = IM, IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
OutputsModel <- RunModel(InputsModel = IM,
RunOptions = RunOptions, Param = PSDini, FUN_MOD = RunModel_GR4J)
test_that("Warm start should give same result as warmed model", {
InputsModel$BasinAreas <- rep(BasinInfo$BasinArea, 3)
InputsModel$Qupstream <- cbind(InputsModel$Qupstream, InputsModel$Qupstream)
InputsModel$LengthHydro <- c(1000, 1500)
ParamSD[1] <- ParamSD[1] / 2 # 2 time step delay
Ind_Run1 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-12-31"))
Ind_Run2 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-12-31"))
# 1990-1991
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
OutputsModel <- RunModel(InputsModel = InputsModel,
RunOptions = RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
# 1990
RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run1))
OutputsModel1 <- RunModel(InputsModel = InputsModel,
RunOptions = RunOptions1, Param = ParamSD, FUN_MOD = RunModel_GR4J)
# Warm start 1991
RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run2,
InputsModel = IM, IndPeriod_Run = Ind_Run2,
IndPeriod_WarmUp = 0L,
IniStates = OutputsModel1$StateEnd)
OutputsModel2 <- RunModel(InputsModel = InputsModel,
RunOptions = RunOptions2, Param = ParamSD, FUN_MOD = RunModel_GR4J)
OutputsModel2 <- RunModel(InputsModel = IM,
RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
# Compare 1991 Qsim from warm started and from 1990-1991
names(OutputsModel2$Qsim) <- NULL
expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730])
})
\ No newline at end of file
})
test_that("Error on Wrong length of iniState$SD", {
OutputsModel1$StateEnd$SD[[1]] <- c(1,1)
RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = IM, IndPeriod_Run = Ind_Run2,
IndPeriod_WarmUp = 0L,
IniStates = OutputsModel1$StateEnd)
expect_error(RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
)
})
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