Commit 92a63f96 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Version 0.8.1.0

parent 56d606a1
Package: airGR
Type: Package
Title: Modelling tools used at Irstea-HBAN (France), including GR4J,
GR5J, GR6J and CemaNeige
Version: 0.7.4
Date: 2014-11-01
Title: Modelling tools used at Irstea-HBAN (France), including GR4J and
CemaNeige
Version: 0.8.1.0
Date: 2015-10-27
Author: Laurent CORON
Maintainer: Laurent CORON <laurent.coron@irstea.fr>, Olivier DELAIGUE
Maintainer: Laurent CORON, Olivier DELAIGUE
<olivier.delaigue@irstea.fr>
Depends: R (>= 3.0.1)
Description: This package brings into R the hydrological modelling tools used
at Irstea-HBAN (France). The package includes several conceptual
rainfall-runoff models and the associated functions for their calibration
and evaluation,including GR4J, GR5J, GR6J and CemaNeige. Use help(airGR)
for package description.
Description: This package brings into R the hydrological modelling tools used at Irstea-HBAN (France).
The package includes several conceptual rainfall-runoff models and the associated functions
for their calibration and evaluation (GR4H, GR4J, GR5J, GR6J, GR2M, GR1A and CemaNeige).
Use help(airGR) for package description.
License: GPL-2
Packaged: 2014-11-25 21:51:31 UTC; H61970
Built: R 3.0.2; x86_64-w64-mingw32; 2015-10-27 18:41:54 UTC; windows
Archs: i386, x64
BasinInfo Data sample: characteristics of a fictional
catchment (L0123001, L0123002 or L0123003)
BasinObs Data sample: time series of observations of a
fictional catchment (L0123001, L0123002 or
L0123003)
Calibration Calibration algorithm which minimises an error
criterion on the model outputs using the
provided functions
Calibration_HBAN Calibration algorithm which minimises the error
criterion using the Irstea-HBAN procedure
Calibration_optim Calibration algorithm which minimises the error
criterion using the stats::optim function
CreateCalibOptions Creation of the CalibOptions object required to
the Calibration functions
CreateInputsCrit Creation of the InputsCrit object required to
the ErrorCrit functions
CreateInputsModel Creation of the InputsModel object required to
the RunModel functions
CreateRunOptions Creation of the RunOptions object required to
the RunModel functions
DataAltiExtrapolation_HBAN
Altitudinal extrapolation of precipitation and
temperature series
ErrorCrit Error criterion using the provided function
ErrorCrit_KGE Error criterion based on the KGE formula
ErrorCrit_KGE2 Error criterion based on the KGE' formula
ErrorCrit_NSE Error criterion based on the NSE formula
ErrorCrit_RMSE Error criterion based on the RMSE
PEdaily_Oudin Computation of daily series of potential
evapotranspiration with Oudin's formula
RunModel Run with the provided hydrological model
function
RunModel_CemaNeige Run with the CemaNeige snow module
RunModel_CemaNeigeGR4J
Run with the CemaNeigeGR4J hydrological model
RunModel_CemaNeigeGR5J
Run with the CemaNeigeGR5J hydrological model
RunModel_CemaNeigeGR6J
Run with the CemaNeigeGR6J hydrological model
RunModel_GR1A Run with the GR1A hydrological model
RunModel_GR2M Run with the GR2M hydrological model
RunModel_GR4H Run with the GR4H hydrological model
RunModel_GR4J Run with the GR4J hydrological model
RunModel_GR5J Run with the GR5J hydrological model
RunModel_GR6J Run with the GR6J hydrological model
SeriesAggreg Conversion of time series to another time-step
(aggregation only)
TransfoParam Transformation of the parameters using the
provided function
TransfoParam_CemaNeige
Transformation of the parameters from the
CemaNeige module
TransfoParam_GR1A Transformation of the parameters from the GR1A
model
TransfoParam_GR2M Transformation of the parameters from the GR2M
model
TransfoParam_GR4H Transformation of the parameters from the GR4H
model
TransfoParam_GR4J Transformation of the parameters from the GR4J
model
TransfoParam_GR5J Transformation of the parameters from the GR5J
model
TransfoParam_GR6J Transformation of the parameters from the GR6J
model
airGR Modelling tools used at Irstea-HBAN (France),
including GR4J and CemaNeige
plot_OutputsModel Default preview of model outputs
745f1ff9a1a987be74e26d327e1001f3 *DESCRIPTION
8cfd5d38abb52c25f9de6da3b5a49dbd *INDEX
4ef88d7fe2a815ae9d537a6a0ce475e2 *Meta/Rd.rds
32a1c5de93e3b6254dbd86b07ba073ba *Meta/data.rds
85e37b304c759576eef3e1715799d846 *Meta/hsearch.rds
13463ac75ee802be8ab7942c547edf89 *Meta/links.rds
db1f9a1a649ce0c9ed5f0bf1321337a4 *Meta/nsInfo.rds
886528f90ab482d4b30d2586da61f7e7 *Meta/package.rds
ef20fcb07a98bb3a13d1dfbce9b66ab0 *NAMESPACE
ebf0fc819595d631b8bf280c4b049940 *R/airGR
44c5327b5161fa3578c84656f9e872c6 *R/airGR.rdb
af79b71953440b5fb2bb5968d1c7c4a5 *R/airGR.rdx
9a4212f4316f102accff6f1b737b591f *data/L0123001.rda
1940776e833cda1019508134b87c65f0 *data/L0123002.rda
517b0e945c588adcb1846db0b138cd10 *data/L0123003.rda
98c6e43c002546b43228d42f834515f2 *help/AnIndex
7d2aecc63e95084cf48a00bb02d63f5c *help/airGR.rdb
42aefbf56becfa3987a4b3e94f8c032d *help/airGR.rdx
ddae52805808b1d0202eea1457a18f12 *help/aliases.rds
06e18f0d21774860b6fe45293af4f4fb *help/paths.rds
e5b5791fbc25b39ec22b8509c01f4d69 *html/00Index.html
444535b9cb76ddff1bab1e1865a3fb14 *html/R.css
c27696c4a2f0ae0ccc8b7cd8ad62b66e *libs/i386/airGR.dll
991f7b24c923dfc3bd5b7420ae1a6cc7 *libs/x64/airGR.dll
File added
File added
File added
File added
# Generated by roxygen2 (4.0.1): do not edit by hand
# Generated by roxygen2 (4.1.1): do not edit by hand
export(Calibration)
export(Calibration_HBAN)
......@@ -19,11 +19,18 @@ export(RunModel_CemaNeige)
export(RunModel_CemaNeigeGR4J)
export(RunModel_CemaNeigeGR5J)
export(RunModel_CemaNeigeGR6J)
export(RunModel_GR1A)
export(RunModel_GR2M)
export(RunModel_GR4H)
export(RunModel_GR4J)
export(RunModel_GR5J)
export(RunModel_GR6J)
export(SeriesAggreg)
export(TransfoParam)
export(TransfoParam_CemaNeige)
export(TransfoParam_GR1A)
export(TransfoParam_GR2M)
export(TransfoParam_GR4H)
export(TransfoParam_GR4J)
export(TransfoParam_GR5J)
export(TransfoParam_GR6J)
......
#' @name BasinInfo
#' @docType data
#' @title Data sample: characteristics of a fictional catchment (L0123001, L0123002 or L0123003)
#' @description
#' R-object containing the code, station's name, area and hypsometric curve of the catchment.
#' @encoding UTF-8
#' @format
#' List named 'BasinInfo' containing
#' \itemize{
#' \item two strings: catchment's code and station's name
#' \item one float: catchment's area in km2
#' \item one numeric vector: catchment's hypsometric curve (min, quantiles 01 to 99 and max) in metres
#' }
#' @examples
#' require(airGR)
#' data(L0123001)
#' str(BasinInfo)
NULL
#' @name BasinObs
#' @docType data
#' @title Data sample: time series of observations of a fictional catchment (L0123001, L0123002 or L0123003)
#' @description
#' R-object containing the times series of precipitation, temperature, potential evapotranspiration and discharges. \cr
#' Times series for L0123001 or L0123002 are at the daily time-step for use with daily models such as GR4J, GR5J, GR6J, CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J.
#' Times series for L0123003 are at the hourly time-step for use with hourly models such as GR4H.
#' @encoding UTF-8
#' @format
#' Data frame named 'BasinObs' containing
#' \itemize{
#' \item one POSIXlt vector: time series dates in the POSIXlt format
#' \item five numeric vectors: time series of catchment average precipitation [mm], catchment average air temperature [degC], catchment average potential evapotranspiration [mm], outlet discharge [l/s], outlet discharge [mm]
#' }
#' @examples
#' require(airGR)
#' data(L0123001)
#' str(BasinObs)
NULL
#*************************************************************************************************
#' 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_HBAN}}, \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_HBAN, Calibration_optim), default=Calibration_HBAN
#' @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_HBAN}} or \code{\link{Calibration_optim}}
#**************************************************************************************************
Calibration <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_CALIB=Calibration_HBAN,FUN_TRANSFO=NULL,quiet=FALSE){
return( FUN_CALIB(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO,quiet=quiet) )
}
#*************************************************************************************************
#' 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_HBAN.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_HBAN <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL,quiet=FALSE){
##_____Arguments_check_____________________________________________________________________
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,"HBAN")==FALSE){ stop("CalibOptions must be of class 'HBAN' if Calibration_HBAN 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); }
}
##_variables_initialisation
ParamFinalR <- NULL; ParamFinalT <- NULL; CritFinal <- NULL;
NRuns <- 0; NIter <- 0;
if("StartParamDistrib" %in% names(CalibOptions)){ PrefilteringType <- 2; } else { PrefilteringType <- 1; }
if(PrefilteringType==1){ NParam <- ncol(CalibOptions$StartParamList); }
if(PrefilteringType==2){ NParam <- ncol(CalibOptions$StartParamDistrib); }
if(NParam>20){ stop("Calibration_HBAN can handle a maximum of 20 parameters \n"); return(NULL); }
HistParamR <- matrix(NA,nrow=500*NParam,ncol=NParam);
HistParamT <- matrix(NA,nrow=500*NParam,ncol=NParam);
HistCrit <- matrix(NA,nrow=500*NParam,ncol=1);
CritName <- NULL;
CritBestValue <- NULL;
Multiplier <- NULL;
CritOptim <- +1E100;
##_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
##_____Parameter_Grid_Screening____________________________________________________________
##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
ProposeCandidatesGrid <- function(DistribParam){
##Managing_matrix_sizes
Nvalmax <- nrow(DistribParam);
NParam <- ncol(DistribParam);
##we_add_columns_to_MatDistrib_until_it_has_20_columns
DistribParam2 <- matrix(NA,nrow=Nvalmax,ncol=20);
DistribParam2[1:Nvalmax,1:NParam] <- DistribParam;
##we_check_the_number_of_values_to_test_for_each_param
NbDistrib <- rep(1,20);
for(iC in 1:20){ NbDistrib[iC] <- max( 1 , Nvalmax-sum(is.na(DistribParam2[,iC])) ); }
##Loop_on_the_various_values_to_test ###(if 4 param and 3 values for each => 3^4 sets)
##NB_we_always_do_20_loops ###which_is_here_the_max_number_of_param_that_can_be_optimised
VECT <- NULL;
for(iL01 in 1:NbDistrib[01]){ for(iL02 in 1:NbDistrib[02]){ for(iL03 in 1:NbDistrib[03]){ for(iL04 in 1:NbDistrib[04]){ for(iL05 in 1:NbDistrib[05]){
for(iL06 in 1:NbDistrib[06]){ for(iL07 in 1:NbDistrib[07]){ for(iL08 in 1:NbDistrib[08]){ for(iL09 in 1:NbDistrib[09]){ for(iL10 in 1:NbDistrib[10]){
for(iL11 in 1:NbDistrib[11]){ for(iL12 in 1:NbDistrib[12]){ for(iL13 in 1:NbDistrib[13]){ for(iL14 in 1:NbDistrib[14]){ for(iL15 in 1:NbDistrib[15]){
for(iL16 in 1:NbDistrib[16]){ for(iL17 in 1:NbDistrib[17]){ for(iL18 in 1:NbDistrib[18]){ for(iL19 in 1:NbDistrib[19]){ for(iL20 in 1:NbDistrib[20]){
VECT <- c(VECT,
DistribParam2[iL01,01],DistribParam2[iL02,02],DistribParam2[iL03,03],DistribParam2[iL04,04],DistribParam2[iL05,05],
DistribParam2[iL06,06],DistribParam2[iL07,07],DistribParam2[iL08,08],DistribParam2[iL09,09],DistribParam2[iL10,10],
DistribParam2[iL11,11],DistribParam2[iL12,12],DistribParam2[iL13,13],DistribParam2[iL14,14],DistribParam2[iL15,15],
DistribParam2[iL16,16],DistribParam2[iL17,17],DistribParam2[iL18,18],DistribParam2[iL19,19],DistribParam2[iL20,20]);
} } } } }
} } } } }
} } } } }
} } } } }
MAT <- matrix(VECT,ncol=20,byrow=TRUE)[,1:NParam];
if(is.matrix(MAT)==FALSE){ MAT <- cbind(MAT); }
Output <- NULL;
Output$NewCandidates <- MAT;
return(Output);
}
##Creation_of_new_candidates_______________________________________________
if(PrefilteringType==1){ CandidatesParamR <- CalibOptions$StartParamList; }
if(PrefilteringType==2){ DistribParamR <- CalibOptions$StartParamDistrib; DistribParamR[,!CalibOptions$OptimParam] <- NA; CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates; }
##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam]; return(x); });
if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); }
##Loop_to_test_the_various_candidates______________________________________
iNewOptim <- 0;
Ncandidates <- nrow(CandidatesParamR);
if(!quiet & 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){
for(k in c(2,4,6,8)){ if(iNew==round(k/10*Ncandidates)){ cat(paste(" ",10*k,"%",sep="")); } }
}
##Model_run
Param <- CandidatesParamR[iNew,];
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);
if(!is.na(OutputsCrit$CritValue)){ if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier;
iNewOptim <- iNew;
} }
##Storage_of_crit_info
if(is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)){
CritName <- OutputsCrit$CritName;
CritBestValue <- OutputsCrit$CritBestValue;
Multiplier <- OutputsCrit$Multiplier;
}
}
if(!quiet & Ncandidates>1){ cat(" 100%) \n"); }
##End_of_first_step_Parameter_Screening____________________________________
ParamStartR <- CandidatesParamR[iNewOptim,]; if(!is.matrix(ParamStartR)){ ParamStartR <- matrix(ParamStartR,nrow=1); }
ParamStartT <- FUN_TRANSFO(ParamStartR,"RT");
CritStart <- CritOptim;
NRuns <- NRuns+nrow(CandidatesParamR);
if(!quiet){
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=""));
cat(paste("\t Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritStart*Multiplier,format="f",digits=4),"\n",sep=""));
}
##Results_archiving________________________________________________________
HistParamR[1,] <- ParamStartR;
HistParamT[1,] <- ParamStartT;
HistCrit[1,] <- CritStart;
##_____Steepest_Descent_Local_Search_______________________________________________________
##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
ProposeCandidatesLoc <- function(NewParamOptimT,OldParamOptimT,RangesT,OptimParam,Pace){
##Format_checking
if(nrow(NewParamOptimT)!=1 | nrow(OldParamOptimT)!=1){ stop("each input set must be a matrix of one single line \n"); return(NULL); }
if(ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT)!=length(OptimParam)){ stop("each input set must have the same number of values \n"); return(NULL); }
##Proposal_of_new_parameter_sets ###(local search providing 2*NParam-1 new sets)
NParam <- ncol(NewParamOptimT);
VECT <- NULL;
for(I in 1:NParam){
##We_check_that_the_current_parameter_should_indeed_be_optimised
if(OptimParam[I]==TRUE){
for(J in 1:2){
Sign <- 2*J-3; #Sign can be equal to -1 or +1
##We_define_the_new_potential_candidate
Add <- TRUE;
PotentialCandidateT <- NewParamOptimT;
PotentialCandidateT[1,I] <- NewParamOptimT[I]+Sign*Pace;
##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
if(PotentialCandidateT[1,I]<RangesT[1,I]){ PotentialCandidateT[1,I] <- RangesT[1,I]; }
if(PotentialCandidateT[1,I]>RangesT[2,I]){ PotentialCandidateT[1,I] <- RangesT[2,I]; }
##We_check_the_set_is_not_outside_the_range_of_possible_values
if( NewParamOptimT[I]==RangesT[1,I] & Sign<0 ){ Add <- FALSE; }
if( NewParamOptimT[I]==RangesT[2,I] & Sign>0 ){ Add <- FALSE; }
##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
if(identical(PotentialCandidateT,OldParamOptimT)){ Add <- FALSE; }
##We_add_the_candidate_to_our_list
if(Add==TRUE){ VECT <- c(VECT,PotentialCandidateT); }
}
}
}
Output <- NULL;
Output$NewCandidatesT <- matrix(VECT,ncol=NParam,byrow=TRUE);
return(Output);
}
##Initialisation_of_variables
if(!quiet){
cat("\t Steepest-descent local search in progress \n");
}
Pace <- 0.64;
PaceDiag <- rep(0,NParam);
CLG <- 0.7^(1/NParam);
Compt <- 0;
CritOptim <- CritStart;
##Conversion_of_real_parameter_values
RangesR <- CalibOptions$SearchRanges;
RangesT <- FUN_TRANSFO(RangesR,"RT");
NewParamOptimT <- ParamStartT;
OldParamOptimT <- ParamStartT;
##START_LOOP_ITER_________________________________________________________
for(ITER in 1:(100*NParam)){
##Exit_loop_when_Pace_becomes_too_small___________________________________
if(Pace<0.01){ break; }
##Creation_of_new_candidates______________________________________________
CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT,OldParamOptimT,RangesT,CalibOptions$OptimParam,Pace)$NewCandidatesT;
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR");
##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam]; return(x); });
if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); }
##Loop_to_test_the_various_candidates_____________________________________
iNewOptim <- 0;
for(iNew in 1:nrow(CandidatesParamR)){
##Model_run
Param <- CandidatesParamR[iNew,];
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);
if(!is.na(OutputsCrit$CritValue)){ if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier;
iNewOptim <- iNew;
} }
}
NRuns <- NRuns+nrow(CandidatesParamR);
##When_a_progress_has_been_achieved_______________________________________
if(iNewOptim!=0){
##We_store_the_optimal_set
OldParamOptimT <- NewParamOptimT;
NewParamOptimT <- matrix(CandidatesParamT[iNewOptim,1:NParam],nrow=1);
Compt <- Compt+1;
##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
if(Compt>2*NParam){
Pace <- Pace*2;
Compt <- 0;
}
##We_update_PaceDiag
VectPace <- NewParamOptimT-OldParamOptimT;
for(iC in 1:NParam){ if(CalibOptions$OptimParam[iC]==TRUE){
if(VectPace[iC]!=0){ PaceDiag[iC] <- CLG*PaceDiag[iC]+(1-CLG)*VectPace[iC]; }
if(VectPace[iC]==0){ PaceDiag[iC] <- CLG*PaceDiag[iC]; }
} }
} else {
##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
Pace <- Pace/2;
Compt <- 0;
}
##Test_of_an_additional_candidate_using_diagonal_progress_________________
if(ITER>4*NParam){
NRuns <- NRuns+1;
iNewOptim <- 0; iNew <- 1;
CandidatesParamT <- NewParamOptimT+PaceDiag; if(!is.matrix(CandidatesParamT)){ CandidatesParamT <- matrix(CandidatesParamT,nrow=1); }
##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
for(iC in 1:NParam){ if(CalibOptions$OptimParam[iC]==TRUE){
if(CandidatesParamT[iNew,iC]<RangesT[1,iC]){ CandidatesParamT[iNew,iC] <- RangesT[1,iC]; }
if(CandidatesParamT[iNew,iC]>RangesT[2,iC]){ CandidatesParamT[iNew,iC] <- RangesT[2,iC]; }
} }
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR");
##Model_run
Param <- CandidatesParamR[iNew,];
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);
if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier;
iNewOptim <- iNew;
}
##When_a_progress_has_been_achieved
if(iNewOptim!=0){
OldParamOptimT <- NewParamOptimT;
NewParamOptimT <- matrix(CandidatesParamT[iNewOptim,1:NParam],nrow=1);
}
}
##Results_archiving_______________________________________________________
NewParamOptimR <- FUN_TRANSFO(NewParamOptimT,"TR");
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="")); }
} ##END_LOOP_ITER_________________________________________________________