diff --git a/R/CreateInputsCrit_Lavenne.R b/R/CreateInputsCrit_Lavenne.R index f9bed6d1ac15d31ed5eb627821663e9dfc6c2a0a..7aae4bf6b54773b394b78e552b7e6c97507314ad 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 7c0f6eaf9a65be2ce3fe270a0dab82f6cd22bb34..63bd8b3ab3e01344e671dea56ed01ac95cf03066 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 675201e8dd746d59e6e667e1e5f6cccdc5ebc039..eccadcbb7fe3ba4f5f3dac3514188d72da7da224 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 1454498504572eb938fff0da3d199b81bca4af64..18a0ad1603ebec7146eb6e3bd6e5fc89a492a0d5 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 db34200d5bf8ba1c2a177215323162c319320643..d23cc8e061070307219924d82aa8e0d5297c7d23 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