Merge branch '134-crash-with-hysteresis-in-cemaneige' into 'dev'

Resolve "Crash with Hysteresis in Cemaneige"

Closes #135 and #134

See merge request !66

R/CreateCalibOptions.R

@@ -23,17 +23,21 @@ CreateCalibOptions <- function(x, ...) {
 #' @export
 CreateCalibOptions.InputsModel <- function(x,
                                            ...) {
-  if (!exists("FUN_MOD") && !is.null(x$FUN_MOD)) {
-    airGR::CreateCalibOptions(
-      FUN_MOD = x$FUN_MOD,
-      IsSD = !is.null(x$Qupstream) & x$FUN_MOD != "RunModel_Lag",
-      ...
-    )
-  } else {
-    airGR::CreateCalibOptions(
-      ...
-    )
+  dots <- list(...)
+  # Add FUN_MOD in parameters if carried by InputsModel
+  if (!"FUN_MOD" %in% names(dots)) {
+    if(!is.null(x$FUN_MOD)) {
+      dots$FUN_MOD <- x$FUN_MOD
+    } else {
+      stop(" The parameter `FUN_MOD` must be defined")
+    }
   }
+  # Automatically define IsSD for intermediate basin GR models
+  dots$IsSD = !is.null(x$Qupstream) & dots$FUN_MOD != "RunModel_Lag"
+  # Add IsHyst in parameters if carried by InputsModel
+  if (!is.null(x$model$IsHyst)) dots$IsHyst <- x$model$IsHyst
+  # Call airGR function
+  do.call(airGR::CreateCalibOptions, dots)
 }
 
 #' @rdname CreateCalibOptions

R/CreateInputsModel.GRiwrm.R

@@ -35,6 +35,8 @@
 #' @param NLayers (optional) named [vector] of [numeric] integer giving the number
 #'        of elevation layers requested [-], required to create CemaNeige module
 #'        inputs, default=5
+#' @param IsHyst [logical] boolean indicating if the hysteresis version of
+#'        CemaNeige is used. See details of [airGR::CreateRunOptions].
 #' @param ... used for compatibility with S3 methods
 #'
 #' @details Meteorological data are needed for the nodes of the network that
@@ -70,7 +72,8 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
                                      PrecipScale = TRUE,
                                      TempMean = NULL, TempMin = NULL,
                                      TempMax = NULL, ZInputs = NULL,
-                                     HypsoData = NULL, NLayers = 5, ...) {
+                                     HypsoData = NULL, NLayers = 5,
+                                     IsHyst = FALSE, ...) {
 
   # Check and format inputs
   varNames <- c("Precip", "PotEvap", "TempMean", "Qobs", "Qmin",
@@ -212,7 +215,8 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
                                  HypsoData = getInputBV(HypsoData, id),
                                  NLayers = getInputBV(NLayers, id, 5),
                                  Qobs = Qobs,
-                                 Qmin = getInputBV(Qmin, id)
+                                 Qmin = getInputBV(Qmin, id),
+                                 IsHyst = IsHyst
                                  )
   }
   attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
@@ -247,7 +251,7 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
 #'
 #' @return \emph{InputsModel} object for one.
 #' @noRd
-CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
+CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin, IsHyst) {
   np <- getNodeProperties(id, griwrm)
 
   if (np$Diversion) {
@@ -326,7 +330,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
 
   # Add the model function
   InputsModel$FUN_MOD <- FUN_MOD
-  featModel <- .GetFeatModel(InputsModel)
+  featModel <- .GetFeatModel(InputsModel, IsHyst)
   InputsModel$isUngauged <- griwrm$model[griwrm$id == id] == "Ungauged"
   InputsModel$gaugedId <- griwrm$donor[griwrm$id == id]
   InputsModel$hasUngaugedNodes <- hasUngaugedNodes(id, griwrm)
@@ -334,8 +338,9 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
     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)
+      hasX4 = grepl("RunModel_(CemaNeige|)GR[456][HJ]", FUN_MOD),
+      iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4),
+      IsHyst = featModel$IsHyst
     )
   InputsModel$hasDiversion <- np$Diversion
   InputsModel$isReservoir <- np$Reservoir
