ErrorCrit_KGE.R

#*****************************************************************************************************************
#' Function which computes an error criterion based on the KGE formula proposed by Gupta et al. (2009).
#'
#' 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 KGE formula
#' @author Laurent Coron (June 2014)
#' @references
#'   Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F. (2009), 
#'       Decomposition of the mean squared error and NSE performance criteria: Implications
#'       for improving hydrological modelling, Journal of Hydrology, 377(1-2), 80-91, doi:10.1016/j.jhydrol.2009.08.003. \cr
#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \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{$SubCritValues  }   \tab   [numeric] values of the sub-criteria \cr
#'          \emph{$SubCritNames   }   \tab   [character] names of the sub-criteria \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_KGE <- 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 <- "KGE[Q]"      ; }
  if(InputsCrit$transfo=="sqrt"){ CritName <- "KGE[sqrt(Q)]"; }
  if(InputsCrit$transfo=="log" ){ CritName <- "KGE[log(Q)]" ; }
  if(InputsCrit$transfo=="inv" ){ CritName <- "KGE[1/Q]"    ; }
  if(InputsCrit$transfo=="sort"){ CritName <- "KGE[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 ;
  if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
  if(sum(!TS_ignore)<365 & !quiet){ warning("\t criterion computed on less than 365 time-steps \n"); }  
  ##Other_variables_preparation