From f8dcff40054680985932e3a6a119338f422e0965 Mon Sep 17 00:00:00 2001
From: Delaigue Olivier <olivier.delaigue@irstea.priv>
Date: Fri, 22 Mar 2019 08:25:08 +0100
Subject: [PATCH] v1.2.11.0 NEW: add an IsHyst argument in CreatCalibOptions to
 use hysteresis #5252

---
 DESCRIPTION               |  2 +-
 NEWS.rmd                  |  5 ++-
 R/CreateCalibOptions.R    | 74 ++++++++++++++++++++++++++-------------
 man/CreateCalibOptions.Rd | 10 +++++-
 4 files changed, 63 insertions(+), 28 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index e7326238..5988d279 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.2.10.0
+Version: 1.2.11.0
 Date: 2019-03-22
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NEWS.rmd b/NEWS.rmd
index 58f87544..20bfb41e 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -13,7 +13,7 @@ output:
 
 
 
-### 1.2.10.0 Release Notes (2019-03-22) 
+### 1.2.11.0 Release Notes (2019-03-22) 
 
 
 
@@ -50,6 +50,9 @@ output:
 - <code>CreateRunOptions()</code> now presents a <code>warnings</code> argument to replace the verbose action (the <code>verbose</code> argument is kept to print messages).
 
 
