ErrorCrit_RMSE_regimeMean.R

#*****************************************************************************************************************
#' Function which computes an error criterion based on the root mean square error (RMSE) of the flow regime.
#'
#' Root mean square error on river regime (i.e. interannual mean daily discharges). \cr
#' 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  ErrorCrit_RMSE_RegimeMean
#' @author Laurent Coron (June 2014)
#' @seealso \code{\link{ErrorCrit_RMSE}}
#' @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_RMSE_regimeMean <- 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      <- "RMSE_RegimeMean[Q]";
  CritValue     <- NA;
  CritBestValue <- +0;
  Multiplier    <- +1; ### WARNING must be equal to -1 or +1 only


##Data_preparation_______________________________
  TS_ignore <- is.na(InputsCrit$Qobs) | is.na(OutputsModel$Qsim) | !InputsCrit$BoolCrit ;
  VarObs <- InputsCrit$Qobs  ;   VarObs[TS_ignore] <- NA; 
  VarSim <- OutputsModel$Qsim;   VarSim[TS_ignore] <- NA;
  if(sum(!TS_ignore)<365 & !quiet){ warning("\t criterion computed on less than 365 time-steps \n"); }  
  if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
    VarObs <- VarObs + InputsCrit$epsilon;
    VarSim <- VarSim + InputsCrit$epsilon;
  } }
  AggregData <- aggregate(cbind(VarObs,VarSim),by=list(as.numeric(format(OutputsModel$DatesR,format="%m%d%H%M%S"))),FUN=mean,na.rm=T);
  colnames(AggregData) <- c("AggregDates","Qobs","Qsim");
  MyRollMean2 <- function(x,n){ return(filter(c(tail(x,n%/%2),x,x[1:(n%/%2)]),rep(1/n,n),sides=2)[(n%/%2+1):(length(x)+n%/%2)]); }
  TimeStep  <- difftime(tail(OutputsModel$DatesR,1),tail(OutputsModel$DatesR,2),units="secs")[[1]];
  Window  <- 31*(24*60*60/TimeStep);
  VarObs2 <- MyRollMean2(AggregData$Qobs,Window);
  VarSim2 <- MyRollMean2(AggregData$Qsim,Window);


##ErrorCrit______________________________________
    Crit <- sqrt(sum((VarSim2-VarObs2)^2,na.rm=TRUE)/sum(!is.na(VarObs2)));
    if(is.numeric(Crit) & is.finite(Crit)){ CritValue <- Crit; }


##Output_________________________________________
  OutputsCrit <- list(CritValue,CritName,CritBestValue,Multiplier,which(TS_ignore));
  names(OutputsCrit) <- c("CritValue","CritName","CritBestValue","Multiplier","Ind_notcomputed");
  return(OutputsCrit);