From 569431360f96d103a8ab1fa2cc3e84d37e4579bd Mon Sep 17 00:00:00 2001 From: Delaigue Olivier <olivier.delaigue@irstea.priv> Date: Thu, 21 Feb 2019 15:41:01 +0100 Subject: [PATCH] v1.1.3.9 CLEAN: indent code of ErrorCrit_* fun --- DESCRIPTION | 2 +- NEWS.rmd | 2 +- R/ErrorCrit_KGE.R | 408 +++++++++++++++++++++--------------------- R/ErrorCrit_KGE2.R | 434 ++++++++++++++++++++++----------------------- R/ErrorCrit_NSE.R | 282 ++++++++++++++--------------- R/ErrorCrit_RMSE.R | 282 ++++++++++++++--------------- 6 files changed, 705 insertions(+), 705 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6dec702c..6e5ec2fe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: airGR Type: Package Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling -Version: 1.1.3.8 +Version: 1.1.3.9 Date: 2019-02-21 Authors@R: c( person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")), diff --git a/NEWS.rmd b/NEWS.rmd index 79d50a96..96d8dec0 100644 --- a/NEWS.rmd +++ b/NEWS.rmd @@ -13,7 +13,7 @@ output: -### 1.1.3.8 Release Notes (2019-02-21) +### 1.1.3.9 Release Notes (2019-02-21) diff --git a/R/ErrorCrit_KGE.R b/R/ErrorCrit_KGE.R index fccc6158..a30ac999 100644 --- a/R/ErrorCrit_KGE.R +++ b/R/ErrorCrit_KGE.R @@ -1,211 +1,211 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) { - ##Arguments_check________________________________ - if (!inherits(InputsCrit, "InputsCrit")) { - stop("InputsCrit must be of class 'InputsCrit' \n") - return(NULL) - } - if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { - stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' or 'Compo' with KGE") - return(NULL) - } - if (!inherits(OutputsModel, "OutputsModel")) { - stop("OutputsModel must be of class 'OutputsModel' \n") - return(NULL) - } - - - ##Initialisation_________________________________ - CritName <- NA - CritVar <- InputsCrit$varObs - if (InputsCrit$transfo == "") { - CritName <- "KGE[CritVar]" - } - if (InputsCrit$transfo == "sqrt") { - CritName <- "KGE[sqrt(CritVar)]" - } - if (InputsCrit$transfo == "log") { - CritName <- "KGE[log(CritVar)]" - } - if (InputsCrit$transfo == "inv") { - CritName <- "KGE[1/CritVar]" - } - if (InputsCrit$transfo == "sort") { - CritName <- "KGE[sort(CritVar)]" - } - CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) - CritValue <- NA - CritBestValue <- +1 - Multiplier <- -1 - ### must be equal to -1 or +1 only - - - ##Data_preparation_______________________________ - VarObs <- InputsCrit$obs - VarObs[!InputsCrit$BoolCrit] <- NA - VarSim <- OutputsModel$Qsim - VarSim[!InputsCrit$BoolCrit] <- NA - - ##Data_transformation - if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { - if (any(VarObs %in% 0)) { - warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - if (any(VarSim %in% 0)) { - warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - } - if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { - VarObs <- VarObs + InputsCrit$epsilon - VarSim <- VarSim + InputsCrit$epsilon - } - if (InputsCrit$transfo == "sqrt") { - VarObs <- sqrt(VarObs) - VarSim <- sqrt(VarSim) - } - if (InputsCrit$transfo == "log") { - VarObs <- log(VarObs) - VarSim <- log(VarSim) - VarSim[VarSim < -1e100] <- NA - } - if (InputsCrit$transfo == "inv") { - VarObs <- 1 / VarObs - VarSim <- 1 / VarSim - VarSim[abs(VarSim) > 1e+100] <- NA - } - if (InputsCrit$transfo == "sort") { - VarSim[is.na(VarObs)] <- NA - VarSim <- sort(VarSim, na.last = TRUE) - VarObs <- sort(VarObs, na.last = TRUE) - InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE) - } - - ##TS_ignore - TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit - Ind_TS_ignore <- which(TS_ignore) - if (length(Ind_TS_ignore) == 0) { - Ind_TS_ignore <- NULL - } - if (sum(!TS_ignore) == 0) { - OutputsCrit <- list(NA) - names(OutputsCrit) <- c("CritValue") - return(OutputsCrit) - } - if (sum(!TS_ignore) == 1) { - OutputsCrit <- list(NA) - names(OutputsCrit) <- c("CritValue") - return(OutputsCrit) - } ### to avoid a problem in standard deviation computation - if (inherits(OutputsModel, "hourly")) { - WarningTS <- 365 - } - if (inherits(OutputsModel, "daily")) { - WarningTS <- 365 - } - if (inherits(OutputsModel, "monthly")) { - WarningTS <- 12 - } - if (inherits(OutputsModel, "yearly")) { - WarningTS <- 3 - } - if (sum(!TS_ignore) < WarningTS & warnings) { - warning("\t criterion computed on less than ", WarningTS, " time-steps ") - } - - ##Other_variables_preparation - meanVarObs <- mean(VarObs[!TS_ignore]) - meanVarSim <- mean(VarSim[!TS_ignore]) - - iCrit <- 0 - SubCritPrint <- NULL - SubCritNames <- NULL - SubCritValues <- NULL - - - ##SubErrorCrit_____KGE_rPearson__________________ - iCrit <- iCrit + 1 - SubCritPrint[iCrit] <- paste(CritName, " cor(sim, obs, \"pearson\") =", sep = "") - SubCritValues[iCrit] <- NA - SubCritNames[iCrit] <- "r" - - Numer <- sum((VarObs[!TS_ignore] - meanVarObs) * (VarSim[!TS_ignore] - meanVarSim)) - Deno1 <- sqrt(sum((VarObs[!TS_ignore] - meanVarObs) ^ 2)) - Deno2 <- sqrt(sum((VarSim[!TS_ignore] - meanVarSim) ^ 2)) - - if (Numer == 0) { - if (Deno1 == 0 & Deno2 == 0) { - Crit <- 1 - } else { - Crit <- 0 - } - } else { - Crit <- Numer / (Deno1 * Deno2) - } - if (is.numeric(Crit) & is.finite(Crit)) { - SubCritValues[iCrit] <- Crit - } - - - ##SubErrorCrit_____KGE_alpha_____________________ - iCrit <- iCrit + 1 - SubCritPrint[iCrit] <- paste(CritName, " sd(sim)/sd(obs) =", sep = "") - SubCritValues[iCrit] <- NA - SubCritNames[iCrit] <- "alpha" - - Numer <- sd(VarSim[!TS_ignore]) - Denom <- sd(VarObs[!TS_ignore]) - - if (Numer == 0 & Denom == 0) { - Crit <- 1 - } else { - Crit <- Numer / Denom - } - if (is.numeric(Crit) & is.finite(Crit)) { - SubCritValues[iCrit] <- Crit - } - - - ##SubErrorCrit_____KGE_beta______________________ - iCrit <- iCrit + 1 - SubCritPrint[iCrit] <- paste(CritName, " mean(sim)/mean(obs) =", sep = "") - SubCritValues[iCrit] <- NA - SubCritNames[iCrit] <- "beta" - - if (meanVarSim == 0 & meanVarObs == 0) { + ##Arguments_check________________________________ + if (!inherits(InputsCrit, "InputsCrit")) { + stop("InputsCrit must be of class 'InputsCrit' \n") + return(NULL) + } + if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { + stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' or 'Compo' with KGE") + return(NULL) + } + if (!inherits(OutputsModel, "OutputsModel")) { + stop("OutputsModel must be of class 'OutputsModel' \n") + return(NULL) + } + + + ##Initialisation_________________________________ + CritName <- NA + CritVar <- InputsCrit$varObs + if (InputsCrit$transfo == "") { + CritName <- "KGE[CritVar]" + } + if (InputsCrit$transfo == "sqrt") { + CritName <- "KGE[sqrt(CritVar)]" + } + if (InputsCrit$transfo == "log") { + CritName <- "KGE[log(CritVar)]" + } + if (InputsCrit$transfo == "inv") { + CritName <- "KGE[1/CritVar]" + } + if (InputsCrit$transfo == "sort") { + CritName <- "KGE[sort(CritVar)]" + } + CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) + CritValue <- NA + CritBestValue <- +1 + Multiplier <- -1 + ### must be equal to -1 or +1 only + + + ##Data_preparation_______________________________ + VarObs <- InputsCrit$obs + VarObs[!InputsCrit$BoolCrit] <- NA + VarSim <- OutputsModel$Qsim + VarSim[!InputsCrit$BoolCrit] <- NA + + ##Data_transformation + if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { + if (any(VarObs %in% 0)) { + warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + if (any(VarSim %in% 0)) { + warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + } + if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { + VarObs <- VarObs + InputsCrit$epsilon + VarSim <- VarSim + InputsCrit$epsilon + } + if (InputsCrit$transfo == "sqrt") { + VarObs <- sqrt(VarObs) + VarSim <- sqrt(VarSim) + } + if (InputsCrit$transfo == "log") { + VarObs <- log(VarObs) + VarSim <- log(VarSim) + VarSim[VarSim < -1e100] <- NA + } + if (InputsCrit$transfo == "inv") { + VarObs <- 1 / VarObs + VarSim <- 1 / VarSim + VarSim[abs(VarSim) > 1e+100] <- NA + } + if (InputsCrit$transfo == "sort") { + VarSim[is.na(VarObs)] <- NA + VarSim <- sort(VarSim, na.last = TRUE) + VarObs <- sort(VarObs, na.last = TRUE) + InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE) + } + + ##TS_ignore + TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit + Ind_TS_ignore <- which(TS_ignore) + if (length(Ind_TS_ignore) == 0) { + Ind_TS_ignore <- NULL + } + if (sum(!TS_ignore) == 0) { + OutputsCrit <- list(NA) + names(OutputsCrit) <- c("CritValue") + return(OutputsCrit) + } + if (sum(!TS_ignore) == 1) { + OutputsCrit <- list(NA) + names(OutputsCrit) <- c("CritValue") + return(OutputsCrit) + } ### to avoid a problem in standard deviation computation + if (inherits(OutputsModel, "hourly")) { + WarningTS <- 365 + } + if (inherits(OutputsModel, "daily")) { + WarningTS <- 365 + } + if (inherits(OutputsModel, "monthly")) { + WarningTS <- 12 + } + if (inherits(OutputsModel, "yearly")) { + WarningTS <- 3 + } + if (sum(!TS_ignore) < WarningTS & warnings) { + warning("\t criterion computed on less than ", WarningTS, " time-steps ") + } + + ##Other_variables_preparation + meanVarObs <- mean(VarObs[!TS_ignore]) + meanVarSim <- mean(VarSim[!TS_ignore]) + + iCrit <- 0 + SubCritPrint <- NULL + SubCritNames <- NULL + SubCritValues <- NULL + + + ##SubErrorCrit_____KGE_rPearson__________________ + iCrit <- iCrit + 1 + SubCritPrint[iCrit] <- paste(CritName, " cor(sim, obs, \"pearson\") =", sep = "") + SubCritValues[iCrit] <- NA + SubCritNames[iCrit] <- "r" + + Numer <- sum((VarObs[!TS_ignore] - meanVarObs) * (VarSim[!TS_ignore] - meanVarSim)) + Deno1 <- sqrt(sum((VarObs[!TS_ignore] - meanVarObs) ^ 2)) + Deno2 <- sqrt(sum((VarSim[!TS_ignore] - meanVarSim) ^ 2)) + + if (Numer == 0) { + if (Deno1 == 0 & Deno2 == 0) { Crit <- 1 } else { - Crit <- meanVarSim / meanVarObs - } - if (is.numeric(Crit) & is.finite(Crit)) { - SubCritValues[iCrit] <- Crit - } - - - ##ErrorCrit______________________________________ - if (sum(is.na(SubCritValues)) == 0) { - CritValue <- (1 - sqrt((SubCritValues[1] - 1)^2 + (SubCritValues[2] - 1)^2 + (SubCritValues[3] - 1)^2)) + Crit <- 0 } - - - ##Verbose______________________________________ - if (verbose) { - message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) - message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " ")) - } - - - ##Output_________________________________________ - OutputsCrit <- list(CritValue = CritValue, - CritName = CritName, - SubCritValues = SubCritValues, - SubCritNames = SubCritNames, - CritBestValue = CritBestValue, - Multiplier = Multiplier, - Ind_notcomputed = Ind_TS_ignore - ) - - class(OutputsCrit) <- c("KGE", "ErrorCrit") - return(OutputsCrit) - + } else { + Crit <- Numer / (Deno1 * Deno2) } + if (is.numeric(Crit) & is.finite(Crit)) { + SubCritValues[iCrit] <- Crit + } + + + ##SubErrorCrit_____KGE_alpha_____________________ + iCrit <- iCrit + 1 + SubCritPrint[iCrit] <- paste(CritName, " sd(sim)/sd(obs) =", sep = "") + SubCritValues[iCrit] <- NA + SubCritNames[iCrit] <- "alpha" + + Numer <- sd(VarSim[!TS_ignore]) + Denom <- sd(VarObs[!TS_ignore]) + + if (Numer == 0 & Denom == 0) { + Crit <- 1 + } else { + Crit <- Numer / Denom + } + if (is.numeric(Crit) & is.finite(Crit)) { + SubCritValues[iCrit] <- Crit + } + + + ##SubErrorCrit_____KGE_beta______________________ + iCrit <- iCrit + 1 + SubCritPrint[iCrit] <- paste(CritName, " mean(sim)/mean(obs) =", sep = "") + SubCritValues[iCrit] <- NA + SubCritNames[iCrit] <- "beta" + + if (meanVarSim == 0 & meanVarObs == 0) { + Crit <- 1 + } else { + Crit <- meanVarSim / meanVarObs + } + if (is.numeric(Crit) & is.finite(Crit)) { + SubCritValues[iCrit] <- Crit + } + + + ##ErrorCrit______________________________________ + if (sum(is.na(SubCritValues)) == 0) { + CritValue <- (1 - sqrt((SubCritValues[1] - 1)^2 + (SubCritValues[2] - 1)^2 + (SubCritValues[3] - 1)^2)) + } + + + ##Verbose______________________________________ + if (verbose) { + message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) + message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " ")) + } + + + ##Output_________________________________________ + OutputsCrit <- list(CritValue = CritValue, + CritName = CritName, + SubCritValues = SubCritValues, + SubCritNames = SubCritNames, + CritBestValue = CritBestValue, + Multiplier = Multiplier, + Ind_notcomputed = Ind_TS_ignore + ) + + class(OutputsCrit) <- c("KGE", "ErrorCrit") + return(OutputsCrit) + +} diff --git a/R/ErrorCrit_KGE2.R b/R/ErrorCrit_KGE2.R index 96237ea2..5e08b2f0 100644 --- a/R/ErrorCrit_KGE2.R +++ b/R/ErrorCrit_KGE2.R @@ -1,228 +1,228 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) { - ##Arguments_check________________________________ - if (!inherits(InputsCrit, "InputsCrit")) { - stop("InputsCrit must be of class 'InputsCrit' \n") - return(NULL) - } - if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { - stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' or 'Compo' with KGE'") - return(NULL) - } - if (!inherits(OutputsModel, "OutputsModel")) { - stop("OutputsModel must be of class 'OutputsModel' \n") - return(NULL) - } - - - ##Initialisation_________________________________ - CritName <- NA - CritVar <- InputsCrit$varObs - if (InputsCrit$transfo == "") { - CritName <- "KGE'[CritVar]" - } - if (InputsCrit$transfo == "sqrt") { - CritName <- "KGE'[sqrt(CritVar)]" - } - if (InputsCrit$transfo == "log") { - CritName <- "KGE'[log(CritVar)]" - } - if (InputsCrit$transfo == "inv") { - CritName <- "KGE'[1/CritVar]" - } - if (InputsCrit$transfo == "sort") { - CritName <- "KGE'[sort(CritVar)]" - } - CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) - CritValue <- NA - CritBestValue <- +1 - Multiplier <- -1 - ### must be equal to -1 or +1 only - - - ##Data_preparation_______________________________ - VarObs <- InputsCrit$obs - VarObs[!InputsCrit$BoolCrit] <- NA - VarSim <- OutputsModel$Qsim - VarSim[!InputsCrit$BoolCrit] <- NA - - ##Data_transformation - if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { - if (any(VarObs %in% 0)) { - warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - if (any(VarSim %in% 0)) { - warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - } - if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { - VarObs <- VarObs + InputsCrit$epsilon - VarSim <- VarSim + InputsCrit$epsilon - } - if (InputsCrit$transfo == "sqrt") { - VarObs <- sqrt(VarObs) - VarSim <- sqrt(VarSim) - } - if (InputsCrit$transfo == "log") { - VarObs <- log(VarObs) - VarSim <- log(VarSim) - VarSim[VarSim < -1e100] <- NA - } - if (InputsCrit$transfo == "inv") { - VarObs <- 1 / VarObs - VarSim <- 1 / VarSim - VarSim[abs(VarSim) > 1e+100] <- NA - } - if (InputsCrit$transfo == "sort") { - VarSim[is.na(VarObs)] <- NA - VarSim <- sort(VarSim, na.last = TRUE) - VarObs <- sort(VarObs, na.last = TRUE) - InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE) - } - - ##TS_ignore - TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit - Ind_TS_ignore <- which(TS_ignore) - if (length(Ind_TS_ignore) == 0) { - Ind_TS_ignore <- NULL - } - if (sum(!TS_ignore) == 0) { - OutputsCrit <- list(NA) - names(OutputsCrit) <- c("CritValue") - return(OutputsCrit) - } - if (sum(!TS_ignore) == 1) { - OutputsCrit <- list(NA) - names(OutputsCrit) <- c("CritValue") - return(OutputsCrit) - } ### to avoid a problem in standard deviation computation - if (inherits(OutputsModel, "hourly")) { - WarningTS <- 365 - } - if (inherits(OutputsModel, "daily")) { - WarningTS <- 365 - } - if (inherits(OutputsModel, "monthly")) { - WarningTS <- 12 - } - if (inherits(OutputsModel, "yearly")) { - WarningTS <- 3 - } - if (sum(!TS_ignore) < WarningTS & warnings) { - warning("\t criterion computed on less than ", WarningTS, " time-steps") - } - - ##Other_variables_preparation - meanVarObs <- mean(VarObs[!TS_ignore]) - meanVarSim <- mean(VarSim[!TS_ignore]) - - iCrit <- 0 - SubCritPrint <- NULL - SubCritNames <- NULL - SubCritValues <- NULL - - - ##SubErrorCrit_____KGE_rPearson__________________ - iCrit <- iCrit + 1 - SubCritPrint[iCrit] <- paste(CritName, " cor(sim, obs, \"pearson\") =", sep = "") - SubCritValues[iCrit] <- NA - SubCritNames[iCrit] <- "r" - - Numer <- sum((VarObs[!TS_ignore] - meanVarObs) * (VarSim[!TS_ignore] - meanVarSim)) - Deno1 <- sqrt(sum((VarObs[!TS_ignore] - meanVarObs)^2)) - Deno2 <- sqrt(sum((VarSim[!TS_ignore] - meanVarSim)^2)) - - if (Numer == 0) { - if (Deno1 == 0 & Deno2 == 0) { - Crit <- 1 - } else { - Crit <- 0 - } - } else { - Crit <- Numer / (Deno1 * Deno2) - } - if (is.numeric(Crit) & is.finite(Crit)) { - SubCritValues[iCrit] <- Crit - } - - - ##SubErrorCrit_____KGE_gamma______________________ - iCrit <- iCrit + 1 - SubCritPrint[iCrit] <- paste(CritName, " cv(sim)/cv(obs) =", sep = "") - SubCritValues[iCrit] <- NA - SubCritNames[iCrit] <- "gamma" - - if (meanVarSim == 0) { - if (sd(VarSim[!TS_ignore]) == 0) { - CVsim <- 1 - } else { - CVsim <- 99999 - } - } else { - CVsim <- sd(VarSim[!TS_ignore]) / meanVarSim - - } - if (meanVarObs == 0) { - if (sd(VarObs[!TS_ignore]) == 0) { - CVobs <- 1 - } else { - CVobs <- 99999 - } - } else { - CVobs <- sd(VarObs[!TS_ignore]) / meanVarObs - } - if (CVsim == 0 & - CVobs == 0) { + ##Arguments_check________________________________ + if (!inherits(InputsCrit, "InputsCrit")) { + stop("InputsCrit must be of class 'InputsCrit' \n") + return(NULL) + } + if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { + stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' or 'Compo' with KGE'") + return(NULL) + } + if (!inherits(OutputsModel, "OutputsModel")) { + stop("OutputsModel must be of class 'OutputsModel' \n") + return(NULL) + } + + + ##Initialisation_________________________________ + CritName <- NA + CritVar <- InputsCrit$varObs + if (InputsCrit$transfo == "") { + CritName <- "KGE'[CritVar]" + } + if (InputsCrit$transfo == "sqrt") { + CritName <- "KGE'[sqrt(CritVar)]" + } + if (InputsCrit$transfo == "log") { + CritName <- "KGE'[log(CritVar)]" + } + if (InputsCrit$transfo == "inv") { + CritName <- "KGE'[1/CritVar]" + } + if (InputsCrit$transfo == "sort") { + CritName <- "KGE'[sort(CritVar)]" + } + CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) + CritValue <- NA + CritBestValue <- +1 + Multiplier <- -1 + ### must be equal to -1 or +1 only + + + ##Data_preparation_______________________________ + VarObs <- InputsCrit$obs + VarObs[!InputsCrit$BoolCrit] <- NA + VarSim <- OutputsModel$Qsim + VarSim[!InputsCrit$BoolCrit] <- NA + + ##Data_transformation + if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { + if (any(VarObs %in% 0)) { + warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + if (any(VarSim %in% 0)) { + warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + } + if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { + VarObs <- VarObs + InputsCrit$epsilon + VarSim <- VarSim + InputsCrit$epsilon + } + if (InputsCrit$transfo == "sqrt") { + VarObs <- sqrt(VarObs) + VarSim <- sqrt(VarSim) + } + if (InputsCrit$transfo == "log") { + VarObs <- log(VarObs) + VarSim <- log(VarSim) + VarSim[VarSim < -1e100] <- NA + } + if (InputsCrit$transfo == "inv") { + VarObs <- 1 / VarObs + VarSim <- 1 / VarSim + VarSim[abs(VarSim) > 1e+100] <- NA + } + if (InputsCrit$transfo == "sort") { + VarSim[is.na(VarObs)] <- NA + VarSim <- sort(VarSim, na.last = TRUE) + VarObs <- sort(VarObs, na.last = TRUE) + InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE) + } + + ##TS_ignore + TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit + Ind_TS_ignore <- which(TS_ignore) + if (length(Ind_TS_ignore) == 0) { + Ind_TS_ignore <- NULL + } + if (sum(!TS_ignore) == 0) { + OutputsCrit <- list(NA) + names(OutputsCrit) <- c("CritValue") + return(OutputsCrit) + } + if (sum(!TS_ignore) == 1) { + OutputsCrit <- list(NA) + names(OutputsCrit) <- c("CritValue") + return(OutputsCrit) + } ### to avoid a problem in standard deviation computation + if (inherits(OutputsModel, "hourly")) { + WarningTS <- 365 + } + if (inherits(OutputsModel, "daily")) { + WarningTS <- 365 + } + if (inherits(OutputsModel, "monthly")) { + WarningTS <- 12 + } + if (inherits(OutputsModel, "yearly")) { + WarningTS <- 3 + } + if (sum(!TS_ignore) < WarningTS & warnings) { + warning("\t criterion computed on less than ", WarningTS, " time-steps") + } + + ##Other_variables_preparation + meanVarObs <- mean(VarObs[!TS_ignore]) + meanVarSim <- mean(VarSim[!TS_ignore]) + + iCrit <- 0 + SubCritPrint <- NULL + SubCritNames <- NULL + SubCritValues <- NULL + + + ##SubErrorCrit_____KGE_rPearson__________________ + iCrit <- iCrit + 1 + SubCritPrint[iCrit] <- paste(CritName, " cor(sim, obs, \"pearson\") =", sep = "") + SubCritValues[iCrit] <- NA + SubCritNames[iCrit] <- "r" + + Numer <- sum((VarObs[!TS_ignore] - meanVarObs) * (VarSim[!TS_ignore] - meanVarSim)) + Deno1 <- sqrt(sum((VarObs[!TS_ignore] - meanVarObs)^2)) + Deno2 <- sqrt(sum((VarSim[!TS_ignore] - meanVarSim)^2)) + + if (Numer == 0) { + if (Deno1 == 0 & Deno2 == 0) { Crit <- 1 } else { - Crit <- CVsim / CVobs + Crit <- 0 } - if (is.numeric(Crit) & is.finite(Crit)) { - SubCritValues[iCrit] <- Crit - } - - - ##SubErrorCrit_____KGE_beta______________________ - iCrit <- iCrit + 1 - SubCritPrint[iCrit] <- paste(CritName, " mean(sim)/mean(obs) =", sep = "") - SubCritValues[iCrit] <- NA - SubCritNames[iCrit] <- "beta" - - if (meanVarSim == 0 & meanVarObs == 0) { - Crit <- 1 + } else { + Crit <- Numer / (Deno1 * Deno2) + } + if (is.numeric(Crit) & is.finite(Crit)) { + SubCritValues[iCrit] <- Crit + } + + + ##SubErrorCrit_____KGE_gamma______________________ + iCrit <- iCrit + 1 + SubCritPrint[iCrit] <- paste(CritName, " cv(sim)/cv(obs) =", sep = "") + SubCritValues[iCrit] <- NA + SubCritNames[iCrit] <- "gamma" + + if (meanVarSim == 0) { + if (sd(VarSim[!TS_ignore]) == 0) { + CVsim <- 1 } else { - Crit <- meanVarSim / meanVarObs - } - if (is.numeric(Crit) & is.finite(Crit)) { - SubCritValues[iCrit] <- Crit - } - - - ##ErrorCrit______________________________________ - if (sum(is.na(SubCritValues)) == 0) { - CritValue <- (1 - sqrt((SubCritValues[1] - 1)^2 + (SubCritValues[2] - 1)^2 + (SubCritValues[3] - 1)^2)) + CVsim <- 99999 } + } else { + CVsim <- sd(VarSim[!TS_ignore]) / meanVarSim - - ##Verbose______________________________________ - if (verbose) { - message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) - message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " ")) + } + if (meanVarObs == 0) { + if (sd(VarObs[!TS_ignore]) == 0) { + CVobs <- 1 + } else { + CVobs <- 99999 } - - - ##Output_________________________________________ - OutputsCrit <- list(CritValue = CritValue, - CritName = CritName, - SubCritValues = SubCritValues, - SubCritNames = SubCritNames, - CritBestValue = CritBestValue, - Multiplier = Multiplier, - Ind_notcomputed = Ind_TS_ignore - ) - - class(OutputsCrit) <- c("KGE2", "ErrorCrit") - return(OutputsCrit) - + } else { + CVobs <- sd(VarObs[!TS_ignore]) / meanVarObs } + if (CVsim == 0 & + CVobs == 0) { + Crit <- 1 + } else { + Crit <- CVsim / CVobs + } + if (is.numeric(Crit) & is.finite(Crit)) { + SubCritValues[iCrit] <- Crit + } + + + ##SubErrorCrit_____KGE_beta______________________ + iCrit <- iCrit + 1 + SubCritPrint[iCrit] <- paste(CritName, " mean(sim)/mean(obs) =", sep = "") + SubCritValues[iCrit] <- NA + SubCritNames[iCrit] <- "beta" + + if (meanVarSim == 0 & meanVarObs == 0) { + Crit <- 1 + } else { + Crit <- meanVarSim / meanVarObs + } + if (is.numeric(Crit) & is.finite(Crit)) { + SubCritValues[iCrit] <- Crit + } + + + ##ErrorCrit______________________________________ + if (sum(is.na(SubCritValues)) == 0) { + CritValue <- (1 - sqrt((SubCritValues[1] - 1)^2 + (SubCritValues[2] - 1)^2 + (SubCritValues[3] - 1)^2)) + } + + + ##Verbose______________________________________ + if (verbose) { + message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) + message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " ")) + } + + + ##Output_________________________________________ + OutputsCrit <- list(CritValue = CritValue, + CritName = CritName, + SubCritValues = SubCritValues, + SubCritNames = SubCritNames, + CritBestValue = CritBestValue, + Multiplier = Multiplier, + Ind_notcomputed = Ind_TS_ignore + ) + + class(OutputsCrit) <- c("KGE2", "ErrorCrit") + return(OutputsCrit) + +} diff --git a/R/ErrorCrit_NSE.R b/R/ErrorCrit_NSE.R index f14eaa99..57e1e5c3 100644 --- a/R/ErrorCrit_NSE.R +++ b/R/ErrorCrit_NSE.R @@ -1,146 +1,146 @@ ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) { - ##Arguments_check________________________________ - if (!inherits(InputsCrit, "InputsCrit")) { - stop("InputsCrit must be of class 'InputsCrit' \n") - return(NULL) - } - if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { - stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' or 'Compo' with NSE") - return(NULL) - } - if (!inherits(OutputsModel, "OutputsModel")) { - stop("OutputsModel must be of class 'OutputsModel' \n") - return(NULL) - } - - - ##Initialisation_________________________________ - CritName <- NA - CritVar <- InputsCrit$varObs - if (InputsCrit$transfo == "") { - CritName <- "NSE[CritVar]" - } - if (InputsCrit$transfo == "sqrt") { - CritName <- "NSE[sqrt(CritVar)]" - } - if (InputsCrit$transfo == "log") { - CritName <- "NSE[log(CritVar)]" - } - if (InputsCrit$transfo == "inv") { - CritName <- "NSE[1/CritVar]" - } - if (InputsCrit$transfo == "sort") { - CritName <- "NSE[sort(CritVar)]" - } - CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) - CritValue <- NA - CritBestValue <- +1 - Multiplier <- -1 - ### must be equal to -1 or +1 only - - - ##Data_preparation_______________________________ - VarObs <- InputsCrit$obs - VarObs[!InputsCrit$BoolCrit] <- NA - VarSim <- OutputsModel$Qsim - VarSim[!InputsCrit$BoolCrit] <- NA - - ##Data_transformation - if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { - if (any(VarObs %in% 0)) { - warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - if (any(VarSim %in% 0)) { - warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - } - if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { - VarObs <- VarObs + InputsCrit$epsilon - VarSim <- VarSim + InputsCrit$epsilon - } - if (InputsCrit$transfo == "sqrt") { - VarObs <- sqrt(VarObs) - VarSim <- sqrt(VarSim) - } - if (InputsCrit$transfo == "log") { - VarObs <- log(VarObs) - VarSim <- log(VarSim) - VarSim[VarSim < -1e100] <- NA - } - if (InputsCrit$transfo == "inv") { - VarObs <- 1 / VarObs - VarSim <- 1 / VarSim - VarSim[abs(VarSim) > 1e+100] <- NA - } - if (InputsCrit$transfo == "sort") { - VarSim[is.na(VarObs)] <- NA - VarSim <- sort(VarSim, na.last = TRUE) - VarObs <- sort(VarObs, na.last = TRUE) - InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE) - } - - ##TS_ignore - TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit - Ind_TS_ignore <- which(TS_ignore) - if (length(Ind_TS_ignore) == 0) { - Ind_TS_ignore <- NULL - } - if (sum(!TS_ignore) == 0) { - OutputsCrit <- list(NA) - names(OutputsCrit) <- c("CritValue") - return(OutputsCrit) - } - if (inherits(OutputsModel, "hourly")) { - WarningTS <- 365 - } - if (inherits(OutputsModel, "daily")) { - WarningTS <- 365 - - } - if (inherits(OutputsModel, "monthly")) { - WarningTS <- 12 - } - if (inherits(OutputsModel, "yearly")) { - WarningTS <- 3 - } - if (sum(!TS_ignore) < WarningTS & warnings) { - warning("\t criterion computed on less than ", WarningTS, " time-steps") - } - - ##Other_variables_preparation - meanVarObs <- mean(VarObs[!TS_ignore]) - meanVarSim <- mean(VarSim[!TS_ignore]) - - - ##ErrorCrit______________________________________ - Emod <- sum((VarSim[!TS_ignore] - VarObs[!TS_ignore])^2) - Eref <- sum((VarObs[!TS_ignore] - mean(VarObs[!TS_ignore]))^2) - - if (Emod == 0 & Eref == 0) { - Crit <- 0 - } else { - Crit <- (1 - Emod / Eref) - } - if (is.numeric(Crit) & is.finite(Crit)) { - CritValue <- Crit - } - - - ##Verbose______________________________________ - if (verbose) { - message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue), "\n") - } - - - ##Output_________________________________________ - OutputsCrit <- list(CritValue = CritValue, - CritName = CritName, - CritBestValue = CritBestValue, - Multiplier = Multiplier, - Ind_notcomputed = Ind_TS_ignore) - - class(OutputsCrit) <- c("NSE", "ErrorCrit") + ##Arguments_check________________________________ + if (!inherits(InputsCrit, "InputsCrit")) { + stop("InputsCrit must be of class 'InputsCrit' \n") + return(NULL) + } + if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { + stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' or 'Compo' with NSE") + return(NULL) + } + if (!inherits(OutputsModel, "OutputsModel")) { + stop("OutputsModel must be of class 'OutputsModel' \n") + return(NULL) + } + + + ##Initialisation_________________________________ + CritName <- NA + CritVar <- InputsCrit$varObs + if (InputsCrit$transfo == "") { + CritName <- "NSE[CritVar]" + } + if (InputsCrit$transfo == "sqrt") { + CritName <- "NSE[sqrt(CritVar)]" + } + if (InputsCrit$transfo == "log") { + CritName <- "NSE[log(CritVar)]" + } + if (InputsCrit$transfo == "inv") { + CritName <- "NSE[1/CritVar]" + } + if (InputsCrit$transfo == "sort") { + CritName <- "NSE[sort(CritVar)]" + } + CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) + CritValue <- NA + CritBestValue <- +1 + Multiplier <- -1 + ### must be equal to -1 or +1 only + + + ##Data_preparation_______________________________ + VarObs <- InputsCrit$obs + VarObs[!InputsCrit$BoolCrit] <- NA + VarSim <- OutputsModel$Qsim + VarSim[!InputsCrit$BoolCrit] <- NA + + ##Data_transformation + if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { + if (any(VarObs %in% 0)) { + warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + if (any(VarSim %in% 0)) { + warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + } + if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { + VarObs <- VarObs + InputsCrit$epsilon + VarSim <- VarSim + InputsCrit$epsilon + } + if (InputsCrit$transfo == "sqrt") { + VarObs <- sqrt(VarObs) + VarSim <- sqrt(VarSim) + } + if (InputsCrit$transfo == "log") { + VarObs <- log(VarObs) + VarSim <- log(VarSim) + VarSim[VarSim < -1e100] <- NA + } + if (InputsCrit$transfo == "inv") { + VarObs <- 1 / VarObs + VarSim <- 1 / VarSim + VarSim[abs(VarSim) > 1e+100] <- NA + } + if (InputsCrit$transfo == "sort") { + VarSim[is.na(VarObs)] <- NA + VarSim <- sort(VarSim, na.last = TRUE) + VarObs <- sort(VarObs, na.last = TRUE) + InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE) + } + + ##TS_ignore + TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit + Ind_TS_ignore <- which(TS_ignore) + if (length(Ind_TS_ignore) == 0) { + Ind_TS_ignore <- NULL + } + if (sum(!TS_ignore) == 0) { + OutputsCrit <- list(NA) + names(OutputsCrit) <- c("CritValue") return(OutputsCrit) + } + if (inherits(OutputsModel, "hourly")) { + WarningTS <- 365 + } + if (inherits(OutputsModel, "daily")) { + WarningTS <- 365 + + } + if (inherits(OutputsModel, "monthly")) { + WarningTS <- 12 + } + if (inherits(OutputsModel, "yearly")) { + WarningTS <- 3 + } + if (sum(!TS_ignore) < WarningTS & warnings) { + warning("\t criterion computed on less than ", WarningTS, " time-steps") + } + + ##Other_variables_preparation + meanVarObs <- mean(VarObs[!TS_ignore]) + meanVarSim <- mean(VarSim[!TS_ignore]) + + + ##ErrorCrit______________________________________ + Emod <- sum((VarSim[!TS_ignore] - VarObs[!TS_ignore])^2) + Eref <- sum((VarObs[!TS_ignore] - mean(VarObs[!TS_ignore]))^2) + + if (Emod == 0 & Eref == 0) { + Crit <- 0 + } else { + Crit <- (1 - Emod / Eref) + } + if (is.numeric(Crit) & is.finite(Crit)) { + CritValue <- Crit + } + + + ##Verbose______________________________________ + if (verbose) { + message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue), "\n") + } + + + ##Output_________________________________________ + OutputsCrit <- list(CritValue = CritValue, + CritName = CritName, + CritBestValue = CritBestValue, + Multiplier = Multiplier, + Ind_notcomputed = Ind_TS_ignore) + + class(OutputsCrit) <- c("NSE", "ErrorCrit") + return(OutputsCrit) } diff --git a/R/ErrorCrit_RMSE.R b/R/ErrorCrit_RMSE.R index e665f3ee..258e9224 100644 --- a/R/ErrorCrit_RMSE.R +++ b/R/ErrorCrit_RMSE.R @@ -1,144 +1,144 @@ ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) { - - - ##Arguments_check________________________________ - if (!inherits(InputsCrit, "InputsCrit")) { - stop("InputsCrit must be of class 'InputsCrit' \n") - return(NULL) - } - if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { - stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' with RMSE") - return(NULL) - } - if (!inherits(OutputsModel, "OutputsModel")) { - stop("OutputsModel must be of class 'OutputsModel' \n") - return(NULL) - } - - - ##Initialisation_________________________________ - CritName <- NA - CritVar <- InputsCrit$varObs - if (InputsCrit$transfo == "") { - CritName <- "RMSE[CritVar]" - } - if (InputsCrit$transfo == "sqrt") { - CritName <- "RMSE[sqrt(CritVar)]" - } - if (InputsCrit$transfo == "log") { - CritName <- "RMSE[log(CritVar)]" - } - if (InputsCrit$transfo == "inv") { - CritName <- "RMSE[1/CritVar]" - } - if (InputsCrit$transfo == "sort") { - CritName <- "RMSE[sort(CritVar)]" - } - CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) - CritValue <- NA - CritBestValue <- +1 - Multiplier <- +1 - ### must be equal to -1 or +1 only - - - ##Data_preparation_______________________________ - VarObs <- InputsCrit$obs - VarObs[!InputsCrit$BoolCrit] <- NA - VarSim <- OutputsModel$Qsim - VarSim[!InputsCrit$BoolCrit] <- NA - - ##Data_transformation - if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { - if (any(VarObs %in% 0)) { - warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - if (any(VarSim %in% 0)) { - warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") - } - } - if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { - VarObs <- VarObs + InputsCrit$epsilon - VarSim <- VarSim + InputsCrit$epsilon - } - if (InputsCrit$transfo == "sqrt") { - VarObs <- sqrt(VarObs) - VarSim <- sqrt(VarSim) - } - if (InputsCrit$transfo == "log") { - VarObs <- log(VarObs) - VarSim <- log(VarSim) - VarSim[VarSim < -1e100] <- NA - } - if (InputsCrit$transfo == "inv") { - VarObs <- 1 / VarObs - VarSim <- 1 / VarSim - VarSim[abs(VarSim) > 1e+100] <- NA - } - if (InputsCrit$transfo == "sort") { - VarSim[is.na(VarObs)] <- NA - VarSim <- sort(VarSim, na.last = TRUE) - VarObs <- sort(VarObs, na.last = TRUE) - InputsCrit$BoolCrit <- - sort(InputsCrit$BoolCrit, decreasing = TRUE) - } - - ##TS_ignore - TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit - Ind_TS_ignore <- which(TS_ignore) - if (length(Ind_TS_ignore) == 0) { - Ind_TS_ignore <- NULL - } - if (sum(!TS_ignore) == 0) { - OutputsCrit <- list(NA) - names(OutputsCrit) <- c("CritValue") - return(OutputsCrit) - } - if (inherits(OutputsModel, "hourly")) { - WarningTS <- 365 - } - if (inherits(OutputsModel, "daily")) { - WarningTS <- 365 - } - if (inherits(OutputsModel, "monthly")) { - WarningTS <- 12 - } - if (inherits(OutputsModel, "yearly")) { - WarningTS <- 3 - } - if (sum(!TS_ignore) < WarningTS & warnings) { - warning("\t criterion computed on less than ", WarningTS, " time-steps") - } - - - ##ErrorCrit______________________________________ - Numer <- sum((VarSim - VarObs)^2, na.rm = TRUE) - Denom <- sum(!is.na(VarObs)) - - if (Numer == 0) { - Crit <- 0 - } else { - Crit <- sqrt(Numer / Denom) - } - if (is.numeric(Crit) & is.finite(Crit)) { - CritValue <- Crit - } - - - ##Verbose______________________________________ - if (verbose) { - message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue), "\n") - } - - - ##Output_________________________________________ - OutputsCrit <- list(CritValue = CritValue, - CritName = CritName, - CritBestValue = CritBestValue, - Multiplier = Multiplier, - Ind_notcomputed = Ind_TS_ignore - ) - - class(OutputsCrit) <- c("RMSE", "ErrorCrit") + + + ##Arguments_check________________________________ + if (!inherits(InputsCrit, "InputsCrit")) { + stop("InputsCrit must be of class 'InputsCrit' \n") + return(NULL) + } + if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { + stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' with RMSE") + return(NULL) + } + if (!inherits(OutputsModel, "OutputsModel")) { + stop("OutputsModel must be of class 'OutputsModel' \n") + return(NULL) + } + + + ##Initialisation_________________________________ + CritName <- NA + CritVar <- InputsCrit$varObs + if (InputsCrit$transfo == "") { + CritName <- "RMSE[CritVar]" + } + if (InputsCrit$transfo == "sqrt") { + CritName <- "RMSE[sqrt(CritVar)]" + } + if (InputsCrit$transfo == "log") { + CritName <- "RMSE[log(CritVar)]" + } + if (InputsCrit$transfo == "inv") { + CritName <- "RMSE[1/CritVar]" + } + if (InputsCrit$transfo == "sort") { + CritName <- "RMSE[sort(CritVar)]" + } + CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) + CritValue <- NA + CritBestValue <- +1 + Multiplier <- +1 + ### must be equal to -1 or +1 only + + + ##Data_preparation_______________________________ + VarObs <- InputsCrit$obs + VarObs[!InputsCrit$BoolCrit] <- NA + VarSim <- OutputsModel$Qsim + VarSim[!InputsCrit$BoolCrit] <- NA + + ##Data_transformation + if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) { + if (any(VarObs %in% 0)) { + warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + if (any(VarSim %in% 0)) { + warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL") + } + } + if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) { + VarObs <- VarObs + InputsCrit$epsilon + VarSim <- VarSim + InputsCrit$epsilon + } + if (InputsCrit$transfo == "sqrt") { + VarObs <- sqrt(VarObs) + VarSim <- sqrt(VarSim) + } + if (InputsCrit$transfo == "log") { + VarObs <- log(VarObs) + VarSim <- log(VarSim) + VarSim[VarSim < -1e100] <- NA + } + if (InputsCrit$transfo == "inv") { + VarObs <- 1 / VarObs + VarSim <- 1 / VarSim + VarSim[abs(VarSim) > 1e+100] <- NA + } + if (InputsCrit$transfo == "sort") { + VarSim[is.na(VarObs)] <- NA + VarSim <- sort(VarSim, na.last = TRUE) + VarObs <- sort(VarObs, na.last = TRUE) + InputsCrit$BoolCrit <- + sort(InputsCrit$BoolCrit, decreasing = TRUE) + } + + ##TS_ignore + TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit + Ind_TS_ignore <- which(TS_ignore) + if (length(Ind_TS_ignore) == 0) { + Ind_TS_ignore <- NULL + } + if (sum(!TS_ignore) == 0) { + OutputsCrit <- list(NA) + names(OutputsCrit) <- c("CritValue") return(OutputsCrit) - } + if (inherits(OutputsModel, "hourly")) { + WarningTS <- 365 + } + if (inherits(OutputsModel, "daily")) { + WarningTS <- 365 + } + if (inherits(OutputsModel, "monthly")) { + WarningTS <- 12 + } + if (inherits(OutputsModel, "yearly")) { + WarningTS <- 3 + } + if (sum(!TS_ignore) < WarningTS & warnings) { + warning("\t criterion computed on less than ", WarningTS, " time-steps") + } + + + ##ErrorCrit______________________________________ + Numer <- sum((VarSim - VarObs)^2, na.rm = TRUE) + Denom <- sum(!is.na(VarObs)) + + if (Numer == 0) { + Crit <- 0 + } else { + Crit <- sqrt(Numer / Denom) + } + if (is.numeric(Crit) & is.finite(Crit)) { + CritValue <- Crit + } + + + ##Verbose______________________________________ + if (verbose) { + message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue), "\n") + } + + + ##Output_________________________________________ + OutputsCrit <- list(CritValue = CritValue, + CritName = CritName, + CritBestValue = CritBestValue, + Multiplier = Multiplier, + Ind_notcomputed = Ind_TS_ignore + ) + + class(OutputsCrit) <- c("RMSE", "ErrorCrit") + return(OutputsCrit) + +} -- GitLab