fix: test with 2 controllers that cancel each other out works!

- fix issue with IniStates in RunModel_Lag (See HYCAR-Hydro/airgr#102 and HYCAR-Hydro/airgr#104)
- add GRiwrm attributes in GRiwrmInputsModel object and griwrm variable in Supervisor environment
- debug RunModel.Supervisor, doSupervision and  functions called by doSupervision

Refs #19

R/CreateInputsModel.GRiwrm.R

@@ -11,8 +11,8 @@
 #' @export
 CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
 
-  InputsModel <- CreateEmptyGRiwrmInputsModel()
-  Qobs[is.na(Qobs)] <- -99 # airGRCreateInputsModel doesn't accept NA values
+  InputsModel <- CreateEmptyGRiwrmInputsModel(x)
+  Qobs[is.na(Qobs)] <- -99 # airGR::CreateInputsModel doesn't accept NA values
 
   for(id in getNodeRanking(x)) {
     message("CreateInputsModel.GRiwrm: Treating sub-basin ", id, "...")
@@ -27,9 +27,10 @@ CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
 #' Create an empty InputsModel object for **airGRiwrm** nodes
 #'
 #' @return \emph{GRiwrmInputsModel} empty object
-CreateEmptyGRiwrmInputsModel <- function() {
+CreateEmptyGRiwrmInputsModel <- function(griwrm) {
   InputsModel <- list()
   class(InputsModel) <- c("GRiwrmInputsModel", class(InputsModel))
+  attr(InputsModel, "GRiwrm") <- griwrm
   return(InputsModel)
 }
 

R/CreateSupervisor.R

@@ -25,7 +25,7 @@
 #' sv <- CreateSupervisor(InputsModel)
 CreateSupervisor <- function(InputsModel) {
   if(!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateInputsModel.GRiwrm)")
-  # Create Supervisor environment in the parent of GlobalEnv
+  # Create Supervisor environment from the parent of GlobalEnv
   e <- new.env(parent = parent.env(globalenv()))
   class(e) <- c("Supervisor", class(e))
 
@@ -35,8 +35,9 @@ CreateSupervisor <- function(InputsModel) {
   # Add pointer to itself in order to assign variable from function environment
   e$supervisor <- e
 
-  # Copy of the InputsModel
+  # Copy of InputsModel, griwrm and prepare OutputsModel
   e$InputsModel <- InputsModel
+  e$griwrm <- attr(InputsModel, "GRiwrm")
   e$OutputsModel <- list()
 
   # Controller list

R/RunModel.Supervisor.R

@@ -9,9 +9,7 @@
 #' @export
 RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
 
-  x$ts.index0 <- sapply(RunOptions, function(x) {
-    x$IndPeriod_Run[1] - 1
-  })
+  x$ts.index0 <- RunOptions[[1]]$IndPeriod_Run[1] - 1
 
   # Run runoff model for each sub-basin
   x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
@@ -20,6 +18,13 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
                 Param = Param[[IM$id]])
     })
   class(x$OutputsModel) <- append(class(x$OutputsModel), "GRiwrmOutputsModel")
+  # Copy simulated pure runoff flows (no SD nodes) to Qupstream in downstream SD nodes
+  for(id in getNoSD_Ids(x$InputsModel)) {
+    downId <- x$InputsModel[[id]]$down
+    x$InputsModel[[downId]]$Qupstream[RunOptions[[downId]]$IndPeriod_Run, id] <-
+      x$OutputsModel[[id]]$Qsim
+  }
+
   # Save Qsim for step by step simulation
   Qsim <- lapply(x$OutputsModel, function(OM) {
     OM$Qsim
@@ -34,19 +39,17 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   # Loop over time steps
   for(iTS in RunOptions[[1]]$IndPeriod_Run) {
     # Run regulation on the whole basin for the current time step
-    x$ts.index <- iTS
+    x$ts.index <- iTS - x$ts.index0
     x$ts.date <- x$InputsModel[[1]]$DatesR[iTS]
-    doSupervision(x)
+    # Regulation occurs from second time step
+    if(iTS > RunOptions[[1]]$IndPeriod_Run[1]) {
+      doSupervision(x)
+      message("Supervision done")
+    }
 
     # Loop over sub-basin using SD model
     for(id in getSD_Ids(x$InputsModel)) {
 
-      # Update InputsModel$Qupstream with simulated upstream flows
-      for(i in which(x$InputsModel[[id]]$UpstreamIsRunoff)) {
-        x$InputsModel[[id]]$Qupstream[iTS, i] <-
-          x$OutputsModel[[x$InputsModel[[id]]$UpstreamNodes[i]]]$Qsim[iTS - x$ts.index0]
-      }
-
       # Run the SD model for the sub-basin and one time step
       RunOptions[[id]]$IndPeriod_Run <- iTS
       RunOptions[[id]]$IniStates <- unlist(x$OutputsModel[[id]]$StateEnd)
@@ -54,9 +57,15 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
         x$InputsModel[[id]],
         RunOptions = RunOptions[[id]],
         Param = Param[[id]],
-        QsimDown = Qsim[[id]][iTS - x$ts.index0]
+        QsimDown = Qsim[[id]][x$ts.index]
       )
-      Qsim[[id]][iTS - x$ts.index0] <- x$OutputsModel[[id]]$Qsim
+      # Storing Qsim in the data.frame Qsim
+      Qsim[[id]][x$ts.index] <- x$OutputsModel[[id]]$Qsim
+      # Routing Qsim to the downstream node
+      if(!is.na(x$InputsModel[[id]]$down)) {
+        x$InputsModel[[x$InputsModel[[id]]$down]]$Qupstream[iTS, i] <-
+          x$OutputsModel[[id]]$Qsim
+      }
     }
   }
   for(id in getSD_Ids(x$InputsModel)) {

R/RunModel_Lag.R

@@ -78,20 +78,17 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
   #message("Initstates: ",paste(IniStates, collapse = ", "))
 
   for (upstream_basin in seq_len(NbUpBasins)) {
-    Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
+    Qupstream <- c(IniStates[[upstream_basin]],
+                   InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin])
     if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
       # Upstream flow with area needs to be converted to m3 by time step
       Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
     }
-    Qupstream1 <- c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])], Qupstream[1:(LengthTs - floor(PT[upstream_basin]))])
-    Qupstream2 <- IniStates[[upstream_basin]]
-    if(LengthTs - floor(PT[upstream_basin]) - 1 > 0) Qupstream2 <- c(Qupstream2, Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)])
-    #message("Qupstream1: ", paste(Qupstream1, collapse = ", "))
-    #message("Qupstream2: ", paste(Qupstream2, collapse = ", "))
-
-    OutputsModel$Qsim <- OutputsModel$Qsim +
-                         Qupstream1 * HUTRANS[1, upstream_basin] +
-                         Qupstream2 * HUTRANS[2, upstream_basin]
+    #message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
+    OutputsModel$Qsim <-
+      OutputsModel$Qsim +
+      Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
+      Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
   }
   # Warning for negative flows
   if (any(OutputsModel$Qsim < 0)) {
@@ -105,12 +102,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
   if ("StateEnd" %in% RunOptions$Outputs_Sim) {
     OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
       LengthTs <- tail(RunOptions$IndPeriod_Run,1)
-      Qupstream <- InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
-      if (!is.na(InputsModel$BasinAreas[x])) {
-        # Upstream flow with area needs to be converted to m3 by time step
-        Qupstream <- Qupstream * InputsModel$BasinAreas[x] * 1e3
-      }
-      return(Qupstream)
+      InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
     })
     #message("StateEnd: ",paste(OutputsModel$StateEnd$SD, collapse = ", "))
   }

