Commit 83a4b5eb authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Retrait de l'algorithme de calibration "optim"

parent 7dfe1d05
......@@ -10,7 +10,6 @@ useDynLib(airGR)
#####################################
export(Calibration)
export(Calibration_Michel)
export(Calibration_optim)
export(CreateCalibOptions)
export(CreateInputsCrit)
export(CreateInputsModel)
......
#*************************************************************************************************
#' 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_Michel}},
#' \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_GR4H )){ FUN_TRANSFO <- TransfoParam_GR4H ; }
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_GR2M )){ FUN_TRANSFO <- TransfoParam_GR2M ; }
if(identical(FUN_MOD,RunModel_GR1A )){ FUN_TRANSFO <- TransfoParam_GR1A ; }
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); }
}
##_RunModelAndCrit
RunModelAndCrit <- function(par,InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO){
ParamT <- NA*CalibOptions$FixedParam;
ParamT[CalibOptions$OptimParam] <- par;
Param <- FUN_TRANSFO(ParamIn=ParamT,Direction="TR");
Param[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam];
OutputsModel <- FUN_MOD(InputsModel=InputsModel,RunOptions=RunOptions,Param=Param);
OutputsCrit <- FUN_CRIT(InputsCrit=InputsCrit,OutputsModel=OutputsModel);
return(OutputsCrit$CritValue*OutputsCrit$Multiplier);
}
##_temporary_change_of_Outputs_Sim
RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal; ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
##_screenPrint
if(!quiet){
cat(paste("\t Calibration in progress (function optim from the stats package) \n",sep=""));
}
##_lower_and_upper_limit_values (transformed)
RangesR <- CalibOptions$SearchRanges;
RangesT <- FUN_TRANSFO(RangesR,"RT");
lower <- RangesT[1,CalibOptions$OptimParam];
upper <- RangesT[2,CalibOptions$OptimParam];
##_starting_values (transformed)
ParamStartT <- FUN_TRANSFO(CalibOptions$StartParam,"RT");
par_start <- ParamStartT[CalibOptions$OptimParam];
##_calibration
RESULT <- optim(par=par_start,fn=RunModelAndCrit,gr=NULL,
InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO, ## arguments for the RunModelAndCrit function (other than par)
method="L-BFGS-B",lower=lower,upper=upper,control=list(),hessian=FALSE)
##_outputs_preparation
ParamFinalT <- NA*ParamStartT;
ParamFinalT[CalibOptions$OptimParam] <- RESULT$par;
ParamFinalR <- FUN_TRANSFO(ParamFinalT,"TR");
ParamFinalR[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam];
CritFinal <- RESULT$value;
##_storage_of_crit_info
OutputsModel <- FUN_MOD(InputsModel=InputsModel,RunOptions=RunOptions,Param=ParamFinalR);
OutputsCrit <- FUN_CRIT(InputsCrit=InputsCrit,OutputsModel=OutputsModel);
CritName <- OutputsCrit$CritName;
CritBestValue <- OutputsCrit$CritBestValue;
Multiplier <- OutputsCrit$Multiplier;
##_screenPrint
if(!quiet){
if(RESULT$convergence==0){
cat(paste("\t Calibration completed: \n",sep=""));
cat(paste("\t Param = ",paste(formatC(ParamFinalR,format="f",width=8,digits=3),collapse=" , "),"\n",sep=""));
cat(paste("\t Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritFinal*Multiplier,format="f",digits=4),"\n",sep=""));
} else {
cat(paste("\t Calibration failed: \n",sep=""));
cat(paste("\t ",RESULT$message,sep=""));
}
}
##_function_output
OutputsCalib <- list(as.double(ParamFinalR),CritFinal*Multiplier,as.integer(RESULT$counts[1]),CritName,CritBestValue);
names(OutputsCalib) <- c("ParamFinalR","CritFinal","NRuns","CritName","CritBestValue");
class(OutputsCalib) <- c("OutputsCalib","optim");
return(OutputsCalib);
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/Calibration_optim.R
\encoding{UTF-8}
\name{Calibration_optim}
\alias{Calibration_optim}
\title{Calibration algorithm which optimises the error criterion selected as objective function using the stats::optim function}
\usage{
Calibration_optim(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD,
FUN_CRIT, FUN_TRANSFO = NULL, quiet = FALSE)
}
\arguments{
\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
\item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
\item{InputsCrit}{[object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details}
\item{CalibOptions}{[object of class \emph{CalibOptions}] see \code{\link{CreateCalibOptions}} for details}
\item{FUN_MOD}{[function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)}
\item{FUN_CRIT}{[function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)}
\item{FUN_TRANSFO}{(optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package, FUN_TRANSFO is automatically defined}
\item{quiet}{(optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE}
}
\value{
[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 selected as objective function 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 selected as objective function \cr
\emph{$CritBestValue} \tab [numeric] theoretical best criterion value \cr
}
}
\description{
Calibration algorithm which optimises the error criterion selected as objective function. \cr
\cr
The algorithm is based on the \code{\link{optim}} function from the "stats" R-package
(using method="L-BFGS-B", i.e. a local optimization quasi-Newton method).
}
\details{
To optimise the exploration of the parameter space, transformation functions are used to convert
the model parameters. This is done using the TransfoParam functions.
}
\examples{
## loading catchment data
library(airGR)
data(L0123001)
## preparation of InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## calibration period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format="\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
## preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_optim)
## calibration
OutputsCalib <- Calibration_optim(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE)
## simulation
Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
## results preview
plot_OutputsModel(OutputsModel = OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
cat(paste(" Crit ", OutputsCrit$CritName, " ",
round(OutputsCrit$CritValue, 4), "\\n", sep = ""))
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
cat(paste(" Crit ", OutputsCrit$CritName, " ",
round(OutputsCrit$CritValue, 4), "\\n", sep = ""))
}
\author{
Laurent Coron (August 2013)
}
\seealso{
\code{\link{Calibration}}, \code{\link{Calibration_Michel}},
\code{\link{RunModel_GR4J}}, \code{\link{TransfoParam_GR4J}}, \code{\link{ErrorCrit_RMSE}},
\code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}},
\code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment