feat: add CreateInputsCrit_DeLavenne

Refs #111

NAMESPACE

@@ -27,6 +27,7 @@ export(CreateCalibOptions)
 export(CreateErrorCrit_GAPX)
 export(CreateIniStates)
 export(CreateInputsCrit)
+export(CreateInputsCrit_DeLavenne)
 export(CreateInputsModel)
 export(CreateRunOptions)
 export(DataAltiExtrapolation_Valery)

R/CreateErrorCrit_GAPX.R

@@ -18,7 +18,7 @@ CreateErrorCrit_GAPX <- function(FUN_TRANSFO) {
       ParamOpt <- EC$VarSim[!EC$TS_ignore]
 
       ## ErrorCrit
-      Crit <- 1 - sum(((ParamApr - ParamOpt) / 20)^2)^0.5
+      Crit <- 1 - sum(((ParamApr - ParamOpt) / 40)^2)^0.5
 
       if (is.numeric(Crit) & is.finite(Crit)) {
         CritValue <- Crit

R/CreateInputsCrit_DeLavenne.R

+CreateInputsCrit_DeLavenne <- function(FUN_CRIT = ErrorCrit_KGE,
+                                      InputsModel,
+                                      RunOptions,
+                                      Obs,
+                                      VarObs = "Q",
+                                      AprParamR,
+                                      AprCrit = 1,
+                                      k = 0.15,
+                                      BoolCrit = NULL,
+                                      transfo = "",
+                                      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")
+  }
+
+
+  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))
+}

man/CreateErrorCrit_GAPX.Rd

@@ -66,7 +66,11 @@ AprParamR <- c(X1 = 157, X2 = 0.8, X3 = 100, X4 = 1.5)
 AprParamT <- TransfoParam_GR4J(AprParamR, "RT")
 
 ## Single efficiency criterion: GAPX with a priori parameters
-InputsCrit <- CreateInputsCrit(ErrorCrit_GAPX, InputsModel, RunOptions, Obs = AprParamT, VarObs = "ParamT")
+InputsCrit <- CreateInputsCrit(ErrorCrit_GAPX,
+                               InputsModel,
+                               RunOptions,
+                               Obs = AprParamT,
+                               VarObs = "ParamT")
 ErrorCrit <- ErrorCrit_GAPX(InputsCrit, OutputsModel)
 str(ErrorCrit)
 }

man/CreateInputsCrit_DeLavenne.Rd

