fix(calibration): correct handle of node with diversion

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,

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

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{

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))
+})

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)
+})