feat: Add initial states handling on RunModel_Lag

Refs #78

R/RunModel_Lag.R

 RunModel_Lag <- function(InputsModel, RunOptions, Param) {
 
   NParam <- 1
-  
+
   ##Arguments_check
   if (!inherits(InputsModel, "InputsModel")) {
     stop("'InputsModel' must be of class 'InputsModel'")
-  }  
+  }
   if (!inherits(InputsModel, "SD")) {
     stop("'InputsModel' must be of class 'SD'")
-  }  
+  }
   if (!inherits(RunOptions, "RunOptions")) {
     stop("'RunOptions' must be of class 'RunOptions'")
   }
@@ -30,7 +30,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
   
   OutputsModel <- InputsModel$OutputsModel
   OutputsModel$QsimDown <- OutputsModel$Qsim
-  
+
   if (inherits(InputsModel, "daily")) {
     TimeStep <- 60 * 60 * 24
   }
@@ -45,7 +45,21 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
   NbUpBasins <- length(InputsModel$LengthHydro)
   LengthTs <- length(OutputsModel$QsimDown)
   OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
-
+  
+  IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
+  if(length(IniSD) > 0) {
+    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+      iStart <- 1
+      if(x > 1) {
+        iStart <- iStart + sum(floor(PT[1:x-1]) + 1)
+      }
+      IniSD[iStart:(iStart + PT[x])]
+    })
+  } else {
+    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+      rep(0, floor(PT[x] + 1))})
+  }
+  
   for (upstream_basin in seq_len(NbUpBasins)) {
     Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
     if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
@@ -53,10 +67,10 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
       Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
     }
     OutputsModel$Qsim <- OutputsModel$Qsim +
-      c(rep(0, floor(PT[upstream_basin])),
+      c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])],
         Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
       HUTRANS[1, upstream_basin] +
-      c(rep(0, floor(PT[upstream_basin] + 1)),
+      c(IniStates[[upstream_basin]],
         Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
       HUTRANS[2, upstream_basin]
   }
@@ -67,5 +81,11 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
   }
   # Convert back Qsim to mm
   OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+
+  if("StateEnd" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
+      Qupstream[(LengthTs - floor(PT[upstream_basin])):LengthTs]})
+  }
+
   return(OutputsModel)
 }
\ No newline at end of file

tests/testthat/test-RunModel_LAG.R

@@ -50,9 +50,9 @@ InputsModel <- CreateInputsModel(
 Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
                which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
 
-RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                InputsModel = InputsModel,
-                               IndPeriod_Run = Ind_Run)
+                               IndPeriod_Run = Ind_Run))
 
 test_that("InputsModel parameter should contain an OutputsModel key", {
   expect_error(
@@ -178,3 +178,35 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
   
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })
+
+test_that("Warm start should give same result as warmed model", {
+  InputsModel$BasinAreas <- rep(BasinInfo$BasinArea, 3)
+  InputsModel$Qupstream <- cbind(InputsModel$Qupstream, InputsModel$Qupstream)
+  InputsModel$LengthHydro <- c(1000, 1500)
+  
+  ParamSD[1] <- ParamSD[1] / 2 # 2 time step delay
+  Ind_Run1 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"), 
+                  which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-12-31"))
+  Ind_Run2 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-01-01"), 
+                  which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-12-31"))
+  # 1990-1991
+  RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                                  InputsModel = InputsModel, IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
+  OutputsModel <- RunModel(InputsModel = InputsModel,
+                                RunOptions = RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
+  # 1990
+  RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Run1))
+  OutputsModel1 <- RunModel(InputsModel = InputsModel,
+                                 RunOptions = RunOptions1, Param = ParamSD, FUN_MOD = RunModel_GR4J)
+  # Warm start 1991
+  RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                  InputsModel = InputsModel, IndPeriod_Run = Ind_Run2, 
+                                  IndPeriod_WarmUp = 0L,
+                                  IniStates = OutputsModel1$StateEnd)
+  OutputsModel2 <- RunModel(InputsModel = InputsModel,
+                                 RunOptions = RunOptions2, Param = ParamSD, FUN_MOD = RunModel_GR4J)
+  # Compare 1991 Qsim from warm started and from 1990-1991
+  names(OutputsModel2$Qsim) <- NULL
+  expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730])
+})
\ No newline at end of file