@@ -446,7 +451,7 @@ hasUngaugedNodes <- function(id, griwrm) {
 #' function to extract model features partially copied from airGR:::.GetFeatModel
 #' @importFrom utils tail
 #' @noRd
-.GetFeatModel <- function(InputsModel) {
+.GetFeatModel <- function(InputsModel, IsHyst) {
   path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
   FeatMod <- read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE)
   NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR",
@@ -461,6 +466,13 @@ hasUngaugedNodes <- function(id, griwrm) {
   if (FeatMod$IsSD) {
     FeatMod$NbParam <- FeatMod$NbParam + 1
   }
+  FeatMod$IsHyst <- FALSE
+  if (grepl("CemaNeige", FeatMod$NameMod)) {
+    FeatMod$IsHyst <- IsHyst
+    if (IsHyst) {
+      FeatMod$NbParam <- FeatMod$NbParam + 2
+    }
+  }
   return(FeatMod)
 }
 

R/CreateRunOptions.R

@@ -43,6 +43,8 @@ CreateRunOptions.InputsModel <- function(x, ...) {
       stop(" The parameter `FUN_MOD` must be defined")
     }
   }
+  # Add IsHyst in parameters if carried by InputsModel
+  if (!is.null(x$model$IsHyst)) dots$IsHyst <- x$model$IsHyst
 
   # Temporary fix waiting for resolution of HYCAR-Hydro/airgr#167
   if (identical(match.fun(dots$FUN_MOD), RunModel_Lag)) {

man/CreateInputsModel.GRiwrm.Rd

@@ -18,6 +18,7 @@
   ZInputs = NULL,
   HypsoData = NULL,
   NLayers = 5,
+  IsHyst = FALSE,
   ...
 )
 }
@@ -70,6 +71,9 @@ if not defined a single elevation is used for CemaNeige}
 of elevation layers requested \link{-}, required to create CemaNeige module
 inputs, default=5}
 
+\item{IsHyst}{\link{logical} boolean indicating if the hysteresis version of
+CemaNeige is used. See details of \link[airGR:CreateRunOptions]{airGR::CreateRunOptions}.}
+
 \item{...}{used for compatibility with S3 methods}
 }
 \value{

tests/testthat/helper_1_RunModel.R

@@ -15,7 +15,8 @@ setupRunModel <-
            runRunOptions = TRUE,
            runRunModel = TRUE,
            griwrm = NULL,
-           Qobs2 = NULL) {
+           Qobs2 = NULL,
+           IsHyst = FALSE) {
 
     data(Severn)
 
@@ -43,12 +44,20 @@ setupRunModel <-
     # Convert meteo data to SD (remove upstream areas)
     Precip <- ConvertMeteoSD(griwrm, PrecipTot)
     PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+    if (IsHyst) {
+      TempMean <- PotEvap+5 # Fake temperatures
+    } else {
+      TempMean <- NULL
+    }
 
     # set up inputs
     if (!runInputsModel)
       return(environment())
     InputsModel <-
-      suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs = Qobs2))
+      suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap,
+                                         TempMean = TempMean,
+                                         Qobs = Qobs2,
+                                         IsHyst = IsHyst))
 
     # RunOptions
     if (!runRunOptions)

tests/testthat/helper_Calibration.R

 runCalibration <- function(nodes = NULL,
                            Qobs2 = NULL,
                            InputsCrit = NULL,
+                           CalibOptions = NULL,
                            FUN_CRIT = ErrorCrit_KGE2,
-                           runRunModel = FALSE) {
+                           runRunModel = FALSE,
+                           IsHyst = FALSE) {
   if (is.null(nodes)) {
     griwrm <- NULL
   } else {
     griwrm <- CreateGRiwrm(nodes)
   }
-  e <- setupRunModel(griwrm = griwrm, runRunModel = runRunModel, Qobs2 = Qobs2)
+  e <- setupRunModel(griwrm = griwrm,
+                     runRunModel = runRunModel,
+                     Qobs2 = Qobs2,
+                     IsHyst = IsHyst)
   for(x in ls(e)) assign(x, get(x, e))
   rm(e)
   np <- getAllNodesProperties(griwrm)
@@ -22,7 +27,8 @@ runCalibration <- function(nodes = NULL,
     )
   }
 