+\encoding{UTF-8}
+
+
+\name{CreateInputsCrit_DeLavenne}
+\alias{CreateInputsCrit_DeLavenne}
+
+
+\title{Creation of the InputsCrit object for De Lavenne Criterion}
+
+
+\description{
+Creation of the \code{InputsCrit} object required to the \code{\link{ErrorCrit}} function. This function defines a composed criterion based on the formula proposed by de Lavenne et al. (2019).
+}
+
+
+\usage{
+CreateInputsCrit_DeLavenne(FUN_CRIT = ErrorCrit_KGE,
+                           InputsModel,
+                           RunOptions,
+                           Obs,
+                           VarObs = "Q",
+                           AprParamR,
+                           AprCrit = 1,
+                           k = 0.15,
+                           BoolCrit = NULL,
+                           transfo = "",
+                           epsilon = NULL)
+}
+
+
+\arguments{
+\item{FUN_CRIT}{[function] error criterion function (e.g. \code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_NSE}})}. Default \code{\link{ErrorCrit_KGE}}
+
+\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
+
+\item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
+
+\item{Obs}{[numeric (atomic or list)] series of observed variable ([mm/time step] for discharge or SWE, [-] for SCA)}
+
+\item{VarObs}{(optional) [character (atomic or list)] names of the observed variable (\code{"Q"} by default, or one of \code{"SCA"}, \code{"SWE"}])}
+
+\item{AprParamR}{[numeric] a priori parameter set from a donor catchment}
+
+\item{AprCrit}{(optional) [numeric] performance criterion of the donor catchment (1 by default)}
+
+\item{k}{(optional) [numeric] coefficient used for the weighted average between \code{FUN_CRIT} and the gap between optimised parameter set and a priori parameter set calculated with the funciton produced by \code{\link{CreateErrorCrit_GAPX}}}
+
+\item{BoolCrit}{(optional) [boolean] boolean (the same length as \code{Obs}) giving the time steps to consider in the computation (all time steps are considered by default. See details)}
+
+\item{transfo}{(optional) [character] name of the transformation applied to the variables (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"}, \code{"boxcox"} or a numeric value for power transformation for \code{FUN_CRIT}. See details of \code{\link{CreateInputsCrit}}}
+
+\item{epsilon}{(optional) [numeric] small value to add to all observations and simulations for \code{FUN_CRIT} when \code{"log"} or \code{"inv"} transformations are used [same unit as \code{Obs}]. See details of \code{\link{CreateInputsCrit}}}
+}
+
+
+\value{
+
+\code{CreateInputsCrit_DeLavenne} returns an object of class \emph{Compo} that is a list of lists such as the one described in \code{\link{CreateInputsCrit}}.
+
+Items \code{Weights} of the criteria are respectively equal to \code{k} and \code{k * max(0,AprCrit)}.
+
+To calculate the De Lavenne criterion, it is necessary to use the \code{ErrorCrit} function as for any other composed criterion.
+}
+
+
+\details{
+
+The parameters \code{FUN_CRIT}, \code{Obs}, \code{VarObs}, \code{BoolCrit}, \code{transfo}, and \code{epsilon} must be use as they would be used for \code{\link{CreateInputsCrit}} in the case of a single criterion.
+
+\code{CreateInputsCrit_DeLavenne} creates a composed criterion in respect with Equations 1 and 2 of de Lavenne et al. (2019).
+
+## --- Period of calculation
+
+Criteria can be calculated over discontinuous periods (i.e. only over winter periods, or when observed discharge is below a certain threshold). To do so, please indicate in \code{Bool_Crit} which indices must be used for calcullation. Discontinuous periods are allowed in the \code{Bool_Crit} argument.
+
+## --- Transformations
+
+Transformations are simple functions applied to the observed and simulated variables used in order to change their distribution. Transformations are often used in hydrology for modifying the weight put on errors made for high flows or low flows. The following transformations are available: \cr \cr
+\itemize{
+  \item \code{""}: no transformation is used (default case)
+  \item \code{"sqrt"}: squared root transformation
+  \item \code{"log"}: logarithmic transformation (see below regarding the specific case of KGE or KGE2)
+  \item \code{"inv"}: inverse transformation
+  \item \code{"sort"}: sort transformation (the simulated and observed variables are sorted from lowest to highest)
+  \item \code{"boxcox"}: Box-Cox transformation (see below for details)
+  \item numeric: power transformation (see below for details)
+}
+We do not advise computing KGE or KGE' with log-transformation as it might be wrongly influenced by discharge values close to 0 or 1 and the criterion value is dependent on the discharge unit. See Santos et al. (2018) for more details and alternative solutions (see the references list below). \cr \cr
+In order to make sure that KGE and KGE2 remain dimensionless and are not impacted by zero values, the Box-Cox transformation (\code{transfo = "boxcox"}) uses the formulation given in Equation 10 of Santos et al. (2018). Lambda is set to 0.25 accordingly. \cr \cr
+The syntax of the power transformation allows a numeric or a string of characters. For example for a squared transformation, the following can be used: \code{transfo = 2}, \code{transfo = "2"} or \code{transfo = "^2"}. Negative values are allowed. Fraction values are not allowed (e.g., \code{"-1/2"} must instead be written \code{"-0.5"}).\cr \cr
+
+## --- The epsilon value
+
+The epsilon value is useful when \code{"log"} or \code{"inv"} transformations are used (to avoid calculation of the inverse or of the logarithm of zero). If an epsilon value is provided, then it is added to the observed and simulated variable time series at each time step and before the application of a transformation. The epsilon value has no effect when the \code{"boxcox"} transformation is used. The impact of this value and a recommendation about the epsilon value to use (usually one hundredth of average observation) are discussed in Pushpalatha et al. (2012) for NSE and in Santos et al. (2018) for KGE and KGE'. \cr \cr
+
+## --- Single, multiple or composite criteria calculation
+
+Users can set the following arguments as atomic or list: \code{FUN_CRIT}, \code{Obs}, \code{VarObs}, \code{BoolCrit}, \code{transfo}, \code{Weights}. If the list format is chosen, all the lists must have the same length. \cr
+Calculation of a single criterion (e.g. NSE computed on discharge) is prepared by providing to \code{CreateInputsCrit} arguments atomics only. \cr
+Calculation of multiple criteria (e.g. NSE computed on discharge and RMSE computed on discharge) is prepared by providing to \code{CreateInputsCrit} arguments lists except for \code{Weights} that must be set as \code{NULL}. \cr
+Calculation of a composite criterion (e.g. the average between NSE computed on discharge and NSE computed on log of discharge) is prepared by providing to \code{CreateInputsCrit} arguments lists including \code{Weights}. \cr
+\code{\link{ErrorCrit_RMSE}} cannot be used in a composite criterion since it is not a unitless value.
+}
+
+
+\examples{
+library(airGR)
+
+## loading catchment data
+data(L0123001)
+
+## preparation of the InputsModel object
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+                                 Precip = BasinObs$P, PotEvap = BasinObs$E)
+
+## calibration period selection
+Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-01"),
+               which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
+
+## preparation of RunOptions object
+RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+                               IndPeriod_Run = Ind_Run)
+
+## simulation
+Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
+
+## single efficiency criterion: Nash-Sutcliffe Efficiency
+InputsCritSingle <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE,
+                                     InputsModel = InputsModel, RunOptions = RunOptions,
+                                     Obs = list(BasinObs$Qmm[Ind_Run]),
+                                     VarObs = "Q", transfo = "",
+                                     Weights = NULL)
+str(InputsCritSingle)
+invisible(ErrorCrit(InputsCrit = InputsCritSingle, OutputsModel = OutputsModel))
+
+## 2 efficiency criteria: RMSE and Nash-Sutcliffe Efficiency
+InputsCritMulti <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_RMSE, ErrorCrit_NSE),
+                                    InputsModel = InputsModel, RunOptions = RunOptions,
+                                    Obs = list(BasinObs$Qmm[Ind_Run], BasinObs$Qmm[Ind_Run]),
+                                    VarObs = list("Q", "Q"), transfo = list("", "sqrt"),
+                                    Weights = NULL)
+str(InputsCritMulti)
+invisible(ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel))
+
+## efficiency composite criterion: Nash-Sutcliffe Efficiency mixing
+##                                 both raw and log-transformed flows
+InputsCritCompo <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_NSE, ErrorCrit_NSE),
+                                    InputsModel = InputsModel, RunOptions = RunOptions,
+                                    Obs = list(BasinObs$Qmm[Ind_Run], BasinObs$Qmm[Ind_Run]),
+                                    VarObs = list("Q", "Q"), transfo = list("", "log"),
+                                    Weights = list(0.4, 0.6))
+str(InputsCritCompo)
+invisible(ErrorCrit(InputsCrit = InputsCritCompo, OutputsModel = OutputsModel))
+}
+
+
+\author{
+David Dorchies
+}
+
+
+\references{
+De Lavenne, A., Andréassian, V., Thirel, G., Ramos, M.-H., Perrin, C., 2019. A Regularization Approach to Improve the Sequential Calibration of a Semidistributed Hydrological Model. Water Resources Research 55, 8821–8839. \url{https://doi.org/10.1029/2018WR024266}
+}
+
+
+\seealso{
+\code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateCalibOptions}}, \code{\link{ErrorCrit}}
+}
+

