Commit 87e2d48c authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Merge branch '78-states-handling-on-lag-model' into 'dev'

Resolve "States handling on Lag model"

Closes #78

See merge request !23
parents ab7dc6e7 3be3e839
Pipeline #18835 passed with stages
in 35 minutes and 1 second
......@@ -3,17 +3,18 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
SD = NULL,
verbose = TRUE) {
ObjectClass <- NULL
UH1n <- 20L
UH2n <- UH1n * 2L
nameFUN_MOD <- as.character(substitute(FUN_MOD))
FUN_MOD <- match.fun(FUN_MOD)
## check FUN_MOD
BOOL <- FALSE
if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H)) {
......@@ -56,7 +57,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & IsIntStore) {
stop("'IsIntStore' cannot be TRUE if GR5H is not used in 'FUN_MOD'")
}
## check InputsModel
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
......@@ -68,13 +69,13 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
!inherits(InputsModel, "CemaNeige")) {
stop("'InputsModel' must be of class 'CemaNeige'")
}
## check states
if (any(eTGCemaNeigeLayers > 0)) {
stop("Positive values are not allowed for 'eTGCemaNeigeLayers'")
}
}
if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
if (is.null(ExpStore)) {
stop("'RunModel_*GR6J' need an 'ExpStore' value")
......@@ -85,7 +86,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
}
ExpStore <- Inf
}
if (identical(FUN_MOD, RunModel_GR2M)) {
if (!is.null(UH1)) {
if (verbose) {
......@@ -100,20 +101,20 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
UH2 <- rep(Inf, UH2n)
}
}
if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !is.null(UH1)) {
if (verbose) {
warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
}
UH1 <- rep(Inf, UH1n)
}
}
if ((!identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & !is.null(IntStore)) {
if (verbose) {
warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", nameFUN_MOD))
}
IntStore <- Inf
}
if ("CemaNeige" %in% ObjectClass & ! "GR" %in% ObjectClass) {
if (!is.null(ProdStore)) {
if (verbose) {
......@@ -170,7 +171,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
}
GthrCemaNeigeLayers <- Inf
GlocmaxCemaNeigeLayers <- Inf
GlocmaxCemaNeigeLayers <- Inf
}
if(!"CemaNeige" %in% ObjectClass &
(!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers) | !is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
......@@ -180,24 +181,24 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
GCemaNeigeLayers <- Inf
eTGCemaNeigeLayers <- Inf
GthrCemaNeigeLayers <- Inf
GlocmaxCemaNeigeLayers <- Inf
GlocmaxCemaNeigeLayers <- Inf
}
## set states
if("CemaNeige" %in% ObjectClass) {
NLayers <- length(InputsModel$LayerPrecip)
} else {
NLayers <- 1
}
## manage NULL values
if (is.null(ExpStore)) {
ExpStore <- Inf
ExpStore <- Inf
}
if (is.null(IntStore)) {
IntStore <- Inf
IntStore <- Inf
}
if (is.null(UH1)) {
if ("hourly" %in% ObjectClass) {
......@@ -232,16 +233,16 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
}
if (any(is.infinite(GlocmaxCemaNeigeLayers))) {
GlocmaxCemaNeigeLayers <- rep(Inf, NLayers)
}
}
# check negative values
if (any(ProdStore < 0) | any(RoutStore < 0) | any(IntStore < 0) |
any(UH1 < 0) | any(UH2 < 0) |
any(GCemaNeigeLayers < 0)) {
stop("Negative values are not allowed for any of 'ProdStore', 'RoutStore', 'IntStore', 'UH1', 'UH2', 'GCemaNeigeLayers'")
}
## check length
if (!is.numeric(ProdStore) || length(ProdStore) != 1L) {
stop("'ProdStore' must be numeric of length one")
......@@ -281,7 +282,23 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
}
}
# SD model state handling
if(!is.null(SD)) {
if(!inherits(InputsModel, "SD")) {
stop("'SD' argument provided and 'InputsModel' is not of class 'SD'")
}
if(!is.list(SD)) {
stop("'SD' argument must be a list")
}
lapply(SD, function(x) {
if(!is.numeric(x)) stop("Each item of 'SD' list argument must be numeric")
})
if(length(SD) != length(InputsModel$LengthHydro)) {
stop("Number of items of 'SD' list argument must be the same as the number of upstream connections",
sprintf(" (%i required, found %i)", length(InputsModel$LengthHydro), length(SD)))
}
}
## format output
IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore, Int = IntStore),
......@@ -291,7 +308,11 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
IniStatesNA <- unlist(IniStates)
IniStatesNA[is.infinite(IniStatesNA)] <- NA
IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
if(!is.null(SD)) {
IniStatesNA$SD <- SD
}
class(IniStatesNA) <- c("IniStates", ObjectClass)
if(IsHyst) {
class(IniStatesNA) <- c(class(IniStatesNA), "hysteresis")
......@@ -299,8 +320,8 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if(IsIntStore) {
class(IniStatesNA) <- c(class(IniStatesNA), "interception")
}
return(IniStatesNA)
}
RunModel_Lag <- function(InputsModel, RunOptions, Param) {
NParam <- 1
##Arguments_check
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
}
if (!inherits(InputsModel, "SD")) {
stop("'InputsModel' must be of class 'SD'")
}
}
if (!inherits(RunOptions, "RunOptions")) {
stop("'RunOptions' must be of class 'RunOptions'")
}
......@@ -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,19 +50,46 @@ 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 (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)
}
IniSD[iStart:(iStart + PT[x])]
})
} else {
IniStates <- lapply(seq_len(NbUpBasins), function(x) {
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(rep(0, floor(PT[upstream_basin])),
c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])],
Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
HUTRANS[1, upstream_basin] +
c(rep(0, floor(PT[upstream_basin] + 1)),
c(IniStates[[upstream_basin]],
Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
HUTRANS[2, upstream_basin]
}
......@@ -67,5 +100,12 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
}
# Convert back Qsim to mm
OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
Qupstream[(LengthTs - floor(PT[x])):LengthTs]
})
}
return(OutputsModel)
}
\ No newline at end of file
}
......@@ -19,7 +19,7 @@ CreateIniStates(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = FALSE,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
verbose = TRUE)
SD = NULL, verbose = TRUE)
}
......@@ -52,6 +52,8 @@ CreateIniStates(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = FALSE,
\item{GlocmaxCemaNeigeLayers}{(optional) [numeric] local melt threshold for hysteresis [mm], possibly used to create the CemaNeige model initial state in case the Linear Hysteresis version is used}
\item{SD}{(optional) [list] of [numeric] states of delayed upstream flows for semi-distributed models, the nature of the state and the unit depend on the model and the unit of the upstream flow}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
}
......
context("CreateRunOptions")
test_that("Warm start of GR4J should give same result as warmed model", {
data(L0123001)
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
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_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
# 1990
RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run1))
OutputsModel1 <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions1, Param = Param)
# Warm start 1991
RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run2,
IndPeriod_WarmUp = 0L,
IniStates = OutputsModel1$StateEnd)
OutputsModel2 <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions2, Param = Param)
# Compare 1991 Qsim from warm started and from 1990-1991
expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730])
})
context("CreateIniStates on SD model")
data(L0123001)
test_that("Error: SD argument provided on non-SD 'InputsModel'", {
InputsModel <-
CreateInputsModel(
FUN_MOD = RunModel_GR4J,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E
)
expect_error(
IniStates <-
CreateIniStates(
FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
ProdStore = 0,
RoutStore = 0,
ExpStore = NULL,
UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
SD = list(rep(0, 10))
),
regexp = "'SD' argument provided and"
)
})
BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
# Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
Qupstream <-
floor((sin((
seq_along(BasinObs$Qmm) / 365 * 2 * 3.14
)) + 1) * mean(BasinObs$Qmm, na.rm = TRUE))
InputsModel <- CreateInputsModel(
FUN_MOD = RunModel_GR4J,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E,
Qupstream = matrix(Qupstream, ncol = 1),
LengthHydro = 1000,
BasinAreas = BasinAreas
)
test_that("Error: Non-list 'SD' argument", {
expect_error(
IniStates <-
CreateIniStates(
FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
ProdStore = 0,
RoutStore = 0,
ExpStore = NULL,
UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
SD = rep(0, 10)
),
regexp = "'SD' argument must be a list"
)
})
test_that("Error: Non-numeric items in 'SD' list argument", {
lapply(list(list(list(rep(0, 10))), list(toto = NULL)),
function(x) {
expect_error(
IniStates <-
CreateIniStates(
FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
ProdStore = 0,
RoutStore = 0,
ExpStore = NULL,
UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
SD = x
),
regexp = "Each item of 'SD' list argument must be numeric"
)
})
})
test_that("Error: Number of items not equal to number of upstream connections", {
lapply(list(list(), list(rep(0, 10), rep(0, 10))),
function(x) {
expect_error(
IniStates <-
CreateIniStates(
FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
ProdStore = 0,
RoutStore = 0,
ExpStore = NULL,
UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
SD = x
),
regexp = "list argument must be the same as the number of upstream"
)
})
})
......@@ -50,13 +50,13 @@ InputsModel <- CreateInputsModel(
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run)
IndPeriod_Run = Ind_Run))
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,13 +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", {
# Warm start 1991
RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = IM, IndPeriod_Run = Ind_Run2,
IndPeriod_WarmUp = 0L,
IniStates = OutputsModel1$StateEnd)
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])
})
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