Source

Target

Showing with 2357 additions and 2015 deletions
+2357 -2015
This diff is collapsed.
This diff is collapsed.
ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT, warnings = TRUE, verbose = TRUE){
return( FUN_CRIT(InputsCrit,OutputsModel, warnings = warnings, verbose = verbose) )
ErrorCrit <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## ---------- Arguments check
if (!inherits(InputsCrit, "InputsCrit")) {
stop("InputsCrit must be of class 'InputsCrit'")
}
if (!inherits(OutputsModel, "OutputsModel")) {
stop("OutputsModel must be of class 'OutputsModel'")
}
## ---------- Criterion computation
## ----- Single criterion
if (inherits(InputsCrit, "Single")) {
FUN_CRIT <- match.fun(InputsCrit$FUN_CRIT)
OutputsCrit <- FUN_CRIT(InputsCrit = InputsCrit,
OutputsModel = OutputsModel,
warnings = warnings,
verbose = verbose)
}
## ----- Multiple criteria or Composite criterion
if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
listOutputsCrit <- lapply(InputsCrit, FUN = function(iInputsCrit) {
FUN_CRIT <- match.fun(iInputsCrit$FUN_CRIT)
FUN_CRIT(InputsCrit = iInputsCrit,
OutputsModel = 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"]]))
listweights <- listweights / sum(listweights)
if (inherits(InputsCrit, "Compo")) {
CritValue <- sum(listValCrit * 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: sum(", msgForm, ")\n")
}
} else {
OutputsCrit <- listOutputsCrit
class(OutputsCrit) <- c("Multi", "ErrorCrit")
}
}
return(OutputsCrit)
}
ErrorCrit_KGE <- 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 <- "KGE[Q]" ; }
if(InputsCrit$transfo=="sqrt"){ CritName <- "KGE[sqrt(Q)]"; }
if(InputsCrit$transfo=="log" ){ CritName <- "KGE[log(Q)]" ; }
if(InputsCrit$transfo=="inv" ){ CritName <- "KGE[1/Q]" ; }
if(InputsCrit$transfo=="sort"){ CritName <- "KGE[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)
ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
##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){ 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 ) );
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
SubCritValues <- rep(NA, 3)
SubCritNames <- c("r", "alpha", "beta")
SubCritPrint <- rep(NA, 3)
if (EC$CritCompute) {
## Other variables preparation
meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
## SubErrorCrit KGE rPearson
SubCritPrint[1L] <- paste0(EC$CritName, " cor(sim, obs, \"pearson\") =")
Numer <- sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) * (EC$VarSim[!EC$TS_ignore] - meanVarSim))
Deno1 <- sqrt(sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) ^ 2))
Deno2 <- sqrt(sum((EC$VarSim[!EC$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[1L] <- Crit
}
## SubErrorCrit KGE alpha
SubCritPrint[2L] <- paste0(EC$CritName, " sd(sim)/sd(obs) =")
Numer <- sd(EC$VarSim[!EC$TS_ignore])
Denom <- sd(EC$VarObs[!EC$TS_ignore])
if (Numer == 0 & Denom == 0) {
Crit <- 1
} else {
Crit <- Numer / Denom
}
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[2L] <- Crit
}
## SubErrorCrit KGE beta
SubCritPrint[3L] <- paste0(EC$CritName, " mean(sim)/mean(obs) =")
if (meanVarSim == 0 & meanVarObs == 0) {
Crit <- 1
} else {
Crit <- meanVarSim / meanVarObs
}
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[3L] <- Crit
}
## ErrorCrit
if (sum(is.na(SubCritValues)) == 0) {
CritValue <- (1 - sqrt((SubCritValues[1L] - 1)^2 + (SubCritValues[2L] - 1)^2 + (SubCritValues[3L] - 1)^2))
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " "))
}
}
##Verbose______________________________________
if(verbose) {
message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue))
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " "))
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
SubCritValues = SubCritValues,
SubCritNames = SubCritNames,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
##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)
}
class(ErrorCrit_KGE) <- c("FUN_CRIT", class(ErrorCrit_KGE))
This diff is collapsed.
This diff is collapsed.
ErrorCrit_RMSE <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){
ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
##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); }
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "RMSE", OutputsModel = OutputsModel, warnings = warnings)
##Initialisation_________________________________
CritName <- NA;
if(InputsCrit$transfo=="" ){ CritName <- "RMSE[Q]" ; }
if(InputsCrit$transfo=="sqrt"){ CritName <- "RMSE[sqrt(Q)]"; }
if(InputsCrit$transfo=="log" ){ CritName <- "RMSE[log(Q)]" ; }
if(InputsCrit$transfo=="inv" ){ CritName <- "RMSE[1/Q]" ; }
if(InputsCrit$transfo=="sort"){ CritName <- "RMSE[sort(Q)]"; }
CritValue <- NA
CritValue <- NA;
CritBestValue <- +1;
Multiplier <- +1; ### must be equal to -1 or +1 only
if (EC$CritCompute) {
## ErrorCrit
Numer <- sum((EC$VarSim - EC$VarObs)^2, na.rm = TRUE)
Denom <- sum(!is.na(EC$VarObs))
if (Numer == 0) {
Crit <- 0
} else {
Crit <- sqrt(Numer / Denom)
}
if (is.numeric(Crit) & is.finite(Crit)) {
CritValue <- Crit
}
##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)
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
}
}
##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))
}
##Output_________________________________________
OutputsCrit <- list(CritValue = CritValue, CritName = CritName, CritBestValue = CritBestValue,
Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore)
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("RMSE", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_RMSE) <- c("FUN_CRIT", class(ErrorCrit_RMSE))
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.