diff --git a/R/Calibration.InputsModel.R b/R/Calibration.InputsModel.R index 3bc19d76b4f0180fd79c69e812ed79bc9b5bc92e..c05ad7c9a0674a8471f5205fb263d45b76313110 100644 --- a/R/Calibration.InputsModel.R +++ b/R/Calibration.InputsModel.R @@ -11,9 +11,15 @@ Calibration.InputsModel <- function(InputsModel, class(OC) <- c("OutputsCalib", class(OC)) return(OC) } else { - airGR::Calibration(InputsModel, - CalibOptions = CalibOptions, - FUN_MOD = InputsModel$FUN_MOD, ...) + if (!is.null(InputsModel$hasDiversion) && InputsModel$hasDiversion) { + class(InputsModel) <- setdiff(class(InputsModel), "SD") + FUN_MOD = RunModel.InputsModel + } else { + FUN_MOD = InputsModel$FUN_MOD + } + airGR::Calibration(InputsModel, + CalibOptions = CalibOptions, + FUN_MOD = FUN_MOD, ...) } } else { airGR::Calibration(InputsModel, diff --git a/R/RunModel.InputsModel.R b/R/RunModel.InputsModel.R index eddeb513921474741018ee8f74566293e32c76c6..f36e4b0e4a08f524678c1c1107cdba6b2aa888fb 100644 --- a/R/RunModel.InputsModel.R +++ b/R/RunModel.InputsModel.R @@ -11,7 +11,24 @@ #' @inherit airGR::RunModel return return #' #' @export -RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) { +RunModel.InputsModel <- function(x = NULL, + RunOptions, + Param, + FUN_MOD = NULL, + InputsModel = NULL, + ...) { + if (is.null(x)) { + if (!is.null(InputsModel)) { + x <- InputsModel + } else { + stop("`x` or `InputsModel` must be defined") + } + } else { + if (!is.null(InputsModel)) { + stop("`x` and `InputsModel` can't be defined at the same time") + } + } + if(is.null(FUN_MOD)) { if (x$isReservoir) { FUN_MOD <- "RunModel_Reservoir" @@ -19,6 +36,7 @@ RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) { FUN_MOD <- x$FUN_MOD } } + FUN_MOD <- match.fun(FUN_MOD) if (identical(FUN_MOD, RunModel_Lag)) { QcontribDown <- list( @@ -31,11 +49,15 @@ RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) { x$BasinAreas[length(x$BasinAreas)] <- 1 OutputsModel <- RunModel_Lag(x, RunOptions, Param, QcontribDown) OutputsModel$DatesR <- x$DatesR[RunOptions$IndPeriod_Run] - } else if((inherits(x, "GR") & !inherits(x, "SD")) | identical(FUN_MOD, RunModel_Reservoir)) { + } else if((inherits(x, "GR") & is.null(x$UpstreamNodes)) | identical(FUN_MOD, RunModel_Reservoir)) { # Upstream basins and Reservoir are launch directly OutputsModel <- FUN_MOD(x, RunOptions, Param) } else { # Intermediate basins (other than reservoir) are launch with SD capabilities + if (!is.null(x$UpstreamNodes) & !inherits(x, "SD")) { + # For calibration of node with diversion + class(x) <- c(class(x), "SD") + } OutputsModel <- airGR::RunModel(x, RunOptions, Param, FUN_MOD) } OutputsModel$RunOptions$TimeStep <- RunOptions$FeatFUN_MOD$TimeStep @@ -54,7 +76,6 @@ RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) { return(OutputsModel) } - #' Model the diversion of a flow from an existing modeled node #' #' On a Diversion node, this function is called after `airGR::RunModel` to diff --git a/man/RunModel.InputsModel.Rd b/man/RunModel.InputsModel.Rd index 7f9e11c71343aae58f065f2feeeb6db7055c0ee8..7b4887d9aadddbbe3a6a3e2415c8de1f00903461 100644 --- a/man/RunModel.InputsModel.Rd +++ b/man/RunModel.InputsModel.Rd @@ -4,7 +4,7 @@ \alias{RunModel.InputsModel} \title{Wrapper for \link[airGR:RunModel]{airGR::RunModel} for one sub-basin} \usage{ -\method{RunModel}{InputsModel}(x, RunOptions, Param, FUN_MOD = NULL, ...) +\method{RunModel}{InputsModel}(x = NULL, RunOptions, Param, FUN_MOD = NULL, InputsModel = NULL, ...) } \arguments{ \item{x}{[object of class \emph{InputsModel}] see \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} for details} @@ -15,6 +15,8 @@ \item{FUN_MOD}{[function] hydrological model function (e.g. \code{\link[airGR]{RunModel_GR4J}}, \code{\link[airGR]{RunModel_CemaNeigeGR4J}})} +\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link[airGR]{CreateInputsModel}} for details} + \item{...}{Further arguments for compatibility with S3 methods} } \value{ diff --git a/tests/testthat/test-Calibration.R b/tests/testthat/test-Calibration.R index 8ff1ccc36d38d22aa65b353dcdb3a0da720e2948..4c86ae251adf00d397842c2cd4ed3cd1678c52e5 100644 --- a/tests/testthat/test-Calibration.R +++ b/tests/testthat/test-Calibration.R @@ -42,19 +42,20 @@ for(x in ls(e)) assign(x, get(x, e)) CalibOptions <- CreateCalibOptions(InputsModel) -test_that("Calibrated parameters remains unchanged", { - InputsCrit <- CreateInputsCrit( - InputsModel = InputsModel, - RunOptions = RunOptions, - Obs = Qobs[IndPeriod_Run,] - ) +InputsCrit <- CreateInputsCrit( + InputsModel = InputsModel, + RunOptions = RunOptions, + Obs = Qobs[IndPeriod_Run,] +) + +OC <- Calibration( + InputsModel = InputsModel, + RunOptions = RunOptions, + InputsCrit = InputsCrit, + CalibOptions = CalibOptions +) - OC <- Calibration( - InputsModel = InputsModel, - RunOptions = RunOptions, - InputsCrit = InputsCrit, - CalibOptions = CalibOptions - ) +test_that("Calibrated parameters remains unchanged", { ParamFinalR <- lapply(OC, "[[", "ParamFinalR") @@ -99,6 +100,8 @@ test_that("Calibration with regularization is OK", { }) }) +skip_on_cran() + test_that("Calibration with Diversion works", { n_div <- rbind(nodes, data.frame(id = "54029", down = "54002", length = 50, area = NA, model = "Diversion")) @@ -124,3 +127,63 @@ test_that("Calibration with Diversion works", { ) expect_length(OC$`54002`$ParamFinalR, 5) }) + +test_that("Derivation and normal connection should return same calibration", { + n_2ol <- nodes[nodes$id %in% c("54095", "54001"), ] + n_2ol[n_2ol$id %in% c("54095", "54001"), c("down", "length")] <- c(NA, NA) + meteoIds <- n_2ol$id + n_2ol$area[n_2ol$id == "54001"] <- + n_2ol$area[n_2ol$id == "54001"] - n_2ol$area[n_2ol$id == "54095"] + n_2ol <- rbind(n_2ol, + data.frame(id = "54095", down = "54001", length = 42, area = NA, model = "Diversion"), + data.frame(id = "upstream", down = "54095", length = 0, area = NA, model = NA)) + g_2ol <- CreateGRiwrm(n_2ol) + + # Add upstream flow on 54095 that is removed by the Diversion + # and derive previously simulated flow in order to get the same Qsim as before + Qinf = matrix(0, nrow = length(DatesR), ncol = 2) + Qinf[IndPeriod_Run, 1] <- OM_GriwrmInputs[["54095"]]$Qsim_m3 + Qinf[IndPeriod_Run, 2] <- - OM_GriwrmInputs[["54095"]]$Qsim_m3 + Qinf[IndPeriod_WarmUp, 1] <- OM_GriwrmInputs[["54095"]]$RunOptions$WarmUpQsim_m3 + Qinf[IndPeriod_WarmUp, 2] <- - OM_GriwrmInputs[["54095"]]$RunOptions$WarmUpQsim_m3 + + colnames(Qinf) <- c("upstream", "54095") + + Qmin = matrix(0, nrow = length(DatesR), ncol = 1) + colnames(Qmin) <- "54095" + + IM_2ol <- CreateInputsModel(g_2ol, + DatesR, + Precip[, meteoIds], + PotEvap[, meteoIds], + Qobs = Qinf, + Qmin = Qmin) + + # Copy area of upstream node to downstream node in order to get + # correct conversion of Qsim in mm + IM_2ol[["54001"]]$BasinAreas[1] <- tail(IM_2ol[["54095"]]$BasinAreas, 1) + + RO_2ol <- setupRunOptions(IM_2ol)$RunOptions + IC_2ol <- CreateInputsCrit( + InputsModel = IM_2ol, + RunOptions = RO_2ol, + Obs = Qobs[IndPeriod_Run,], + ) + CO_2ol <- CreateCalibOptions(IM_2ol) + CO_2ol[["54095"]]$FixedParam[1] <- 1 + OC_2ol <- Calibration( + InputsModel = IM_2ol, + RunOptions = RO_2ol, + InputsCrit = IC_2ol, + CalibOptions = CO_2ol + ) + ParamRef <- ParamMichel[names(IM_2ol)] + ParamRef[["54095"]] <- c(1, ParamRef[["54095"]]) + ParamFinalR <- lapply(OC_2ol, "[[", "ParamFinalR") + lapply(names(ParamFinalR), function(id) expect_equal(OC_2ol[[id]]$CritFinal, + OC[[id]]$CritFinal, + tolerance = 1E-5)) + #Excepted parameter #2 of GR4J all others are equal (precision 3/1000) + lapply(names(ParamFinalR), function(id) + expect_equal(ParamFinalR[[!!id]][-3] / ParamRef[[!!id]][-3], rep(1, 4), tolerance = 3E-3)) +}) diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R index 344a5a88c2e72cc446cd1f3d50f7d99c3099d252..611caec457960c86c2967db8ffeaeac093da49bf 100644 --- a/tests/testthat/test-RunModel.R +++ b/tests/testthat/test-RunModel.R @@ -150,3 +150,76 @@ test_that("RunModel_Lag should work", { expect_s3_class(OM, "GRiwrmOutputsModel") expect_true(all(!is.na(attr(OM, "Qm3s")))) }) + +test_that("Upstream node - equal Diversion should return same results", { + n_2ol <- nodes[nodes$id %in% c("54095", "54001"), ] + n_2ol[n_2ol$id %in% c("54095", "54001"), c("down", "length")] <- c(NA, NA) + meteoIds <- n_2ol$id + n_2ol$area[n_2ol$id == "54001"] <- + n_2ol$area[n_2ol$id == "54001"] - n_2ol$area[n_2ol$id == "54095"] + n_2ol <- rbind(n_2ol, + data.frame(id = "54095", down = "54001", length = 42, area = NA, model = "Diversion"), + data.frame(id = "upstream", down = "54095", length = 0, area = NA, model = NA)) + g_2ol <- CreateGRiwrm(n_2ol) + + # Add upstream flow on 54095 that is removed by the Diversion + # and derive previously simulated flow in order to get the same Qsim as before + Qinf = matrix(0, nrow = length(DatesR), ncol = 2) + Qinf[IndPeriod_Run, 1] <- OM_GriwrmInputs[["54095"]]$Qsim_m3 + Qinf[IndPeriod_Run, 2] <- - OM_GriwrmInputs[["54095"]]$Qsim_m3 + Qinf[IndPeriod_WarmUp, 1] <- OM_GriwrmInputs[["54095"]]$RunOptions$WarmUpQsim_m3 + Qinf[IndPeriod_WarmUp, 2] <- - OM_GriwrmInputs[["54095"]]$RunOptions$WarmUpQsim_m3 + + colnames(Qinf) <- c("upstream", "54095") + + Qmin = matrix(0, nrow = length(DatesR), ncol = 1) + colnames(Qmin) <- "54095" + + IM_2ol <- CreateInputsModel(g_2ol, + DatesR, + Precip[, meteoIds], + PotEvap[, meteoIds], + Qobs = Qinf, + Qmin = Qmin) + + RO_2ol <- setupRunOptions(IM_2ol)$RunOptions + P_2ol <- ParamMichel[names(IM_2ol)] + P_2ol[["54095"]] <- c(1, P_2ol[["54095"]]) + OM_2ol <- RunModel(IM_2ol, RO_2ol, P_2ol) + + # Is the diversion correctly taken into account? + expect_equal( + OM_2ol[["54095"]]$Qdiv_m3, + OM_GriwrmInputs[["54095"]]$Qsim_m3 + ) + + # Is 54001 InputsModel correctly updated? + IM_54001_div <- UpdateQsimUpstream(IM_2ol[["54001"]], + RO_2ol[["54001"]], + OM_2ol) + IM_54001 <- UpdateQsimUpstream(InputsModel[["54001"]], + RunOptions[["54001"]], + OM_GriwrmInputs) + expect_equal( + IM_54001_div$Qupstream, + IM_54001$Qupstream + ) + + # All simulated flows with or without div must be equal + sapply(names(IM_2ol), function(id) { + expect_equal(OM_2ol[[!!id]]$Qsimdown, OM_GriwrmInputs[[!!id]]$Qsimdown) + expect_equal(OM_2ol[[!!id]]$Qsim_m3, OM_GriwrmInputs[[!!id]]$Qsim_m3) + }) + + id <- "54095" + IM_54095_div <- IM_2ol[[id]] + class(IM_54095_div) <- setdiff(class(IM_54095_div), "SD") + + OM_airGR <- airGR::RunModel(IM_54095_div, + RunOptions = RO_2ol[[id]], + Param = P_2ol[[id]], + FUN_MOD = RunModel.InputsModel) + + expect_equal(OM_airGR$Qsimdown, OM_GriwrmInputs[[!!id]]$Qsimdown) + expect_equal(OM_airGR$Qsim_m3, OM_GriwrmInputs[[!!id]]$Qsim_m3) +})