From dd7c0dbb4a18b8148aa706c939139d135856631e Mon Sep 17 00:00:00 2001
From: Thirel Guillaume <guillaume.thirel@inrae.fr>
Date: Thu, 30 Apr 2020 15:04:45 +0200
Subject: [PATCH] feat: Modification of Michel's calibration for SD models
Refs #34
---
NAMESPACE | 1 +
R/Calibration_Michel.R | 107 +++++++++++++++++++++++++++++------------
R/CreateCalibOptions.R | 98 ++++++++++++++++++++++++++++++++-----
R/TransfoParam_LAG.R | 39 +++++++++++++++
4 files changed, 203 insertions(+), 42 deletions(-)
create mode 100644 R/TransfoParam_LAG.R
diff --git a/NAMESPACE b/NAMESPACE
index 17665070..858dd293 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -56,6 +56,7 @@ export(TransfoParam_GR5H)
export(TransfoParam_GR4J)
export(TransfoParam_GR5J)
export(TransfoParam_GR6J)
+export(TransfoParam_LAG)
export(plot.OutputsModel)
exportPattern(".FortranOutputs")
exportPattern(".ErrorCrit")
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 759630ea..dc813358 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -57,27 +57,31 @@ Calibration_Michel <- function(InputsModel,
##_check_FUN_TRANSFO
if (is.null(FUN_TRANSFO)) {
- if (identical(FUN_MOD, RunModel_GR4H )) {
- FUN_TRANSFO <- TransfoParam_GR4H
+ if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
+ FUN1 <- TransfoParam_GR4H
}
- if (identical(FUN_MOD, RunModel_GR5H )) {
- FUN_TRANSFO <- TransfoParam_GR5H
+ if (identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
+ FUN1 <- TransfoParam_GR5H
}
- if (identical(FUN_MOD, RunModel_GR4J )) {
- FUN_TRANSFO <- TransfoParam_GR4J
+ if (identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
+ FUN1 <- TransfoParam_GR4J
}
- if (identical(FUN_MOD, RunModel_GR5J )) {
- FUN_TRANSFO <- TransfoParam_GR5J
+ if (identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
+ FUN1 <- TransfoParam_GR5J
}
- if (identical(FUN_MOD, RunModel_GR6J )) {
- FUN_TRANSFO <- TransfoParam_GR6J
+ if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
+ FUN1 <- TransfoParam_GR6J
}
- if (identical(FUN_MOD, RunModel_GR2M )) {
- FUN_TRANSFO <- TransfoParam_GR2M
+ if (identical(FUN_MOD, RunModel_GR2M )) {
+ FUN1 <- TransfoParam_GR2M
}
- if (identical(FUN_MOD, RunModel_GR1A )) {
- FUN_TRANSFO <- TransfoParam_GR1A
+ if (identical(FUN_MOD, RunModel_GR1A )) {
+ FUN1 <- TransfoParam_GR1A
}
+ ##_set_FUN_LAG
+ if (inherits(CalibOptions, "SD")) {
+ FUN_LAG <- TransfoParam_LAG
+ }
if (identical(FUN_MOD, RunModel_CemaNeige )) {
if (inherits(CalibOptions, "hysteresis")) {
FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
@@ -85,29 +89,35 @@ Calibration_Michel <- function(InputsModel,
FUN_TRANSFO <- TransfoParam_CemaNeige
}
}
+ if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H) |
+ identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_GR6J)) {
+ if (!inherits(CalibOptions, "SD")) {
+ FUN_TRANSFO <- FUN1
+ } else {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ ParamOut[, 1:(NParam-1)] <- FUN1(ParamIn[, 1:(NParam-1)], Direction)
+ ParamOut[, NParam ] <- FUN_LAG(ParamIn[, NParam ], Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
+ }
if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H) |
identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
- FUN1 <- TransfoParam_GR4H
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
- FUN1 <- TransfoParam_GR5H
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
- FUN1 <- TransfoParam_GR4J
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
- FUN1 <- TransfoParam_GR5J
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- FUN1 <- TransfoParam_GR6J
- }
if (inherits(CalibOptions, "hysteresis")) {
FUN2 <- TransfoParam_CemaNeigeHyst
} else {
FUN2 <- TransfoParam_CemaNeige
}
- if (inherits(CalibOptions, "hysteresis")) {
+ if (inherits(CalibOptions, "hysteresis") & !inherits(CalibOptions, "SD")) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
@@ -122,7 +132,8 @@ Calibration_Michel <- function(InputsModel,
}
return(ParamOut)
}
- } else {
+ }
+ if (!inherits(CalibOptions, "hysteresis") & !inherits(CalibOptions, "SD")) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
@@ -138,6 +149,40 @@ Calibration_Michel <- function(InputsModel,
return(ParamOut)
}
}
+ if (inherits(CalibOptions, "hysteresis") & inherits(CalibOptions, "SD")) {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ ParamOut[, 1:(NParam-5)] <- FUN1( ParamIn[, 1:(NParam-5)], Direction)
+ ParamOut[, (NParam-4):(NParam-1)] <- FUN2( ParamIn[, (NParam-4):(NParam-1)], Direction)
+ ParamOut[, NParam ] <- FUN_LAG(ParamIn[, NParam ], Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
+ if (!inherits(CalibOptions, "hysteresis") & inherits(CalibOptions, "SD")) {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ ParamOut[, 1:(NParam-3)] <- FUN1( ParamIn[, 1:(NParam-3)], Direction)
+ ParamOut[, (NParam-2):(NParam-1)] <- FUN2( ParamIn[, (NParam-2):(NParam-1)], Direction)
+ ParamOut[, NParam ] <- FUN_LAG(ParamIn[, NParam ], Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
}
if (is.null(FUN_TRANSFO)) {
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 56dedd96..b3475e84 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -2,6 +2,7 @@ CreateCalibOptions <- function(FUN_MOD,
FUN_CALIB = Calibration_Michel,
FUN_TRANSFO = NULL,
IsHyst = FALSE,
+ IsSD = FALSE,
FixedParam = NULL,
SearchRanges = NULL,
StartParamList = NULL,
@@ -17,6 +18,9 @@ CreateCalibOptions <- function(FUN_MOD,
if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
stop("'IsHyst' must be a logical of length 1")
}
+ if (!is.logical(IsSD) | length(IsSD) != 1L) {
+ stop("'IsSD' must be a logical of length 1")
+ }
##check_FUN_MOD
BOOL <- FALSE
@@ -75,6 +79,9 @@ CreateCalibOptions <- function(FUN_MOD,
if (IsHyst) {
ObjectClass <- c(ObjectClass, "hysteresis")
}
+ if (IsSD) {
+ ObjectClass <- c(ObjectClass, "SD")
+ }
if (!BOOL) {
stop("incorrect 'FUN_MOD' for use in 'CreateCalibOptions'")
return(NULL)
@@ -139,11 +146,32 @@ CreateCalibOptions <- function(FUN_MOD,
} else {
FUN2 <- 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")) > 0) {
- FUN_TRANSFO <- FUN1
+ if (!IsSD) {
+ FUN_TRANSFO <- FUN1
+ } else {
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (!Bool) {
+ ParamIn <- rbind(ParamIn)
+ }
+ ParamOut <- NA * ParamIn
+ NParam <- ncol(ParamIn)
+ ParamOut[, 1:(NParam-1)] <- FUN1(ParamIn[, 1:(NParam-1)], Direction)
+ ParamOut[, NParam ] <- FUN_LAG(ParamIn[, NParam ], Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
} else {
- if (IsHyst) {
+ if (IsHyst & !IsSD) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
@@ -151,14 +179,15 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam - 4)] <- FUN1(ParamIn[, 1:(NParam - 4)], Direction)
+ ParamOut[, 1:(NParam - 4) ] <- FUN1(ParamIn[, 1:(NParam - 4) ], Direction)
ParamOut[, (NParam - 3):NParam] <- FUN2(ParamIn[, (NParam - 3):NParam], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
- } else {
+ }
+ if (!IsHyst & !IsSD) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
@@ -169,7 +198,7 @@ CreateCalibOptions <- function(FUN_MOD,
if (NParam <= 3) {
ParamOut[, 1:(NParam - 2)] <- FUN1(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
} else {
- ParamOut[, 1:(NParam - 2)] <- FUN1(ParamIn[, 1:(NParam - 2)], Direction)
+ ParamOut[, 1:(NParam - 2)] <- FUN1( ParamIn[, 1:(NParam - 2)], Direction)
}
ParamOut[, (NParam - 1):NParam] <- FUN2(ParamIn[, (NParam - 1):NParam], Direction)
if (!Bool) {
@@ -178,6 +207,44 @@ CreateCalibOptions <- function(FUN_MOD,
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[, 1:(NParam - 5) ] <- FUN1( ParamIn[, 1:(NParam - 5) ], Direction)
+ ParamOut[, (NParam - 4):(NParam - 1)] <- FUN2( ParamIn[, (NParam - 4):(NParam - 1)], Direction)
+ ParamOut[, NParam ] <- FUN_LAG(ParamIn[, 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 - 3)] <- FUN1(cbind(ParamIn[, 1:(NParam - 3)]), Direction)
+ } else {
+ ParamOut[, 1:(NParam - 3)] <- FUN1( ParamIn[, 1:(NParam - 3)], Direction)
+ }
+ ParamOut[, (NParam - 2):(NParam - 1)] <- FUN2( ParamIn[, (NParam - 2):(NParam - 1)], Direction)
+ ParamOut[, NParam ] <- FUN_LAG(ParamIn[, NParam], Direction)
+ if (!Bool) {
+ ParamOut <- ParamOut[1, ]
+ }
+ return(ParamOut)
+ }
+ }
}
}
if (is.null(FUN_TRANSFO)) {
@@ -228,6 +295,9 @@ CreateCalibOptions <- function(FUN_MOD,
if (IsHyst) {
NParam <- NParam + 2
}
+ if (IsSD) {
+ NParam <- NParam + 1
+ }
##check_FixedParam
if (is.null(FixedParam)) {
@@ -276,22 +346,22 @@ CreateCalibOptions <- function(FUN_MOD,
if ("GR4H" %in% ObjectClass) {
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)
+ +6.01, -0.50, +5.14, -8.87), ncol = 4, byrow = TRUE)
}
if (("GR5H" %in% ObjectClass) & ("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);
+ +4.29, +0.16, +5.39, -7.39, +3.33), ncol=5, byrow = TRUE);
}
if (("GR5H" %in% ObjectClass) & !("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);
+ +4.01, -0.04, +5.43, -7.53, -5.33), ncol=5, byrow = TRUE);
}
if ("GR4J" %in% ObjectClass) {
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)
+ +6.07, -0.02, +4.42, -8.06), ncol = 4, byrow = TRUE)
}
if ("GR5J" %in% ObjectClass) {
ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45,
@@ -329,12 +399,12 @@ CreateCalibOptions <- function(FUN_MOD,
if (("CemaNeigeGR5H" %in% ObjectClass) & ("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);
+ +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)) {
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);
+ +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE);
}
if ("CemaNeigeGR4J" %in% ObjectClass) {
ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
@@ -358,6 +428,12 @@ CreateCalibOptions <- function(FUN_MOD,
+7.00, +7.00), ncol = 2, byrow = TRUE)
ParamT <- cbind(ParamT, ParamTHyst)
}
+ if (IsSD) {
+ ParamTSD <- matrix(c(+1.25,
+ +2.50,
+ +5.00), ncol = 1, byrow = TRUE)
+ ParamT <- cbind(ParamT, ParamTSD)
+ }
StartParamList <- NULL
StartParamDistrib <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
diff --git a/R/TransfoParam_LAG.R b/R/TransfoParam_LAG.R
new file mode 100644
index 00000000..c484bf78
--- /dev/null
+++ b/R/TransfoParam_LAG.R
@@ -0,0 +1,39 @@
+TransfoParam_LAG <- function(ParamIn, Direction) {
+
+ ## number of model parameters
+ NParam <- 1L
+
+
+ ## check arguments
+ isVecParamIn <- is.vector(ParamIn)
+ if (isVecParamIn) {
+ ParamIn <- matrix(ParamIn, nrow = 1)
+ }
+ if (!inherits(ParamIn, "matrix")) {
+ stop("'ParamIn' must be of class 'matrix'")
+ }
+ if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
+ stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
+ }
+ if (ncol(ParamIn) != NParam) {
+ stop(sprintf("the LAG model requires %i parameters", NParam))
+ }
+
+
+ ## transformation
+ if (Direction == "TR") {
+ ParamOut <- 20 * (ParamIn + 10) / 20.0
+ }
+ if (Direction == "RT") {
+ ParamOut <- ParamIn * 20.0 / 20 - 10
+ }
+
+
+ ## check output
+ if (isVecParamIn) {
+ ParamOut <- as.vector(ParamOut)
+ }
+
+ return(ParamOut)
+
+}
--
GitLab