diff --git a/DESCRIPTION b/DESCRIPTION index da6556c119c3217ef736dfc15d8320ebb12dc19f..ca9bb5c5d543fcdd5f46cde168ffe5a9fae6e98b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: airGR Type: Package Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling -Version: 1.1.0.0 -Date: 2018-10-17 +Version: 1.1.1.0 +Date: 2018-10-22 Authors@R: c( 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")), diff --git a/R/ErrorCrit.R b/R/ErrorCrit.R index 06890369d2651c596b770be7df97f2f54ef14c62..7bbd69fb1df2c54a63decfd64da50bc41aa5e0de 100644 --- a/R/ErrorCrit.R +++ b/R/ErrorCrit.R @@ -1,4 +1,75 @@ -ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT, warnings = TRUE, verbose = TRUE){ - return( FUN_CRIT(InputsCrit,OutputsModel, warnings = warnings, verbose = verbose) ) +ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE) { + + ## ---------- 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) + } diff --git a/man/ErrorCrit.Rd b/man/ErrorCrit.Rd index 47e8ce17396ea3c118f0418f8708f5275d067026..26c7bb3f846e4735854f46508462c6b8379825d5 100644 --- a/man/ErrorCrit.Rd +++ b/man/ErrorCrit.Rd @@ -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{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{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{ 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{ library(airGR) @@ -40,65 +62,55 @@ library(airGR) ## loading catchment data data(L0123001) -## preparation of the InputsModel object + InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR, 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"), which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y")=="31/12/1999")) -## preparation of the RunOptions object -RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, - InputsModel = InputsModel, IndPeriod_Run = Ind_Run) +## preparation of RunOptions object +RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, + IndPeriod_Run = Ind_Run) ## simulation -Param <- c(734.568, -0.840, 109.809, 1.971) -OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions, - Param = Param, FUN = RunModel_GR4J) - -## 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) - -## efficiency criterion: Nash-Sutcliffe Efficiency on log-transformed flows -transfo <- "log" -InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, - RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run], - transfo = transfo) -OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel) - -## efficiency criterion: Nash-Sutcliffe Efficiency above a threshold (q75\%) -BoolCrit <- rep(TRUE, length(BasinObs$Qmm[Ind_Run])) -BoolCrit[BasinObs$Qmm[Ind_Run]<quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)] <- FALSE -InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, - RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run], - BoolCrit = BoolCrit) -OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel) - -## 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) - -## 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) +Param <- c(257.238, 1.012, 88.235, 2.208) +OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param) + +## single efficiency criterion: Nash-Sutcliffe Efficiency +InputsCritSingle <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, + InputsModel = InputsModel, RunOptions = RunOptions, + obs = list(BasinObs$Qmm[Ind_Run]), + varObs = "Qobs", transfo = "", + weight = NULL) +str(ErrorCrit(InputsCrit = InputsCritSingle, OutputsModel = OutputsModel)) + +## 2 efficiency critera: RMSE and the Nash-Sutcliffe Efficiency +InputsCritMulti <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_RMSE, ErrorCrit_NSE), + InputsModel = InputsModel, RunOptions = RunOptions, + obs = list(BasinObs$Qmm[Ind_Run], BasinObs$Qmm[Ind_Run]), + varObs = list("Qobs", "Qobs"), transfo = list("", "sqrt"), + weight = NULL) +str(ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel)) + +## efficiency composite criterion: Nash-Sutcliffe Efficiency mixing +## both raw and log-transformed flows +InputsCritCompo <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_NSE, ErrorCrit_NSE), + InputsModel = InputsModel, RunOptions = RunOptions, + obs = list(BasinObs$Qmm[Ind_Run], BasinObs$Qmm[Ind_Run]), + varObs = list("Qobs", "Qobs"), transfo = list("", "log"), + weight = list(0.4, 0.6)) +str(ErrorCrit(InputsCrit = InputsCritCompo, OutputsModel = OutputsModel)) } \author{ -Laurent Coron +Olivier Delaigue } \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}} }