R/utils.R

@@ -14,35 +14,53 @@ getSD_Ids <- function(InputsModel) {
   names(InputsModel)[bSDs]
 }
 
+#' Id of sub-basins not using SD model
+#'
+#' @param InputsModel `GRiwrmInputsModel` object
+#'
+#' @return [character] IDs of the sub-basins not using SD model
+#'
+getNoSD_Ids <- function(InputsModel) {
+  if (!inherits(InputsModel, "GRiwrmInputsModel")) {
+    stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
+  }
+  bSDs <- sapply(InputsModel, function (IM) {
+    !inherits(IM, "SD")
+  })
+  names(InputsModel)[bSDs]
+}
+
 
 #' Retrieve data in the model for the current time steps
 #'
 #' This function should be call inside a Supervisor
 #'
 #' @param loc location of the data
-#' @param supervisor `Supervisor` (See [CreateSupervisor])
+#' @param sv a `Supervisor` (See [CreateSupervisor])
 #'
 #' @return [numeric] retrieved data at the location
-getDataFromLocation <- function(loc, supervisor) {
-  if (grep("\\[[0-9]+\\]$", loc)) {
+getDataFromLocation <- function(loc, sv) {
+  if (length(grep("\\[[0-9]+\\]$", loc)) > 0) {
     stop("Reaching output of other controller is not implemented yet")
   } else {
-    supervisor$OutputsModel[[loc]]$Qsim[supervisor$ts.index - 1]
+    sv$InputsModel[[sv$griwrm$down[sv$griwrm$id == loc]]]$Qupstream[sv$ts.index0 + sv$ts.index - 1, loc]
   }
 }
 
 
 #' Write data to model input for the current time step
 #'
-#' @param control [list] A row of the `U` [data.frame] from a `Controller`
-#' @param supervisor `Supervisor` (See [CreateSupervisor])
+#' @param control [vector] A row of the `U` [data.frame] from a `Controller`
+#' @param sv `Supervisor` (See [CreateSupervisor])
 #'
 #' @return [NULL]
-setDataToLocation <- function(control, supervisor) {
-  node <- supervisor$InputsModel[[control$loc]]$down
+setDataToLocation <- function(control, sv) {
+  message("setDataToLocation[", control[1], "] <- ", control[2])
+  node <- sv$griwrm$down[sv$griwrm$id == control[1]]
   # ! Qupstream contains warm up period and run period => the index is shifted
-  supervisor$InputsModel[[node]]$Qupstream[supervisor$ts.index0[node] + supervisor$ts.index, control$loc] <-
-    control$v
+  sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index, control[1]] <-
+    as.numeric(control[2])
+  message("setDataToLocation[", control[1], "] <- ", control[2])
 }
 
 
@@ -55,12 +73,12 @@ doSupervision <- function(supervisor) {
   for (id in names(supervisor$controllers)) {
     # Read Y from locations in the model
     supervisor$controllers[[id]]$Y$v <-
-      sapply(supervisor$controllers[[id]]$Y$loc, getDataFromLocation, supervisor = supervisor)
+      sapply(supervisor$controllers[[id]]$Y$loc, getDataFromLocation, sv = supervisor)
     # Run logic
     supervisor$controllers[[id]]$U$v <-
       sapply(supervisor$controllers[[id]]$Y$v, supervisor$controllers[[id]]$FUN)
     # Write U to locations in the model
-    sapply(supervisor$controllers[[id]]$U, setDataToLocation, supervisor = supervisor)
+    apply(supervisor$controllers[[id]]$U, 1, setDataToLocation, sv = supervisor)
   }
   return()
 }

tests/testthat/test-RunModel.R

 context("RunModel.Supervisor")
 
+# Load data
+data(Severn)
+
+# Network configuration
+nodes <- Severn$BasinsInfo[c(1,2,5), c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes$distance_downstream <- nodes$distance_downstream * 1000 # Conversion km -> m
+nodes$model <- NA
+nodes$model[1] <- "RunModel_GR4J"
+griwrm <- GRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+
+# InputsModel
+DatesR <- Severn$BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$peti}))
+Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
+InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+
+# RunOptions
+nTS <- 365
+IndPeriod_Run <- seq(
+  length(InputsModel[[1]]$DatesR) - nTS + 1,
+  length(InputsModel[[1]]$DatesR)
+)
+IndPeriod_WarmUp = seq(IndPeriod_Run[1]-366,IndPeriod_Run[1]-1)
+RunOptions <- CreateRunOptions(
+  InputsModel = InputsModel,
+  IndPeriod_WarmUp = IndPeriod_WarmUp,
+  IndPeriod_Run = IndPeriod_Run
+)
+
+# RunModel.GRiwrmInputsModel
+Param <- list("54057" = c(0.727,  175.493,   -0.082,    0.029,    4.654))
+OM_GriwrmInputs <- RunModel(
+  InputsModel,
+  RunOptions = RunOptions,
+  Param = Param
+)
+
 test_that("RunModelSupervisor with no regulation should returns same results as RunModel.GRiwrmInputsModel", {
-  data(Severn)
-  # Network configuration
-  nodes <- Severn$BasinsInfo[c(1,2,5), c("gauge_id", "downstream_id", "distance_downstream", "area")]
-  nodes$distance_downstream <- nodes$distance_downstream * 1000 # Conversion km -> m
-  nodes$model <- NA
-  nodes$model[1] <- "RunModel_GR4J"
-  griwrm <- GRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
-  # InputsModel
-  DatesR <- Severn$BasinsObs[[1]]$DatesR
-  PrecipTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$precipitation}))
-  PotEvapTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$peti}))
-  Precip <- ConvertMeteoSD(griwrm, PrecipTot)
-  PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
-  Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
-  InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
-  # RunOptions
-  nTS <- 365
-  IndPeriod_Run <- seq(
-    length(InputsModel[[1]]$DatesR) - nTS + 1,
-    length(InputsModel[[1]]$DatesR)
-  )
-  IndPeriod_WarmUp = seq(IndPeriod_Run[1]-366,IndPeriod_Run[1]-1)
-  RunOptions <- CreateRunOptions(
-    InputsModel = InputsModel,
-    IndPeriod_WarmUp = IndPeriod_WarmUp,
-    IndPeriod_Run = IndPeriod_Run
-  )
-  Param <- list("54057" = c(0.727,  175.493,   -0.082,    0.029,    4.654))
-  OM_GriwrmInputs <- RunModel(
-    InputsModel,
+  sv <- CreateSupervisor(InputsModel)
+  OM_Supervisor <- RunModel(
+    sv,
     RunOptions = RunOptions,
     Param = Param
   )
-  supervisor <- CreateSupervisor(InputsModel)
+  expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})
+
+test_that("RunModelSupervisor with two regulations that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
+  # Add 2 nodes to the network
+  griwrm2 <- rbind(griwrm,
+                  data.frame(
+                    id = c("R1", "R2"),
+                    down = "54057",
+                    length = 100000,
+                    area = NA,
+                    model = NA
+                  ))
+  # Add Qobs for the 2 new nodes
+  Qobs2 <- cbind(Qobs, matrix(data = rep(0, 2*nrow(Qobs)), ncol = 2))
+  colnames(Qobs2) <- c(colnames(Qobs2)[1:6], "R1", "R2")
+  InputsModel <- CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qobs2)
+  sv <- CreateSupervisor(InputsModel)
+  # Function to withdraw half of the measured flow
+  fWithdrawal <- function(y) { message("fWithdrawal(", y, ")"); -y/2 }
+  # Function to release half of the the measured flow
+  fRelease <- function(y) { message("fWithdrawal(", y, ")"); y/2 }
+  # Controller that withdraw half of the flow measured at node "54002" at location "R1"
+  createController(sv, "Withdrawal", Y = c("54002"), U = c("R1"), FUN = fWithdrawal)
+  # Controller that release half of the flow measured at node "54002" at location "R2"
+  createController(sv, "Release", Y = c("54002"), U = c("R2"), FUN = fRelease)
   OM_Supervisor <- RunModel(
-    supervisor,
+    sv,
     RunOptions = RunOptions,
     Param = Param
   )
   expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
 })
+