R/ErrorCrit.R

-#*****************************************************************************************************************
-#' Function which computes an error criterion with the provided function.
-#*****************************************************************************************************************
-#' @title  Error criterion using the provided function
-#' @author Laurent Coron (June 2014)
-#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}
-#' @example tests/example_ErrorCrit.R
-#' @useDynLib airGR
-#' @encoding UTF-8
-#' @export
-#_FunctionInputs__________________________________________________________________________________________________
-#' @param  InputsCrit      [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
-#' @param  OutputsModel    [object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details
-#' @param  FUN_CRIT        [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
-#' @param  quiet           (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
-#_FunctionOutputs_________________________________________________________________________________________________
-#' @return  [list] list containing the function outputs, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details
-#*****************************************************************************************************************'
-ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT,quiet=FALSE){
-    return( FUN_CRIT(InputsCrit,OutputsModel,quiet=quiet) )
+ErrorCrit <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
+
+
+  ## ---------- Arguments check
+
+  if (!inherits(InputsCrit, "InputsCrit")) {
+    stop("InputsCrit must be of class 'InputsCrit'")
+  }
+  if (!inherits(OutputsModel, "OutputsModel")) {
+    stop("OutputsModel must be of class 'OutputsModel'")
+  }
+
+
+  ## ---------- Criterion computation
+
+  ## ----- Single criterion
+  if (inherits(InputsCrit, "Single")) {
+    FUN_CRIT <- match.fun(InputsCrit$FUN_CRIT)
+    OutputsCrit <- FUN_CRIT(InputsCrit = InputsCrit,
+                            OutputsModel = OutputsModel,
+                            warnings = warnings,
+                            verbose = verbose)
+  }
+
+
+  ## ----- Multiple criteria or Composite criterion
+
+  if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
+    listOutputsCrit <- lapply(InputsCrit, FUN = function(iInputsCrit) {
+      FUN_CRIT <- match.fun(iInputsCrit$FUN_CRIT)
+      FUN_CRIT(InputsCrit = iInputsCrit,
+               OutputsModel = OutputsModel,
+               warnings = warnings,
+               verbose = verbose)
+    })
+
+    listValCrit  <- sapply(listOutputsCrit, function(x) x[["CritValue"]])
+    listNameCrit <- sapply(listOutputsCrit, function(x) x[["CritName"]])
+    listweights  <- unlist(lapply(InputsCrit, function(x) x[["Weights"]]))
+    listweights  <- listweights / sum(listweights)
+
+    if (inherits(InputsCrit, "Compo")) {
+      CritValue <- sum(listValCrit * listweights)
+      OutputsCritCompo <- list(MultiCritValues  = listValCrit,
+                               MultiCritNames   = listNameCrit,
+                               MultiCritWeights = listweights)
+      OutputsCrit <- list(CritValue       = CritValue,
+                          CritName        = "Composite",
+                          CritBestValue   = +1,
+                          Multiplier      = -1,
+                          Ind_notcomputed = NULL,
+                          CritCompo       = OutputsCritCompo,
+                          MultiCrit       = listOutputsCrit)
+      class(OutputsCrit) <- c("Compo", "ErrorCrit")
+      if (verbose) {
+        message("------------------------------------\n")
+        message("Crit. Composite = ", sprintf("%.4f", CritValue))
+        msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
+        msgForm <- unlist(strsplit(msgForm, split = ","))
+        msgFormSep <- rep(c(",", ",", ",\n\t\t    "), times = ceiling(length(msgForm)/3))[1:length(msgForm)]
+        msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
+        msgForm <- gsub("\\,\\\n\\\t\\\t    $|\\,$", "", msgForm)
+        message("\tFormula: sum(", msgForm, ")\n")
+      }
+    } else {
+      OutputsCrit <- listOutputsCrit
+      class(OutputsCrit) <- c("Multi", "ErrorCrit")
+    }
+
+  }
+
+  return(OutputsCrit)
+
 }
 

R/ErrorCrit_NSE.R

