Commit ddf30bac authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Le mode bavard est désormais nommé verbose = TRUE (remplace quiet = FALSE)

parent 78936193
#*************************************************************************************************
#' Calibration algorithm which minimises the error criterion using the provided functions. \cr
#*************************************************************************************************
#' @title Calibration algorithm which minimises an error criterion on the model outputs using the provided functions
#' @author Laurent Coron (June 2014)
#' @seealso \code{\link{Calibration_Michel}}, \code{\link{Calibration_optim}},
#' \code{\link{RunModel}}, \code{\link{ErrorCrit}}, \code{\link{TransfoParam}},
#' \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}},
#' \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
#' @example tests/example_Calibration.R
#' @export
#' @encoding UTF-8
#_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_CALIB (optional) [function] calibration algorithm function (e.g. Calibration_Michel, Calibration_optim), default=Calibration_Michel
#' @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] see \code{\link{Calibration_Michel}} or \code{\link{Calibration_optim}}
#**************************************************************************************************
Calibration <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_CALIB=Calibration_Michel,FUN_TRANSFO=NULL,quiet=FALSE){
return( FUN_CALIB(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO,quiet=quiet) )
Calibration <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_CALIB=Calibration_Michel,FUN_TRANSFO=NULL, verbose = TRUE){
return( FUN_CALIB(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO, verbose = verbose) )
}
#*************************************************************************************************
#' Calibration algorithm which minimises the error criterion. \cr
#' \cr
#' The algorithm is based on a local search procedure.
#' First, a screening is performed using either a rough predefined grid or a list of parameter sets
#' and then a simple steepest descent local search algorithm is performed.
#'
#' A screening is first performed either from a rough predefined grid (considering various initial
#' values for each paramete) or from a list of initial parameter sets. \cr
#' The best set identified in this screening is then used as a starting point for the steepest
#' descent local search algorithm. \cr
#' For this search, the parameters are used in a transformed version, to obtain uniform
#' variation ranges (and thus a similar pace), while the true ranges might be quite different. \cr
#' At each iteration, we start from a parameter set of NParam values (NParam being the number of
#' free parameters of the chosen hydrological model) and we determine the 2*NParam-1 new candidates
#' by changing one by one the different parameters (+/- pace). \cr
#' All these candidates are tested and the best one kept to be the starting point for the next
#' iteration. At the end of each iteration, the pace is either increased or decreased to adapt
#' the progression speed. A diagonal progress can occasionally be done. \cr
#' The calibration algorithm stops when the pace becomes too small. \cr
#'
#' 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 Irstea-HBAN procedure
#' @author Laurent Coron (August 2013)
#' @references
#' Michel, C. (1991),
#' Hydrologie appliquée aux petits bassins ruraux, Hydrology handout (in French), Cemagref, Antony, France.
#' @example tests/example_Calibration_Michel.R
#' @seealso \code{\link{Calibration}}, \code{\link{Calibration_optim}},
#' \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{$NIter } \tab [numeric] number of iterations during the calibration \cr
#' \emph{$NRuns } \tab [numeric] number of model runs done during the calibration \cr
#' \emph{$HistParamR } \tab [numeric] table showing the progression steps in the search for optimal set: parameter values \cr
#' \emph{$HistCrit } \tab [numeric] table showing the progression steps in the search for optimal set: criterion values \cr
#' \emph{$MatBoolCrit } \tab [boolean] table giving the requested and actual time steps when the model is calibrated \cr
#' \emph{$CritName } \tab [character] name of the calibration criterion \cr
#' \emph{$CritBestValue} \tab [numeric] theoretical best criterion value \cr
#' }
#**************************************************************************************************
Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL,quiet=FALSE){
Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL, verbose = TRUE){
##_____Arguments_check_____________________________________________________________________
......@@ -162,13 +103,13 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
##Loop_to_test_the_various_candidates______________________________________
iNewOptim <- 0;
Ncandidates <- nrow(CandidatesParamR);
if(!quiet & Ncandidates>1){
if(verbose & Ncandidates>1){
if(PrefilteringType==1){ cat(paste("\t List-Screening in progress (",sep="")); }
if(PrefilteringType==2){ cat(paste("\t Grid-Screening in progress (",sep="")); }
cat("0%");
}
for(iNew in 1:nrow(CandidatesParamR)){
if(!quiet & Ncandidates>1){
if(verbose & Ncandidates>1){
for(k in c(2,4,6,8)){ if(iNew==round(k/10*Ncandidates)){ cat(paste(" ",10*k,"%",sep="")); } }
}
##Model_run
......@@ -187,7 +128,7 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
Multiplier <- OutputsCrit$Multiplier;
}
}
if(!quiet & Ncandidates>1){ cat(" 100%) \n"); }
if(verbose & Ncandidates>1){ cat(" 100%) \n"); }
##End_of_first_step_Parameter_Screening____________________________________
......@@ -195,7 +136,7 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
ParamStartT <- FUN_TRANSFO(ParamStartR,"RT");
CritStart <- CritOptim;
NRuns <- NRuns+nrow(CandidatesParamR);
if(!quiet){
if(verbose){
if(Ncandidates> 1){ cat(paste("\t Screening completed (",NRuns," runs): \n",sep="")); }
if(Ncandidates==1){ cat(paste("\t Starting point for steepest-descent local search: \n",sep="")); }
cat(paste("\t Param = ",paste(formatC(ParamStartR,format="f",width=8,digits=3),collapse=" , "),"\n",sep=""));
......@@ -249,7 +190,7 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
##Initialisation_of_variables
if(!quiet){
if(verbose){
cat("\t Steepest-descent local search in progress \n");
}
Pace <- 0.64;
......@@ -355,7 +296,7 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
HistParamR[ITER+1,] <- NewParamOptimR;
HistParamT[ITER+1,] <- NewParamOptimT;
HistCrit[ITER+1,] <- CritOptim;
### if(!quiet){ cat(paste("\t Iter ",formatC(ITER,format="d",width=3)," Crit ",formatC(CritOptim,format="f",digits=4)," Pace ",formatC(Pace,format="f",digits=4),"\n",sep="")); }
### if(verbose){ cat(paste("\t Iter ",formatC(ITER,format="d",width=3)," Crit ",formatC(CritOptim,format="f",digits=4)," Pace ",formatC(Pace,format="f",digits=4),"\n",sep="")); }
......@@ -364,7 +305,7 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
##Case_when_the_starting_parameter_set_remains_the_best_solution__________
if(CritOptim==CritStart & !quiet){
if(CritOptim==CritStart & verbose){
cat("\t No progress achieved \n");
}
......@@ -373,7 +314,7 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
ParamFinalT <- NewParamOptimT;
CritFinal <- CritOptim;
NIter <- 1+ITER;
if(!quiet){
if(verbose){
cat(paste("\t Calibration completed (",NIter," iterations, ",NRuns," runs): \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=""));
......
#*************************************************************************************************
#' Creation of the InputsCrit object required to the ErrorCrit functions.
#'
#' Users wanting to use FUN_CRIT functions that are not included in
#' the package must create their own InputsCrit object accordingly.
#*************************************************************************************************
#' @title Creation of the InputsCrit object required to the ErrorCrit functions
#' @author Laurent Coron (June 2014)
#' @seealso \code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateCalibOptions}}
#' @example tests/example_ErrorCrit.R
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________
#' @param FUN_CRIT [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
#' @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 Qobs [numeric] series of observed discharges [mm]
#' @param BoolCrit (optional) [boolean] boolean giving the time steps to consider in the computation (all time steps are consider by default)
#' @param transfo (optional) [character] name of the transformation (e.g. "", "sqrt", "log", "inv", "sort")
#' @param Ind_zeroes (optional) [numeric] indices of the time-steps where zeroes are observed
#' @param epsilon (optional) [numeric] epsilon to add to all Qobs and Qsim if \emph{$Ind_zeroes} is not empty
#_FunctionOutputs_________________________________________________________________________________
#' @return [list] object of class \emph{InputsCrit} containing the data required to evaluate the model outputs; it can include the following:
#' \tabular{ll}{
#' \emph{$BoolCrit } \tab [boolean] boolean giving the time steps to consider in the computation \cr
#' \emph{$Qobs } \tab [numeric] series of observed discharges [mm] \cr
#' \emph{$transfo } \tab [character] name of the transformation (e.g. "", "sqrt", "log", "inv", "sort") \cr
#' \emph{$Ind_zeroes} \tab [numeric] indices of the time-steps where zeroes are observed \cr
#' \emph{$epsilon } \tab [numeric] epsilon to add to all Qobs and Qsim if \emph{$Ind_zeroes} is not empty \cr
#' }
#**************************************************************************************************
CreateInputsCrit <- function(FUN_CRIT,InputsModel,RunOptions,Qobs,BoolCrit=NULL,transfo="",Ind_zeroes=NULL,epsilon=NULL){
ObjectClass <- NULL;
......
#*************************************************************************************************
#' Creation of the InputsModel object required to the RunModel functions.
#'
#' Users wanting to use FUN_MOD functions that are not included in
#' the package must create their own InputsModel object accordingly.
#*************************************************************************************************
#' @title Creation of the InputsModel object required to the RunModel functions
#' @author Laurent Coron (June 2014)
#' @seealso \code{\link{RunModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}, \code{\link{DataAltiExtrapolation_Valery}}
#' @example tests/example_RunModel.R
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________
#' @param FUN_MOD [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
#' @param DatesR [POSIXlt] vector of dates required to create the GR model and CemaNeige module inputs
#' @param Precip [numeric] time series of total precipitation (catchment average) [mm], required to create the GR model and CemaNeige module inputs
#' @param PotEvap [numeric] time series of potential evapotranspiration (catchment average) [mm], required to create the GR model inputs
#' @param TempMean (optional) [numeric] time series of mean air temperature [degC], required to create the CemaNeige module inputs
#' @param TempMin (optional) [numeric] time series of min air temperature [degC], possibly used to create the CemaNeige module inputs
#' @param TempMax (optional) [numeric] time series of max air temperature [degC], possibly used to create the CemaNeige module inputs
#' @param ZInputs (optional) [numeric] real giving the mean elevation of the Precip and Temp series (before extrapolation) [m]
#' @param HypsoData (optional) [numeric] vector of 101 reals: min, q01 to q99 and max of catchment elevation distribution [m], required to create the GR model inputs, if not defined a single elevation is used for CemaNeige
#' @param NLayers (optional) [numeric] integer giving the number of elevation layers requested [-], required to create the GR model inputs, default=5
#' @param quiet (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
#_FunctionOutputs_________________________________________________________________________________
#' @return [list] object of class \emph{InputsModel} containing the data required to evaluate the model outputs; it can include the following:
#' \tabular{ll}{
#' \emph{$DatesR } \tab [POSIXlt] vector of dates \cr
#' \emph{$Precip } \tab [numeric] time series of total precipitation (catchment average) [mm] \cr
#' \emph{$PotEvap } \tab [numeric] time series of potential evapotranspiration (catchment average) [mm], \cr\tab defined if FUN_MOD includes GR4H, GR4J, GR5J, GR6J, GR2M or GR1A \cr \cr
#' \emph{$LayerPrecip } \tab [list] list of time series of precipitation (layer average) [mm], \cr\tab defined if FUN_MOD includes CemaNeige \cr \cr
#' \emph{$LayerTempMean } \tab [list] list of time series of mean air temperature (layer average) [degC], \cr\tab defined if FUN_MOD includes CemaNeige \cr \cr
#' \emph{$LayerFracSolidPrecip} \tab [list] list of time series of solid precip. fract. (layer average) [-], \cr\tab defined if FUN_MOD includes CemaNeige \cr \cr
#' }
#**************************************************************************************************
CreateInputsModel <- function(FUN_MOD,DatesR,Precip,PotEvap=NULL,TempMean=NULL,TempMin=NULL,TempMax=NULL,ZInputs=NULL,HypsoData=NULL,NLayers=5,quiet=FALSE){
CreateInputsModel <- function(FUN_MOD,DatesR,Precip,PotEvap=NULL,TempMean=NULL,TempMin=NULL,TempMax=NULL,ZInputs=NULL,HypsoData=NULL,NLayers=5, verbose = TRUE){
ObjectClass <- NULL;
......@@ -109,11 +74,11 @@ CreateInputsModel <- function(FUN_MOD,DatesR,Precip,PotEvap=NULL,TempMean=NULL,T
if(is.na(ZInputs) | !is.numeric(ZInputs)){ stop("\t ZInputs must be a single numeric value if not null \n"); return(NULL); }
}
if(is.null(HypsoData)){
if(!quiet){ warning("\t HypsoData is missing => a single layer is used and no extrapolation is made \n"); }
if(verbose){ warning("\t HypsoData is missing => a single layer is used and no extrapolation is made \n"); }
HypsoData <- as.numeric(rep(NA,101)); ZInputs <- as.numeric(NA); NLayers <- as.integer(1);
}
if(is.null(ZInputs)){
if(!quiet & !identical(HypsoData,as.numeric(rep(NA,101)))){ warning("\t ZInputs is missing => HypsoData[51] is used \n"); }
if(verbose & !identical(HypsoData,as.numeric(rep(NA,101)))){ warning("\t ZInputs is missing => HypsoData[51] is used \n"); }
ZInputs <- HypsoData[51];
}
}
......@@ -122,15 +87,15 @@ CreateInputsModel <- function(FUN_MOD,DatesR,Precip,PotEvap=NULL,TempMean=NULL,T
##check_NA_values
BOOL_NA <- rep(FALSE,length(DatesR));
if("GR" %in% ObjectClass){
BOOL_NA_TMP <- (Precip < 0) | is.na(Precip ); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < 0 or NA values detected in Precip series \n"); } }
BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < 0 or NA values detected in PotEvap series \n"); } }
BOOL_NA_TMP <- (Precip < 0) | is.na(Precip ); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < 0 or NA values detected in Precip series \n"); } }
BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < 0 or NA values detected in PotEvap series \n"); } }
}
if("CemaNeige" %in% ObjectClass){
BOOL_NA_TMP <- (Precip < 0 ) | is.na(Precip ); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < 0 or NA values detected in Precip series \n"); } }
BOOL_NA_TMP <- (TempMean<(-150)) | is.na(TempMean); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < -150) or NA values detected in TempMean series \n"); } }
BOOL_NA_TMP <- (Precip < 0 ) | is.na(Precip ); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < 0 or NA values detected in Precip series \n"); } }
BOOL_NA_TMP <- (TempMean<(-150)) | is.na(TempMean); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < -150) or NA values detected in TempMean series \n"); } }
if(!is.null(TempMin) & !is.null(TempMax)){
BOOL_NA_TMP <- (TempMin<(-150)) | is.na(TempMin); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < -150) or NA values detected in TempMin series \n"); } }
BOOL_NA_TMP <- (TempMax<(-150)) | is.na(TempMax); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(!quiet){ warning("\t Values < -150) or NA values detected in TempMax series \n"); } } }
BOOL_NA_TMP <- (TempMin<(-150)) | is.na(TempMin); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < -150) or NA values detected in TempMin series \n"); } }
BOOL_NA_TMP <- (TempMax<(-150)) | is.na(TempMax); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < -150) or NA values detected in TempMax series \n"); } } }
}
if(sum(BOOL_NA)!=0){
WTxt <- NULL;
......@@ -144,14 +109,14 @@ CreateInputsModel <- function(FUN_MOD,DatesR,Precip,PotEvap=NULL,TempMean=NULL,T
DatesR <- DatesR[Select];
WTxt <- paste(WTxt,"\t -> data were trunced to keep the most recent available time-steps \n",sep="");
WTxt <- paste(WTxt,"\t -> ",length(Select)," time-steps were kept \n",sep="");
if(!is.null(WTxt) & !quiet){ warning(WTxt); }
if(!is.null(WTxt) & verbose){ warning(WTxt); }
}
##DataAltiExtrapolation_Valery
if("CemaNeige" %in% ObjectClass){
RESULT <- DataAltiExtrapolation_Valery(DatesR=DatesR,Precip=Precip,TempMean=TempMean,TempMin=TempMin,TempMax=TempMax,ZInputs=ZInputs,HypsoData=HypsoData,NLayers=NLayers,quiet=quiet);
if(!quiet){ if(NLayers==1){ cat(paste("\t Input series were successfully created on 1 elevation layer for use by CemaNeige \n",sep=""));
RESULT <- DataAltiExtrapolation_Valery(DatesR=DatesR,Precip=Precip,TempMean=TempMean,TempMin=TempMin,TempMax=TempMax,ZInputs=ZInputs,HypsoData=HypsoData,NLayers=NLayers, verbose = verbose);
if(verbose){ if(NLayers==1){ cat(paste("\t Input series were successfully created on 1 elevation layer for use by CemaNeige \n",sep=""));
} else { cat(paste("\t Input series were successfully created on ",NLayers," elevation layers for use by CemaNeige \n",sep="")); } }
}
......
CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod_Run,IniStates=NULL,IniResLevels=NULL,
Outputs_Cal=NULL,Outputs_Sim="all",RunSnowModule=TRUE,MeanAnSolidPrecip=NULL,quiet=FALSE){
Outputs_Cal=NULL,Outputs_Sim="all",RunSnowModule=TRUE,MeanAnSolidPrecip=NULL, verbose = TRUE){
ObjectClass <- NULL;
......@@ -91,11 +91,11 @@ CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod
if((IndPeriod_Run[1]-1)!=tail(IndPeriod_WarmUp,1) & !identical(IndPeriod_WarmUp,as.integer(0))){
WTxt <- paste(WTxt,"\t Model warm-up period is not directly before the model run period \n",sep=""); }
}
if(!is.null(WTxt) & !quiet){ warning(WTxt); }
if(!is.null(WTxt) & verbose){ warning(WTxt); }
##check_IniStates_and_IniResLevels
if(is.null(IniStates) & is.null(IniResLevels) & !quiet){
if(is.null(IniStates) & is.null(IniResLevels) & verbose){
warning("\t Model states initialisation not defined -> default configuration used \n"); }
if("CemaNeige" %in% ObjectClass){ NLayers <- length(InputsModel$LayerPrecip); } else { NLayers <- 0; }
NState <- NULL;
......@@ -190,7 +190,7 @@ CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod
if(inherits(InputsModel,"yearly" )){ Factor <- 1; }
if(is.null(Factor)){ stop("InputsModel must be of class 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
MeanAnSolidPrecip <- rep(mean(SolidPrecip)*Factor,NLayers); ### default value: same Gseuil for all layers
if(!quiet){ warning("\t MeanAnSolidPrecip not defined -> it was automatically set to c(",paste(round(MeanAnSolidPrecip),collapse=","),") \n"); }
if(verbose){ warning("\t MeanAnSolidPrecip not defined -> it was automatically set to c(",paste(round(MeanAnSolidPrecip),collapse=","),") \n"); }
}
if("CemaNeige" %in% ObjectClass & !is.null(MeanAnSolidPrecip)){
if(!is.vector( MeanAnSolidPrecip) ){ stop(paste("MeanAnSolidPrecip must be a vector of numeric values \n",sep="")); return(NULL); }
......@@ -205,13 +205,13 @@ CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod
WTxt <- NULL;
WTxt <- paste(WTxt,"\t PliqAndMelt was not defined in Outputs_Cal but is needed to feed the hydrological model with the snow module outputs \n",sep="");
WTxt <- paste(WTxt,"\t -> it was automatically added \n",sep="");
if(!is.null(WTxt) & !quiet){ warning(WTxt); }
if(!is.null(WTxt) & verbose){ warning(WTxt); }
Outputs_Cal <- c(Outputs_Cal,"PliqAndMelt"); }
if("PliqAndMelt" %in% Outputs_Sim == FALSE & "all" %in% Outputs_Sim == FALSE){
WTxt <- NULL;
WTxt <- paste(WTxt,"\t PliqAndMelt was not defined in Outputs_Sim but is needed to feed the hydrological model with the snow module outputs \n",sep="");
WTxt <- paste(WTxt,"\t -> it was automatically added \n",sep="");
if(!is.null(WTxt) & !quiet){ warning(WTxt); }
if(!is.null(WTxt) & verbose){ warning(WTxt); }
Outputs_Sim <- c(Outputs_Sim,"PliqAndMelt"); }
}
......
#*****************************************************************************************************************
#' Function which extrapolates the precipitation and air temperature series for different elevation layers (method from Valery, 2010).
#'
#' Elevation layers of equal surface are created the 101 elevation quantiles (\emph{HypsoData})
#' and the number requested elevation layers (\emph{NLayers}). \cr
#' Forcing data (precipitation and air temperature) are extrapolated using gradients from Valery (2010).
#' (e.g. gradP=0.0004 [m-1] for France and gradT=0.434 [degreC/100m] for January, 1st). \cr
#' This function is used by the \emph{CreateInputsModel} function. \cr
#*****************************************************************************************************************
#' @title Altitudinal extrapolation of precipitation and temperature series
#' @author Laurent Coron, Pierre Brigode (June 2014)
#' @references
#' Turcotte, R., L.-G. Fortin, V. Fortin, J.-P. Fortin and J.-P. Villeneuve (2007),
#' Operational analysis of the spatial distribution and the temporal evolution of the snowpack water equivalent
#' in southern Quebec, Canada, Nordic Hydrology, 38(3), 211, doi:10.2166/nh.2007.009. \cr
#' Valéry, A. (2010), Modélisation précipitations-débit sous influence nivale ? : Elaboration d'un module neige
#' et évaluation sur 380 bassins versants, PhD thesis (in french), AgroParisTech, Paris, France. \cr
#' USACE (1956), Snow Hydrology, pp. 437, U.S. Army Coprs of Engineers (USACE) North Pacific Division, Portland, Oregon, USA.
#' @seealso \code{\link{CreateInputsModel}}, \code{\link{RunModel_CemaNeigeGR4J}}
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________________________
#' @param DatesR [POSIXlt] vector of dates
#' @param Precip [numeric] time series of daily total precipitation (catchment average) [mm]
#' @param TempMean [numeric] time series of daily mean air temperature [degC]
#' @param TempMin (optional) [numeric] time series of daily min air temperature [degC]
#' @param TempMax (optional) [numeric] time series of daily max air temperature [degC]
#' @param ZInputs [numeric] real giving the mean elevation of the Precip and Temp series (before extrapolation) [m]
#' @param HypsoData [numeric] vector of 101 reals: min, q01 to q99 and max of catchment elevation distribution [m]
#' @param NLayers [numeric] integer giving the number of elevation layers requested [-]
#' @param quiet (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
#_FunctionOutputs_________________________________________________________________________________________________
#' @return list containing the extrapolated series of precip. and air temp. on each elevation layer
#' \tabular{ll}{
#' \emph{$LayerPrecip } \tab [list] list of time series of daily precipitation (layer average) [mm] \cr
#' \emph{$LayerTempMean } \tab [list] list of time series of daily mean air temperature (layer average) [degC] \cr
#' \emph{$LayerTempMin } \tab [list] list of time series of daily min air temperature (layer average) [degC] \cr
#' \emph{$LayerTempMax } \tab [list] list of time series of daily max air temperature (layer average) [degC] \cr
#' \emph{$LayerFracSolidPrecip} \tab [list] list of time series of daily solid precip. fract. (layer average) [-] \cr
#' \emph{$ZLayers } \tab [numeric] vector of median elevation for each layer \cr
#' }
#*****************************************************************************************************************
DataAltiExtrapolation_Valery <- function(DatesR,Precip,TempMean,TempMin=NULL,TempMax=NULL,ZInputs,HypsoData,NLayers,quiet=FALSE){
DataAltiExtrapolation_Valery <- function(DatesR,Precip,TempMean,TempMin=NULL,TempMax=NULL,ZInputs,HypsoData,NLayers, verbose = TRUE){
##Altitudinal_gradient_functions_______________________________________________________________
......
#*****************************************************************************************************************
#' Function which computes an error criterion with the provided function.
#*****************************************************************************************************************
#' @title Error criterion using the provided function
#' @author Laurent Coron (June 2014)
#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}
#' @example tests/example_ErrorCrit.R
#' @useDynLib airGR
#' @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 FUN_CRIT [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
#' @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, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details
#*****************************************************************************************************************'
ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT,quiet=FALSE){
return( FUN_CRIT(InputsCrit,OutputsModel,quiet=quiet) )
ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT, verbose = TRUE){
return( FUN_CRIT(InputsCrit,OutputsModel, verbose = verbose) )
}
#*****************************************************************************************************************
#' 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){
ErrorCrit_KGE <- function(InputsCrit,OutputsModel, verbose = TRUE){
##Arguments_check________________________________
......@@ -72,7 +39,7 @@ ErrorCrit_KGE <- function(InputsCrit,OutputsModel,quiet=FALSE){
if(inherits(OutputsModel,"daily" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"monthly")){ WarningTS <- 12; }
if(inherits(OutputsModel,"yearly" )){ WarningTS <- 3; }
if(sum(!TS_ignore)<WarningTS & !quiet){ warning(paste("\t criterion computed on less than ",WarningTS," time-steps \n",sep="")); }
if(sum(!TS_ignore)<WarningTS & verbose){ warning(paste("\t criterion computed on less than ",WarningTS," time-steps \n",sep="")); }
##Other_variables_preparation
meanVarObs <- mean(VarObs[!TS_ignore]);
meanVarSim <- mean(VarSim[!TS_ignore]);
......
#*****************************************************************************************************************
#' Function which computes an error criterion based on the KGE' formula proposed by Kling et al. (2012).
#'
#' 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
#' Kling, H., Fuchs, M. and Paulin, M. (2012),
#' Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios,
#' Journal of Hydrology, 424-425, 264-277, doi:10.1016/j.jhydrol.2012.01.011.
#' @seealso \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}
#' @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_KGE2 <- function(InputsCrit,OutputsModel,quiet=FALSE){
ErrorCrit_KGE2 <- function(InputsCrit,OutputsModel, verbose = TRUE){
##Arguments_check________________________________
......@@ -75,7 +39,7 @@ ErrorCrit_KGE2 <- function(InputsCrit,OutputsModel,quiet=FALSE){
if(inherits(OutputsModel,"daily" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"monthly")){ WarningTS <- 12; }
if(inherits(OutputsModel,"yearly" )){ WarningTS <- 3; }
if(sum(!TS_ignore)<WarningTS & !quiet){ warning(paste("\t criterion computed on less than ",WarningTS," time-steps \n",sep="")); }
if(sum(!TS_ignore)<WarningTS & verbose){ warning(paste("\t criterion computed on less than ",WarningTS," time-steps \n",sep="")); }
##Other_variables_preparation
meanVarObs <- mean(VarObs[!TS_ignore]);
meanVarSim <- mean(VarSim[!TS_ignore]);
......
#*****************************************************************************************************************
#' 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){
ErrorCrit_NSE <- function(InputsCrit,OutputsModel, verbose = TRUE){
##Arguments_check________________________________
......@@ -69,7 +38,7 @@ ErrorCrit_NSE <- function(InputsCrit,OutputsModel,quiet=FALSE){
if(inherits(OutputsModel,"daily" )){ WarningTS <- 365; }
if(inherits(OutputsModel,"monthly")){ WarningTS <- 12; }
if(inherits(OutputsModel,"yearly" )){ WarningTS <- 3; }
if(sum(!TS_ignore)<WarningTS & !quiet){ warning(paste("\t criterion computed on less than ",WarningTS," time-steps \n",sep="")); }
if(sum(!TS_ignore)<WarningTS & verbose){ warning(paste("\t criterion computed on less than ",WarningTS," time-steps \n",sep="")); }
##Other_variables_preparation
meanVarObs <- mean(VarObs[!TS_ignore]);
meanVarSim <- mean(VarSim[!TS_ignore]);
......
#*****************************************************************************************************************
#' Function which computes an error criterion based on the root mean square error (RMSE).
#'
#' 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 RMSE
#' @author Laurent Coron (June 2014)
#' @seealso \code{\link{ErrorCrit_NSE}}, \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
#' }