feat(RunModel_Lag): handling of warm-up simulated flows

- Implementation and tests

Refs #123

R/RunModel_Lag.R

@@ -69,9 +69,21 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
       IniSD[iStart:(iStart + PT[x])]
     })
   } else {
-    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
-      rep(0, floor(PT[x] + 1))
-    })
+    IniStates <- lapply(
+      seq_len(NbUpBasins),
+      function(iUpBasins) {
+        iWarmUp <- seq.int(
+          max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
+          max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)])
+        )
+        ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
+        if(length(ini) != floor(PT[iUpBasins] + 1)) {
+          # If warm-up period is not enough long complete beginning with first value
+          ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
+        }
+        return(ini)
+      }
+    )
   }
   # message("Initstates: ", paste(IniStates, collapse = ", "))
 

tests/testthat/test-RunModel_Lag.R

@@ -165,7 +165,7 @@ Qm3GR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e3
 test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
   OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
   Qm3SdSim <- OutputsSD$Qsim_m3
-  Qm3UpstLagObs <- c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * InputsModel$BasinAreas[1] * 1e3
+  Qm3UpstLagObs <- Qupstream[Ind_Run - 1] * InputsModel$BasinAreas[1] * 1e3
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })
 
@@ -174,7 +174,7 @@ test_that("1 input with lag of 0.5 time step delay out gives an output delayed o
                         Param = c(InputsModel$LengthHydro * 1e3 / (12 * 3600), Param),
                         FUN_MOD = RunModel_GR4J)
   Qm3SdSim <- OutputsSD$Qsim_m3
-  Qm3UpstLagObs <- (Qupstream[Ind_Run] + c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])) / 2 * InputsModel$BasinAreas[1] * 1e3
+  Qm3UpstLagObs <- (Qupstream[Ind_Run] + Qupstream[Ind_Run - 1]) / 2 * InputsModel$BasinAreas[1] * 1e3
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })
 
@@ -233,7 +233,7 @@ InputsCrit <- CreateInputsCrit(
   InputsModel = InputsModel,
   RunOptions = RunOptions,
   VarObs = "Q",
-  Obs = (c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
+  Obs = (Qupstream[Ind_Run - 1] * BasinAreas[1L] +
            BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas)
 )
 
@@ -283,7 +283,7 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
   expect_false(any(is.na(OutputsSD$Qsim)))
 
   Qm3SdSim <- OutputsSD$Qsim_m3
-  Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
+  Qm3UpstLagObs <- InputsModel$Qupstream[Ind_Run - 1]
 
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })
@@ -338,6 +338,21 @@ test_that("Error on Wrong length of iniState$SD", {
                                   InputsModel = IM, IndPeriod_Run = Ind_Run2,
                                   IndPeriod_WarmUp = 0L,
                                   IniStates = OutputsModel1$StateEnd)
-  expect_error(RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
+  expect_error(
+    RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
   )
 })
+
+test_that("First Qupstream time steps must be repeated if warm-up period is too short", {
+  IM2 <- IM[2558:3557]
+  IM2$BasinAreas[3] <- 0
+  IM2$Qupstream <- matrix(rep(1:1000, 2), ncol = 2)
+  RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                  InputsModel = IM2, IndPeriod_Run = seq_len(1000),
+                                  IndPeriod_WarmUp = 0L)
+  OM2 <- RunModel(InputsModel = IM2,
+                  RunOptions = RunOptions2,
+                  Param = PSDini,
+                  FUN_MOD = RunModel_GR4J)
+  expect_equal(OM2$Qsim_m3[1:3], rep(2,3))
+})