Merge branch '54-calibration-with-parameter-regularisation' into 'dev'

Resolve "Calibration with parameter regularisation"

Closes #54

See merge request !24

NAMESPACE

@@ -31,6 +31,9 @@ export(CreateInputsModel)
 export(CreateRunOptions)
 export(CreateSupervisor)
 export(RunModel)
+export(getNoSD_Ids)
+export(getSD_Ids)
+export(isNodeDownstream)
 import(airGR)
 importFrom(grDevices,rainbow)
 importFrom(graphics,matplot)

R/Calibration.GRiwrmInputsModel.R

@@ -8,11 +8,28 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
                                           useUpstreamQsim = TRUE,
                                           ...) {
 
+  # Argument checks
+
+  # We invoke the mandatory arguments here for avoiding
+  # a messy error message on "get(x)" if an argument is missing
+  InputsModel
+  RunOptions
+  InputsCrit
+  CalibOptions
+
+  # Checking argument classes
+  vars2check <- c("InputsModel", "RunOptions", "InputsCrit", "CalibOptions")
+  lapply(vars2check, function(x) {
+    if (!inherits(get(x), paste0("GRiwrm", x))) {
+      stop(sprintf("'%1$s' must be of class GRiwrm%1$s, type '?Create%1$s' for help", x))
+    }
+  })
+
   OutputsCalib <- list()
-  class(OutputsCalib) <- append(class(OutputsCalib), "GRiwrmOutputsCalib")
+  class(OutputsCalib) <- append("GRiwrmOutputsCalib", class(OutputsCalib))
 
   OutputsModel <- list()
-  class(OutputsModel) <- append(class(OutputsModel), "GRiwrmOutputsModel")
+  class(OutputsModel) <- append("GRiwrmOutputsModel", class(OutputsModel))
 
   for(IM in InputsModel) {
     message("Calibration.GRiwrmInputsModel: Treating sub-basin ", IM$id, "...")
@@ -22,10 +39,16 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
       IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]], OutputsModel)
     }
 
+    if (inherits(InputsCrit[[IM$id]], "InputsCritLavenneFunction")) {
+      IC <- getInputsCrit_Lavenne(IM$id, OutputsModel, InputsCrit)
+    } else {
+      IC <- InputsCrit[[IM$id]]
+    }
+
     OutputsCalib[[IM$id]] <- Calibration(
       InputsModel = IM,
       RunOptions = RunOptions[[IM$id]],
-      InputsCrit = InputsCrit[[IM$id]],
+      InputsCrit = IC,
       CalibOptions = CalibOptions[[IM$id]],
       ...
     )
@@ -44,3 +67,30 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
   return(OutputsCalib)
 
 }
+
+#' Create InputsCrit for De Lavenne regularisation
+#'
+#' Internal function that run [airGR::CreateInputsCrit_Lavenne] on-the-fly with a priori upstream
+#' sub-catchment parameters grabbed during network calibration process.
+#'
+#' @param id [character] the id of the current sub-catchment
+#' @param OutputsModel \[GRiwrmOutputsModel\] object with simulation results of upstream sub-catchments run with calibrated parameters
+#' @param InputsCrit \[InputsCritLavenneFunction\] object internally created by [CreateInputsCrit.GRiwrmInputsModel]
+#'
+#' @return \[InputsCrit\] object with De Lavenne regularisation
+#' @noRd
+#'
+getInputsCrit_Lavenne <- function(id, OutputsModel, InputsCrit) {
+  if (!inherits(InputsCrit[[id]], "InputsCritLavenneFunction")) {
+    stop("'InputsCrit[[id]]' must be of class InputsCritLavenneFunction")
+  }
+  AprioriId <- attr(InputsCrit[[id]], "AprioriId")
+  Lavenne_FUN <- attr(InputsCrit[[id]], "Lavenne_FUN")
+  AprParamR <- OutputsModel[[AprioriId]]$Param
+  if(!inherits(OutputsModel[[AprioriId]], "SD")) {
+    # Add neutral velocity parameter for upstream catchment
+    AprParamR <- c(NA, AprParamR)
+  }
+  AprCrit <- ErrorCrit(InputsCrit[[AprioriId]], OutputsModel[[AprioriId]])$CritValue
+  return(Lavenne_FUN(AprParamR, AprCrit))
+}

