Merge branch '123-add-warm-up-period-simulated-flows-in-outputsmodel' into 'dev'

Resolve "Add warm-up period simulated flows in OutputsModel"

Closes #123

See merge request !49

.Rbuildignore

@@ -10,3 +10,4 @@
 ^Rplots\.pdf$
 ^ci$
 ^data-raw$
+^revdep$

.gitignore

@@ -12,6 +12,8 @@ packrat/lib*/
 *.pdf
 !man/figures/*.pdf
 
+# revdep
+/revdep/
 
 ######################################################################################################
 ### Generic .gitignore for R (source: https://github.com/github/gitignore/blob/master/R.gitignore) ###

.regressionignore

@@ -33,3 +33,33 @@ RunModel_GR6J RunOptions
 RunModel_Lag RunOptions
 RunModel RunOptions
 SeriesAggreg RunOptions
+Calibration OutputsModel
+Calibration_Michel OutputsModel
+CreateCalibOptions OutputsModel
+CreateIniStates OutputsModel
+CreateInputsCrit OutputsModel
+CreateInputsModel OutputsModel
+CreateRunOptions OutputsModel
+ErrorCrit OutputsModel
+ErrorCrit_KGE OutputsModel
+ErrorCrit_KGE2 OutputsModel
+ErrorCrit_NSE OutputsModel
+ErrorCrit_RMSE OutputsModel
+Imax OutputsModel
+RunModel OutputsModel
+RunModel_CemaNeige OutputsModel
+RunModel_CemaNeigeGR4J OutputsModel
+RunModel_CemaNeigeGR5J OutputsModel
+RunModel_CemaNeigeGR6J OutputsModel
+RunModel_GR1A OutputsModel
+RunModel_GR2M OutputsModel
+RunModel_GR4H OutputsModel
+RunModel_GR4J OutputsModel
+RunModel_GR5H OutputsModel
+RunModel_GR5J OutputsModel
+RunModel_GR6J OutputsModel
+RunModel_Lag OutputsModel
+SeriesAggreg OutputsModel
+Param_Sets_GR4J OutputsModel_Val
+RunModel_Lag OutputsModelDown
+SeriesAggreg SimulatedMonthlyRegime

R/CreateRunOptions.R

@@ -295,31 +295,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
   ##check_Outputs_Cal_and_Sim
 
   ##Outputs_all
-  Outputs_all <- NULL
-  if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4H")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR5H")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4J")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR5J")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR6J")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR2M)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR2M")$GR)
-  }
-  if (identical(FUN_MOD, RunModel_GR1A)) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR1A")$GR)
-  }
-  if ("CemaNeige" %in% ObjectClass) {
-    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = NULL, isCN = TRUE)$CN)
-  }
+  Outputs_all <- c("DatesR", unlist(FortranOutputs), "WarmUpQsim", "StateEnd")
 
   ##check_Outputs_Sim
   if (!is.vector(Outputs_Sim)) {
@@ -332,9 +308,9 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
     stop("'Outputs_Sim' must not contain NA")
   }
   if ("all" %in% Outputs_Sim) {
-    Outputs_Sim <- c("DatesR", Outputs_all, "StateEnd")
+    Outputs_Sim <- Outputs_all
   }
-  Test <- which(!Outputs_Sim %in% c("DatesR", Outputs_all, "StateEnd"))
+  Test <- which(!Outputs_Sim %in% Outputs_all)
   if (length(Test) != 0) {
     stop(paste0( "'Outputs_Sim' is incorrectly defined: ",
                  paste(Outputs_Sim[Test], collapse = ", "), " not found"))
@@ -366,10 +342,9 @@ CreateRunOptions <- function(FUN_MOD, InputsModel,
     }
   }
   if ("all" %in% Outputs_Cal) {
-    Outputs_Cal <- c("DatesR", Outputs_all, "StateEnd")
-
+    Outputs_Cal <- Outputs_all
   }
-  Test <- which(!Outputs_Cal %in% c("DatesR", Outputs_all, "StateEnd"))
+  Test <- which(!Outputs_Cal %in% Outputs_all)
   if (length(Test) != 0) {
     stop(paste0("'Outputs_Cal' is incorrectly defined: ",
                 paste(Outputs_Cal[Test], collapse = ", "), " not found"))

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(
+          from = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
+          to   = 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 = ", "))
 

R/SeriesAggreg.OutputsModel.R

@@ -2,6 +2,6 @@ SeriesAggreg.OutputsModel <- function(x, Format, ...) {
   SeriesAggreg.list(x,
                     Format,
                     ConvertFun = NA,
-                    except = "StateEnd",
+                    except = c("WarmUpQsim", "StateEnd"),
                     ...)
 }

R/UtilsRunModel.R

@@ -10,23 +10,19 @@
 #' @noRd
 #'
 .GetOutputsModelGR <- function(InputsModel,
-                                RunOptions,
-                                RESULTS,
-                                LInputSeries,
-                                CemaNeigeLayers = NULL) {
+                               RunOptions,
+                               RESULTS,
+                               LInputSeries,
+                               CemaNeigeLayers = NULL) {
 
   IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
   FortranOutputs <- RunOptions$FortranOutputs$GR
 
-  if ("all" %in% RunOptions$Outputs_Sim) {
-    IndOutputs <- as.integer(1:length(FortranOutputs))
-  } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
-  }
+  IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
 
   OutputsModel <- list()
 
-  if("DatesR" %in% RunOptions$Outputs_Sim) {
+  if ("DatesR" %in% RunOptions$Outputs_Sim) {
     OutputsModel$DatesR <- InputsModel$DatesR[RunOptions$IndPeriod_Run]
   }
 
@@ -36,11 +32,17 @@
   OutputsModel <- c(OutputsModel,
                     lapply(seqOutputs, function(i) RESULTS$Outputs[IndPeriod2, i]))
 
-  if(!is.null(CemaNeigeLayers)) {
+  if (!is.null(CemaNeigeLayers)) {
     OutputsModel$CemaNeigeLayers <- CemaNeigeLayers
   }
 
-  if("StateEnd" %in% RunOptions$Outputs_Sim) {
+  if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$WarmUpQsim <- RESULTS$Outputs[seq_len(length(RunOptions$IndPeriod_WarmUp)),
+                                               which(FortranOutputs == "Qsim")]
+    class(OutputsModel$WarmUpQsim) <- c("WarmUpOutputsModelItem", class(OutputsModel$WarmUpQsim))
+  }
+
+  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
     OutputsModel$StateEnd <- RESULTS$StateEnd
   }
 

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,14 +174,14 @@ 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)
 })
 
 test_that("Qupstream in different units should return the same result", {
   OutputsSD_mm <- RunModel(InputsModel, RunOptions,
-                        Param = ParamSD,
-                        FUN_MOD = RunModel_GR4J)
+                           Param = ParamSD,
+                           FUN_MOD = RunModel_GR4J)
   InputsModel_m3 <- CreateInputsModel(
     FUN_MOD = RunModel_GR4J,
     DatesR = BasinObs$DatesR,
@@ -208,8 +208,8 @@ test_that("Qupstream in different units should return the same result", {
     QupstrUnit = "m3/s"
   )
   OutputsSD_m3s <- RunModel(InputsModel_m3s, RunOptions,
-                           Param = ParamSD,
-                           FUN_MOD = RunModel_GR4J)
+                            Param = ParamSD,
+                            FUN_MOD = RunModel_GR4J)
   expect_equal(OutputsSD_mm$Qsim, OutputsSD_m3s$Qsim)
 
   InputsModel_ls <- CreateInputsModel(
@@ -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))
+})

vignettes/V05_sd_model.Rmd

 ---
-title: "Simulating a reservoir with semi-distributed GR4J model"
+title: "Simulated vs observed upstream flows in calibration of semi-distributed GR4J model"
 author: "David Dorchies"
 bibliography: V00_airgr_ref.bib
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteEngine{knitr::rmarkdown}
-  %\VignetteIndexEntry{Simulating a reservoir with semi-distributed GR4J model}
+  %\VignetteIndexEntry{Simulated vs observed upstream flows in calibration of semi-distributed GR4J model}
   %\VignetteEncoding{UTF-8}
 ---
 
@@ -124,10 +124,13 @@ OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
                                         FUN_MOD = RunModel_GR4J)
 ```
 
-To run the complete model, we should substitute the observed upstream flow by the simulated one:
+To run the complete model, we should substitute the observed upstream flow by the simulated one for the entire period of simulation (warm-up + run):
 
 ```{r}
 InputsModelDown2 <- InputsModelDown1
+# Simulated flow during warm-up period 
+InputsModelDown2$Qupstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$WarmUpQsim
+# Simulated flow during run period
 InputsModelDown2$Qupstream[Ind_Run] <- OutputsModelUp$Qsim
 ```