-#*****************************************************************************************************************
-#' Function which computes an error criterion based on the NSE formula proposed by Nash & Sutcliffe (1970).
-#'
-#' In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows 
-#' the use of the function for model calibration: the product CritValue*Multiplier is the criterion to be minimised 
-#' (e.g. Multiplier=+1 for RMSE, Multiplier=-1 for NSE).
-#*****************************************************************************************************************
-#' @title  Error criterion based on the NSE formula
-#' @author Laurent Coron (June 2014)
-#' @references
-#'   Nash, J.E. and Sutcliffe, J.V. (1970), 
-#'       River flow forecasting through conceptual models part 1.
-#'       A discussion of principles, Journal of Hydrology, 10(3), 282-290, doi:10.1016/0022-1694(70)90255-6. \cr
-#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_KGE2}}
-#' @examples ## see example of the ErrorCrit function
-#' @encoding UTF-8
-#' @export
-#_FunctionInputs__________________________________________________________________________________________________
-#' @param  InputsCrit      [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
-#' @param  OutputsModel    [object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details
-#' @param  quiet           (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
-#_FunctionOutputs_________________________________________________________________________________________________
-#' @return  [list] list containing the function outputs organised as follows:
-#'          \tabular{ll}{
-#'          \emph{$CritValue      }   \tab   [numeric] value of the criterion \cr
-#'          \emph{$CritName       }   \tab   [character] name of the criterion \cr
-#'          \emph{$CritBestValue  }   \tab   [numeric] theoretical best criterion value \cr
-#'          \emph{$Multiplier     }   \tab   [numeric] integer indicating whether the criterion is indeed an error (+1) or an efficiency (-1) \cr
-#'          \emph{$Ind_notcomputed}   \tab   [numeric] indices of the time-steps where InputsCrit$BoolCrit=FALSE or no data is available \cr
-#'          }
-#*****************************************************************************************************************
-ErrorCrit_NSE <- function(InputsCrit,OutputsModel,quiet=FALSE){
+ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
 
+  ## Arguments check
+  if (!inherits(OutputsModel, "OutputsModel")) {
+    stop("'OutputsModel' must be of class 'OutputsModel'")
+  }
 
-##Arguments_check________________________________
-  if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }  
-  if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }  
+  EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "NSE", OutputsModel = OutputsModel, warnings = warnings)
 
+  CritValue <- NA
 
-##Initialisation_________________________________
-  CritName <- NA;
-  if(InputsCrit$transfo==""    ){ CritName <- "NSE[Q]"      ; }
-  if(InputsCrit$transfo=="sqrt"){ CritName <- "NSE[sqrt(Q)]"; }
-  if(InputsCrit$transfo=="log" ){ CritName <- "NSE[log(Q)]" ; }
-  if(InputsCrit$transfo=="inv" ){ CritName <- "NSE[1/Q]"    ; }
-  if(InputsCrit$transfo=="sort"){ CritName <- "NSE[sort(Q)]"; }
-  CritValue       <- NA;
-  CritBestValue   <- +1;
-  Multiplier      <- -1; ### must be equal to -1 or +1 only
+  if (EC$CritCompute) {
+    ## ErrorCrit
+    Emod <- sum((EC$VarSim[!EC$TS_ignore] - EC$VarObs[!EC$TS_ignore])^2)
+    Eref <- sum((EC$VarObs[!EC$TS_ignore] - mean(EC$VarObs[!EC$TS_ignore]))^2)
 
+    if (Emod == 0 & Eref == 0) {
+      Crit <- 0
+    } else {
+      Crit <- (1 - Emod / Eref)
+    }
+    if (is.numeric(Crit) & is.finite(Crit)) {
+      CritValue <- Crit
+    }
 
-##Data_preparation_______________________________
-  VarObs <- InputsCrit$Qobs  ; VarObs[!InputsCrit$BoolCrit] <- NA; 
-  VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA;  
-  ##Data_transformation
-  if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
-    VarObs <- VarObs + InputsCrit$epsilon;
-    VarSim <- VarSim + InputsCrit$epsilon;
-  } }
-  if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
-  if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim      < -1E100] <- NA; }
-  if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs    ; VarSim <- 1/VarSim    ; VarSim[abs(VarSim) > 1E+100] <- NA; }
-  if(InputsCrit$transfo=="sort"){ VarObs <- sort(VarObs); VarSim <- sort(VarSim); }
-  ##TS_ignore
-  TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ;
-  Ind_TS_ignore <- which(TS_ignore); if(length(Ind_TS_ignore)==0){ Ind_TS_ignore <- NULL; }
-  if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
-  if(inherits(OutputsModel,"hourly" )){ WarningTS <- 365; }
-  if(inherits(OutputsModel,"daily"  )){ WarningTS <- 365; }
-  if(inherits(OutputsModel,"monthly")){ WarningTS <-  12; }
-  if(inherits(OutputsModel,"yearly" )){ WarningTS <-   3; }
-  if(sum(!TS_ignore)<WarningTS & !quiet){ warning(paste("\t criterion computed on less than ",WarningTS," time-steps \n",sep="")); }
-  ##Other_variables_preparation
-  meanVarObs <- mean(VarObs[!TS_ignore]);
-  meanVarSim <- mean(VarSim[!TS_ignore]);
+    ## Verbose
+    if (verbose) {
+      message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
+    }
+  }
 
 
-##ErrorCrit______________________________________
-  Emod <- sum((VarSim[!TS_ignore]-VarObs[!TS_ignore])^2);
-  Eref <- sum((VarObs[!TS_ignore]-mean(VarObs[!TS_ignore]))^2);
-  if(Emod==0 & Eref==0){ Crit <- 0; } else { Crit <- (1-Emod/Eref); }
-  if(is.numeric(Crit) & is.finite(Crit)){ CritValue <- Crit; }
+  ## Output
+  OutputsCrit <- list(CritValue       = CritValue,
+                      CritName        = EC$CritName,
+                      CritBestValue   = EC$CritBestValue,
+                      Multiplier      = EC$Multiplier,
+                      Ind_notcomputed = EC$Ind_TS_ignore)
 
