test(sd): Add calibration test

- Check if calibrated params from simulation are equal to initial parameters
- typo corrections

Refs #34

tests/testthat/test-RunModel_LAG.R

@@ -4,82 +4,131 @@ data(L0123001)
 
 
 test_that("'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'", {
-    expect_error(
-        InputsModel <- CreateInputsModel(
-            FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
-            Precip = BasinObs$P, PotEvap = BasinObs$E,
-            QobsUpstr = matrix(BasinObs$Qmm * 2,ncol=1),
-            LengthHydro = matrix(c(1),nrow=1),
-            BasinAreas = matrix(c(1),nrow=1)
-        ),
-        regexp = "'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'"
-    )
+  expect_error(
+    InputsModel <- CreateInputsModel(
+      FUN_MOD = RunModel_GR4J,
+      DatesR = BasinObs$DatesR,
+      Precip = BasinObs$P,
+      PotEvap = BasinObs$E,
+      QobsUpstr = matrix(BasinObs$Qmm, ncol = 1),
+      LengthHydro = matrix(c(1), nrow = 1),
+      BasinAreas = matrix(c(1), nrow = 1)
+    ),
+    regexp = "'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'"
+  )
 })
 
 test_that("'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'", {
-    expect_error(
-        InputsModel <- CreateInputsModel(
-            FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
-            Precip = BasinObs$P, PotEvap = BasinObs$E,
-            QobsUpstr = matrix(BasinObs$Qmm, ncol = 1),
-            LengthHydro = matrix(c(1), nrow = 1),
-            BasinAreas = matrix(c(1, 2),nrow = 1)
-        ),
-        regexp = "'QobsUpstr' cannot contain any NA value"
-    )
+  expect_error(
+    InputsModel <- CreateInputsModel(
+      FUN_MOD = RunModel_GR4J,
+      DatesR = BasinObs$DatesR,
+      Precip = BasinObs$P,
+      PotEvap = BasinObs$E,
+      QobsUpstr = matrix(BasinObs$Qmm, ncol = 1),
+      LengthHydro = matrix(c(1), nrow = 1),
+      BasinAreas = matrix(c(1, 2), nrow = 1)
+    ),
+    regexp = "'QobsUpstr' cannot contain any NA value"
+  )
 })
 
 QobsUpstr = BasinObs$Qmm
 QobsUpstr[is.na(QobsUpstr)] = mean(QobsUpstr, na.rm = TRUE)
 
 InputsModel <- CreateInputsModel(
-    FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
-    Precip = BasinObs$P, PotEvap = BasinObs$E,
-    QobsUpstr = matrix(QobsUpstr,ncol=1),
-    LengthHydro = matrix(c(1),nrow=1),
-    BasinAreas = matrix(c(BasinInfo$BasinArea,BasinInfo$BasinArea),nrow=1)
+  FUN_MOD = RunModel_GR4J,
+  DatesR = BasinObs$DatesR,
+  Precip = BasinObs$P,
+  PotEvap = BasinObs$E,
+  QobsUpstr = matrix(QobsUpstr, ncol = 1),
+  LengthHydro = matrix(c(1), nrow = 1),
+  BasinAreas = matrix(c(
+    BasinInfo$BasinArea * 2, BasinInfo$BasinArea
+  ), nrow = 1)
 )
 
-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")
-)
+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,
-    InputsModel = InputsModel, IndPeriod_Run = Ind_Run
-)
+RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                               InputsModel = InputsModel,
+                               IndPeriod_Run = Ind_Run)
 
 Param = c(257.237556, 1.012237, 88.234673, 2.207958) #  From vignettes/V01_get_started
 
-OutputsGR4JOnly <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
+OutputsGR4JOnly <-
+  RunModel_GR4J(InputsModel = InputsModel,
+                RunOptions = RunOptions,
+                Param = Param)
 
 
 test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
-  UpstBasinArea = InputsModel$BasinAreas[1,1]
-  InputsModel$BasinAreas[1,1] <<- 0
-  OutputsSD <- RunModel(InputsModel, RunOptions, Param = c(Param, 1), FUN_MOD = RunModel_GR4J)
+  UpstBasinArea = InputsModel$BasinAreas[1, 1]
+  InputsModel$BasinAreas[1, 1] <<- 0
+  OutputsSD <-
+    RunModel(InputsModel,
+             RunOptions,
+             Param = c(Param, 1),
+             FUN_MOD = RunModel_GR4J)
   expect_equal(OutputsGR4JOnly$Qsim, OutputsSD$Qsim)
-  InputsModel$BasinAreas[1,1] <<- UpstBasinArea
+  InputsModel$BasinAreas[1, 1] <<- UpstBasinArea
 })
 
-test_that("Downstream basin with nil area and nul upstream length should return same Qdown as Qupstream alone", {
-  InputsModel <- CreateInputsModel(
-    FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
-    Precip = BasinObs$P, PotEvap = BasinObs$E,
-    QobsUpstr = matrix(QobsUpstr,ncol=1),
-    LengthHydro = matrix(c(0),nrow=1),
-    BasinAreas = matrix(c(BasinInfo$BasinArea,0),nrow=1)
-  )
-  OutputsSD <- RunModel(InputsModel, RunOptions, Param = c(Param, 1), FUN_MOD = RunModel_GR4J)
-  expect_equal(OutputsSD$Qsim, QobsUpstr[Ind_Run])
-})
+test_that(
+  "Downstream basin with nil area and nul upstream length should return same Qdown as Qupstream alone",
+  {
+    InputsModelZeroDown <- CreateInputsModel(
+      FUN_MOD = RunModel_GR4J,
+      DatesR = BasinObs$DatesR,
+      Precip = BasinObs$P,
+      PotEvap = BasinObs$E,
+      QobsUpstr = matrix(QobsUpstr, ncol = 1),
+      LengthHydro = matrix(c(0), nrow = 1),
+      BasinAreas = matrix(c(BasinInfo$BasinArea, 0), nrow = 1)
+    )
+    OutputsSD <-
+      RunModel(InputsModelZeroDown,
+               RunOptions,
+               Param = c(Param, 1),
+               FUN_MOD = RunModel_GR4J)
+    expect_equal(OutputsSD$Qsim, QobsUpstr[Ind_Run])
+  }
+)
+
+ParamSD = c(Param, InputsModel$LengthHydro * 1000 / (24 * 60 * 60))
+OutputsSD <-
+  RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
 
 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
-    ParamSD = c(Param, InputsModel$LengthHydro * 1000 / (24 * 60 * 60))
-    OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
-    QlsSdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
-    QlsUpstLagObs <- c(0, QobsUpstr[Ind_Run[1:(length(Ind_Run)-1)] + 1]) * InputsModel$BasinAreas[1] * 1E6 / 86400
-    expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
+  QlsGR4Only <-
+    OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
+  QlsSdSim <-
+    OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
+  QlsUpstLagObs <-
+    c(0, QobsUpstr[Ind_Run[1:(length(Ind_Run) - 1)] + 1]) * 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]
+  )
+  CalibOptions <- CreateCalibOptions(
+    FUN_MOD = RunModel_GR4J, 
+    FUN_CALIB = Calibration_Michel, 
+    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)
 })