diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 938d200bee886d9916de746355d2fb2788b32956..a9474e78e38c6722e4a6b86eaedc66d477ee9ecc 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -35,11 +35,11 @@ default:
extends: .check
scheduled_tests_patched:
- only:
- refs:
- - dev
- - master
- - schedules
+ # only:
+ # refs:
+ # - dev
+ # - master
+ # - schedules
variables:
R_VERSION: "patched"
extends: .scheduled_tests
@@ -108,10 +108,10 @@ check_as_cran_oldrel:
revdepcheck_patched:
stage: revdepcheck
- only:
- refs:
- - tags
- - schedules
+ # only:
+ # refs:
+ # - tags
+ # - schedules
variables:
R_VERSION: "patched"
script:
diff --git a/.regressionignore b/.regressionignore
index f7eeb15b1d582fbb8eb0c291683b03a4467b6029..44aced51723929d2b68509f61dd2e192ef30f2d9 100644
--- a/.regressionignore
+++ b/.regressionignore
@@ -1,65 +1,26 @@
# .test-regression.ignore contains the list of topic/variables produces by
# documentation examples that should be ignore in the regression test
# The format of this file is: 5 lines of comments followed by one line by
-# ignored variable : [Topic]<SPACE>[Variable].
+# ignored variable : [Topic]<SPACE>[Variable] or *<SPACE>[Variable] for every variable whatever the topic
# Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel
-Calibration_Michel RunOptions
-Calibration RunOptions
-CreateCalibOptions RunOptions
-CreateIniStates RunOptions
-CreateInputsCrit RunOptions
-CreateInputsModel RunOptions
-CreateRunOptions RunOptions
-ErrorCrit_KGE RunOptions
-ErrorCrit_KGE2 RunOptions
-ErrorCrit_NSE RunOptions
-ErrorCrit_RMSE RunOptions
-ErrorCrit RunOptions
-Imax RunOptions
+* RunOptions$Outputs_Sim
+* RunOptions$Outputs_Cal
+* RunOptions$FeatFUN_MOD
+* RunOptions$FortranOutputs
+RunModel_CemaNeige RunOptions$Outputs_Cal
Param_Sets_GR4J RunOptions_Cal
Param_Sets_GR4J RunOptions_Val
-RunModel_CemaNeige RunOptions
-RunModel_CemaNeigeGR4J RunOptions
-RunModel_CemaNeigeGR5J RunOptions
-RunModel_CemaNeigeGR6J RunOptions
-RunModel_GR1A RunOptions
-RunModel_GR2M RunOptions
-RunModel_GR4H RunOptions
-RunModel_GR4J RunOptions
-RunModel_GR5H RunOptions
-RunModel_GR5J RunOptions
-RunModel_GR6J RunOptions
-RunModel_Lag RunOptions
-RunModel RunOptions
-SeriesAggreg RunOptions
-Calibration OutputsModel
-Calibration_Michel OutputsModel
-CreateCalibOptions OutputsModel
-CreateIniStates OutputsModel
-CreateInputsCrit OutputsModel
-CreateInputsModel OutputsModel
-CreateRunOptions OutputsModel
-ErrorCrit OutputsModel
-ErrorCrit_KGE OutputsModel
-ErrorCrit_KGE2 OutputsModel
-ErrorCrit_NSE OutputsModel
-ErrorCrit_RMSE OutputsModel
-Imax OutputsModel
-RunModel OutputsModel
-RunModel_CemaNeige OutputsModel
-RunModel_CemaNeigeGR4J OutputsModel
-RunModel_CemaNeigeGR5J OutputsModel
-RunModel_CemaNeigeGR6J OutputsModel
-RunModel_GR1A OutputsModel
-RunModel_GR2M OutputsModel
-RunModel_GR4H OutputsModel
-RunModel_GR4J OutputsModel
-RunModel_GR5H OutputsModel
-RunModel_GR5J OutputsModel
-RunModel_GR6J OutputsModel
-RunModel_Lag OutputsModel
-SeriesAggreg OutputsModel
+* OutputsModel$Param
+* OutputsModel$WarmUpQsim
+* OutputsModel$StateEnd
Param_Sets_GR4J OutputsModel_Val
+RunModel_Lag InputsModel
+RunModel_Lag OutputsModel
RunModel_Lag OutputsModelDown
SeriesAggreg SimulatedMonthlyRegime
+* InputsCrit$FUN_CRIT
+* InputsCritSingle$FUN_CRIT
+* InputsCritCompo
+* InputsCritMulti
+Param_Sets_GR4J InputsCrit_Val$FUN_CRIT
diff --git a/DESCRIPTION b/DESCRIPTION
index a6074eef6714be39fc082d8ad56b648b53dae38d..2d31e79c54e9669fbb9557f1fbaa7cc25b679f06 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -38,3 +38,4 @@ BugReports: https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues
NeedsCompilation: yes
Encoding: UTF-8
VignetteBuilder: knitr
+RoxygenNote: 7.1.1
diff --git a/NAMESPACE b/NAMESPACE
index 9cf8783bf95a9607e34baa47f3b10cf0602cbd38..138db937ab679c768065c17552e32dd0effebf62 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,8 +24,10 @@ S3method(SeriesAggreg, OutputsModel)
export(Calibration)
export(Calibration_Michel)
export(CreateCalibOptions)
+export(CreateErrorCrit_GAPX)
export(CreateIniStates)
export(CreateInputsCrit)
+export(CreateInputsCrit_Lavenne)
export(CreateInputsModel)
export(CreateRunOptions)
export(DataAltiExtrapolation_Valery)
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 1ad475be2811e1dd7e47042362adabe31044c103..462cd26185a8d7b89f69269b36e8a40389d0359e 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -404,7 +404,7 @@ Calibration_Michel <- function(InputsModel,
listNameCrit <- OutputsCrit$CritCompo$MultiCritNames
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)]
+ 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, ")")
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index f6d5ffb60c6ebb770d9fb0a3b0439d9df2885ad4..400b4d4d7dbcfddcb33ab68f19421a78ebb43b29 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -21,67 +21,15 @@ CreateCalibOptions <- function(FUN_MOD,
if (!is.logical(IsSD) | length(IsSD) != 1L) {
stop("'IsSD' must be a logical of length 1")
}
+
## check FUN_MOD
- BOOL <- FALSE
+ FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD)
+ FeatFUN_MOD$IsHyst <- IsHyst
+ FeatFUN_MOD$IsSD <- IsSD
+ ObjectClass <- FeatFUN_MOD$Class
- if (identical(FUN_MOD, RunModel_GR4H)) {
- ObjectClass <- c(ObjectClass, "GR4H")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR5H)) {
- ObjectClass <- c(ObjectClass, "GR5H")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR4J)) {
- ObjectClass <- c(ObjectClass, "GR4J")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR5J)) {
- ObjectClass <- c(ObjectClass, "GR5J")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR6J)) {
- ObjectClass <- c(ObjectClass, "GR6J")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR2M)) {
- ObjectClass <- c(ObjectClass, "GR2M")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR1A)) {
- ObjectClass <- c(ObjectClass, "GR1A")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeige)) {
- ObjectClass <- c(ObjectClass, "CemaNeige")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
- ObjectClass <- c(ObjectClass, "CemaNeigeGR4H")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
- ObjectClass <- c(ObjectClass, "CemaNeigeGR5H")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
- ObjectClass <- c(ObjectClass, "CemaNeigeGR4J")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
- ObjectClass <- c(ObjectClass, "CemaNeigeGR5J")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- ObjectClass <- c(ObjectClass, "CemaNeigeGR6J")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_Lag)) {
- ObjectClass <- c(ObjectClass, "Lag")
- if (IsSD) {
+ if (identical(FUN_MOD, RunModel_Lag) && IsSD) {
stop("RunModel_Lag should not be used with 'isSD=TRUE'")
- }
- BOOL <- TRUE
}
if (IsHyst) {
ObjectClass <- c(ObjectClass, "hysteresis")
@@ -89,10 +37,6 @@ CreateCalibOptions <- function(FUN_MOD,
if (IsSD) {
ObjectClass <- c(ObjectClass, "SD")
}
- if (!BOOL) {
- stop("incorrect 'FUN_MOD' for use in 'CreateCalibOptions'")
- return(NULL)
- }
## check FUN_CALIB
BOOL <- FALSE
@@ -109,202 +53,11 @@ CreateCalibOptions <- function(FUN_MOD,
## check FUN_TRANSFO
if (is.null(FUN_TRANSFO)) {
- ## set FUN1
- if (identical(FUN_MOD, RunModel_GR4H) |
- identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
- FUN_GR <- TransfoParam_GR4H
- }
- if (identical(FUN_MOD, RunModel_GR5H) |
- identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
- FUN_GR <- TransfoParam_GR5H
- }
- if (identical(FUN_MOD, RunModel_GR4J) |
- identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
- FUN_GR <- TransfoParam_GR4J
- }
- if (identical(FUN_MOD, RunModel_GR5J) |
- identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
- FUN_GR <- TransfoParam_GR5J
- }
- if (identical(FUN_MOD, RunModel_GR6J) |
- identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- FUN_GR <- TransfoParam_GR6J
- }
- if (identical(FUN_MOD, RunModel_GR2M)) {
- FUN_GR <- TransfoParam_GR2M
- }
- if (identical(FUN_MOD, RunModel_GR1A)) {
- FUN_GR <- TransfoParam_GR1A
- }
- if (identical(FUN_MOD, RunModel_CemaNeige)) {
- if (IsHyst) {
- FUN_GR <- TransfoParam_CemaNeigeHyst
- } else {
- FUN_GR <- TransfoParam_CemaNeige
- }
- }
- if (identical(FUN_MOD, RunModel_Lag)) {
- FUN_GR <- TransfoParam_Lag
- }
- if (is.null(FUN_GR)) {
- stop("'FUN_GR' was not found")
- return(NULL)
- }
- ## set FUN2
- if (IsHyst) {
- FUN_SNOW <- TransfoParam_CemaNeigeHyst
- } else {
- FUN_SNOW <- TransfoParam_CemaNeige
- }
- ## set FUN_LAG
- if (IsSD) {
- FUN_LAG <- TransfoParam_Lag
- }
- ## set FUN_TRANSFO
- if (sum(ObjectClass %in% c("GR4H", "GR5H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige", "Lag")) > 0) {
- if (!IsSD) {
- FUN_TRANSFO <- FUN_GR
- } else {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 2:NParam] <- FUN_GR(ParamIn[, 2:NParam], Direction)
- ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- } else {
- if (IsHyst & !IsSD) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam - 4) ] <- FUN_GR(ParamIn[, 1:(NParam - 4)], Direction)
- ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- if (!IsHyst & !IsSD) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- if (NParam <= 3) {
- ParamOut[, 1:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
- } else {
- ParamOut[, 1:(NParam - 2)] <- FUN_GR(ParamIn[, 1:(NParam - 2)], Direction)
- }
- ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- if (IsHyst & IsSD) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 2:(NParam - 4) ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
- ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
- ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- if (!IsHyst & IsSD) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- if (NParam <= 3) {
- ParamOut[, 2:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
- } else {
- ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)], Direction)
- }
- ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
- ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- }
- }
- if (is.null(FUN_TRANSFO)) {
- stop("'FUN_TRANSFO' was not found")
- return(NULL)
+ FUN_TRANSFO <- .FunTransfo(FeatFUN_MOD)
}
## NParam
- if ("GR4H" %in% ObjectClass) {
- NParam <- 4
- }
- if ("GR5H" %in% ObjectClass) {
- NParam <- 5
- }
- if ("GR4J" %in% ObjectClass) {
- NParam <- 4
- }
- if ("GR5J" %in% ObjectClass) {
- NParam <- 5
- }
- if ("GR6J" %in% ObjectClass) {
- NParam <- 6
- }
- if ("GR2M" %in% ObjectClass) {
- NParam <- 2
- }
- if ("GR1A" %in% ObjectClass) {
- NParam <- 1
- }
- if ("CemaNeige" %in% ObjectClass) {
- NParam <- 2
- }
- if ("CemaNeigeGR4H" %in% ObjectClass) {
- NParam <- 6
- }
- if ("CemaNeigeGR5H" %in% ObjectClass) {
- NParam <- 7
- }
- if ("CemaNeigeGR4J" %in% ObjectClass) {
- NParam <- 6
- }
- if ("CemaNeigeGR5J" %in% ObjectClass) {
- NParam <- 7
- }
- if ("CemaNeigeGR6J" %in% ObjectClass) {
- NParam <- 8
- }
- if ("Lag" %in% ObjectClass) {
- NParam <- 1
- }
+ NParam <- FeatFUN_MOD$NbParam
if (IsHyst) {
NParam <- NParam + 2
@@ -357,80 +110,80 @@ CreateCalibOptions <- function(FUN_MOD,
## check StartParamList and StartParamDistrib default values
if (("HBAN" %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib))) {
- if ("GR4H" %in% ObjectClass) {
+ if ("GR4H" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69,
+5.58, -0.85, +4.74, -9.47,
+6.01, -0.50, +5.14, -8.87), ncol = 4, byrow = TRUE)
}
- if (("GR5H" %in% ObjectClass) & ("interception" %in% ObjectClass)) {
+ if (("GR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34,
+3.74, -0.41, +4.78, -8.94, -3.33,
+4.29, +0.16, +5.39, -7.39, +3.33), ncol = 5, byrow = TRUE)
}
- if (("GR5H" %in% ObjectClass) & !("interception" %in% ObjectClass)) {
+ if (("GR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49,
+3.62, -0.19, +4.80, -9.00, -6.31,
+4.01, -0.04, +5.43, -7.53, -5.33), ncol = 5, byrow = TRUE)
}
- if ("GR4J" %in% ObjectClass) {
+ if ("GR4J" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05,
+5.51, -0.61, +3.74, -8.51,
+6.07, -0.02, +4.42, -8.06), ncol = 4, byrow = TRUE)
}
- if ("GR5J" %in% ObjectClass) {
+ if ("GR5J" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45,
+5.55, -0.46, +3.75, -9.09, -4.69,
+6.10, -0.11, +4.43, -8.60, -0.66), ncol = 5, byrow = TRUE)
}
- if ("GR6J" %in% ObjectClass) {
+ if ("GR6J" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
+3.90, -0.50, +4.10, -8.70, +0.10, +4.00,
+4.50, +0.50, +5.00, -8.10, +1.10, +5.00), ncol = 6, byrow = TRUE)
}
- if ("GR2M" %in% ObjectClass) {
+ if ("GR2M" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+5.03, -7.15,
+5.22, -6.74,
+5.85, -6.37), ncol = 2, byrow = TRUE)
}
- if ("GR1A" %in% ObjectClass) {
+ if ("GR1A" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(-1.69,
-0.38,
+1.39), ncol = 1, byrow = TRUE)
}
- if ("CemaNeige" %in% ObjectClass) {
+ if ("CemaNeige" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(-9.96, +6.63,
-9.14, +6.90,
+4.10, +7.21), ncol = 2, byrow = TRUE)
}
- if ("CemaNeigeGR4H" %in% ObjectClass) {
+ if ("CemaNeigeGR4H" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69, -9.96, +6.63,
+5.58, -0.85, +4.74, -9.47, -9.14, +6.90,
+6.01, -0.50, +5.14, -8.87, +4.10, +7.21), ncol = 6, byrow = TRUE)
}
- if (("CemaNeigeGR5H" %in% ObjectClass) & ("interception" %in% ObjectClass)) {
+ if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34, -9.96, +6.63,
+3.74, -0.41, +4.78, -8.94, -3.33, -9.14, +6.90,
+4.29, +0.16, +5.39, -7.39, +3.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
}
- if (("CemaNeigeGR5H" %in% ObjectClass) & !("interception" %in% ObjectClass)) {
+ if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49, -9.96, +6.63,
+3.62, -0.19, +4.80, -9.00, -6.31, -9.14, +6.90,
+4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
}
- if ("CemaNeigeGR4J" %in% ObjectClass) {
+ if ("CemaNeigeGR4J" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
+5.51, -0.61, +3.74, -8.51, -9.14, +6.90,
+6.07, -0.02, +4.42, -8.06, +4.10, +7.21), ncol = 6, byrow = TRUE)
}
- if ("CemaNeigeGR5J" %in% ObjectClass) {
+ if ("CemaNeigeGR5J" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45, -9.96, +6.63,
+5.55, -0.46, +3.75, -9.09, -4.69, -9.14, +6.90,
+6.10, -0.11, +4.43, -8.60, -0.66, +4.10, +7.21), ncol = 7, byrow = TRUE)
}
- if ("CemaNeigeGR6J" %in% ObjectClass) {
+ if ("CemaNeigeGR6J" == FeatFUN_MOD$CodeMod) {
ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -9.96, +6.63,
+3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -9.14, +6.90,
+4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = 8, byrow = TRUE)
diff --git a/R/CreateErrorCrit_GAPX.R b/R/CreateErrorCrit_GAPX.R
new file mode 100644
index 0000000000000000000000000000000000000000..33b424fa2c0634e6008488c8738b130a1f273619
--- /dev/null
+++ b/R/CreateErrorCrit_GAPX.R
@@ -0,0 +1,47 @@
+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)
+}
diff --git a/R/CreateIniStates.R b/R/CreateIniStates.R
index 9fe48b137b3a9abbda09836d93f5907f1ec1fd39..84e6b53cfe56717803c9439f7f566750a21e1b95 100644
--- a/R/CreateIniStates.R
+++ b/R/CreateIniStates.R
@@ -12,45 +12,11 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
UH1n <- 20L
UH2n <- UH1n * 2L
- nameFUN_MOD <- as.character(substitute(FUN_MOD))
FUN_MOD <- match.fun(FUN_MOD)
- ## check FUN_MOD
- BOOL <- FALSE
- if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H)) {
- ObjectClass <- c(ObjectClass, "GR", "hourly")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR4J) |
- identical(FUN_MOD, RunModel_GR5J) |
- identical(FUN_MOD, RunModel_GR6J)) {
- ObjectClass <- c(ObjectClass, "GR", "daily")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR2M)) {
- ObjectClass <- c(ObjectClass, "GR", "monthly")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR1A)) {
- stop("'RunModel_GR1A' does not require 'IniStates' object")
- }
- if (identical(FUN_MOD, RunModel_CemaNeige)) {
- ObjectClass <- c(ObjectClass, "CemaNeige", "daily")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
- ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "hourly")
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
- identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
- identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "daily")
- BOOL <- TRUE
- }
- if (!BOOL) {
- stop("incorrect 'FUN_MOD' for use in 'CreateIniStates'")
- }
+ FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, DatesR = InputsModel$DatesR)
+ ObjectClass <- FeatFUN_MOD$Class
+
if (!"CemaNeige" %in% ObjectClass & IsHyst) {
stop("'IsHyst' cannot be TRUE if CemaNeige is not used in 'FUN_MOD'")
}
@@ -82,7 +48,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
}
} else if (!is.null(ExpStore)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'ExpStore'. Value set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'ExpStore'. Value set to NA", FeatFUN_MOD$NameFunMod))
}
ExpStore <- Inf
}
@@ -90,13 +56,13 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if (identical(FUN_MOD, RunModel_GR2M)) {
if (!is.null(UH1)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
UH1 <- rep(Inf, UH1n)
}
if (!is.null(UH2)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'UH2'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
UH2 <- rep(Inf, UH2n)
}
@@ -104,13 +70,13 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !is.null(UH1)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
UH1 <- rep(Inf, UH1n)
}
if ((!identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & !is.null(IntStore)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
IntStore <- Inf
}
@@ -118,57 +84,57 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if ("CemaNeige" %in% ObjectClass & ! "GR" %in% ObjectClass) {
if (!is.null(ProdStore)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'ProdStore'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'ProdStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
}
ProdStore <- Inf
if (!is.null(RoutStore)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'RoutStore'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'RoutStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
}
RoutStore <- Inf
if (!is.null(ExpStore)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'ExpStore'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'ExpStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
}
ExpStore <- Inf
if (!is.null(IntStore)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
}
IntStore <- Inf
if (!is.null(UH1)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
}
UH1 <- rep(Inf, UH1n)
if (!is.null(UH2)) {
if (verbose) {
- warning(sprintf("'%s' does not require 'UH2'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
}
UH2 <- rep(Inf, UH2n)
}
if (IsIntStore & is.null(IntStore)) {
- stop(sprintf("'%s' need values for 'IntStore'", nameFUN_MOD))
+ stop(sprintf("'%s' need values for 'IntStore'", FeatFUN_MOD$NameFunMod))
}
if ("CemaNeige" %in% ObjectClass & !IsHyst &
(is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers))) {
- stop(sprintf("'%s' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'", nameFUN_MOD))
+ stop(sprintf("'%s' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'", FeatFUN_MOD$NameFunMod))
}
if ("CemaNeige" %in% ObjectClass & IsHyst &
(is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers) |
is.null(GthrCemaNeigeLayers) | is.null(GlocmaxCemaNeigeLayers))) {
- stop(sprintf("'%s' need values for 'GCemaNeigeLayers', 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'", nameFUN_MOD))
+ stop(sprintf("'%s' need values for 'GCemaNeigeLayers', 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'", FeatFUN_MOD$NameFunMod))
}
if ("CemaNeige" %in% ObjectClass & !IsHyst &
(!is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
if (verbose) {
- warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
GthrCemaNeigeLayers <- Inf
GlocmaxCemaNeigeLayers <- Inf
@@ -176,7 +142,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if (!"CemaNeige" %in% ObjectClass &
(!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers) | !is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
if (verbose) {
- warning(sprintf("'%s' does not require 'GCemaNeigeLayers' 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
+ warning(sprintf("'%s' does not require 'GCemaNeigeLayers' 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
GCemaNeigeLayers <- Inf
eTGCemaNeigeLayers <- Inf
diff --git a/R/CreateInputsCrit.R b/R/CreateInputsCrit.R
index 4107b328b4f37834caacd541adbc1adf2829f5b4..049e89a037ab60bcec6fc11b082bd324459e8143 100644
--- a/R/CreateInputsCrit.R
+++ b/R/CreateInputsCrit.R
@@ -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,8 +191,15 @@ CreateInputsCrit <- function(FUN_CRIT,
}
## check 'Obs'
- if (!is.vector(iListArgs2$Obs) | length(iListArgs2$Obs) != LLL | !is.numeric(iListArgs2$Obs)) {
- stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
+ 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", L2), call. = FALSE)
}
## check 'BoolCrit'
@@ -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
}
diff --git a/R/CreateInputsCrit_Lavenne.R b/R/CreateInputsCrit_Lavenne.R
new file mode 100644
index 0000000000000000000000000000000000000000..f9bed6d1ac15d31ed5eb627821663e9dfc6c2a0a
--- /dev/null
+++ b/R/CreateInputsCrit_Lavenne.R
@@ -0,0 +1,52 @@
+CreateInputsCrit_Lavenne <- function(FUN_CRIT = ErrorCrit_KGE,
+ InputsModel,
+ RunOptions,
+ Obs,
+ VarObs = "Q",
+ AprParamR,
+ AprCrit = 1,
+ k = 0.15,
+ BoolCrit = NULL,
+ transfo = "sqrt",
+ epsilon = NULL) {
+
+ # Check parameters
+ if (!is.numeric(AprCrit) || length(AprCrit) != 1 || AprCrit > 1) {
+ stop("'AprCrit' must be a numeric of length 1 with a maximum value of 1")
+ }
+ if (!is.numeric(k) || length(k) != 1 || k < 0 || k > 1) {
+ stop("'k' must be a numeric of length 1 with a value between 0 and 1")
+ }
+ if (!is.null(BoolCrit) && !is.logical(BoolCrit)) {
+ stop("'BoolCrit must be logical")
+ }
+ if (!is.character(transfo)) {
+ stop("'transfo' must be character")
+ }
+ if (!is.null(epsilon) && !is.numeric(epsilon)) {
+ stop("'epsilon' must be numeric")
+ }
+ if (!is.numeric(AprParamR) || length(AprParamR) != RunOptions$FeatFUN_MOD$NbParam) {
+ stop("'AprParamR' must be a numeric vector of length ",
+ RunOptions$FeatFUN_MOD$NbParam)
+ }
+
+
+ FUN_TRANSFO <- .FunTransfo(RunOptions$FeatFUN_MOD)
+
+ AprParamT <- FUN_TRANSFO(AprParamR, "RT")
+
+ ErrorCrit_GAPX <- CreateErrorCrit_GAPX(FUN_TRANSFO)
+
+
+
+ CreateInputsCrit(FUN_CRIT = list(FUN_CRIT, ErrorCrit_GAPX),
+ InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Obs = list(Obs, AprParamT),
+ VarObs = c("Q", "ParamT"),
+ Weights = c(1 - k, k * max(0, AprCrit)),
+ BoolCrit = list(BoolCrit, NULL),
+ transfo = list(transfo, ""),
+ epsilon = list(epsilon, NULL))
+}
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index 8821935b9b971faf543f3a5f47735376a5a1edb2..04fe0c818e49977fbbef106844ece8e996d1c1e4 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -32,9 +32,18 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
if (IsIntStore) {
ObjectClass <- c(ObjectClass, "interception")
}
- if (IsHyst) {
- ObjectClass <- c(ObjectClass, "hysteresis")
- FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 2
+ if ("CemaNeige" %in% FeatFUN_MOD$Class) {
+ FeatFUN_MOD$IsHyst <- IsHyst
+ if (IsHyst) {
+ ObjectClass <- c(ObjectClass, "hysteresis")
+ FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 2
+ }
+ }
+
+ ## SD model
+ FeatFUN_MOD$IsSD <- inherits(InputsModel, "SD")
+ if (FeatFUN_MOD$IsSD) {
+ FeatFUN_MOD$NbParam <- FeatFUN_MOD$NbParam + 1
}
if (!"CemaNeige" %in% ObjectClass & "hysteresis" %in% ObjectClass) {
@@ -295,7 +304,10 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
##check_Outputs_Cal_and_Sim
##Outputs_all
- Outputs_all <- c("DatesR", unlist(FortranOutputs), "WarmUpQsim", "StateEnd")
+ Outputs_all <- c("DatesR", unlist(FortranOutputs), "WarmUpQsim", "StateEnd", "Param")
+ if (FeatFUN_MOD$IsSD) {
+ Outputs_all <- c(Outputs_all, "QsimDown", "Qsim_m3")
+ }
##check_Outputs_Sim
if (!is.vector(Outputs_Sim)) {
@@ -321,15 +333,14 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
##check_Outputs_Cal
if (is.null(Outputs_Cal)) {
if ("GR" %in% ObjectClass) {
- Outputs_Cal <- c("Qsim")
+ Outputs_Cal <- c("Qsim", "Param")
+ if ("CemaNeige" %in% ObjectClass) {
+ Outputs_Cal <- c("PliqAndMelt", Outputs_Cal)
+ }
}
if ("CemaNeige" %in% ObjectClass) {
Outputs_Cal <- c("all")
}
- if ("GR" %in% ObjectClass &
- "CemaNeige" %in% ObjectClass) {
- Outputs_Cal <- c("PliqAndMelt", "Qsim")
- }
} else {
if (!is.vector(Outputs_Cal)) {
stop("'Outputs_Cal' must be a vector of characters")
diff --git a/R/ErrorCrit_KGE.R b/R/ErrorCrit_KGE.R
index 008e64e9f2298cb70d612a6e099b1d385c1d62b3..cd5d53f73d25212433fdbede8542afebb60a2de9 100644
--- a/R/ErrorCrit_KGE.R
+++ b/R/ErrorCrit_KGE.R
@@ -1,29 +1,29 @@
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
@@ -31,18 +31,18 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
Crit <- 0
}
} else {
- Crit <- Numer / (Deno1 * Deno2)
+ 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 {
@@ -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))
+
diff --git a/R/ErrorCrit_KGE2.R b/R/ErrorCrit_KGE2.R
index a9313d90a5e4fe17ce020e4811b1aae1a536a8d0..bba5632391f21b4efe406369a5f68f6569f9034d 100644
--- a/R/ErrorCrit_KGE2.R
+++ b/R/ErrorCrit_KGE2.R
@@ -1,29 +1,29 @@
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)
-
+
+ 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
@@ -31,15 +31,15 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
Crit <- 0
}
} else {
- Crit <- Numer / (Deno1 * Deno2)
+ Crit <- Numer / (Deno1 * Deno2)
}
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))
diff --git a/R/ErrorCrit_NSE.R b/R/ErrorCrit_NSE.R
index 4a7a30587fcf770493b17967c0a95403e11532d8..0f66c2d3c71ba7094f256762e4dc54e90f83464b 100644
--- a/R/ErrorCrit_NSE.R
+++ b/R/ErrorCrit_NSE.R
@@ -1,23 +1,23 @@
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)
-
+
+ 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))
diff --git a/R/ErrorCrit_RMSE.R b/R/ErrorCrit_RMSE.R
index 6cab33bc57147355c0603f6cfd86daf8904c1bee..845ff0a50f3bcbc6e8f75b8ec5846a07ef91df45 100644
--- a/R/ErrorCrit_RMSE.R
+++ b/R/ErrorCrit_RMSE.R
@@ -1,19 +1,19 @@
ErrorCrit_RMSE <- 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 = "RMSE", OutputsModel = OutputsModel, warnings = warnings)
-
+
+ EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "RMSE", OutputsModel = OutputsModel, warnings = warnings)
+
CritValue <- NA
-
+
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 {
@@ -22,22 +22,24 @@ ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
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("RMSE", "ErrorCrit")
return(OutputsCrit)
-
+
}
+
+class(ErrorCrit_RMSE) <- c("FUN_CRIT", class(ErrorCrit_RMSE))
diff --git a/R/RunModel.R b/R/RunModel.R
index c21b803cfb83b02de22136b46a0ad95f38d94424..d68bf066453a2fdad3d907ff366b6a63e675a737 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -5,6 +5,7 @@ RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, ...) {
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
# Lag model take one parameter at the beginning of the vector
iFirstParamRunOffModel <- 2
+ RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - 1
} else {
# All parameters
iFirstParamRunOffModel <- 1
diff --git a/R/RunModel_CemaNeigeGR4H.R b/R/RunModel_CemaNeigeGR4H.R
index 51318c1faa81aea3bfad90c75b2eb1771b65a6b3..7bba7814f6ebac4e39ef288900fdddf27a8e8c1d 100644
--- a/R/RunModel_CemaNeigeGR4H.R
+++ b/R/RunModel_CemaNeigeGR4H.R
@@ -167,5 +167,6 @@ RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
RunOptions,
RESULTS,
LInputSeries,
+ Param,
CemaNeigeLayers)
}
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index 34f15d5ea3c3a546e30a5de3b857c4484e315194..28eceb59e732aba01819e2e7f58b4f6416711c44 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -166,5 +166,6 @@ RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
RunOptions,
RESULTS,
LInputSeries,
+ Param,
CemaNeigeLayers)
}
diff --git a/R/RunModel_CemaNeigeGR5H.R b/R/RunModel_CemaNeigeGR5H.R
index a78946b618549ee741ee2c06daf58000ea147505..8cb99dc0d13c4c7a23be6d39a901a4658ded93d6 100644
--- a/R/RunModel_CemaNeigeGR5H.R
+++ b/R/RunModel_CemaNeigeGR5H.R
@@ -177,5 +177,6 @@ RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
RunOptions,
RESULTS,
LInputSeries,
+ Param,
CemaNeigeLayers)
}
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index cdc4faaeae4263dd5a54f3559f67ec7a35b4b899..1b42b04fcc9d7d220d48947621a4539f8570ce21 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -164,5 +164,6 @@ RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
RunOptions,
RESULTS,
LInputSeries,
+ Param,
CemaNeigeLayers)
}
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index 6199e4178d8e02f7ffd8930eee623105e8f3b51b..dc011ce480f668ebc20f866354c7f50597007c8e 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -169,6 +169,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
RunOptions,
RESULTS,
LInputSeries,
+ Param,
CemaNeigeLayers)
}
diff --git a/R/RunModel_GR1A.R b/R/RunModel_GR1A.R
index 87ff6aceb6d53f9b745f2542c2ed85b29f223b2f..b4ddd4f5a784b4b34ca7989c1e8ddc4ba64aa2b6 100644
--- a/R/RunModel_GR1A.R
+++ b/R/RunModel_GR1A.R
@@ -47,5 +47,6 @@ RunModel_GR1A <- function(InputsModel, RunOptions, Param) {
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
- LInputSeries)
+ LInputSeries,
+ Param)
}
diff --git a/R/RunModel_GR2M.R b/R/RunModel_GR2M.R
index dd2f696d5cc30f0dd5e98f8f7a4198f4964aca6c..6b89b73dafe44c16ee065b69bddc3edbf317c5b6 100644
--- a/R/RunModel_GR2M.R
+++ b/R/RunModel_GR2M.R
@@ -69,5 +69,6 @@ RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
- LInputSeries)
+ LInputSeries,
+ Param)
}
diff --git a/R/RunModel_GR4H.R b/R/RunModel_GR4H.R
index f70f3a1d0f165c51ef451e44d56ca50f83a685b6..01caa6d78f781517f1baa63a5638cb2db8f7fd55 100644
--- a/R/RunModel_GR4H.R
+++ b/R/RunModel_GR4H.R
@@ -74,5 +74,6 @@ RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
- LInputSeries)
+ LInputSeries,
+ Param)
}
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index fab8517e6383fbcfcf0244bd285d547ead1c3ea2..88f129d3c789df37c013b4779bf9d6dc64cd4c50 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -69,5 +69,6 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
- LInputSeries)
+ LInputSeries,
+ Param)
}
diff --git a/R/RunModel_GR5H.R b/R/RunModel_GR5H.R
index 1d0139741fbcd7cda8bd567ee388ba7ef9805580..6c5fb4323b477b1316d51feede865e405c23206f 100644
--- a/R/RunModel_GR5H.R
+++ b/R/RunModel_GR5H.R
@@ -88,5 +88,6 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
- LInputSeries)
+ LInputSeries,
+ Param)
}
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index 73e17c6cb32116164131153a9e6decf7f49ce39b..f5d2bdc34560ac16076e50a5af4dd985d99c7921 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -74,5 +74,6 @@ RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
- LInputSeries)
+ LInputSeries,
+ Param)
}
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index 659a0b5cbece1c88aa883fc3a5d4247fccd0f98e..8175af9564636f3c2392451918877fd7faf1dcbe 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -79,5 +79,6 @@ RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
- LInputSeries)
+ LInputSeries,
+ Param)
}
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 7889967059b1646cc2e8c1d368897f8f3a28fdad..63c99bcbde3e487a1caad606c40dee665b59a068 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -1,7 +1,7 @@
RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
NParam <- 1
- ##Arguments_check
+ ## argument check
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
@@ -21,35 +21,49 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
if (is.null(QcontribDown$Qsim)) {
stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
}
- OutputsModel <- QcontribDown
- OutputsModel$QsimDown <- OutputsModel$Qsim
+ if (length(QcontribDown$Qsim) != length(RunOptions$IndPeriod_Run)) {
+ stop("Time series Qsim in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
+ }
+ if (!is.null(QcontribDown$WarmUpQsim) &&
+ length(QcontribDown$WarmUpQsim) != length(RunOptions$IndPeriod_WarmUp) &&
+ RunOptions$IndPeriod_WarmUp != 0L) {
+ stop("Time series WarmUpQsim in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_WarmUp'")
+ }
} else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
- OutputsModel <- list()
- class(OutputsModel) <- c("OutputsModel", class(OutputsModel))
- OutputsModel$QsimDown <- QcontribDown
+ if (length(QcontribDown) != length(RunOptions$IndPeriod_Run)) {
+ stop("'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
+ }
} else {
stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
}
- if (length(OutputsModel$QsimDown) != length(RunOptions$IndPeriod_Run)) {
- stop("Time series in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
+
+ # data set up
+ NbUpBasins <- length(InputsModel$LengthHydro)
+ if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
+ RunOptions$IndPeriod_WarmUp <- NULL
}
+ IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
+ LInputSeries <- as.integer(length(IndPeriod1))
+ IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
- if (inherits(InputsModel, "hourly")) {
- TimeStep <- 60 * 60
- } else if (inherits(InputsModel, "daily")) {
- TimeStep <- 60 * 60 * 24
- } else {
- stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
+ if (inherits(QcontribDown, "OutputsModel")) {
+ OutputsModel <- QcontribDown
+ if (is.null(OutputsModel$WarmUpQsim)) {
+ OutputsModel$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
+ }
+ QsimDown <- c(OutputsModel$WarmUpQsim, OutputsModel$Qsim)
+ } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
+ OutputsModel <- list()
+ class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
+ QsimDown <- c(rep(NA, length(RunOptions$IndPeriod_WarmUp)),
+ QcontribDown)
}
- # propagation time from upstream meshes to outlet
- PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / TimeStep
+ ## propagation time from upstream meshes to outlet
+ PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / RunOptions$FeatFUN_MOD$TimeStep
HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
- NbUpBasins <- length(InputsModel$LengthHydro)
- LengthTs <- length(OutputsModel$QsimDown)
- OutputsModel$Qsim_m3 <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
-
+ ## set up initial states
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
if (length(IniSD) > 0) {
if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
@@ -66,15 +80,15 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
if (x > 1) {
iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
}
- IniSD[iStart:(iStart + PT[x])]
+ as.vector(IniSD[iStart:(iStart + PT[x])])
})
} else {
IniStates <- lapply(
seq_len(NbUpBasins),
function(iUpBasins) {
iWarmUp <- seq.int(
- from = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
- to = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)])
+ from = max(1, IndPeriod1[1] - floor(PT[iUpBasins]) - 1),
+ to = max(1, IndPeriod1[1] - 1)
)
ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
if (length(ini) != floor(PT[iUpBasins] + 1)) {
@@ -85,19 +99,37 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
}
)
}
- # message("Initstates: ", paste(IniStates, collapse = ", "))
+ # message("IniStates: ", paste(IniStates, collapse = ", "))
+
+ ## Lag model computation
+ Qsim_m3 <- QsimDown *
+ InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
for (upstream_basin in seq_len(NbUpBasins)) {
Qupstream <- c(IniStates[[upstream_basin]],
- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin])
+ InputsModel$Qupstream[IndPeriod1, upstream_basin])
# message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
- OutputsModel$Qsim_m3 <- OutputsModel$Qsim_m3 +
- Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
- Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
+ Qsim_m3 <- Qsim_m3 +
+ Qupstream[2:(1 + LInputSeries)] * HUTRANS[1, upstream_basin] +
+ Qupstream[1:LInputSeries] * HUTRANS[2, upstream_basin]
+ }
+
+ ## OutputsModel
+
+ if ("Qsim_m3" %in% RunOptions$Outputs_Sim) {
+ OutputsModel$Qsim_m3 <- Qsim_m3[IndPeriod2]
+ }
+
+ if ("Qsim" %in% RunOptions$Outputs_Sim) {
+ # Convert back Qsim to mm
+ OutputsModel$Qsim <- Qsim_m3[IndPeriod2] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+ # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
+ }
+
+ if ("QsimDown" %in% RunOptions$Outputs_Sim) {
+ # Convert back Qsim to mm
+ OutputsModel$QsimDown <- QsimDown[IndPeriod2]
}
- # Convert back Qsim to mm
- OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
- # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
# Warning for negative flows or NAs only in extended outputs
if (length(RunOptions$Outputs_Sim) > 2) {
@@ -112,12 +144,26 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
}
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
- OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
+ SD <- lapply(seq(NbUpBasins), function(x) {
lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
})
+ if(is.null(OutputsModel$StateEnd)) {
+ OutputsModel$StateEnd <- CreateIniStates(RunModel_Lag, InputsModel, SD = SD)
+ } else {
+ OutputsModel$StateEnd$SD <- SD
+ }
# message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
}
+ if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+ OutputsModel$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+ }
+
+ if ("Param" %in% RunOptions$Outputs_Sim) {
+ OutputsModel$Param <- c(Param, OutputsModel$Param)
+ }
+
+ class(OutputsModel) <- c(class(OutputsModel), "SD")
return(OutputsModel)
}
diff --git a/R/SeriesAggreg.OutputsModel.R b/R/SeriesAggreg.OutputsModel.R
index d065182bab683414545511a75f78aa09e6f7c58d..b714aeb185e6f5adf384d75e826624d1d7ae04aa 100644
--- a/R/SeriesAggreg.OutputsModel.R
+++ b/R/SeriesAggreg.OutputsModel.R
@@ -2,6 +2,6 @@ SeriesAggreg.OutputsModel <- function(x, Format, ...) {
SeriesAggreg.list(x,
Format,
ConvertFun = NA,
- except = c("WarmUpQsim", "StateEnd"),
+ except = c("WarmUpQsim", "StateEnd", "Param"),
...)
}
diff --git a/R/Utils.R b/R/Utils.R
index d5e14535fc3839bd4564168b99a56f40d24a39b2..19892f3ebc110cf4e5a590361ca795e263fb82aa 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -62,8 +62,10 @@
res$TimeStep <- res$TimeStep * 3600
res$TimeStepMean <- as.integer(res$TimeStepMean * 3600)
res$Class <- c(res$TimeUnit, res$Class)
+ res$CodeModHydro <- res$CodeMod
if (grepl("CemaNeige", res$NameFunMod)) {
res$Class <- c(res$Class, "CemaNeige")
+ res$CodeModHydro <- gsub("CemaNeige", "", res$CodeMod)
}
res$Class <- res$Class[!is.na(res$Class)]
if (!is.null(DatesR)) {
@@ -71,7 +73,6 @@
stop("the time step of the model inputs must be ", res$TimeUnit)
}
}
- res$CodeModHydro <- gsub("CemaNeige", "", res$CodeMod)
return(res)
}
}
diff --git a/R/UtilsCalibOptions.R b/R/UtilsCalibOptions.R
new file mode 100644
index 0000000000000000000000000000000000000000..aa4aaf5f4eeb8ea4418fbd2f86b2d8b53e17efae
--- /dev/null
+++ b/R/UtilsCalibOptions.R
@@ -0,0 +1,132 @@
+.FunTransfo <- function(FeatFUN_MOD) {
+
+ IsHyst <- FeatFUN_MOD$IsHyst
+ IsSD <- FeatFUN_MOD$IsSD
+
+ ## set FUN_GR
+ if (FeatFUN_MOD$NameFunMod == "Cemaneige") {
+ if (IsHyst) {
+ FUN_GR <- TransfoParam_CemaNeigeHyst
+ } else {
+ FUN_GR <- TransfoParam_CemaNeige
+ }
+ } else {
+ # fatal error if the TransfoParam function does not exist
+ FUN_GR <- match.fun(sprintf("TransfoParam_%s", FeatFUN_MOD$CodeModHydro))
+ }
+
+ ## set FUN_SNOW
+ if ("CemaNeige" %in% FeatFUN_MOD$Class) {
+ if (IsHyst) {
+ FUN_SNOW <- TransfoParam_CemaNeigeHyst
+ } else {
+ FUN_SNOW <- TransfoParam_CemaNeige
+ }
+ }
+
+ ## set FUN_LAG
+ if (IsSD) {
+ FUN_LAG <- TransfoParam_Lag
+ }
+
+ ## set FUN_TRANSFO
+ if (! "CemaNeige" %in% FeatFUN_MOD$Class) {
+ if (!IsSD) {
+ FUN_TRANSFO <- FUN_GR
+ } else {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ ParamOut[, 2:NParam] <- FUN_GR(ParamIn[, 2:NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
+ } else {
+ if (IsHyst & !IsSD) {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ ParamOut[, 1:(NParam - 4) ] <- FUN_GR(ParamIn[, 1:(NParam - 4)], Direction)
+ ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
+ if (!IsHyst & !IsSD) {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ if (NParam <= 3) {
+ ParamOut[, 1:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
+ } else {
+ ParamOut[, 1:(NParam - 2)] <- FUN_GR(ParamIn[, 1:(NParam - 2)], Direction)
+ }
+ ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
+ if (IsHyst & IsSD) {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ ParamOut[, 2:(NParam - 4) ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
+ ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
+ if (!IsHyst & IsSD) {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ if (NParam <= 3) {
+ ParamOut[, 2:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
+ } else {
+ ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)], Direction)
+ }
+ ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
+ }
+ if (is.null(FUN_TRANSFO)) {
+ stop("'FUN_TRANSFO' was not found")
+ }
+ return(FUN_TRANSFO)
+}
diff --git a/R/UtilsErrorCrit.R b/R/UtilsErrorCrit.R
index 4e200bb28ec42c6378c4ab75ab7cfb59868a66ca..8c8609f003934bdfa7f798fff18f30cdde989e30 100644
--- a/R/UtilsErrorCrit.R
+++ b/R/UtilsErrorCrit.R
@@ -9,11 +9,7 @@
stop("'InputsCrit' must be of class 'InputsCrit'", call. = FALSE)
}
if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
- if (Crit == "RMSE") {
- stop("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' with RMSE", call. = FALSE)
- } else {
- stop(paste0("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' or 'Compo' with ", Crit), call. = FALSE)
- }
+ stop(paste0("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' or 'Compo' with ", Crit), call. = FALSE)
}
@@ -39,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
}
@@ -48,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
@@ -115,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)
+ }
}
diff --git a/R/UtilsRunModel.R b/R/UtilsRunModel.R
index 096aff84a47a57e68253a14b75b3d3a9b56115e3..846f018214f2b4ac21805d727b3974e6be1896ac 100644
--- a/R/UtilsRunModel.R
+++ b/R/UtilsRunModel.R
@@ -4,6 +4,7 @@
#' @param RunOptions output of [CreateRunOptions]
#' @param RESULTS outputs of [.Fortran]
#' @param LInputSeries number of time steps of warm-up + run periods
+#' @param Param [numeric] vector of model parameters
#' @param CemaNeigeLayers outputs of Cemaneige pre-process
#'
#' @return OutputsModel object
@@ -13,6 +14,7 @@
RunOptions,
RESULTS,
LInputSeries,
+ Param,
CemaNeigeLayers = NULL) {
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
@@ -46,6 +48,10 @@
OutputsModel$StateEnd <- RESULTS$StateEnd
}
+ if ("Param" %in% RunOptions$Outputs_Sim) {
+ OutputsModel$Param <- Param
+ }
+
class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
return(OutputsModel)
diff --git a/inst/modelsFeatures/FeatModelsGR.csv b/inst/modelsFeatures/FeatModelsGR.csv
index 21359f84b61a9b1b3e6b83822e20ad5449a66132..38d9e5ca7bfab6eef3e8932e1f3c9ebadd5a220a 100644
--- a/inst/modelsFeatures/FeatModelsGR.csv
+++ b/inst/modelsFeatures/FeatModelsGR.csv
@@ -12,3 +12,4 @@ CemaNeigeGR5J;CemaNeigeGR5J;7;daily;NA;GR;airGR
CemaNeigeGR6J;CemaNeigeGR6J;8;daily;NA;GR;airGR
CemaNeigeGR4H;CemaNeigeGR4H;6;hourly;NA;GR;airGR
CemaNeigeGR5H;CemaNeigeGR5H;7;hourly;NA;GR;airGR
+Lag;Lag;1;NA;NA;SD;airGR
diff --git a/man/CreateErrorCrit_GAPX.Rd b/man/CreateErrorCrit_GAPX.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..487e42f540c040e75f02991b2df525c0f336d698
--- /dev/null
+++ b/man/CreateErrorCrit_GAPX.Rd
@@ -0,0 +1,93 @@
+\encoding{UTF-8}
+
+
+\name{CreateErrorCrit_GAPX}
+\alias{CreateErrorCrit_GAPX}
+
+
+\title{Creation of the ErrorCrit_GAPX function}
+
+
+\description{
+Function which creates the function \code{ErrorCrit_GAPX}.
+
+The produced function \code{ErrorCrit_GAPX} allows computing an error criterion based on the GAPX formula proposed by Lavenne et al. (2019).
+}
+
+
+\usage{
+CreateErrorCrit_GAPX(FUN_TRANSFO)
+}
+
+
+\arguments{
+\item{FUN_TRANSFO}{[function] The parameter transformation function used with the model}
+}
+
+
+\value{
+ [function] function \code{ErrorCrit_GAPX} embedding the parameter transformation function used with the model
+}
+
+\details{
+In addition to the criterion value, the function outputs include a multiplier (-1 or +1) that allows
+the use of the function for model calibration: the product \code{CritValue * Multiplier} is the criterion to be minimised
+(\code{Multiplier = -1} for NSE).
+}
+
+\examples{
+library(airGR)
+
+## loading catchment data
+data(L0123001)
+
+## preparation of the InputsModel object
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P, PotEvap = BasinObs$E)
+
+## calibration period selection
+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"))
+
+## preparation of RunOptions object
+RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+ IndPeriod_Run = Ind_Run)
+
+## simulation
+Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions, Param = Param)
+
+## Creation of the ErrorCrit GAPX function
+ErrorCrit_GAPX <- CreateErrorCrit_GAPX(TransfoParam_GR4J)
+
+## The "a priori" parameters for GAPX
+AprParamR <- c(X1 = 157, X2 = 0.8, X3 = 100, X4 = 1.5)
+AprParamT <- TransfoParam_GR4J(AprParamR, "RT")
+
+## Single efficiency criterion: GAPX with a priori parameters
+InputsCrit <- CreateInputsCrit(ErrorCrit_GAPX,
+ InputsModel,
+ RunOptions,
+ Obs = AprParamT,
+ VarObs = "ParamT")
+ErrorCrit <- ErrorCrit_GAPX(InputsCrit, OutputsModel)
+str(ErrorCrit)
+}
+
+
+\author{
+David Dorchies
+}
+
+\references{
+de Lavenne, A., Andréassian, V., Thirel, G., Ramos, M.-H. and Perrin, C. (2019).
+ A Regularization Approach to Improve the Sequential Calibration of a Semidistributed Hydrological Model.
+ Water Resources Research 55, 8821–8839. \doi{10.1029/2018WR024266}
+}
+
+\seealso{
+\code{\link{CreateInputsCrit}}, \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}},
+\code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_KGE2}}
+}
+
diff --git a/man/CreateInputsCrit_Lavenne.Rd b/man/CreateInputsCrit_Lavenne.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..30f5e428b95e6b6c88fded04c3afa30b8d2ff301
--- /dev/null
+++ b/man/CreateInputsCrit_Lavenne.Rd
@@ -0,0 +1,129 @@
+\encoding{UTF-8}
+
+
+\name{CreateInputsCrit_Lavenne}
+\alias{CreateInputsCrit_Lavenne}
+
+
+\title{Creation of the InputsCrit object for Lavenne Criterion}
+
+
+\description{
+Creation of the \code{InputsCrit} object required to the \code{\link{ErrorCrit}} function. This function defines a composite criterion based on the formula proposed by Lavenne et al. (2019).
+}
+
+
+\usage{
+CreateInputsCrit_Lavenne(FUN_CRIT = ErrorCrit_KGE,
+ InputsModel,
+ RunOptions,
+ Obs,
+ VarObs = "Q",
+ AprParamR,
+ AprCrit = 1,
+ k = 0.15,
+ BoolCrit = NULL,
+ transfo = "sqrt",
+ epsilon = NULL)
+}
+
+
+\arguments{
+\item{FUN_CRIT}{[function] error criterion function (e.g. \code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_NSE}}). Default \code{\link{ErrorCrit_KGE}}}
+
+\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
+
+\item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
+
+\item{Obs}{[numeric (atomic or list)] series of observed variable ([mm/time step] for discharge or SWE, [-] for SCA)}
+
+\item{VarObs}{(optional) [character (atomic or list)] names of the observed variable (\code{"Q"} by default, or one of \code{"SCA"}, \code{"SWE"})}
+
+\item{AprParamR}{[numeric] a priori parameter set from a donor catchment}
+
+\item{AprCrit}{(optional) [numeric] performance criterion of the donor catchment (1 by default)}
+
+\item{k}{(optional) [numeric] coefficient used for the weighted average between \code{FUN_CRIT} and the gap between the optimised parameter set and an a priori parameter set calculated with the function produced by \code{\link{CreateErrorCrit_GAPX}}}
+
+\item{BoolCrit}{(optional) [boolean] boolean (the same length as \code{Obs}) giving the time steps to consider in the computation (all time steps are considered by default. See details)}
+
+\item{transfo}{(optional) [character] name of the transformation applied to the variables (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"}, \code{"boxcox"} or a numeric value for power transformation for \code{FUN_CRIT}. Default value is \code{"sqrt"}. See details of \code{\link{CreateInputsCrit}}}
+
+\item{epsilon}{(optional) [numeric] small value to add to all observations and simulations for \code{FUN_CRIT} when \code{"log"} or \code{"inv"} transformations are used [same unit as \code{Obs}]. See details of \code{\link{CreateInputsCrit}}}
+}
+
+
+\value{
+
+[list] object of class InputsCrit containing the data required to evaluate the model outputs (see \code{\link{CreateInputsCrit}} for more details).
+
+\code{CreateInputsCrit_Lavenne} returns an object of class \emph{Compo}.
+
+Items \code{Weights} of the criteria are respectively equal to \code{k} and \code{k * max(0, AprCrit)}.
+
+To calculate the Lavenne criterion, it is necessary to use the \code{ErrorCrit} function as for any other composite criterion.
+}
+
+
+\details{
+
+The parameters \code{FUN_CRIT}, \code{Obs}, \code{VarObs}, \code{BoolCrit}, \code{transfo}, and \code{epsilon} must be used as they would be used for \code{\link{CreateInputsCrit}} in the case of a single criterion.
+
+\code{\link{ErrorCrit_RMSE}} cannot be used in a composite criterion since it is not a unitless value.
+
+
+\code{CreateInputsCrit_Lavenne} creates a composite criterion in respect with Equations 1 and 2 of de Lavenne et al. (2019).
+}
+
+
+\examples{
+library(airGR)
+
+## loading catchment data
+data(L0123001)
+
+## preparation of the InputsModel object
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P, PotEvap = BasinObs$E)
+
+## calibration period selection
+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"))
+
+## preparation of RunOptions object
+RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+ IndPeriod_Run = Ind_Run)
+
+## simulation
+Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions, Param = Param)
+
+## The "a priori" parameters for the Lavenne formula
+AprParamR <- c(X1 = 157, X2 = 0.8, X3 = 100, X4 = 1.5)
+
+## Single efficiency criterion: GAPX with a priori parameters
+IC_DL <- CreateInputsCrit_Lavenne(InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Obs = BasinObs$Qmm[Ind_Run],
+ AprParamR = AprParamR)
+str(ErrorCrit(InputsCrit = IC_DL, OutputsModel = OutputsModel))
+}
+
+
+\author{
+David Dorchies
+}
+
+
+\references{
+de Lavenne, A., Andréassian, V., Thirel, G., Ramos, M.-H. and Perrin, C. (2019).
+ A Regularization Approach to Improve the Sequential Calibration of a Semidistributed Hydrological Model.
+ Water Resources Research 55, 8821–8839. \doi{10.1029/2018WR024266}
+}
+
+
+\seealso{
+\code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateCalibOptions}}, \code{\link{ErrorCrit}}
+}
+
diff --git a/man/RunModel_Lag.Rd b/man/RunModel_Lag.Rd
index c5eebcc043809c6758af67e5c86d8fb534a834f4..ee34a5980f5085480b464f6824c7159ce616cf0d 100644
--- a/man/RunModel_Lag.Rd
+++ b/man/RunModel_Lag.Rd
@@ -51,6 +51,26 @@ The list value contains an extra item named \code{QsimDown} which is a copy of t
library(airGR)
data(L0123001)
+## ---- simulation of the hydrological catchment with GR4J
+
+InputsModelDown <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P, PotEvap = BasinObs$E)
+
+## run period selection
+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"))
+
+## creation of the RunOptions object
+RunOptionsDown <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModelDown, IndPeriod_Run = Ind_Run)
+
+## simulation of the runoff of the catchment with a GR4J model
+Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+OutputsModelDown <- RunModel_GR4J(InputsModel = InputsModelDown,
+ RunOptions = RunOptionsDown, Param = Param)
+
+## ---- specifications of the reservoir
+
## the reservoir withdraws 1 m3/s when it's possible considering the flow observed in the basin
Qupstream <- matrix(-sapply(BasinObs$Qls / 1000 - 1, function(x) {
min(1, max(0, x, na.rm = TRUE))
@@ -67,32 +87,21 @@ BasinAreas <- c(NA, BasinInfo$BasinArea)
## delay time between the reservoir and the catchment outlet is 2 days and the distance is 150 km
LengthHydro <- 150
-InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+## ---- simulation of the basin with the reservoir influence
+
+InputsModelInf <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
Qupstream = Qupstream, LengthHydro = LengthHydro,
BasinAreas = BasinAreas)
-
-
-## ---- simulation of the basin with the reservoir influence
-
-## run period selection
-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"))
-
## creation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
-
-## simulation of the runoff of the catchment with a GR4J model
-Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
-OutputsModelDown <- RunModel_GR4J(InputsModel = InputsModel,
- RunOptions = RunOptions, Param = Param)
+ InputsModel = InputsModelInf, IndPeriod_Run = Ind_Run)
## with a delay of 2 days for 150 km, the flow velocity is 75 km per day
Velocity <- (LengthHydro * 1e3 / 2) / (24 * 60 * 60) ## Conversion km/day -> m/s
## run the lag model which routes precipitation-runoff model and upstream flows
-OutputsModel <- RunModel_Lag(InputsModel = InputsModel,
+OutputsModel <- RunModel_Lag(InputsModel = InputsModelInf,
RunOptions = RunOptions,
Param = Velocity,
QcontribDown = OutputsModelDown)
diff --git a/tests/testthat/regression.R b/tests/testthat/regression.R
index 1732b2c544430557ada828b159717e8e0d8e2179..147b73c500f537fa8ce0d7dba1c1983e1e9e5e05 100644
--- a/tests/testthat/regression.R
+++ b/tests/testthat/regression.R
@@ -1,9 +1,9 @@
context("Compare example outputs with CRAN")
CompareWithStable <- function(refVarFile, testDir, regIgnore) {
- v <- data.frame(topic = basename(dirname(refVarFile)),
+ v <- list(topic = basename(dirname(refVarFile)),
var = gsub("\\.rds$", "", basename(refVarFile)))
- if (is.null(regIgnore) || all(apply(regIgnore, 1, function(x) !all(x == v)))) {
+ if (is.null(regIgnore) || !any(apply(regIgnore, 1, function(x) {v$var == x[2] && (x[1] == "*" || x[1] == v$topic)}))) {
test_that(paste("Compare", v$topic, v$var), {
testVarFile <- paste0(
file.path(testDir, v$topic, v$var),
@@ -13,6 +13,17 @@ CompareWithStable <- function(refVarFile, testDir, regIgnore) {
if (file.exists(testVarFile)) {
testVar <- readRDS(testVarFile)
refVar <- readRDS(refVarFile)
+ if (!is.null(regIgnore)) {
+ regIgnore$mainVar <- gsub("\\$.*$", "", regIgnore[,2])
+ varIgnore <- regIgnore[apply(regIgnore, 1, function(x) {v$var == x[3] && (x[1] == "*" || x[1] == v$topic)}), 2]
+ if (length(varIgnore) > 0) {
+ itemIgnore <- gsub("^.*\\$", "", varIgnore)
+ for(item in itemIgnore) {
+ testVar[[item]] <- NULL
+ refVar[[item]] <- NULL
+ }
+ }
+ }
expect_equivalent(testVar, refVar)
}
})
diff --git a/tests/testthat/test-CreateErrorCrit_GAPX.R b/tests/testthat/test-CreateErrorCrit_GAPX.R
new file mode 100644
index 0000000000000000000000000000000000000000..3b41c84347915e655308a6229ae428df4159fce0
--- /dev/null
+++ b/tests/testthat/test-CreateErrorCrit_GAPX.R
@@ -0,0 +1,43 @@
+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)
+
+ErrorCrit_GAPX <- CreateErrorCrit_GAPX(TransfoParam_GR4J)
+ParamT <- TransfoParam_GR4J(Param, "RT")
+
+test_that("ErrorCrit should return 1 for same parameters", {
+ IC <- CreateInputsCrit(ErrorCrit_GAPX, InputsModel, RunOptions, Obs = ParamT, VarObs = "ParamT")
+ expect_equal(ErrorCrit_GAPX(IC, OutputsModel)$CritValue, 1)
+})
+
+test_that("ErrorCrit should return 1-nbParam^0.5/40 for ParamT shifted by 1", {
+ ParamT <- ParamT + 1
+ IC <- CreateInputsCrit(ErrorCrit_GAPX, InputsModel, RunOptions, Obs = ParamT, VarObs = "ParamT")
+ expect_equal(ErrorCrit_GAPX(IC, OutputsModel)$CritValue,
+ 1 - RunOptions$FeatFUN_MOD$NbParam^0.5 / 20)
+})
+
+test_that("ErrorCrit should return 1-(nbParam-1)^0.5/40 for ParamT shifted by 1 with one NA", {
+ ParamT <- ParamT + 1
+ ParamT[1] <- NA
+ IC <- CreateInputsCrit(ErrorCrit_GAPX, InputsModel, RunOptions, Obs = ParamT, VarObs = "ParamT")
+ expect_equal(suppressWarnings(ErrorCrit_GAPX(IC, OutputsModel)$CritValue),
+ 1 - (RunOptions$FeatFUN_MOD$NbParam - 1)^0.5 / 20)
+ expect_warning(ErrorCrit_GAPX(IC, OutputsModel)$CritValue,
+ regexp = "criterion GAPX computed on less than 4 parameters")
+})
diff --git a/tests/testthat/test-CreateInputsCrit.R b/tests/testthat/test-CreateInputsCrit.R
new file mode 100644
index 0000000000000000000000000000000000000000..bbfb0ae935a948776306690806e9013ce1752f4a
--- /dev/null
+++ b/tests/testthat/test-CreateInputsCrit.R
@@ -0,0 +1,40 @@
+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"
+ )
+})
diff --git a/tests/testthat/test-CreateInputsCrit_Lavenne.R b/tests/testthat/test-CreateInputsCrit_Lavenne.R
new file mode 100644
index 0000000000000000000000000000000000000000..c61f1a0055f362696e10409ca883e2771db78361
--- /dev/null
+++ b/tests/testthat/test-CreateInputsCrit_Lavenne.R
@@ -0,0 +1,51 @@
+context("CreateInputsCrit_Lavenne")
+
+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)
+
+k <- 0.15
+IC_DL <- CreateInputsCrit_Lavenne(FUN_CRIT = ErrorCrit_KGE,
+ InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Obs = BasinObs$Qmm[Ind_Run],
+ AprParamR = Param,
+ k = k)
+
+test_that("should return a composed ErrorCrit", {
+ expect_s3_class(IC_DL, "Compo")
+ expect_length(IC_DL, 2)
+ AprParamT <- TransfoParam_GR4J(Param, "RT")
+ expect_equal(IC_DL$IC2$Obs, AprParamT)
+})
+
+test_that("should return KGE*(1-k)+k with parameters matching a priori parameters", {
+ IC_KGE <- CreateInputsCrit(ErrorCrit_KGE,
+ InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Obs = BasinObs$Qmm[Ind_Run],
+ transfo = "sqrt")
+ expect_equal(ErrorCrit_KGE(IC_KGE, OutputsModel)$CritValue * (1 - k) + k,
+ ErrorCrit(IC_DL, OutputsModel)$CritValue)
+})
+
+test_that("should return proper error if mismatch number of parameters", {
+ expect_error(
+ CreateInputsCrit_Lavenne(FUN_CRIT = ErrorCrit_KGE,
+ InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Obs = BasinObs$Qmm[Ind_Run],
+ AprParamR = Param[-1],
+ k = k),
+ regex = "'AprParamR' must be a numeric vector of length 4"
+ )
+})
diff --git a/tests/testthat/test-CreateiniStates.R b/tests/testthat/test-CreateiniStates.R
index 17ac824be4cdb090bfdca6c7d2ddd3b8dc089a54..3f1cb1c6f0ba7e1b9e427295b5829d8bcdccceff 100644
--- a/tests/testthat/test-CreateiniStates.R
+++ b/tests/testthat/test-CreateiniStates.R
@@ -96,3 +96,9 @@ test_that("Error: Number of items not equal to number of upstream connections",
)
})
})
+
+test_that("FUN = RunModel_lag must work", {
+ IS <- CreateIniStates(RunModel_Lag, InputsModel, SD = list(rep(0, 10)))
+ expect_equal(IS$SD[[1]], rep(0, 10))
+ expect_s3_class(IS, "SD")
+})
diff --git a/tests/testthat/test-RunModel_Lag.R b/tests/testthat/test-RunModel_Lag.R
index 81537b6bda92f61c04007ef532d60207141d7f07..23889967f71f1fef5983d51c86c92931481842fd 100644
--- a/tests/testthat/test-RunModel_Lag.R
+++ b/tests/testthat/test-RunModel_Lag.R
@@ -57,8 +57,10 @@ test_that("QcontribDown parameter should be a numeric vector or an OutputModel o
Param <- c(257.237556, 1.012237, 88.234673, 2.207958) # From vignettes/V01_get_started
+RunOptionsGR4J <- RunOptions
+RunOptionsGR4J$FeatFUN_MOD$NbParam <- 4
OutputsGR4JOnly <- RunModel_GR4J(InputsModel = InputsModel,
- RunOptions = RunOptions,
+ RunOptions = RunOptionsGR4J,
Param = Param)
test_that("QcontribDown should contain a Qsim key", {
@@ -79,6 +81,14 @@ test_that("'QcontribDown$Qim' should have the same lenght as 'RunOptions$IndPeri
)
})
+test_that("OutputsModel must have a item 'QsimDown' equal to GR4J Qsim contribution", {
+ expect_equal(OutputsGR4JOnly$Qsim,
+ RunModel_Lag(InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Param = 1,
+ QcontribDown = OutputsGR4JOnly)$QsimDown)
+})
+
test_that("RunModel(FUN=RunModel_Lag) should give same result as RunModel_Lag", {
QcontribDown <- OutputsGR4JOnly
Output_RunModel_Lag <- RunModel_Lag(InputsModel = InputsModel,
diff --git a/vignettes/V00_airgr_ref.bib b/vignettes/V00_airgr_ref.bib
index 419563ce6b64db4d2cbe0650532437ed573f14de..e618e4fbb7723b5837347c8f7be774935309811f 100644
--- a/vignettes/V00_airgr_ref.bib
+++ b/vignettes/V00_airgr_ref.bib
@@ -178,4 +178,25 @@
year = {2019},
keywords = {airGRcite},
pages = {70--81}
+}
+
+@article{delavenne_regularization_2019,
+ ids = {lavenneRegularizationApproachImprove2019a},
+ title = {A {{Regularization Approach}} to {{Improve}} the {{Sequential Calibration}} of a {{Semidistributed Hydrological Model}}},
+ author = {de Lavenne, Alban and Andréassian, Vazken and Thirel, Guillaume and Ramos, Maria-Helena and Perrin, Charles},
+ date = {2019},
+ journaltitle = {Water Resources Research},
+ volume = {55},
+ pages = {8821--8839},
+ issn = {1944-7973},
+ doi = {10.1029/2018WR024266},
+ url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018WR024266},
+ urldate = {2021-06-07},
+ abstract = {In semidistributed hydrological modeling, sequential calibration usually refers to the calibration of a model by considering not only the flows observed at the outlet of a catchment but also the different gauging points inside the catchment from upstream to downstream. While sequential calibration aims to optimize the performance at these interior gauged points, we show that it generally fails to improve performance at ungauged points. In this paper, we propose a regularization approach for the sequential calibration of semidistributed hydrological models. It consists in adding a priori information on optimal parameter sets for each modeling unit of the semidistributed model. Calibration iterations are then performed by jointly maximizing simulation performance and minimizing drifts from the a priori parameter sets. The combination of these two sources of information is handled by a parameter k to which the method is quite sensitive. The method is applied to 1,305 catchments in France over 30 years. The leave-one-out validation shows that, at locations considered as ungauged, model simulations are significantly improved (over all the catchments, the median KGE criterion is increased from 0.75 to 0.83 and the first quartile from 0.35 to 0.66), while model performance at gauged points is not significantly impacted by the use of the regularization approach. Small catchments benefit most from this calibration strategy. These performances are, however, very similar to the performances obtained with a lumped model based on similar conceptualization.},
+ annotation = {\_eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2018WR024266},
+ file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\QNCGSJXB\\Lavenne et al. - 2019 - A Regularization Approach to Improve the Sequentia.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\5DMRIV8Y\\2018WR024266.html},
+ keywords = {regularization,semidistributed model,stepwise calibration,ungauged basins},
+ langid = {english},
+ number = {11},
+ options = {useprefix=true}
}
\ No newline at end of file
diff --git a/vignettes/V05_sd_model.Rmd b/vignettes/V05_sd_model.Rmd
index 800db1ee9ee77a6487f07122ba15040da61ca3aa..6c17770e9c503ee166a9ecc34c723106e65986e5 100644
--- a/vignettes/V05_sd_model.Rmd
+++ b/vignettes/V05_sd_model.Rmd
@@ -22,16 +22,18 @@ The **airGR** package implements semi-distributed model capabilities using a lag
`RunModel_Lag` documentation gives an example of simulating the influence of a reservoir in a lumped model. Try `example(RunModel_Lag)` to get it.
-In this vignette, we show how to calibrate 2 sub-catchments in series with a semi-distributed model consisting of 2 GR4J models. For doing this we compare two strategies for calibrating the downstream subcatchment:
+In this vignette, we show how to calibrate 2 sub-catchments in series with a semi-distributed model consisting of 2 GR4J models.
+For doing this we compare 3 strategies for calibrating the downstream subcatchment:
- using upstream observed flows
- using upstream simulated flows
+- using upstream simulated flows and parameter regularisation [@delavenne_regularization_2019]
We finally compare these calibrations with a theoretical set of parameters.
+This comparison is based on the Kling-Gupta Efficiency computed on the root-squared discharges as performance criterion.
## Model description
-
```{r, warning=FALSE, include=FALSE}
library(airGR)
options(digits = 3)
@@ -60,6 +62,13 @@ summary(cbind(QObsUp = BasinObs$Qmm, QObsDown))
options(digits = 3)
```
+With a delay of 2 days between the 2 gauging stations, the theoretical Velocity parameter should be equal to:
+
+```{r}
+Velocity <- 100 * 1e3 / (2 * 86400)
+paste("Velocity: ", format(Velocity), "m/s")
+```
+
# Calibration of the upstream subcatchment
The operations are exactly the same as the ones for a GR4J lumped model. So we do exactly the same operations as in the [Get Started](V01_get_started.html) vignette.
@@ -73,9 +82,11 @@ RunOptionsUp <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModelUp,
IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
IniStates = NULL, IniResLevels = NULL)
-InputsCritUp <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModelUp,
+# Error criterion is KGE computed on the root-squared discharges
+InputsCritUp <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelUp,
RunOptions = RunOptionsUp,
- VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run])
+ VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run],
+ transfo = "sqrt")
CalibOptionsUp <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
OutputsCalibUp <- Calibration_Michel(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
InputsCrit = InputsCritUp, CalibOptions = CalibOptionsUp,
@@ -90,9 +101,11 @@ OutputsModelUp <- RunModel_GR4J(InputsModel = InputsModelUp, RunOptions = RunOpt
```
-# Calibration of the downstream subcatchment with upstream flow observations
+# Calibration of the downstream subcatchment
+
+## Creation of the InputsModel objects
-we need to create the `InputsModel` object completed with upstream information:
+we need to create `InputsModel` objects completed with upstream information with upstream observed flow for the calibration of first case and upstream simulated flows for the other cases:
```{r}
InputsModelDown1 <- CreateInputsModel(
@@ -104,16 +117,38 @@ InputsModelDown1 <- CreateInputsModel(
)
```
-And then calibrate the combination of Lag model for upstream flow transfer and GR4J model for the runoff of the downstream subcatchment:
+For using upstream simulated flows, we should concatenate a vector with the simulated flows for the entire period of simulation (warm-up + run):
+
+```{r}
+Qsim_upstream <- rep(NA, length(BasinObs$DatesR))
+# Simulated flow during warm-up period (365 days before run period)
+Qsim_upstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$WarmUpQsim
+# Simulated flow during run period
+Qsim_upstream[Ind_Run] <- OutputsModelUp$Qsim
+
+InputsModelDown2 <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P, PotEvap = BasinObs$E,
+ Qupstream = matrix(Qsim_upstream, ncol = 1), # upstream observed flow
+ LengthHydro = 100, # distance between upstream catchment outlet & the downstream one [km]
+ BasinAreas = c(180, 180) # upstream and downstream areas [km²]
+)
+```
+
+
+## Calibration with upstream flow observations
+
+We calibrate the combination of Lag model for upstream flow transfer and GR4J model for the runoff of the downstream subcatchment:
```{r}
RunOptionsDown <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModelDown1,
IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
IniStates = NULL, IniResLevels = NULL)
-InputsCritDown <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModelDown1,
+InputsCritDown <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelDown1,
RunOptions = RunOptionsDown,
- VarObs = "Q", Obs = QObsDown[Ind_Run])
+ VarObs = "Q", Obs = QObsDown[Ind_Run],
+ transfo = "sqrt")
CalibOptionsDown <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel,
IsSD = TRUE) # specify that it's a SD model
@@ -124,16 +159,6 @@ OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
FUN_MOD = RunModel_GR4J)
```
-To run the complete model, we should substitute the observed upstream flow by the simulated one for the entire period of simulation (warm-up + run):
-
-```{r}
-InputsModelDown2 <- InputsModelDown1
-# Simulated flow during warm-up period
-InputsModelDown2$Qupstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$WarmUpQsim
-# Simulated flow during run period
-InputsModelDown2$Qupstream[Ind_Run] <- OutputsModelUp$Qsim
-```
-
`RunModel` is run in order to automatically combine GR4J and Lag models.
```{r}
@@ -146,11 +171,11 @@ OutputsModelDown1 <- RunModel(InputsModel = InputsModelDown2,
Performance of the model validation is then:
```{r}
-CritDown1 <- ErrorCrit_NSE(InputsCritDown, OutputsModelDown1)
+KGE_down1 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown1)
```
-# Calibration of the downstream subcatchment with upstream simulated flow
+## Calibration with upstream simulated flow
We calibrate the model with the `InputsModel` object previously created for substituting the observed upstream flow with the simulated one:
@@ -163,67 +188,105 @@ OutputsCalibDown2 <- Calibration_Michel(InputsModel = InputsModelDown2,
ParamDown2 <- OutputsCalibDown2$ParamFinalR
```
+## Calibration with upstream simulated flow and parameter regularisation
-# Discussion
+The regularisation follow the method proposed by @delavenne_regularization_2019.
-## Identification of Velocity parameter
+As a priori parameter set, we use the calibrated parameter set of the upstream catchment and the theoretical velocity:
-The theoretical Velocity parameter should be equal to:
+```{r}
+ParamDownTheo <- c(Velocity, OutputsCalibUp$ParamFinalR)
+```
+
+The Lavenne criterion is initialised with the a priori parameter set and the value of the KGE of the upstream basin.
```{r}
-Velocity <- InputsModelDown1$LengthHydro * 1e3 / (2 * 86400)
-paste(format(Velocity), "m/s")
+IC_Lavenne <- CreateInputsCrit_Lavenne(InputsModel = InputsModelDown2,
+ RunOptions = RunOptionsDown,
+ Obs = QObsDown[Ind_Run],
+ AprParamR = ParamDownTheo,
+ AprCrit = OutputsCalibUp$CritFinal)
```
+The Lavenne criterion is used instead of the KGE for calibration with regularisation
+
+```{r}
+OutputsCalibDown3 <- Calibration_Michel(InputsModel = InputsModelDown2,
+ RunOptions = RunOptionsDown,
+ InputsCrit = IC_Lavenne,
+ CalibOptions = CalibOptionsDown,
+ FUN_MOD = RunModel_GR4J)
+```
+
+The KGE is then calculated for performance comparisons:
+
+```{r}
+OutputsModelDown3 <- RunModel(InputsModel = InputsModelDown2,
+ RunOptions = RunOptionsDown,
+ Param = OutputsCalibDown3$ParamFinalR,
+ FUN_MOD = RunModel_GR4J)
+KGE_down3 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown3)
+```
+
+
+# Discussion
+
+## Identification of Velocity parameter
+
Both calibrations overestimate this parameter:
```{r}
mVelocity <- matrix(c(Velocity,
OutputsCalibDown1$ParamFinalR[1],
- OutputsCalibDown2$ParamFinalR[1]),
+ OutputsCalibDown2$ParamFinalR[1],
+ OutputsCalibDown3$ParamFinalR[1]),
ncol = 1,
dimnames = list(c("theoretical",
"calibrated with observed upstream flow",
- "calibrated with simulated upstream flow"),
+ "calibrated with simulated upstream flow",
+ "calibrated with sim upstream flow and regularisation"),
c("Velocity parameter")))
knitr::kable(mVelocity)
```
## Value of the performance criteria with theoretical calibration
-Theoretically, the parameters of the downstream GR4J model should be the same as the upstream one and we know the lag time. So this set of parameter should give a better performance criteria:
+Theoretically, the parameters of the downstream GR4J model should be the same as the upstream one with the velocity as extra parameter:
```{r}
-ParamDownTheo <- c(Velocity, OutputsCalibUp$ParamFinalR)
OutputsModelDownTheo <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = ParamDownTheo,
FUN_MOD = RunModel_GR4J)
-CritDownTheo <- ErrorCrit_NSE(InputsCritDown, OutputsModelDownTheo)
+KGE_downTheo <- ErrorCrit_KGE(InputsCritDown, OutputsModelDownTheo)
```
-
## Parameters and performance of each subcatchment for all calibrations
```{r}
comp <- matrix(c(0, OutputsCalibUp$ParamFinalR,
rep(OutputsCalibDown1$ParamFinalR, 2),
OutputsCalibDown2$ParamFinalR,
+ OutputsCalibDown3$ParamFinalR,
ParamDownTheo),
ncol = 5, byrow = TRUE)
comp <- cbind(comp, c(OutputsCalibUp$CritFinal,
OutputsCalibDown1$CritFinal,
- CritDown1$CritValue,
+ KGE_down1$CritValue,
OutputsCalibDown2$CritFinal,
- CritDownTheo$CritValue))
-colnames(comp) <- c("Velocity", paste0("X", 1:4), "NSE")
+ KGE_down3$CritValue,
+ KGE_downTheo$CritValue))
+colnames(comp) <- c("Velocity", paste0("X", 1:4), "KGE(√Q)")
rownames(comp) <- c("Calibration of the upstream subcatchment",
"Calibration 1 with observed upstream flow",
"Validation 1 with simulated upstream flow",
"Calibration 2 with simulated upstream flow",
+ "Calibration 3 with simulated upstream flow and regularisation",
"Validation theoretical set of parameters")
knitr::kable(comp)
```
-Even if calibration with observed upstream flows gives an improved performance criteria, in validation using simulated upstream flows the result is quite similar as the performance obtained with the calibration with upstream simulated flows. The theoretical set of parameters give also an equivalent performance but still underperforming the calibration 2 one.
+Even if calibration with observed upstream flows gives an improved performance criteria, in validation using simulated upstream flows the result is quite similar as the performance obtained with the calibration with upstream simulated flows. The theoretical set of parameters give also an equivalent performance but still underperforming the calibration 2 one. Regularisation allows to get similar performance as the one for calibration with simulated flows but with the big advantage of having parameters closer to the theoretical ones (Especially for the velocity parameter).
+
+# References