Merge branch '88-regularisation-taking-into-account-x4-transformation' into 'dev'

Resolve "Regularisation: taking into account X4 transformation"

Closes #93 and #88

See merge request !45

DESCRIPTION

@@ -37,5 +37,5 @@ VignetteBuilder:
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.0
+RoxygenNote: 7.2.1
 Language: en-US

R/Calibration.GRiwrmInputsModel.R

@@ -54,7 +54,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
       IM$FUN_MOD <- "RunModel_Ungauged"
       attr(RunOptions[[id]], "GRiwrmRunOptions") <- l$RunOptions
     } else {
-      if(useUpstreamQsim && any(IM$UpstreamIsRunoff)) {
+      if (useUpstreamQsim && any(IM$UpstreamIsRunoff)) {
         # Update InputsModel$Qupstream with simulated upstream flows
         IM <- UpdateQsimUpstream(IM, RunOptions[[id]], OutputsModel)
       }
@@ -73,20 +73,18 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
       g <- attr(IM, "GRiwrm")
       Ids <- g$id[g$donor == id & !is.na(g$model)]
       # Extract the X4 calibrated for the whole intermediate basin
-      PS <- attr(IM[[id]], "ParamSettings")
-      if(PS$hasX4) {
-        X4 <- OutputsCalib[[id]]$ParamFinalR[PS$iX4] # Global parameter
+      if(IM[[id]]$model$hasX4) {
+        X4 <- OutputsCalib[[id]]$ParamFinalR[IM[[id]]$model$iX4] # Global parameter
         subBasinAreas <- calcSubBasinAreas(IM)
       }
       for (uId in Ids) {
         # Add OutputsCalib for ungauged nodes
         OutputsCalib[[uId]] <- OutputsCalib[[id]]
         # Copy parameters and transform X4 relatively to the sub-basin area
-        PS <- attr(IM[[uId]], "ParamSettings")
         OutputsCalib[[uId]]$ParamFinalR <-
-          OutputsCalib[[uId]]$ParamFinalR[PS$Indexes]
-        if(PS$hasX4) {
-          OutputsCalib[[uId]]$ParamFinalR[PS$iX4] <-
+          OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$indexParamUngauged]
+        if(IM[[id]]$model$hasX4) {
+          OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$iX4] <-
             X4 * (subBasinAreas[uId] / sum(subBasinAreas)) ^ 0.3
         }
       }
@@ -133,6 +131,10 @@ getInputsCrit_Lavenne <- function(id, OutputsModel, InputsCrit) {
     # Add default velocity parameter for a priori upstream catchment
     AprParamR <- c(AprCelerity, AprParamR)
   }
+  if (attr(InputsCrit[[id]], "model")$hasX4) {
+    featMod <- attr(InputsCrit[[id]], "model")
+    AprParamR[featMod$iX4] <- AprParamR[featMod$iX4] * featMod$X4Ratio
+  }
   AprCrit <- ErrorCrit(InputsCrit[[AprioriId]], OutputsModel[[AprioriId]])$CritValue
   return(Lavenne_FUN(AprParamR, AprCrit))
 }
@@ -195,14 +197,7 @@ updateParameters4Ungauged <- function(GaugedId,
         rep(FALSE, length(InputsModel[[id]]$UpstreamIsRunoff))
     }
   }
-  # Add extra info for Param processing
-  nbParam <- RunOptions[[GaugedId]]$FeatFUN_MOD$NbParam
-  for (id in names(InputsModel)) {
-    attr(InputsModel[[id]], "ParamSettings") <-
-      list(Indexes = ifelse(inherits(InputsModel[[id]], "SD"), 1, 2):nbParam,
-           hasX4 = grepl("RunModel_GR[456][HJ]", InputsModel[[id]]$FUN_MOD),
-           iX4 = ifelse(inherits(InputsModel[[id]], "SD"), 5, 4))
-  }
+
   # Add class InputsModel for airGR::Calibration checks
   class(InputsModel) <- c("InputsModel", class(InputsModel))
 
