Commit 52e9d535 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

vv1.0.13.3 ErrorCrit_NSE cleaned

Showing with 132 additions and 72 deletions
+132 -72
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.0.13.2 Version: 1.0.13.3
Date: 2018-08-29 Date: 2018-08-29
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")),
......
...@@ -14,7 +14,7 @@ output: ...@@ -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 #### Deprectated and defunct
......
ErrorCrit_NSE <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){ 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______________________________________ ##Arguments_check________________________________
Emod <- sum((VarSim[!TS_ignore]-VarObs[!TS_ignore])^2); if (inherits(InputsCrit, "InputsCrit") == FALSE) {
Eref <- sum((VarObs[!TS_ignore]-mean(VarObs[!TS_ignore]))^2); stop("InputsCrit must be of class 'InputsCrit' \n")
if(Emod==0 & Eref==0){ Crit <- 0; } else { Crit <- (1-Emod/Eref); } return(NULL)
if(is.numeric(Crit) & is.finite(Crit)){ CritValue <- Crit; } }
if (inherits(OutputsModel, "OutputsModel") == FALSE) {
stop("OutputsModel must be of class 'OutputsModel' \n")
##Verbose______________________________________ return(NULL)
if(verbose) { }
message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue))
}
##Initialisation_________________________________
CritName <- NA
##Output_________________________________________ if (InputsCrit$transfo == "") {
OutputsCrit <- list(CritValue = CritValue, CritName = CritName, CritName <- "NSE[Q]"
CritBestValue = CritBestValue, }
Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore) if (InputsCrit$transfo == "sqrt") {
return(OutputsCrit) 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)
} }
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