+- <code>CreateCalibOptions()</code> now presents a <code>IsHyst</code> argument to give the possibility to use the hysteresis with CemaNeige.
+
+
 
 #### Major user-visible changes
 
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 3f24b564..b5b2025f 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -1,6 +1,7 @@
 CreateCalibOptions <- function(FUN_MOD,
                                FUN_CALIB = Calibration_Michel,
                                FUN_TRANSFO = NULL,
+                               IsHyst = FALSE,
                                FixedParam = NULL,
                                SearchRanges = NULL,
                                StartParamList = NULL,
@@ -13,7 +14,9 @@ CreateCalibOptions <- function(FUN_MOD,
     if(!is.null(FUN_TRANSFO)) {
       FUN_TRANSFO <- match.fun(FUN_TRANSFO)
     }
-    
+    if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+      stop("'IsHyst' must be a 'logical' of length 1")
+    }
     ##check_FUN_MOD
     BOOL <- FALSE
     
@@ -57,6 +60,9 @@ CreateCalibOptions <- function(FUN_MOD,
       ObjectClass <- c(ObjectClass, "CemaNeigeGR6J")
       BOOL <- TRUE
     }
+	if (IsHyst) {
+	  ObjectClass <- c(ObjectClass, "hysteresis")
+	}
     if (!BOOL) {
       stop("incorrect FUN_MOD for use in CreateCalibOptions")
       return(NULL)
@@ -100,7 +106,7 @@ CreateCalibOptions <- function(FUN_MOD,
         FUN1 <- TransfoParam_GR1A
       }
 	  if (identical(FUN_MOD, RunModel_CemaNeige)) {
-        if (inherits(FUN_MOD, "hysteresis")) {
+        if (IsHyst) {
           FUN1 <- TransfoParam_CemaNeigeHyst
         } else {
           FUN1 <- TransfoParam_CemaNeige
@@ -111,30 +117,50 @@ CreateCalibOptions <- function(FUN_MOD,
         return(NULL)
       }
       ##_set_FUN2
-      FUN2 <- TransfoParam_CemaNeige
-      
+      if (IsHyst) {
+        FUN2 <- TransfoParam_CemaNeigeHyst
+      } else {
+        FUN2 <- TransfoParam_CemaNeige
+      }		
       ##_set_FUN_TRANSFO
       if (sum(ObjectClass %in% c("GR4H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige")) > 0) {
         FUN_TRANSFO <- FUN1
       } else {
-        FUN_TRANSFO <- function(ParamIn, Direction) {
-          Bool <- is.matrix(ParamIn)
-          if (Bool == FALSE) {
-            ParamIn <- rbind(ParamIn)
+	     if (IsHyst) {
+           FUN_TRANSFO <- function(ParamIn, Direction) {
+             Bool <- is.matrix(ParamIn)
+             if (Bool == FALSE) {
+               ParamIn <- rbind(ParamIn)
+             }
+             ParamOut <- NA * ParamIn
+             NParam   <- ncol(ParamIn)
+             ParamOut[, 1:(NParam - 4)] <- FUN1(ParamIn[, 1:(NParam - 4)], Direction)
+             ParamOut[, (NParam - 3):NParam] <- FUN2(ParamIn[, (NParam - 3):NParam], Direction)
+             if (Bool == FALSE) {
+               ParamOut <- ParamOut[1, ]
+             }
+             return(ParamOut)
           }
-          ParamOut <- NA * ParamIn
-          NParam   <- ncol(ParamIn)
-          if (NParam <= 3) {
-            ParamOut[, 1:(NParam - 2)] <- FUN1(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
-          } else {
-            ParamOut[, 1:(NParam - 2)] <- FUN1(ParamIn[, 1:(NParam - 2)], Direction)
-          }
-          ParamOut[, (NParam - 1):NParam] <- FUN2(ParamIn[, (NParam - 1):NParam], Direction)
-          if (Bool == FALSE) {
-            ParamOut <- ParamOut[1, ]
-          }
-          return(ParamOut)
-        }
+		} else {
+          FUN_TRANSFO <- function(ParamIn, Direction) {
+            Bool <- is.matrix(ParamIn)
+            if (Bool == FALSE) {
+              ParamIn <- rbind(ParamIn)
+            }
+            ParamOut <- NA * ParamIn
+            NParam   <- ncol(ParamIn)
+            if (NParam <= 3) {
+              ParamOut[, 1:(NParam - 2)] <- FUN1(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
+            } else {
+              ParamOut[, 1:(NParam - 2)] <- FUN1(ParamIn[, 1:(NParam - 2)], Direction)
+            }
+            ParamOut[, (NParam - 1):NParam] <- FUN2(ParamIn[, (NParam - 1):NParam], Direction)
+            if (Bool == FALSE) {
+              ParamOut <- ParamOut[1, ]
+            }
+            return(ParamOut)
+		  }
+		}
       }
     }
     if (is.null(FUN_TRANSFO)) {
@@ -173,10 +199,9 @@ CreateCalibOptions <- function(FUN_MOD,
     if ("CemaNeigeGR6J" %in% ObjectClass) {
       NParam <- 8
     }
-	if (inherits(FUN_MOD, "hysteresis")) {
+	if (IsHyst) {
 	  NParam <- NParam + 2
 	}
-
     
     ##check_FixedParam
     if (is.null(FixedParam)) {
@@ -200,7 +225,6 @@ CreateCalibOptions <- function(FUN_MOD,
     if (is.null(SearchRanges)) {
       ParamT <-  matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
                         ncol = NParam, byrow = TRUE)
-      
       SearchRanges <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
       
     } else {
@@ -297,7 +321,7 @@ CreateCalibOptions <- function(FUN_MOD,
                              +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = 8, byrow = TRUE)
         }
 	  # }
-      if (inherits(FUN_MOD, "hysteresis")) {
+      if (IsHyst) {
         ParamTHyst <- matrix(c(-9.08, -6.99,
                                -8.00, -3.20,
                                -6.40, +9.99), ncol = 2, byrow = TRUE)
diff --git a/man/CreateCalibOptions.Rd b/man/CreateCalibOptions.Rd
index 4d4d5a35..cc64c571 100644
--- a/man/CreateCalibOptions.Rd
+++ b/man/CreateCalibOptions.Rd
@@ -10,7 +10,7 @@
 
 \usage{
 CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
-  FUN_TRANSFO = NULL, FixedParam = NULL,
+  FUN_TRANSFO = NULL, IsHyst = FALSE, FixedParam = NULL,
   SearchRanges = NULL, StartParamList = NULL,
   StartParamDistrib = NULL)
 }
@@ -23,6 +23,8 @@ CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
 
 \item{FUN_TRANSFO}{(optional) [function] model parameters transformation function, if the \code{FUN_MOD} used is native in the package, \code{FUN_TRANSFO} is automatically defined}
 
+\item{IsHyst}{[boolean] boolean indicating if the hysteresis version of CemaNeige is used. See details}
+
 \item{FixedParam}{(optional) [numeric] vector giving the values set for the non-optimised parameter values (NParam columns, 1 line)
 \cr Example:
 \tabular{llllll}{
@@ -75,6 +77,12 @@ Creation of the \emph{CalibOptions} object required by the \code{Calibration*} f
 \details{
 Users wanting to use \code{FUN_MOD}, \code{FUN_CALIB} or \code{FUN_TRANSFO} functions that are not included in 
 the package must create their own \code{CalibOptions} object accordingly.
+
+##### CemaNeige version #####
+
+If \code{IsHyst = FALSE}, the original CemaNeige version from Valéry et al. (2014) is used.  \cr
+If \code{IsHyst = TRUE}, the CemaNeige version from Riboust et al. (2019) is used. Compared to the original version, this version of CemaNeige needs two more parameters and it includes a representation of the hysteretic relationship between the Snow Cover Area (SCA) and the Snow Water Equivalent (SWE) in the catchment. The hysteresis included in airGR is the Modified Linear hysteresis (LH*); it is represented on panel b) of Fig. 3 in Riboust et al. (2019). Riboust et al. (2019) advise to use the LH* version of CemaNeige with parameters calibrated using an objective function combining 75 \% of KGE calculated on discharge simulated from a rainfall-runoff model compared to observed discharge and 5 \% of KGE calculated on SCA on 5 CemaNeige elevation bands compared to satellite (e.g. MODIS) SCA (see Eq. (18), Table 3 and Fig. 6). Riboust et al. (2019)'s tests were realized with GR4J as the chosen rainfall-runoff model. \cr \cr
+
 }
 
 
-- 
GitLab