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){


##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); }  


##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


##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; }