diff --git a/DESCRIPTION b/DESCRIPTION index 2f60741d8208c93e91ccf9b6721216b14694eae8..029b5ae56d050ffa073d34b9f3a921e421809b50 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.0.13.2 +Version: 1.0.13.3 Date: 2018-08-29 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 9a0e80b44042f8a674b7faf14644110bd85a20f5..8abc58225e300bbe97535dc53fbe6acf2a207ba2 100644 --- a/NEWS.rmd +++ b/NEWS.rmd @@ -14,7 +14,7 @@ output: -### 1.0.13.2 Release Notes (2018-08-29) +### 1.0.13.3 Release Notes (2018-08-29) #### Deprectated and defunct diff --git a/R/ErrorCrit_NSE.R b/R/ErrorCrit_NSE.R index 742aa202f260d7fac7196cb33e5b032487108ef5..d41e0d2c96ab0bf254b49665bb8e4867000a8eb4 100644 --- a/R/ErrorCrit_NSE.R +++ b/R/ErrorCrit_NSE.R @@ -1,73 +1,133 @@ -ErrorCrit_NSE <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){ - - -##Arguments_check________________________________ - if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); } - if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); } - - -##Initialisation_________________________________ - CritName <- NA; - if(InputsCrit$transfo=="" ){ CritName <- "NSE[Q]" ; } - if(InputsCrit$transfo=="sqrt"){ CritName <- "NSE[sqrt(Q)]"; } - if(InputsCrit$transfo=="log" ){ CritName <- "NSE[log(Q)]" ; } - if(InputsCrit$transfo=="inv" ){ CritName <- "NSE[1/Q]" ; } - if(InputsCrit$transfo=="sort"){ CritName <- "NSE[sort(Q)]"; } - CritValue <- NA; - CritBestValue <- +1; - Multiplier <- -1; ### must be equal to -1 or +1 only - - -##Data_preparation_______________________________ - VarObs <- InputsCrit$Qobs ; VarObs[!InputsCrit$BoolCrit] <- NA; - VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA; - ##Data_transformation - if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){ - 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_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) { -##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)) - } - - -##Output_________________________________________ - OutputsCrit <- list(CritValue = CritValue, CritName = CritName, - CritBestValue = CritBestValue, - Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore) - return(OutputsCrit) - + ##Arguments_check________________________________ + if (inherits(InputsCrit, "InputsCrit") == FALSE) { + stop("InputsCrit must be of class 'InputsCrit' \n") + return(NULL) + } + if (inherits(OutputsModel, "OutputsModel") == FALSE) { + stop("OutputsModel must be of class 'OutputsModel' \n") + return(NULL) + } + + + ##Initialisation_________________________________ + CritName <- NA + if (InputsCrit$transfo == "") { + CritName <- "NSE[Q]" + } + if (InputsCrit$transfo == "sqrt") { + CritName <- "NSE[sqrt(Q)]" + } + if (InputsCrit$transfo == "log") { + CritName <- "NSE[log(Q)]" + } + if (InputsCrit$transfo == "inv") { + CritName <- "NSE[1/Q]" + } + if (InputsCrit$transfo == "sort") { + CritName <- "NSE[sort(Q)]" + } + CritValue <- NA + CritBestValue <- +1 + Multiplier <- -1 + ### must be equal to -1 or +1 only + + + ##Data_preparation_______________________________ + VarObs <- InputsCrit$Qobs + VarObs[!InputsCrit$BoolCrit] <- NA + VarSim <- OutputsModel$Qsim + VarSim[!InputsCrit$BoolCrit] <- NA + + ##Data_transformation + if ("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)) { + if (length(InputsCrit$Ind_zeroes) > 0) { + print(typeof(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)) + } + + + ##Output_________________________________________ + OutputsCrit <- list(CritValue = CritValue, + CritName = CritName, + CritBestValue = CritBestValue, + Multiplier = Multiplier, + Ind_notcomputed = Ind_TS_ignore) + return(OutputsCrit) } - -