From b7230c949a38ffb4cf2729c0d1495ec738e7463e Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@inrae.fr>
Date: Tue, 24 May 2022 20:57:23 +0200
Subject: [PATCH] fix(RunModel_LLR): bug corrections
Version of 24/05/2022
Refs HYCAR-Hydro/airgr#152, HYCAR-Hydro/airgr#153
---
R/CreateInputsCrit_Lavenne.R | 104 +++++++++++++--------------
R/RunModel_LLR.R | 26 ++++---
R/TransfoParam_LLR.R | 2 +-
R/UtilsCalibOptions.R | 12 +++-
inst/modelsFeatures/FeatModelsGR.csv | 2 +-
5 files changed, 79 insertions(+), 67 deletions(-)
diff --git a/R/CreateInputsCrit_Lavenne.R b/R/CreateInputsCrit_Lavenne.R
index f9bed6d1..7aae4bf6 100644
--- a/R/CreateInputsCrit_Lavenne.R
+++ b/R/CreateInputsCrit_Lavenne.R
@@ -1,52 +1,52 @@
-CreateInputsCrit_Lavenne <- function(FUN_CRIT = ErrorCrit_KGE,
- InputsModel,
- RunOptions,
- Obs,
- VarObs = "Q",
- AprParamR,
- AprCrit = 1,
- k = 0.15,
- BoolCrit = NULL,
- transfo = "sqrt",
- epsilon = NULL) {
-
- # Check parameters
- if (!is.numeric(AprCrit) || length(AprCrit) != 1 || AprCrit > 1) {
- stop("'AprCrit' must be a numeric of length 1 with a maximum value of 1")
- }
- if (!is.numeric(k) || length(k) != 1 || k < 0 || k > 1) {
- stop("'k' must be a numeric of length 1 with a value between 0 and 1")
- }
- if (!is.null(BoolCrit) && !is.logical(BoolCrit)) {
- stop("'BoolCrit must be logical")
- }
- if (!is.character(transfo)) {
- stop("'transfo' must be character")
- }
- if (!is.null(epsilon) && !is.numeric(epsilon)) {
- stop("'epsilon' must be numeric")
- }
- if (!is.numeric(AprParamR) || length(AprParamR) != RunOptions$FeatFUN_MOD$NbParam) {
- stop("'AprParamR' must be a numeric vector of length ",
- RunOptions$FeatFUN_MOD$NbParam)
- }
-
-
- FUN_TRANSFO <- .FunTransfo(RunOptions$FeatFUN_MOD)
-
- AprParamT <- FUN_TRANSFO(AprParamR, "RT")
-
- ErrorCrit_GAPX <- CreateErrorCrit_GAPX(FUN_TRANSFO)
-
-
-
- CreateInputsCrit(FUN_CRIT = list(FUN_CRIT, ErrorCrit_GAPX),
- InputsModel = InputsModel,
- RunOptions = RunOptions,
- Obs = list(Obs, AprParamT),
- VarObs = c("Q", "ParamT"),
- Weights = c(1 - k, k * max(0, AprCrit)),
- BoolCrit = list(BoolCrit, NULL),
- transfo = list(transfo, ""),
- epsilon = list(epsilon, NULL))
-}
+CreateInputsCrit_Lavenne <- function(FUN_CRIT = ErrorCrit_KGE,
+ InputsModel,
+ RunOptions,
+ Obs,
+ VarObs = "Q",
+ AprParamR,
+ AprCrit = 1,
+ k = 0.15,
+ BoolCrit = NULL,
+ transfo = "sqrt",
+ epsilon = NULL) {
+
+ # Check parameters
+ if (!is.numeric(AprCrit) || length(AprCrit) != 1 || AprCrit > 1) {
+ stop("'AprCrit' must be a numeric of length 1 with a maximum value of 1")
+ }
+ if (!is.numeric(k) || length(k) != 1 || k < 0 || k > 1) {
+ stop("'k' must be a numeric of length 1 with a value between 0 and 1")
+ }
+ if (!is.null(BoolCrit) && !is.logical(BoolCrit)) {
+ stop("'BoolCrit must be logical")
+ }
+ if (!is.character(transfo)) {
+ stop("'transfo' must be character")
+ }
+ if (!is.null(epsilon) && !is.numeric(epsilon)) {
+ stop("'epsilon' must be numeric")
+ }
+ if (!is.numeric(AprParamR) || length(AprParamR) != RunOptions$FeatFUN_MOD$NbParam) {
+ stop("'AprParamR' must be a numeric vector of length ",
+ RunOptions$FeatFUN_MOD$NbParam)
+ }
+
+
+ FUN_TRANSFO <- .FunTransfo(RunOptions$FeatFUN_MOD, RunOptions$FUN_SD)
+
+ AprParamT <- FUN_TRANSFO(AprParamR, "RT")
+
+ ErrorCrit_GAPX <- CreateErrorCrit_GAPX(FUN_TRANSFO)
+
+
+
+ CreateInputsCrit(FUN_CRIT = list(FUN_CRIT, ErrorCrit_GAPX),
+ InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Obs = list(Obs, AprParamT),
+ VarObs = c("Q", "ParamT"),
+ Weights = c(1 - k, k * max(0, AprCrit)),
+ BoolCrit = list(BoolCrit, NULL),
+ transfo = list(transfo, ""),
+ epsilon = list(epsilon, NULL))
+}
diff --git a/R/RunModel_LLR.R b/R/RunModel_LLR.R
index 7c0f6eaf..63bd8b3a 100644
--- a/R/RunModel_LLR.R
+++ b/R/RunModel_LLR.R
@@ -65,35 +65,39 @@ RunModel_LLR <- function(InputsModel, RunOptions, Param, QcontribDown) {
## Parameters set up
-
TParam <- Param[1L]
KParam <- Param[2L]
## propagation time from upstream meshes to outlet
- PT <- round(InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep)
+ PT <- floor(InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep)
PK <- sqrt(InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * KParam
C0 <- exp(-1/PK)
C1 <- ((PK/1) * (1-C0))-C0
C2 <- 1 - (PK/1) * (1-C0)
-
## Lag model computation
- Qsim_m3 <- matrix(QsimDown *
- InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+ Qsim_m3 <- matrix(0
, nrow = length(IndPeriod1), ncol = NbUpBasins)
-
-
+ QsimDown_input <- matrix(QsimDown *
+ InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3, ncol = 1)
for (upstream_basin in seq_len(NbUpBasins)) {
Qupstream_m3 <- InputsModel$Qupstream[IndPeriod1, upstream_basin]
Qroute <- Qllr <- Qupstream_m3
+ if (is.na(Qroute[1])){
+ Qroute[1] <- 0
+ }
+ Qroute_t1 <- Qroute[1]
for (q_upstream in seq_along(Qupstream_m3)[2:max(seq_along(Qupstream_m3))]){
- Qroute[q_upstream] <- C0 * Qroute[q_upstream - 1] + C1 * Qupstream_m3[q_upstream - 1] + C2 * Qupstream_m3[q_upstream]
+ if (!is.na(Qroute[q_upstream - 1])){
+ Qroute_t1 <- Qroute[q_upstream - 1]
+ }
+ Qroute[q_upstream] <- C0[upstream_basin] * Qroute_t1 + C1[upstream_basin] * Qupstream_m3[q_upstream - 1] + C2[upstream_basin] * Qupstream_m3[q_upstream]
Qllr[q_upstream + PT[upstream_basin]] <- Qroute[q_upstream]
}
- Qsim_m3[, upstream_basin] <- Qllr[IndPeriod1]
+ Qsim_m3[, upstream_basin] <- Qllr[(1 +length(Qllr) - length(IndPeriod1)) :
+ (length(IndPeriod1) + (length(Qllr) - length(IndPeriod1)))]
}
-
-
+ Qsim_m3 <- cbind(QsimDown_input, Qsim_m3)
Qsim_m3 <- rowSums(Qsim_m3)
## OutputsModel
diff --git a/R/TransfoParam_LLR.R b/R/TransfoParam_LLR.R
index 675201e8..eccadcbb 100644
--- a/R/TransfoParam_LLR.R
+++ b/R/TransfoParam_LLR.R
@@ -28,7 +28,7 @@ TransfoParam_LLR <- function(ParamIn, Direction) {
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- ParamIn[, 1] * 20.0 / 20 - 10
- ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
+ ParamOut[, 2] <- ParamIn[, 2] * 20.0 / 20 - 10
}
diff --git a/R/UtilsCalibOptions.R b/R/UtilsCalibOptions.R
index 14544985..18a0ad16 100644
--- a/R/UtilsCalibOptions.R
+++ b/R/UtilsCalibOptions.R
@@ -108,7 +108,11 @@
NParam <- ncol(ParamIn)
ParamOut[, 2:(NParam - 4) ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
- ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
+ if (identical(RunModel_Lag, FUN_SD)){
+ ParamOut[, 1 ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
+ } else if (identical(RunModel_LLR, FUN_SD)){
+ ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
+ }
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -129,7 +133,11 @@
ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)], Direction)
}
ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
- ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
+ if (identical(RunModel_Lag, FUN_SD)){
+ ParamOut[, 1 ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
+ } else if (identical(RunModel_LLR, FUN_SD)){
+ ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
+ }
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
diff --git a/inst/modelsFeatures/FeatModelsGR.csv b/inst/modelsFeatures/FeatModelsGR.csv
index db34200d..d23cc8e0 100644
--- a/inst/modelsFeatures/FeatModelsGR.csv
+++ b/inst/modelsFeatures/FeatModelsGR.csv
@@ -13,4 +13,4 @@ CemaNeigeGR6J;CemaNeigeGR6J;8;daily;NA;GR;airGR
CemaNeigeGR4H;CemaNeigeGR4H;6;hourly;NA;GR;airGR
CemaNeigeGR5H;CemaNeigeGR5H;7;hourly;NA;GR;airGR
Lag;Lag;1;NA;NA;SD;airGR
-LinearLagAndRoute;Lag;2;NA;NA;SD;airGR
+LLR;LLR;2;NA;NA;SD;airGR
--
GitLab