fix(RunModel_LLR): bug corrections

Version of 24/05/2022

Refs HYCAR-Hydro/airgr#152, HYCAR-Hydro/airgr#153

R/CreateInputsCrit_Lavenne.R

-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))
+}

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
 

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
   }
 
 

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, ]
         }

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