@@ -256,10 +251,9 @@ RunModel_Ungauged <- function(InputsModel, RunOptions, Param) {
   SBVI <- sum(calcSubBasinAreas(InputsModel))
   # Compute Param for each sub-basin
   P <- lapply(InputsModel, function(IM) {
-    PS <- attr(IM, "ParamSettings")
-    p <- Param[PS$Indexes]
-    if(PS$hasX4) {
-      p[PS$iX4] <- Param[PS$iX4] * (IM$BasinAreas[length(IM$BasinAreas)] / SBVI) ^ 0.3
+    p <- Param[IM$model$indexParamUngauged]
+    if(IM$model$hasX4) {
+      p[IM$model$iX4] <- Param[IM$model$iX4] * (IM$BasinAreas[length(IM$BasinAreas)] / SBVI) ^ 0.3
     }
     return(p)
   })

R/CreateInputsCrit.GRiwrmInputsModel.R

@@ -61,6 +61,10 @@ CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel,
         stop("'AprioriIds': the node \"", AprioriIds[id],
              "\" is an ungauged upstream node of the node \"", id,"\"")
       }
+      if (!identical(InputsModel[[id]]$FUN_MOD, InputsModel[[AprioriIds[id]]]$FUN_MOD)) {
+        stop("'AprioriIds': the node \"", AprioriIds[id],
+             "\" must use the same hydrological model as the node \"", id,"\"")
+      }
     })
   }
 
@@ -88,6 +92,11 @@ CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel,
         )
       attr(InputsCrit[[IM$id]], "AprioriId") <- AprioriIds[IM$id]
       attr(InputsCrit[[IM$id]], "AprCelerity") <- AprCelerity
+      attr(InputsCrit[[IM$id]], "model") <- IM$model
+      if (IM$model$hasX4) {
+        attr(InputsCrit[[IM$id]], "model")$X4Ratio <-
+          (tail(IM$BasinAreas, 1) / tail(InputsModel[[AprioriIds[IM$id]]]$BasinAreas, 1))^0.3
+      }
       class(InputsCrit[[IM$id]]) <- c("InputsCritLavenneFunction", class(InputsCrit[[IM$id]]))
     }
   }

R/CreateInputsModel.GRiwrm.R

@@ -194,10 +194,13 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) {
 
   # Add the model function
   InputsModel$FUN_MOD <- FUN_MOD
+  featModel <- .GetFeatModel(InputsModel)
   InputsModel$isUngauged <- griwrm$model[griwrm$id == id] == "Ungauged"
   InputsModel$gaugedId <- griwrm$donor[griwrm$id == id]
   InputsModel$hasUngaugedNodes <- hasUngaugedNodes(id, griwrm)
-
+  InputsModel$model <- list(indexParamUngauged = ifelse(inherits(InputsModel, "SD"), 0, 1) + seq.int(featModel$NbParam),
+                            hasX4 = grepl("RunModel_GR[456][HJ]", FUN_MOD),
+                            iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4))
   return(InputsModel)
 }
 
@@ -282,3 +285,24 @@ hasUngaugedNodes <- function(id, griwrm) {
   }
   return(FALSE)
 }
+
+
+#' function to extract model features partially copied from airGR:::.GetFeatModel
+#' @noRd
+.GetFeatModel <- function(InputsModel) {
+  path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
+  FeatMod <- read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE)
+  NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR",
+                       yes  = paste("RunModel", FeatMod$NameMod, sep = "_"),
+                       no   = FeatMod$NameMod)
+  IdMod <- which(sapply(NameFunMod, FUN = function(x) identical(InputsModel$FUN_MOD, x)))
+  if (length(IdMod) < 1) {
+    stop("'FUN_MOD' must be one of ", paste(NameFunMod, collapse = ", "))
+  }
+  FeatMod <- as.list(FeatMod[IdMod, ])
+  FeatMod$IsSD <- inherits(InputsModel, "SD")
+  if (FeatMod$IsSD) {
+    FeatMod$NbParam <- FeatMod$NbParam + 1
+  }
+  return(FeatMod)
+}

tests/testthat/test-CreateInputsCrit.R

@@ -109,6 +109,18 @@ test_that("Lavenne criterion: wrong sub-catchment order should throw error", {
   )
 })
 
+test_that("Lavenne criterion: current node and a priori node must use the same model", {
+  InputsModel[["54032"]]$FUN_MOD <- RunModel_GR6J
+  expect_error(
+    CreateInputsCrit(InputsModel = InputsModel,
+                     RunOptions = RunOptions,
+                     Obs = Qobs[IndPeriod_Run,],
+                     AprioriIds = AprioriIds,
+                     transfo = "sqrt"),
+    regexp = "must use the same hydrological model"
+  )
+})
+
 test_that("Ungauged node as Apriori node should throw an error", {
   nodes$model[nodes$gauge_id == "54001"] <- "Ungauged"
   griwrm <- CreateGRiwrm(