tests/testthat/test-CreateErrorCrit_GAPX.R

@@ -17,10 +17,16 @@ Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
 OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
                               RunOptions = RunOptions, Param = Param)
 
+ErrorCrit_GAPX <- CreateErrorCrit_GAPX(TransfoParam_GR4J)
+ParamT <- TransfoParam_GR4J(Param, "RT")
+
 test_that("ErrorCrit should return 1 for same parameters", {
-  ErrorCrit_GAPX <- CreateErrorCrit_GAPX(TransfoParam_GR4J)
-  ParamT <- TransfoParam_GR4J(Param, "RT")
   IC <- CreateInputsCrit(ErrorCrit_GAPX, InputsModel, RunOptions, Obs = ParamT, VarObs = "ParamT")
   expect_equal(ErrorCrit_GAPX(IC, OutputsModel)$CritValue, 1)
 })
 
+test_that("ErrorCrit should return 1-nbParam^0.5/40 for ParamT shifted by 1", {
+  ParamT <- ParamT + 1
+  IC <- CreateInputsCrit(ErrorCrit_GAPX, InputsModel, RunOptions, Obs = ParamT, VarObs = "ParamT")
+  expect_equal(ErrorCrit_GAPX(IC, OutputsModel)$CritValue, 1 - RunOptions$FeatFUN_MOD$NbParam^0.5 / 40)
+})

tests/testthat/test-CreateInputsCrit_DeLavenne.R

+context("CreateInputsCrit_DeLavenne")
+
+data(L0123001)
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+                                 Precip = BasinObs$P, PotEvap = BasinObs$E)
+Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
+               which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1999-12-31"))
+RunOptions <- suppressWarnings(
+  CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                   InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
+)
+Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
+                              RunOptions = RunOptions, Param = Param)
+
+k <- 0.15
+IC_DL <- CreateInputsCrit_DeLavenne(FUN_CRIT = ErrorCrit_KGE,
+                                InputsModel = InputsModel,
+                                RunOptions = RunOptions,
+                                Obs = BasinObs$Qmm[Ind_Run],
+                                AprParamR = Param,
+                                k = k)
+
+test_that("should return a composed ErrorCrit", {
+  expect_s3_class(IC_DL, "Compo")
+  expect_length(IC_DL, 2)
+  AprParamT <- TransfoParam_GR4J(Param, "RT")
+  expect_equal(IC_DL$IC2$Obs, AprParamT)
+})
+
+test_that("should return KGE*(1-k)+k with parameters matching a priori parameters", {
+  IC_KGE <- CreateInputsCrit(ErrorCrit_KGE,
+                             InputsModel = InputsModel,
+                             RunOptions = RunOptions,
+                             Obs = BasinObs$Qmm[Ind_Run])
+  expect_equal(ErrorCrit_KGE(IC_KGE, OutputsModel)$CritValue * (1 - k) + k,
+               ErrorCrit(IC_DL, OutputsModel)$CritValue)
+})
+