-
-##Output_________________________________________
-  OutputsCrit <- list(CritValue,CritName,CritBestValue,Multiplier,Ind_TS_ignore);
-  names(OutputsCrit) <- c("CritValue","CritName","CritBestValue","Multiplier","Ind_notcomputed");
-  return(OutputsCrit);
+  class(OutputsCrit) <- c("NSE", "ErrorCrit")
+  return(OutputsCrit)
 
 }
 
-
+class(ErrorCrit_NSE) <- c("FUN_CRIT", class(ErrorCrit_NSE))

R/Imax.R

+
+Imax <- function(InputsModel,
+                 IndPeriod_Run,
+                 TestedValues = seq(from = 0.1, to = 3, by = 0.1)) {
+
+
+  ## ---------- check arguments
+
+  ##  InputsModel
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if (!inherits(InputsModel, "hourly")) {
+    stop("'InputsModel' must be of class 'hourly'")
+  }
+
+  ## IndPeriod_Run
+  if (!is.vector(IndPeriod_Run)) {
+    stop("'IndPeriod_Run' must be a vector of numeric values")
+  }
+  if (!inherits(IndPeriod_Run, "integer")) {
+    stop("'IndPeriod_Run' must be of type integer")
+  }
+  if (!identical(as.integer(IndPeriod_Run), IndPeriod_Run[1]:IndPeriod_Run[length(IndPeriod_Run)])) {
+    stop("'IndPeriod_Run' must be a continuous sequence of integers")
+  }
+
+  ## TestedValues
+  if (!(is.numeric(TestedValues))) {
+    stop("'TestedValues' must be 'numeric'")
+  }
+
+
+  ## ---------- hourly inputs aggregation
+
+  ## aggregate data at the daily time step
+  daily_data <- SeriesAggreg(InputsModel[IndPeriod_Run], Format = "%Y%m%d")
+
+
+  ## ---------- calculate interception
+
+  ## calculate total interception of daily GR models on the period
+  cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
+  if (anyNA(cum_daily)) {
+    stop("'IndPeriod_Run' must be set to select 24 hours by day")
+  }
+
+  ## calculate total interception of the GR5H interception store on the period
+  ## and compute difference with daily values
+  differences <- array(NA, c(length(TestedValues)))
+  for (Imax in TestedValues) {
+    C0 <- 0
+    cum_hourly <- 0
+    for (i in IndPeriod_Run) {
+      Ec  <- min(InputsModel$PotEvap[i], InputsModel$Precip[i] + C0)
+      Pth <- max(0, InputsModel$Precip[i] - (Imax-C0) - Ec)
+      C0  <- C0 + InputsModel$Precip[i] - Ec - Pth
+      cum_hourly <- cum_hourly + Ec
+    }
+    differences[which(Imax == TestedValues)] <- abs(cum_hourly - cum_daily)
+  }
+
+  ## return the Imax value that minimises the difference
+  return(TestedValues[which.min(differences)])
+
+}