Commit 755717fc authored by Dorchies David's avatar Dorchies David
Browse files

feat: add ErrorCrit_GAPX

- Modified CreateInputsCrit for handling parameters

Refs #111
parent 4e33703e
......@@ -15,3 +15,8 @@ Param_Sets_GR4J RunOptions_Val
Param_Sets_GR4J OutputsModel_Val
RunModel_Lag OutputsModelDown
SeriesAggreg SimulatedMonthlyRegime
* InputsCrit$FUN_CRIT
* InputsCritSingle$FUN_CRIT
* InputsCritCompo
* InputsCritMulti
Param_Sets_GR4J InputsCrit_Val$FUN_CRIT
......@@ -24,6 +24,7 @@ S3method(SeriesAggreg, OutputsModel)
export(Calibration)
export(Calibration_Michel)
export(CreateCalibOptions)
export(CreateErrorCrit_GAPX)
export(CreateIniStates)
export(CreateInputsCrit)
export(CreateInputsModel)
......
CreateErrorCrit_GAPX <- function(FUN_TRANSFO) {
FUN_CRIT <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
OutputsModel$ParamT <- FUN_TRANSFO(OutputsModel$Param, "RT")
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "GAPX", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
if (EC$CritCompute) {
ParamApr <- EC$VarObs[!EC$TS_ignore]
ParamOpt <- EC$VarSim[!EC$TS_ignore]
## ErrorCrit
Crit <- 1 - sum(((ParamApr - ParamOpt) / 20)^2)^0.5
if (is.numeric(Crit) & is.finite(Crit)) {
CritValue <- Crit
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("GAPX", "ErrorCrit")
return(OutputsCrit)
}
class(FUN_CRIT) <- c("FUN_CRIT", class(FUN_CRIT))
return(FUN_CRIT)
}
......@@ -48,9 +48,22 @@ CreateInputsCrit <- function(FUN_CRIT,
## check 'Obs' and definition of idLayer
vecObs <- unlist(Obs)
if (length(vecObs) %% LLL != 0 | !is.numeric(vecObs)) {
stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
if (!is.numeric(unlist(Obs))) {
stop("'Obs' must be a (list of) vector(s) of numeric values")
}
Obs2 <- Obs
if ("ParamT" %in% VarObs) {
if (is.list(Obs2)) {
Obs2[[which(VarObs == "ParamT")]] <- NULL
} else {
Obs2 <- NULL
}
}
if (!is.null(Obs2)) {
vecObs <- unlist(Obs2)
if (length(vecObs) %% LLL != 0) {
stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
}
}
if (!is.list(Obs)) {
idLayer <- list(1L)
......@@ -154,7 +167,7 @@ CreateInputsCrit <- function(FUN_CRIT,
listArgs2 <- lapply(seq_along(listArgs$FUN_CRIT), function(i) lapply(listArgs, "[[", i))
## preparation of warning messages
inVarObs <- c("Q", "SCA", "SWE")
inVarObs <- c("Q", "SCA", "SWE", "ParamT")
msgVarObs <- "'VarObs' must be a (list of) character vector(s) and one of %s"
msgVarObs <- sprintf(msgVarObs, paste(sapply(inVarObs, shQuote), collapse = ", "))
inTransfo <- c("", "sqrt", "log", "inv", "sort", "boxcox") # pow is not checked by inTransfo, but appears in the warning message and checkef after (see ## check 'transfo')
......@@ -166,9 +179,11 @@ CreateInputsCrit <- function(FUN_CRIT,
InputsCrit <- lapply(listArgs2, function(iListArgs2) {
## define FUN_CRIT as a character string
iListArgs2$FUN_CRIT <- match.fun(iListArgs2$FUN_CRIT)
## check 'FUN_CRIT'
if (!(identical(iListArgs2$FUN_CRIT, ErrorCrit_NSE ) | identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE ) |
identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE2) | identical(iListArgs2$FUN_CRIT, ErrorCrit_RMSE))) {
if (!all(class(iListArgs2$FUN_CRIT) == c("FUN_CRIT", "function"))) {
stop("incorrect 'FUN_CRIT' for use in 'CreateInputsCrit'", call. = FALSE)
}
if (identical(iListArgs2$FUN_CRIT, ErrorCrit_RMSE) & length(listArgs$Weights) > 1 & all(!is.null(unlist(listArgs$Weights)))) {
......@@ -176,7 +191,14 @@ CreateInputsCrit <- function(FUN_CRIT,
}
## check 'Obs'
if (!is.vector(iListArgs2$Obs) | length(iListArgs2$Obs) != LLL | !is.numeric(iListArgs2$Obs)) {
if (iListArgs2$VarObs == "ParamT") {
# Parameter for regularisation
L2 <- RunOptions$FeatFUN_MOD$NbParam
} else {
# Observation time series
L2 <- LLL
}
if (!is.vector(iListArgs2$Obs) | length(iListArgs2$Obs) != L2 | !is.numeric(iListArgs2$Obs)) {
stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
}
......@@ -187,7 +209,7 @@ CreateInputsCrit <- function(FUN_CRIT,
if (!is.logical(iListArgs2$BoolCrit)) {
stop("'BoolCrit' must be a (list of) vector(s) of boolean", call. = FALSE)
}
if (length(iListArgs2$BoolCrit) != LLL) {
if (length(iListArgs2$BoolCrit) != L2) {
stop("'BoolCrit' and the period defined in 'RunOptions' must have the same length", call. = FALSE)
}
......@@ -283,13 +305,6 @@ CreateInputsCrit <- function(FUN_CRIT,
})
names(InputsCrit) <- paste0("IC", seq_along(InputsCrit))
## define FUN_CRIT as a characater string
listErrorCrit <- c("ErrorCrit_KGE", "ErrorCrit_KGE2", "ErrorCrit_NSE", "ErrorCrit_RMSE")
InputsCrit <- lapply(InputsCrit, function(i) {
i$FUN_CRIT <- listErrorCrit[sapply(listErrorCrit, function(j) identical(i$FUN_CRIT, get(j)))]
i
})
listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
inCnVarObs <- c("SCA", "SWE")
if (!"ZLayers" %in% names(InputsModel)) {
......@@ -314,7 +329,7 @@ CreateInputsCrit <- function(FUN_CRIT,
## define idLayer as an index of the layer to use
for (iInCnVarObs in unique(listVarObs)) {
if (iInCnVarObs == "Q") {
if (!iInCnVarObs %in% inCnVarObs) {
for (i in which(listVarObs == iInCnVarObs)) {
InputsCrit[[i]]$idLayer <- NA
}
......
ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
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
......@@ -36,13 +36,13 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
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 {
......@@ -51,10 +51,10 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
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 {
......@@ -63,20 +63,20 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
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 = " "))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
......@@ -85,8 +85,11 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("KGE", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_KGE) <- c("FUN_CRIT", class(ErrorCrit_KGE))
ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE2", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
SubCritValues <- rep(NA, 3)
SubCritNames <- c("r", "gamma", "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
......@@ -36,10 +36,10 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[1L] <- Crit
}
## SubErrorCrit KGE gamma
SubCritPrint[2L] <- paste0(EC$CritName, " cv(sim)/cv(obs) =")
if (meanVarSim == 0) {
if (sd(EC$VarSim[!EC$TS_ignore]) == 0) {
CVsim <- 1
......@@ -48,7 +48,7 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
}
} else {
CVsim <- sd(EC$VarSim[!EC$TS_ignore]) / meanVarSim
}
if (meanVarObs == 0) {
if (sd(EC$VarObs[!EC$TS_ignore]) == 0) {
......@@ -68,10 +68,10 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
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 {
......@@ -80,20 +80,20 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
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 = " "))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
......@@ -102,8 +102,10 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("KGE2", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_KGE2) <- c("FUN_CRIT", class(ErrorCrit_KGE2))
ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "NSE", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
if (EC$CritCompute) {
## Other variables preparation
meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
## ErrorCrit
Emod <- sum((EC$VarSim[!EC$TS_ignore] - EC$VarObs[!EC$TS_ignore])^2)
Eref <- sum((EC$VarObs[!EC$TS_ignore] - mean(EC$VarObs[!EC$TS_ignore]))^2)
if (Emod == 0 & Eref == 0) {
Crit <- 0
} else {
......@@ -26,22 +26,24 @@ ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
if (is.numeric(Crit) & is.finite(Crit)) {
CritValue <- Crit
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("NSE", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_NSE) <- c("FUN_CRIT", class(ErrorCrit_NSE))
......@@ -4,9 +4,6 @@ ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
if (!inherits(InputsCrit, "Single")) {
stop("'ErrorCrit_RMSE' can only be used with 'InputsCrit' of class 'Single'")
}
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "RMSE", OutputsModel = OutputsModel, warnings = warnings)
......@@ -44,3 +41,5 @@ ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
return(OutputsCrit)
}
class(ErrorCrit_RMSE) <- c("FUN_CRIT", class(ErrorCrit_RMSE))
......@@ -35,7 +35,7 @@
CritBestValue <- +1
Multiplier <- +1
}
if (Crit %in% c("NSE", "KGE", "KGE2")) {
if (Crit %in% c("NSE", "KGE", "KGE2", "GAPX")) {
CritBestValue <- +1
Multiplier <- -1
}
......@@ -44,15 +44,14 @@
## Data preparation
VarObs <- InputsCrit$Obs
VarObs[!InputsCrit$BoolCrit] <- NA
if (InputsCrit$VarObs == "Q") {
VarSim <- OutputsModel$Qsim
}
if (InputsCrit$VarObs == "SCA") {
VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "Gratio"))
}
if (InputsCrit$VarObs == "SWE") {
VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "SnowPack"))
}
VarSim <- switch(
InputsCrit$VarObs,
Q = OutputsModel$Qsim,
SCA = rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "Gratio")),
SWE = rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "SnowPack")),
ParamT = OutputsModel$ParamT
)
VarSim[!InputsCrit$BoolCrit] <- NA
......@@ -111,9 +110,16 @@
} else {
CritCompute <- TRUE
}
WarningTS <- 10
if (sum(!TS_ignore) < WarningTS & warnings) {
warning("\t criterion computed on less than ", WarningTS, " time-steps", call. = FALSE)
if (Crit != "GAPX") {
WarningTS <- 10
if (sum(!TS_ignore) < WarningTS & warnings) {
warning("\t criterion computed on less than ", WarningTS, " time-steps", call. = FALSE)
}
} else {
WarningTS <- 4 # For at least daily time step models (GR4J)
if (sum(!TS_ignore) < WarningTS & warnings) {
warning("\t criterion GAPX computed on less than ", WarningTS, " parameters", call. = FALSE)
}
}
......
context("CreateErrorCrit_GAPX")
test_that("Function should return ErrorCrit function", {
expect_equal(class(CreateErrorCrit_GAPX(TransfoParam_GR1A)), c("FUN_CRIT", "function"))
})
data(L0123001)
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1999-12-31"))
RunOptions <- suppressWarnings(
CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
test_that("ErrorCrit should return 1 for same parameters", {
ErrorCrit_GAPX <- CreateErrorCrit_GAPX(TransfoParam_GR4J)
ParamT <- TransfoParam_GR4J(Param, "RT")
IC <- CreateInputsCrit(ErrorCrit_GAPX, InputsModel, RunOptions, Obs = ParamT, VarObs = "ParamT")
expect_equal(ErrorCrit_GAPX(IC, OutputsModel)$CritValue, 1)
})
context("CreateInputsCrit")
data(L0123001)
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1999-12-31"))
RunOptions <- suppressWarnings(
CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
test_that("KGE crit with log transform should return a warning", {
expect_warning(
CreateInputsCrit(
FUN_CRIT = ErrorCrit_KGE,
InputsModel = InputsModel,
RunOptions = RunOptions,
Obs = BasinObs$Qmm[Ind_Run],
transfo = "log"
),
regex = "we do not advise using the KGE with a log transformation on Obs"
)
})
test_that("Composed crit with two identical should return a warning", {
expect_warning(
CreateInputsCrit(
FUN_CRIT = list(ErrorCrit_KGE, ErrorCrit_KGE),
InputsModel = InputsModel,
RunOptions = RunOptions,
Obs = list(BasinObs$Qmm[Ind_Run], BasinObs$Qmm[Ind_Run]),
VarObs = list("Q", "Q")
),
regex = "the criteria list are identical"
)
})
Markdown is supported
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