Commit 7d2b9422 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

v1.1.1.0 NEW: ErrorCrit can compute single, multiple or composite criterion

Showing with 135 additions and 52 deletions
+135 -52
Package: airGR Package: airGR
Type: Package Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.1.0.0 Version: 1.1.1.0
Date: 2018-10-17 Date: 2018-10-22
Authors@R: c( Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")), person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")), person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
......
ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT, warnings = TRUE, verbose = TRUE){ ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE) {
return( FUN_CRIT(InputsCrit,OutputsModel, warnings = warnings, verbose = verbose) )
## ---------- Arguments check
if (!inherits(InputsCrit, "InputsCrit")) {
stop("InputsCrit must be of class 'InputsCrit' \n")
return(NULL)
}
if (!inherits(OutputsModel, "OutputsModel")) {
stop("OutputsModel must be of class 'OutputsModel' \n")
return(NULL)
}
if (!missing(FUN_CRIT)) {
warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object", call. = FALSE)
}
## ---------- Criterion computation
## ----- Single criterion
if (inherits(InputsCrit, "Single")) {
OutputsCrit <- InputsCrit$FUN_CRIT(InputsCrit,
OutputsModel,
warnings = warnings,
verbose = verbose)
}
## ----- Multiple criteria or Composite criterion
if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
listOutputsCrit <- lapply(InputsCrit, FUN = function(iInputsCrit) {
iInputsCrit$FUN_CRIT(iInputsCrit, OutputsModel, warnings = warnings, verbose = verbose)
})
listValCrit <- sapply(listOutputsCrit, function(x) x[["CritValue"]])
listNameCrit <- sapply(listOutputsCrit, function(x) x[["CritName"]])
listweights <- unlist(lapply(InputsCrit, function(x) x[["weights"]]))
if (inherits(InputsCrit, "Compo")) {
CritValue <- sum(listValCrit * (listweights / sum(listweights)))
OutputsCritCompo <- list(MultiCritValues = listValCrit,
MultiCritNames = listNameCrit,
MultiCritWeights = listweights)
OutputsCrit <- list(CritValue = CritValue,
CritName = "Composite",
CritBestValue = +1,
Multiplier = -1,
Ind_notcomputed = NULL,
CritCompo = OutputsCritCompo,
MultiCrit = listOutputsCrit)
class(OutputsCrit) <- c("Compo", "ErrorCrit")
if (verbose) {
message("------------------------------------\n")
message("Crit. Composite = ", sprintf("%.4f", CritValue))
msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
msgForm <- unlist(strsplit(msgForm, split = ","))
msgFormSep <- rep(c(",", ",", ",\n\t\t "), times = ceiling(length(msgForm)/3))[1: length(msgForm)]
msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
msgForm <- gsub("\\,\\\n\\\t\\\t $|\\,$", "", msgForm)
message("\tFormula: mean(", msgForm, ")\n")
}
} else {
OutputsCrit <- listOutputsCrit
class(OutputsCrit) <- c("Multi", "ErrorCrit")
}
}
return(OutputsCrit)
} }
...@@ -18,21 +18,43 @@ ErrorCrit(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE) ...@@ -18,21 +18,43 @@ ErrorCrit(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE)
\item{OutputsModel}{[object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details} \item{OutputsModel}{[object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details}
\item{FUN_CRIT}{[function] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}})} \item{FUN_CRIT}{(deprecated) [function] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}})}
\item{warnings}{(optional) [boolean] boolean indicating if the warning messages are shown, default = \code{TRUE}} \item{warnings}{(optional) [boolean] boolean indicating if the warning messages are shown, default = \code{TRUE}}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}} \item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
} }
\value{
[list] list containing the function outputs, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details
}
\description{ \description{
Function which computes an error criterion with the provided function. Function which computes an error criterion with the provided function.
} }
\value{
If \code{InputsCrit} is of class \emph{Single}:
\tabular{ll}{
[list] containing the \code{ErrorCrit_*} functions outputs, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details \cr
}
If \code{InputsCrit} is of class \emph{Multi}:
\tabular{ll}{
[list] of list containing the \code{ErrorCrit_*} functions outputs, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details \cr
}
If \code{InputsCrit} is of class \emph{Compo}:
\tabular{ll}{
\emph{$CritValue } \tab [numeric] value of the composite criterion \cr
\emph{$CritName } \tab [character] name of the composite 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{$CritCompo$MultiCritValues} \tab [numeric] values of the sub-criteria \cr
\emph{$CritCompo$MultiCritNames} \tab [numeric] names of the sub-criteria \cr
\emph{$CritCompo$MultiCritWeights} \tab [character] weighted values of the sub-criteria \cr
\emph{$MultiCrit } \tab [list] of list containing the \code{ErrorCrit_*} functions outputs, see \code{\link{ErrorCrit_NSE}} or \code{\link{ErrorCrit_KGE}} for details \cr
}
}
\examples{ \examples{
library(airGR) library(airGR)
...@@ -40,65 +62,55 @@ library(airGR) ...@@ -40,65 +62,55 @@ library(airGR)
## loading catchment data ## loading catchment data
data(L0123001) data(L0123001)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR, InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E) Precip = BasinObs$P, PotEvap = BasinObs$E)
## run period selection ## calibration period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y")=="01/01/1990"), Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y")=="01/01/1990"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y")=="31/12/1999")) which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y")=="31/12/1999"))
## preparation of the RunOptions object ## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run) IndPeriod_Run = Ind_Run)
## simulation ## simulation
Param <- c(734.568, -0.840, 109.809, 1.971) Param <- c(257.238, 1.012, 88.235, 2.208)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions, OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
Param = Param, FUN = RunModel_GR4J)
## single efficiency criterion: Nash-Sutcliffe Efficiency
## efficiency criterion: Nash-Sutcliffe Efficiency InputsCritSingle <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE,
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, InputsModel = InputsModel, RunOptions = RunOptions,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run]) obs = list(BasinObs$Qmm[Ind_Run]),
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel) varObs = "Qobs", transfo = "",
weight = NULL)
## efficiency criterion: Nash-Sutcliffe Efficiency on log-transformed flows str(ErrorCrit(InputsCrit = InputsCritSingle, OutputsModel = OutputsModel))
transfo <- "log"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, ## 2 efficiency critera: RMSE and the Nash-Sutcliffe Efficiency
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run], InputsCritMulti <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_RMSE, ErrorCrit_NSE),
transfo = transfo) InputsModel = InputsModel, RunOptions = RunOptions,
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel) obs = list(BasinObs$Qmm[Ind_Run], BasinObs$Qmm[Ind_Run]),
varObs = list("Qobs", "Qobs"), transfo = list("", "sqrt"),
## efficiency criterion: Nash-Sutcliffe Efficiency above a threshold (q75\%) weight = NULL)
BoolCrit <- rep(TRUE, length(BasinObs$Qmm[Ind_Run])) str(ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel))
BoolCrit[BasinObs$Qmm[Ind_Run]<quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)] <- FALSE
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, ## efficiency composite criterion: Nash-Sutcliffe Efficiency mixing
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run], ## both raw and log-transformed flows
BoolCrit = BoolCrit) InputsCritCompo <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_NSE, ErrorCrit_NSE),
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel) InputsModel = InputsModel, RunOptions = RunOptions,
obs = list(BasinObs$Qmm[Ind_Run], BasinObs$Qmm[Ind_Run]),
## efficiency criterion: Kling-Gupta Efficiency varObs = list("Qobs", "Qobs"), transfo = list("", "log"),
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel, weight = list(0.4, 0.6))
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run]) str(ErrorCrit(InputsCrit = InputsCritCompo, OutputsModel = OutputsModel))
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency below a threshold (q10\%) on log-transformed flows
transfo <- "log"
BoolCrit <- rep(TRUE, length(BasinObs$Qmm[Ind_Run]))
BoolCrit[BasinObs$Qmm[Ind_Run]>quantile(BasinObs$Qmm[Ind_Run], 0.10, na.rm = TRUE)] <- FALSE
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit, transfo = transfo)
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
} }
\author{ \author{
Laurent Coron Olivier Delaigue
} }
\seealso{ \seealso{
\code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_KGE2}} \code{\link{CreateInputsCrit}}, \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_KGE2}}
} }
Supports Markdown
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