-  CalibOptions <- CreateCalibOptions(InputsModel)
+  if (is.null(CalibOptions)) CalibOptions <- CreateCalibOptions(InputsModel)
+
   OutputsCalib <- Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions)
   Param <- sapply(OutputsCalib, "[[", "ParamFinalR")
   return(environment())

tests/testthat/test-RunModel_Ungauged.R

@@ -233,3 +233,59 @@ test_that("Donor node with diversion should work", {
   for(x in ls(e)) assign(x, get(x, e))
   expect_equal(OC_ref$`54032`$CritFinal, OutputsCalib$`54032`$CritFinal, tolerance = 1E-3)
 })
+
+test_that("Cemaneige with hysteresis works",  {
+  nodes <- loadSevernNodes()
+  nodes <- nodes[nodes$id %in% c("54057", "54032", "54001"), ]
+  nodes$model <- "RunModel_CemaNeigeGR4J"
+  nodes$model[nodes$id != 54057] <- "Ungauged"
+  griwrm <- CreateGRiwrm(nodes)
+
+  # # The custom ErrorCrit function !!!
+  ErrorCrit_KGE3 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
+    OutputsCritQ <- suppressMessages(
+      ErrorCrit_KGE2(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
+    )
+    InputsCrit$Obs <- InputsCrit$SCA #a adapter
+    OutputsModel$Qsim <- OutputsModel$CemaNeigeLayers[[1]]$Gratio #a adapter
+    OutputsCritSCA <- suppressMessages(
+      ErrorCrit_KGE2(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
+    )
+    OutputsCritQ$CritValue <-
+      (OutputsCritQ$CritValue + OutputsCritSCA$CritValue) / 2
+    OutputsCritQ$CritName <- "(0.5 * KGE2[Q] + 0.5 * KGE2[SCA])"
+    return(OutputsCritQ)
+  }
+  class(ErrorCrit_KGE3) <- c("FUN_CRIT", class(ErrorCrit_KGE3))
+
+  e <- suppressWarnings(
+    setupRunModel(griwrm = griwrm, runRunModel = FALSE, IsHyst = TRUE)
+  )
+  for(x in ls(e)) assign(x, get(x, e))
+
+  expect_true(all(sapply(InputsModel, function(x) x$model$hasX4)))
+
+  np <- getAllNodesProperties(griwrm)
+  InputsCrit <- CreateInputsCrit(
+    InputsModel,
+    FUN_CRIT = ErrorCrit_KGE3,
+    RunOptions = RunOptions,
+    Obs = Qobs[IndPeriod_Run, np$id[np$RunOff & np$calibration == "Gauged"], drop = FALSE],
+  )
+  InputsCrit$`54057`$SCA <- runif(length(IndPeriod_Run)) # Fake SCA
+  CalibOptions <- CreateCalibOptions(InputsModel)
+  CO <- lapply(CalibOptions, function(x) {
+    x$StartParamList <- matrix(
+      c(0.605,  320.596,   -0.042,   37.991,    2.221,    0.705,    6.764,   85.000,    0.850),
+      nrow = 1)
+    x$StartParamDistrib <- NULL
+    x
+  })
+  class(CO) <- class(CalibOptions)
+  e <- suppressWarnings(
+    runCalibration(nodes, InputsCrit = InputsCrit, CalibOptions = CO, IsHyst = TRUE)
+  )
+  for(x in ls(e)) assign(x, get(x, e))
+  expect_equal(sapply(Param, length),
+               c("54057" = 9, "54032" = 9, "54001" = 8))
+})