Calibration_optim.R

#*************************************************************************************************
#' Calibration algorithm which minimises the error criterion. \cr
#' \cr
#' The algorithm is based on the "optim" function from the "stats" R-package 
#' (using method="L-BFGS-B", i.e. a local optimization quasi-Newton method).
#'
#' To optimise the exploration of the parameter space, transformation functions are used to convert
#' the model parameters. This is done using the TransfoParam functions.
#*************************************************************************************************
#' @title  Calibration algorithm which minimises the error criterion using the stats::optim function
#' @author Laurent Coron (August 2013)
#' @example tests/example_Calibration_optim.R
#' @seealso \code{\link{Calibration}}, \code{\link{Calibration_HBAN}},
#'          \code{\link{RunModel_GR4J}}, \code{\link{TransfoParam_GR4J}}, \code{\link{ErrorCrit_RMSE}},
#'          \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, 
#'          \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________
#' @param  InputsModel         [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
#' @param  RunOptions          [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
#' @param  InputsCrit          [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
#' @param  CalibOptions        [object of class \emph{CalibOptions}] see \code{\link{CreateCalibOptions}} for details
#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
#' @param  FUN_CRIT            [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
#' @param  FUN_TRANSFO         (optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package FUN_TRANSFO is automatically defined
#' @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{$ParamFinalR  }  \tab   [numeric] parameter set obtained at the end of the calibration \cr
#'          \emph{$CritFinal    }  \tab   [numeric] error criterion obtained at the end of the calibration \cr
#'          \emph{$Nruns        }  \tab   [numeric] number of model runs done during the calibration \cr
#'          \emph{$CritName     }  \tab   [character] name of the calibration criterion \cr
#'          \emph{$CritBestValue}  \tab   [numeric] theoretical best criterion value \cr
#'          }
#**************************************************************************************************
Calibration_optim <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL,quiet=FALSE){
    

    ##_check_class
    if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
    if(inherits(RunOptions,"RunOptions")==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }  
    if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }  
    if(inherits(CalibOptions,"CalibOptions")==FALSE){ stop("CalibOptions must be of class 'CalibOptions' \n"); return(NULL); }  
    if(inherits(CalibOptions,"optim")==FALSE){ stop("CalibOptions must be of class 'optim' if Calibration_optim is used \n"); return(NULL); }  


    ##_check_FUN_TRANSFO
    if(is.null(FUN_TRANSFO)){
      if(identical(FUN_MOD,RunModel_GR4J         )){ FUN_TRANSFO <- TransfoParam_GR4J     ; }
      if(identical(FUN_MOD,RunModel_GR5J         )){ FUN_TRANSFO <- TransfoParam_GR5J     ; }
      if(identical(FUN_MOD,RunModel_GR6J         )){ FUN_TRANSFO <- TransfoParam_GR6J     ; }
      if(identical(FUN_MOD,RunModel_CemaNeige    )){ FUN_TRANSFO <- TransfoParam_CemaNeige; }
      if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
        if(identical(FUN_MOD,RunModel_CemaNeigeGR4J)){ FUN1 <- TransfoParam_GR4J; FUN2 <- TransfoParam_CemaNeige; }
        if(identical(FUN_MOD,RunModel_CemaNeigeGR5J)){ FUN1 <- TransfoParam_GR5J; FUN2 <- TransfoParam_CemaNeige; }
        if(identical(FUN_MOD,RunModel_CemaNeigeGR6J)){ FUN1 <- TransfoParam_GR6J; FUN2 <- TransfoParam_CemaNeige; }
        FUN_TRANSFO <- function(ParamIn,Direction){
          Bool <- is.matrix(ParamIn);
          if(Bool==FALSE){ ParamIn <- rbind(ParamIn); }
          ParamOut <- NA*ParamIn;
          NParam   <- ncol(ParamIn);
          ParamOut[,         1:(NParam-2)] <- FUN1(ParamIn[,         1:(NParam-2)],Direction);
          ParamOut[,(NParam-1):NParam    ] <- FUN2(ParamIn[,(NParam-1):NParam    ],Direction);
          if(Bool==FALSE){ ParamOut <- ParamOut[1,]; }
          return(ParamOut);
        }
      }
      if(is.null(FUN_TRANSFO)){ stop("FUN_TRANSFO was not found (in Calibration function) \n"); return(NULL);  }