v1.6.2.3 fix: error in lag formulation

- Following correction proposed #34 (comment 22395)
- Correct tests consequently

Refs #34

DESCRIPTION

 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.6.2.2
-Date: 2020-06-05
+Version: 1.6.2.3
+Date: 2020-07-16
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
   person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),

NEWS.md

 ## Release History of the airGR Package
 
+### 1.6.2.2 Release Notes (2020-07-15)
+
+#### Bug fixes
+
+- Major error in lag model implementation
+
+____________________________________________________________________________________
+
+
 ### 1.6.2.1 Release Notes (2020-06-05)
 
 #### New features

R/RunModel.R

@@ -30,10 +30,10 @@ RunModel <- function (InputsModel, RunOptions, Param, FUN_MOD) {
       }
       OutputsModelDown$Qsim <- OutputsModelDown$Qsim +
         c(rep(0, floor(PT[upstream_basin])),
-          Qupstream[(1 + floor(PT[upstream_basin])):LengthTs]) *
+          Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
         HUTRANS[1, upstream_basin] +
         c(rep(0, floor(PT[upstream_basin] + 1)),
-          Qupstream[(2 + floor(PT[upstream_basin])):LengthTs]) *
+          Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
         HUTRANS[2, upstream_basin]
     }
     # Warning for negative flows

tests/testthat/test-RunModel_LAG.R

@@ -2,7 +2,6 @@ context("RunModel_LAG")
 
 data(L0123001)
 
-
 test_that("'BasinAreas' must have one more element than 'LengthHydro'", {
   expect_error(
     InputsModel <- CreateInputsModel(
@@ -18,6 +17,8 @@ test_that("'BasinAreas' must have one more element than 'LengthHydro'", {
   )
 })
 
+BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
+
 test_that("'Qupstream' cannot contain any NA value", {
   expect_error(
     InputsModel <- CreateInputsModel(
@@ -27,14 +28,14 @@ test_that("'Qupstream' cannot contain any NA value", {
       PotEvap = BasinObs$E,
       Qupstream = matrix(BasinObs$Qmm, ncol = 1),
       LengthHydro = 1,
-      BasinAreas = c(1, 2)
+      BasinAreas = BasinAreas
     ),
     regexp = "'Qupstream' cannot contain any NA value"
   )
 })
 
-Qupstream = BasinObs$Qmm
-Qupstream[is.na(Qupstream)] = mean(Qupstream, na.rm = TRUE)
+# Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
+Qupstream = floor((sin((1:length(BasinObs$Qmm)/365*2*3.14))+1)*mean(BasinObs$Qmm, na.rm = T))
 
 InputsModel <- CreateInputsModel(
   FUN_MOD = RunModel_GR4J,
@@ -43,7 +44,7 @@ InputsModel <- CreateInputsModel(
   PotEvap = BasinObs$E,
   Qupstream = matrix(Qupstream, ncol = 1),
   LengthHydro = 1000,
-  BasinAreas = c(BasinInfo$BasinArea * 2, BasinInfo$BasinArea)
+  BasinAreas = BasinAreas
 )
 
 Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
@@ -60,7 +61,6 @@ OutputsGR4JOnly <-
                 RunOptions = RunOptions,
                 Param = Param)
 
-
 test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
   UpstBasinArea = InputsModel$BasinAreas[1]
   InputsModel$BasinAreas[1] <- 0
@@ -88,25 +88,40 @@ test_that(
 
 ParamSD = c(Param, InputsModel$LengthHydro / (24 * 60 * 60)) # Speed corresponding to one time step delay
 
+QlsGR4Only <-
+  OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
+
 test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
-  QlsGR4Only <-
-    OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
   OutputsSD <-
     RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
   QlsSdSim <-
     OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
   QlsUpstLagObs <-
-    c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)] + 1]) * InputsModel$BasinAreas[1] * 1E6 / 86400
+    c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * InputsModel$BasinAreas[1] * 1E6 / 86400
   expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
 })
 
+test_that("1 input with lag of 0.5 time step delay out gives an output delayed of 0.5 time step", {
+  OutputsSD <-
+    RunModel(InputsModel, RunOptions, Param = c(Param, InputsModel$LengthHydro / (12 * 3600)), FUN_MOD = RunModel_GR4J)
+  QlsSdSim <-
+    OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
+  QlsUpstLagObs <-
+    (Qupstream[Ind_Run] + c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]))/2 * InputsModel$BasinAreas[1] * 1E6 / 86400
+  expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
+})
+
+
 test_that("Params from calibration with simulated data should be similar to initial params", {
   InputsCrit <- CreateInputsCrit(
     FUN_CRIT = ErrorCrit_NSE,
     InputsModel = InputsModel,
     RunOptions = RunOptions,
     VarObs = "Q",
-    Obs = BasinObs$Qmm[Ind_Run]
+    Obs = (
+      c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1] +
+      BasinObs$Qmm[Ind_Run] * BasinAreas[2]
+    ) / sum(BasinAreas)
   )
   CalibOptions <- CreateCalibOptions(
     FUN_MOD = RunModel_GR4J,
@@ -114,31 +129,32 @@ test_that("Params from calibration with simulated data should be similar to init
     IsSD = TRUE
   )
   OutputsCalib <- Calibration_Michel(
-      InputsModel = InputsModel,
-      RunOptions = RunOptions,
-      InputsCrit = InputsCrit,
-      CalibOptions = CalibOptions,
-      FUN_MOD = RunModel_GR4J
-    )
-  expect_equal(OutputsCalib$ParamFinalR, ParamSD, tolerance = 1E-3)
+    InputsModel = InputsModel,
+    RunOptions = RunOptions,
+    InputsCrit = InputsCrit,
+    CalibOptions = CalibOptions,
+    FUN_MOD = RunModel_GR4J
+  )
+  expect_equal(OutputsCalib$ParamFinalR[1:4] / ParamSD[1:4], rep(1, 4), tolerance = 1E-2)
+  expect_equal(OutputsCalib$ParamFinalR[5], ParamSD[5], tolerance = 2E-3)
 })
 
 test_that("1 no area input with lag of 1 time step delay out gives an output delayed of one time step converted to mm", {
   Qm3GR4Only <-
-    OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E3
+    OutputsGR4JOnly$Qsim * BasinAreas[2] * 1E3
   # Specify that upstream flow is not related to an area
-  InputsModel$BasinAreas = c(NA, BasinInfo$BasinArea)
-  # Convert upstream flow to Liter/day
-  InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinInfo$BasinArea * 2
+  InputsModel$BasinAreas = c(NA, BasinAreas[2])
+  # Convert upstream flow to m3/day
+  InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1E3
+
   OutputsSD <-
     RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
-  
+
   expect_false(any(is.na(OutputsSD$Qsim)))
-  
-  Qm3SdSim <-
-    OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1E3
+
+  Qm3SdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1E3
   Qm3UpstLagObs <-
-    c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)] + 1])
-  
+    c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
+
   expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
 })