R/CreateCalibOptions.GRiwrmInputsModel.R

@@ -3,6 +3,7 @@
 CreateCalibOptions.GRiwrmInputsModel <- function(InputsModel, ...) {
 
   CalibOptions <- list()
+  class(CalibOptions) <- c("GRiwrmCalibOptions", class(CalibOptions))
 
   for(IM in InputsModel) {
     CalibOptions[[IM$id]] <- CreateCalibOptions.InputsModel(

R/CreateInputsCrit.GRiwrmInputsModel.R

@@ -4,9 +4,56 @@ CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel,
                                                FUN_CRIT = airGR::ErrorCrit_NSE,
                                                RunOptions,
                                                Obs,
+                                               AprioriIds = NULL,
+                                               k = 0.15,
                                                ...) {
+  # Parameter checks
+
+  # We invoke the mandatory arguments here for avoiding
+  # a messy error message on "get(x)" if an argument is missing
+  InputsModel
+  RunOptions
+  Obs
+
+  # Checking argument classes
+  lVars2Check <- list(InputsModel = "GRiwrmInputsModel",
+                      RunOptions = "GRiwrmRunOptions",
+                      Obs = c("matrix", "data.frame"))
+  lapply(names(lVars2Check), function(argName) {
+    b <- sapply(lVars2Check[[argName]], function(argClass) {
+      !inherits(get(argName), argClass)
+    })
+    if (all(b)) {
+      stop(sprintf("'%s' must be of class %s", argName, paste(lVars2Check[[argName]], collapse = " or ")))
+    }
+  })
+
+  if (!is.null(AprioriIds)) {
+    AprioriIds <- unlist(AprioriIds)
+    if (!is.character(AprioriIds) || is.null(names(AprioriIds))) {
+      stop("Argument 'AprioriIds' must be a named list or a named vector of characters")
+    }
+    if (length(unique(names(AprioriIds))) != length(names(AprioriIds))) {
+      stop("Each name of AprioriIds items must be unique: duplicate entry detected")
+    }
+    lapply(names(AprioriIds), function(id) {
+      if (!id %in% names(InputsModel)) {
+        stop("'Each item of names(AprioriIds) must be an id of a simulated node:",
+             " the id \"", id ,"\" is unknown")
+      }
+      if (!AprioriIds[id] %in% names(InputsModel)) {
+        stop("'Each item of AprioriIds must be an id of a simulated node:",
+             " the id \"", id ,"\" is unknown")
+      }
+      if (! isNodeDownstream(InputsModel, AprioriIds[id], id)) {
+        stop("'AprioriIds': the node \"", AprioriIds[id],
+             "\" is not upstream the node \"", id,"\"")
+      }
+    })
+  }
+
   InputsCrit <- list()
-  class(InputsCrit) <- append(class(InputsCrit), "GRiwrmInputsCrit")
+  class(InputsCrit) <- append("GRiwrmInputsCrit", class(InputsCrit))
 
   for(IM in InputsModel) {
     InputsCrit[[IM$id]] <- CreateInputsCrit.InputsModel(
@@ -16,7 +63,34 @@ CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel,
       Obs = Obs[, IM$id],
       ...
     )
+    if (!is.null(AprioriIds) && IM$id %in% names(AprioriIds)) {
+      # De Lavenne regularisation for this sub-catchment
+      attr(InputsCrit[[IM$id]], "Lavenne_FUN") <-
+        CreateLavenneFunction(
+          InputsModel = IM,
+          FUN_CRIT = FUN_CRIT,
+          RunOptions = RunOptions[[IM$id]],
+          Obs = Obs[, IM$id],
+          k = k,
+          ...
+        )
+      attr(InputsCrit[[IM$id]], "AprioriId") <- AprioriIds[IM$id]
+      class(InputsCrit[[IM$id]]) <- c("InputsCritLavenneFunction", class(InputsCrit[[IM$id]]))
+    }
   }
 
   return(InputsCrit)
 }
+
+CreateLavenneFunction <- function(InputsModel, FUN_CRIT, RunOptions, Obs, k, ...) {
+  function(AprParamR, AprCrit) {
+    CreateInputsCrit_Lavenne(FUN_CRIT = FUN_CRIT,
+                               InputsModel = InputsModel,
+                               RunOptions = RunOptions,
+                               Obs = Obs,
+                               AprParamR = AprParamR,
+                               AprCrit = AprCrit,
+                               k = k,
+                               ...)
+  }
+}

R/CreateInputsCrit.R

@@ -6,6 +6,8 @@
 #' @param FUN_CRIT \[function (atomic or list)\] error criterion function (e.g. [airGR::ErrorCrit_RMSE], [airGR::ErrorCrit_NSE])
 #' @param RunOptions object of class \emph{RunOptions} or \emph{GRiwrmRunOptions}, see [CreateRunOptions]
 #' @param Obs [numeric], [matrix] or [data.frame] series of observed flows, see details
+#' @param AprioriIds (optional) named [list] or named [vector] of [character] used for the parameter regularisation (see details)
+#' @param k (optional) [numeric] weight coefficient used in the parameter regularisation (See [airGR::CreateInputsCrit_Lavenne])
 #' @param ... arguments passed to [airGR::CreateInputsCrit], see details
 #'
 #' @details See [airGR::CreateInputsCrit] documentation for a complete list of arguments.
@@ -14,10 +16,15 @@
 #'
 #' With a \emph{GRiwrmInputsModel} object, all arguments are applied on each sub-catchments of the network.
 #'
+#' Parameter regularisation consists of defining a priori parameters which are used in a composed criterion based on the formula proposed by de Lavenne et al. (2019) (See [airGR::CreateInputsCrit_Lavenne]).
+#' The parameter `AprioriIds` allows to define which upstream sub-catchment is used for providing a priori parameters. Its format is as follows: `AprioriIds <- c("Downstream sub-catchment 1" = "A priori upstream sub-catchment 1", ...)` where the quoted strings are the ids of the sub-catchments. See vignettes for more details.
+#'
 #' @return Depending on the class of `InputsModel` argument (respectively `InputsModel` and `GRiwrmInputsModel` object), the returned value is respectively:
 #' - a `InputsCrit` object (See [airGR::CreateInputsCrit])
 #' - a `GRiwrmInputsCrit` object which is a [list] of `InputsCrit` object with one item per modelled sub-catchment
 #'
+#' @references De Lavenne, A., Andréassian, V., Thirel, G., Ramos, M.-H., Perrin, C., 2019. A Regularization Approach to Improve the Sequential Calibration of a Semidistributed Hydrological Model. Water Resources Research 55, 8821–8839. \doi{10.1029/2018WR024266}
+#'
 #' @rdname CreateInputsCrit
 #' @export
 CreateInputsCrit <- function(InputsModel, ...) {

R/Severn.R

-#' Catchment Attributes and Hydro-Meteorological Timeseries for some gauging stations on the Severn River
+#' Catchment attributes and hydro-meteorological timeseries for some gauging stations on the Severn River
 #'
 #' @format a [list] with 2 items:
 #'
 #' - "BasinsInfo" which contains a [data.frame] with Gauging station identifier, name, coordinates (GPS), area (km2), mean elevation (m), station type, flow period start and end, the bank full flow (m3/s), the identifier of the following downstream station and the distance to the following downstream station
 #' - "BasinObs" which contains a [list] with an item by gauging station which contains a [data.frame] with [POSIXct] dates, precipitations (mm/time step), potential evapotranspiration (mm/time step) and measured flows (mm/time step)
 #'
-#' @source \url{https://doi.org/10.5285/8344E4F3-D2EA-44F5-8AFA-86D2987543A9}
+#' @source These data are extracted from the CAMEL-GB dataset.
+#'
+#' Coxon, G.; Addor, N.; Bloomfield, J.P.; Freer, J.; Fry, M.; Hannaford, J.; Howden, N.J.K.; Lane, R.; Lewis, M.; Robinson, E.L.; Wagener, T.; Woods, R. (2020). Catchment attributes and hydro-meteorological timeseries for 671 catchments across Great Britain (CAMELS-GB). NERC Environmental Information Data Centre. (Dataset). \doi{10.5285/8344E4F3-D2EA-44F5-8AFA-86D2987543A9}
 "Severn"

R/utils.R

@@ -3,7 +3,7 @@
 #' @param InputsModel \[`GRiwrmInputsModel` object\]
 #'
 #' @return [character] IDs of the sub-basins using SD model
-#' @noRd
+#' @export
 getSD_Ids <- function(InputsModel) {
   if (!inherits(InputsModel, "GRiwrmInputsModel")) {
     stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
@@ -19,7 +19,7 @@ getSD_Ids <- function(InputsModel) {
 #' @param InputsModel \[`GRiwrmInputsModel` object\]
 #'
 #' @return [character] IDs of the sub-basins not using the SD model
-#' @noRd
+#' @export
 getNoSD_Ids <- function(InputsModel) {
   if (!inherits(InputsModel, "GRiwrmInputsModel")) {
     stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
@@ -152,6 +152,7 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
   return(dfQsim)
 }
 
+
 #' Convert IniStates list into a vector
 #'
 #' @param IniStates see [CreateIniStates]
@@ -164,3 +165,20 @@ serializeIniStates <- function(IniStates) {
   IniStates[is.na(IniStates)] <- 0
   return(IniStates)
 }
+
+
+#' Check if a node is downstream another one
+#'
+#' @param InputsModel \[`GRiwrmInputsModel` object\] see [CreateInputsModel.GRiwrm] for details
+#' @param current_node [character] with the id of the current node
+#' @param down_node [character] with the id of the node for which we want to know if it is downstream `current_node`
+#'
+#' @return [logical] `TRUE` if the node with the id `down_node` is downstream the node with the id `current_node`
+#' @export
+#'
+isNodeDownstream <- function(InputsModel, current_node, down_node) {
+  current_down_node <- InputsModel[[current_node]]$down
+  if (is.na(current_down_node)) return(FALSE)
+  if (current_down_node == down_node) return(TRUE)
+  return(isNodeDownstream(InputsModel, current_down_node, down_node))
+}

man/CreateInputsCrit.Rd

@@ -12,6 +12,8 @@
   FUN_CRIT = airGR::ErrorCrit_NSE,
   RunOptions,
   Obs,
+  AprioriIds = NULL,
+  k = 0.15,
   ...
 )
 
@@ -28,6 +30,10 @@ CreateInputsCrit(InputsModel, ...)
 
 \item{Obs}{\link{numeric}, \link{matrix} or \link{data.frame} series of observed flows, see details}
 
+\item{AprioriIds}{(optional) named \link{list} or named \link{vector} of \link{character} used for the parameter regularisation (see details)}
+
+\item{k}{(optional) \link{numeric} weight coefficient used in the parameter regularisation (See \link[airGR:CreateInputsCrit_Lavenne]{airGR::CreateInputsCrit_Lavenne})}
+
 \item{...}{arguments passed to \link[airGR:CreateInputsCrit]{airGR::CreateInputsCrit}, see details}
 }
 \value{
@@ -46,4 +52,10 @@ See \link[airGR:CreateInputsCrit]{airGR::CreateInputsCrit} documentation for a c
 \code{Obs} argument is equivalent to the same argument in \link[airGR:CreateInputsCrit]{airGR::CreateInputsCrit} except that it must a \link{matrix} or a \link{data.frame} if \code{InputsModel} is a \emph{GRiwrmInputsModel} object. Then, each column of the \link{matrix} or \link{data.frame} represents the observations of one of the simulated node with the name of the columns representing the id of each node.
 
 With a \emph{GRiwrmInputsModel} object, all arguments are applied on each sub-catchments of the network.
+
+Parameter regularisation consists of defining a priori parameters which are used in a composed criterion based on the formula proposed by de Lavenne et al. (2019) (See \link[airGR:CreateInputsCrit_Lavenne]{airGR::CreateInputsCrit_Lavenne}).
+The parameter \code{AprioriIds} allows to define which upstream sub-catchment is used for providing a priori parameters. Its format is as follows: \code{AprioriIds <- c("Downstream sub-catchment 1" = "A priori upstream sub-catchment 1", ...)} where the quoted strings are the ids of the sub-catchments. See vignettes for more details.
+}
+\references{
+De Lavenne, A., Andréassian, V., Thirel, G., Ramos, M.-H., Perrin, C., 2019. A Regularization Approach to Improve the Sequential Calibration of a Semidistributed Hydrological Model. Water Resources Research 55, 8821–8839. \doi{10.1029/2018WR024266}
 }

man/Severn.Rd

@@ -3,7 +3,7 @@
 \docType{data}
 \name{Severn}
 \alias{Severn}
-\title{Catchment Attributes and Hydro-Meteorological Timeseries for some gauging stations on the Severn River}
+\title{Catchment attributes and hydro-meteorological timeseries for some gauging stations on the Severn River}
 \format{
 a \link{list} with 2 items:
 \itemize{
@@ -12,12 +12,14 @@ a \link{list} with 2 items:
 }
 }
 \source{
-\url{https://doi.org/10.5285/8344E4F3-D2EA-44F5-8AFA-86D2987543A9}
+These data are extracted from the CAMEL-GB dataset.
+
+Coxon, G.; Addor, N.; Bloomfield, J.P.; Freer, J.; Fry, M.; Hannaford, J.; Howden, N.J.K.; Lane, R.; Lewis, M.; Robinson, E.L.; Wagener, T.; Woods, R. (2020). Catchment attributes and hydro-meteorological timeseries for 671 catchments across Great Britain (CAMELS-GB). NERC Environmental Information Data Centre. (Dataset). \doi{10.5285/8344E4F3-D2EA-44F5-8AFA-86D2987543A9}
 }
 \usage{
 Severn
 }
 \description{
-Catchment Attributes and Hydro-Meteorological Timeseries for some gauging stations on the Severn River
+Catchment attributes and hydro-meteorological timeseries for some gauging stations on the Severn River
 }
 \keyword{datasets}

man/getNoSD_Ids.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{getNoSD_Ids}
+\alias{getNoSD_Ids}
+\title{Function to obtain the ID of sub-basins not using SD model}
+\usage{
+getNoSD_Ids(InputsModel)
+}
+\arguments{
+\item{InputsModel}{[\code{GRiwrmInputsModel} object]}
+}
+\value{
+\link{character} IDs of the sub-basins not using the SD model
+}
+\description{
+Function to obtain the ID of sub-basins not using SD model
+}

man/getSD_Ids.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{getSD_Ids}
+\alias{getSD_Ids}
+\title{Function to obtain the ID of sub-basins using SD model}
+\usage{
+getSD_Ids(InputsModel)
+}
+\arguments{
+\item{InputsModel}{[\code{GRiwrmInputsModel} object]}
+}
+\value{
+\link{character} IDs of the sub-basins using SD model
+}
+\description{
+Function to obtain the ID of sub-basins using SD model
+}

man/isNodeDownstream.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{isNodeDownstream}
+\alias{isNodeDownstream}
+\title{Check if a node is downstream another one}
+\usage{
+isNodeDownstream(InputsModel, current_node, down_node)
+}
+\arguments{
+\item{InputsModel}{[\code{GRiwrmInputsModel} object] see \link{CreateInputsModel.GRiwrm} for details}
+
+\item{current_node}{\link{character} with the id of the current node}
+
+\item{down_node}{\link{character} with the id of the node for which we want to know if it is downstream \code{current_node}}
+}
+\value{
+\link{logical} \code{TRUE} if the node with the id \code{down_node} is downstream the node with the id \code{current_node}
+}
+\description{
+Check if a node is downstream another one
+}

tests/testthat/helper_RunModel.R

+#' Prepare useful variables for GRiwrm tests
+#'
+#' @return [environment] with the variables (See examples section)
+#' @noRd
+#'
+#' @examples
+#' # data set up
+#' e <- setupRunModel()
+# variables are copied from environment 'e' to the current environment
+# https://stackoverflow.com/questions/9965577/r-copy-move-one-environment-to-another
+#' for(x in ls(e)) assign(x, get(x, e))
+#'
+setupRunModel <- function() {
+  data(Severn)
+
+  # Format observation
+  BasinsObs <- Severn$BasinsObs
+  DatesR <- BasinsObs[[1]]$DatesR
+  PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+  PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+  Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
+
+  # Set network
+  nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+  nodes$distance_downstream <- nodes$distance_downstream
+  nodes$model <- "RunModel_GR4J"
+  griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+
+  # Convert meteo data to SD (remove upstream areas)
+  Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+  PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+
+  # Calibration parameters
+  ParamMichel <- list(
+    `54057` = c(0.779999999999999, 57.9743110789593, -1.23788116619639, 0.960789439152323, 2.47147147147147),
+    `54032` = c(1.37562057772709, 1151.73462496385, -0.379248293750608, 6.2243898378232, 8.23716221550954),
+    `54001` = c(1.03, 24.7790862245877, -1.90430150145153, 21.7584023961971, 1.37837837837838),
+    `54095` = c(256.844150254651, 0.0650458497009288, 57.523675209819, 2.71809513102128),
+    `54002` = c(419.437754485522, 0.12473266292168, 13.0379482833606, 2.12230907892238),
+    `54029` = c(219.203385553954, 0.389211590110934, 48.4242150713452, 2.00300300300301)
+  )
+
+  # set up inputs
+  InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+
+  # RunOptions
+  nTS <- 365
+  IndPeriod_Run <- seq(
+    length(InputsModel[[1]]$DatesR) - nTS + 1,
+    length(InputsModel[[1]]$DatesR)
+  )
+  IndPeriod_WarmUp = seq(IndPeriod_Run[1]-365,IndPeriod_Run[1]-1)
+  RunOptions <- CreateRunOptions(
+    InputsModel = InputsModel,
+    IndPeriod_WarmUp = IndPeriod_WarmUp,
+    IndPeriod_Run = IndPeriod_Run
+  )
+
+  # RunModel.GRiwrmInputsModel
+  OM_GriwrmInputs <- RunModel(
+    InputsModel,
+    RunOptions = RunOptions,
+    Param = ParamMichel
+  )
+  return(environment())
+}

tests/testthat/test-Calibration.R

+# data set up
+e <- setupRunModel()
+# variables are copied from environment 'e' to the current environment
+# https://stackoverflow.com/questions/9965577/r-copy-move-one-environment-to-another
+for(x in ls(e)) assign(x, get(x, e))
+
+context("Calibration.GRiwrmInputsModel")
+
+CalibOptions <- CreateCalibOptions(InputsModel = InputsModel)
+
+test_that("Calibrated parameters remains unchanged", {
+  InputsCrit <- CreateInputsCrit(
+    InputsModel = InputsModel,
+    RunOptions = RunOptions,
+    Obs = Qobs[IndPeriod_Run,]
+  )
+
+  OC <- Calibration(
+    InputsModel = InputsModel,
+    RunOptions = RunOptions,
+    InputsCrit = InputsCrit,
+    CalibOptions = CalibOptions
+  )
+
+  ParamFinalR <- lapply(OC, "[[", "ParamFinalR")
+
+  lapply(names(ParamFinalR), function(id) expect_equal(ParamFinalR[[id]], ParamMichel[[id]]))
+
+})
+
+test_that("Calibration with regularisation is OK", {
+  InputsCrit <- CreateInputsCrit(
+    InputsModel = InputsModel,
+    RunOptions = RunOptions,
+    Obs = Qobs[IndPeriod_Run,],
+    AprioriIds = c(
+      "54057" = "54032",
+      "54032" = "54001",
+      "54001" = "54095"
+    )
+  )
+
+  OC <- Calibration(
+    InputsModel = InputsModel,
+    RunOptions = RunOptions,
+    InputsCrit = InputsCrit,
+    CalibOptions = CalibOptions
+  )
+
+  ParamLavenne <- lapply(OC, "[[", "ParamFinalR")
+  expect_equal(OC[["54095"]]$CritFinal, ErrorCrit(
+    InputsCrit[["54095"]],
+    RunModel(InputsModel, RunOptions, ParamLavenne)[["54095"]]
+  )$CritValue)
+  OM <- RunModel(InputsModel, RunOptions, ParamLavenne)
+  lapply(names(OC), function(id) {
+    expect_gt(
+      ErrorCrit(
+        InputsCrit[[id]],
+        OM[[id]]
+      )$CritValue,
+      0.9
+    )
+  })
+})

tests/testthat/test-CreateInputsCrit.R

+# data set up
+e <- setupRunModel()
+# variables are copied from environment 'e' to the current environment
+# https://stackoverflow.com/questions/9965577/r-copy-move-one-environment-to-another
+for(x in ls(e)) assign(x, get(x, e))
+
+context("CreateInputsCrit.GRiwrmInputsModel")
+
+test_that("Wrong argument class should throw error", {
+  expect_error(CreateInputsCrit(InputsModel = InputsModel[[1]],
+                                RunOptions = RunOptions,
+                                Obs = Qobs[IndPeriod_Run,]))
+  expect_error(CreateInputsCrit.GRiwrmInputsModel(InputsModel = InputsModel[[1]],
+                                RunOptions = RunOptions,
+                                Obs = Qobs[IndPeriod_Run,]),
+               regexp = "GRiwrmInputsModel")
+
+  expect_error(CreateInputsCrit(InputsModel = InputsModel,
+                                RunOptions = RunOptions[[1]],
+                                Obs = Qobs[IndPeriod_Run,]),
+               regexp = "GRiwrmRunOptions")
+  expect_error(CreateInputsCrit(InputsModel = InputsModel,
+                                RunOptions = RunOptions,
+                                Obs = 1),
+               regexp = "matrix or data.frame")
+})
+
+test_that("De Lavenne criterion is OK", {
+  IC <- CreateInputsCrit(InputsModel = InputsModel,
+                         RunOptions = RunOptions,
+                         Obs = Qobs[IndPeriod_Run,],
+                         AprioriIds = c("54057" = "54032", "54032" = "54001", "54001" = "54095"))
+  expect_s3_class(IC[["54057"]], "InputsCritLavenneFunction")
+  Lavenne_FUN <- attr(IC[["54057"]], "Lavenne_FUN")
+  IC57 <- Lavenne_FUN(ParamMichel[["54032"]], 0.9)
+  expect_s3_class(IC57, "InputsCrit")
+  expect_s3_class(IC57, "Compo")
+})
+
+test_that("De Lavenne criterion: wrong sub-catchment order should throw error", {
+  expect_error(
+    CreateInputsCrit(InputsModel = InputsModel,
+                     RunOptions = RunOptions,
+                     Obs = Qobs[IndPeriod_Run,],
+                     AprioriIds = c("54057" = "54032", "54032" = "54001", "54001" = "54029")),
+    regexp = "is not upstream the node"
+  )
+})

tests/testthat/test-RunModel.R

 # data set up
-
-data(Severn)
-
-# Format observation
-BasinsObs <- Severn$BasinsObs
-DatesR <- BasinsObs[[1]]$DatesR
-PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
-PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
-Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
-
-# Set network
-nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
-nodes$distance_downstream <- nodes$distance_downstream
-nodes$model <- "RunModel_GR4J"
-griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
-
-# Convert meteo data to SD (remove upstream areas)
-Precip <- ConvertMeteoSD(griwrm, PrecipTot)
-PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
-
-# Calibration parameters
-ParamMichel <- list(