diff --git a/.regressionignore b/.regressionignore
index a6a2adffef1b4fc3f0fba87f3b0acd671c965838..1dd20c9de2d828cd149545d7fbc7da4326b40af2 100644
--- a/.regressionignore
+++ b/.regressionignore
@@ -17,6 +17,12 @@ RunModel_GR2M NewTabSeries
RunModel_GR2M NewTimeFormat
RunModel_GR2M OutputsModel
RunModel_GR2M RunOptions
+RunModel_GR1A OutputsModel
+Calibration_Michel CalibOptions
+Calibration CalibOptions
+CreateCalibOptions CalibOptions
+
+# New version of the SeriesAggreg function
RunModel_GR2M TabSeries
RunModel_GR2M TimeFormat
SeriesAggreg BasinInfo
diff --git a/DESCRIPTION b/DESCRIPTION
index 67edaf4e9b4f3a67d5767765384f19687f0f6799..32468507acef67c47ed52d79dd8d14ae743a74cd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -21,7 +21,15 @@ Authors@R: c(
person("Audrey", "Valéry", role = c("ctb"))
)
Depends: R (>= 3.1.0)
-Suggests: knitr, rmarkdown, coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, Rmalschains, testthat, imputeTS
+Imports:
+ graphics,
+ grDevices,
+ stats,
+ utils
+Suggests:
+ knitr, rmarkdown,
+ coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, imputeTS, Rmalschains,
+ testthat
Description: Hydrological modelling tools developed at INRAE-Antony (HYCAR Research Unit, France). The package includes several conceptual rainfall-runoff models (GR4H, GR5H, GR4J, GR5J, GR6J, GR2M, GR1A), a snow accumulation and melt model (CemaNeige) and the associated functions for their calibration and evaluation. Use help(airGR) for package description and references.
License: GPL-2
URL: https://hydrogr.github.io/airGR/
diff --git a/NEWS.md b/NEWS.md
index 2517a356f14436c12717fdbae10f4cc0404b15ea..e05beedc8ab353471e66308793fd17e105b16815 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -17,7 +17,7 @@
____________________________________________________________________________________
-### 1.6.3.73 Release Notes (2020-11-24)
+### 1.6.3.78 Release Notes (2021-01-05)
#### New features
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index aaad7c8a6160000a934cd163173bc175f7d65741..cfa85ac63c9c52b86e691caa3f7dfa46b3e4fcba 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -1,29 +1,34 @@
-Calibration_Michel <- function(InputsModel,
- RunOptions,
- InputsCrit,
+Calibration_Michel <- function(InputsModel,
+ RunOptions,
+ InputsCrit,
CalibOptions,
- FUN_MOD,
+ FUN_MOD,
FUN_CRIT, # deprecated
- FUN_TRANSFO = NULL,
+ FUN_TRANSFO = NULL,
verbose = TRUE) {
-
-
+
+
FUN_MOD <- match.fun(FUN_MOD)
if (!missing(FUN_CRIT)) {
FUN_CRIT <- match.fun(FUN_CRIT)
}
+
+ # Handling 'FUN_TRANSFO' from direct argument or provided by 'CaliOptions'
if (!is.null(FUN_TRANSFO)) {
FUN_TRANSFO <- match.fun(FUN_TRANSFO)
+ } else if(!is.null(CalibOptions$FUN_TRANSFO)) {
+ FUN_TRANSFO <- CalibOptions$FUN_TRANSFO
+ } else {
+ stop("'FUN_TRANSFO' is not provided neither as 'FUN_TRANSFO' argument or in 'CaliOptions' argument")
}
-
-
+
##_____Arguments_check_____________________________________________________________________
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
- }
+ }
if (!inherits(RunOptions, "RunOptions")) {
stop("'RunOptions' must be of class 'RunOptions'")
- }
+ }
if (!inherits(InputsCrit, "InputsCrit")) {
stop("'InputsCrit' must be of class 'InputsCrit'")
}
@@ -46,151 +51,15 @@ Calibration_Michel <- function(InputsModel,
}
if (!inherits(CalibOptions, "CalibOptions")) {
stop("'CalibOptions' must be of class 'CalibOptions'")
- }
+ }
if (!inherits(CalibOptions, "HBAN")) {
stop("'CalibOptions' must be of class 'HBAN' if 'Calibration_Michel' is used")
}
if (!missing(FUN_CRIT)) {
warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object")
- }
-
-
- ##_check_FUN_TRANSFO
- if (is.null(FUN_TRANSFO)) {
- if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
- FUN1 <- TransfoParam_GR4H
- }
- if (identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
- FUN1 <- TransfoParam_GR5H
- }
- if (identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
- FUN1 <- TransfoParam_GR4J
- }
- if (identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
- FUN1 <- TransfoParam_GR5J
- }
- if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- FUN1 <- TransfoParam_GR6J
- }
- if (identical(FUN_MOD, RunModel_GR2M )) {
- FUN1 <- TransfoParam_GR2M
- }
- if (identical(FUN_MOD, RunModel_GR1A )) {
- FUN1 <- TransfoParam_GR1A
- }
- ##_set_FUN_LAG
- if (inherits(CalibOptions, "SD")) {
- FUN_LAG <- TransfoParam_Lag
- }
- if (identical(FUN_MOD, RunModel_CemaNeige )) {
- if (inherits(CalibOptions, "hysteresis")) {
- FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
- } else {
- FUN_TRANSFO <- TransfoParam_CemaNeige
- }
- }
- if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H) |
- identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_GR6J)) {
- if (!inherits(CalibOptions, "SD")) {
- FUN_TRANSFO <- FUN1
- } else {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 2:NParam] <- FUN1(ParamIn[, 2:NParam], Direction)
- ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H) |
- identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- if (inherits(CalibOptions, "hysteresis")) {
- FUN2 <- TransfoParam_CemaNeigeHyst
- } else {
- FUN2 <- TransfoParam_CemaNeige
- }
- if (inherits(CalibOptions, "hysteresis") & !inherits(CalibOptions, "SD")) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam-4)] <- FUN1(ParamIn[, 1:(NParam-4)], Direction)
- ParamOut[, (NParam-3):NParam ] <- FUN2(ParamIn[, (NParam-3):NParam ], Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- if (!inherits(CalibOptions, "hysteresis") & !inherits(CalibOptions, "SD")) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam-2)] <- FUN1(ParamIn[, 1:(NParam-2)], Direction)
- ParamOut[, (NParam-1):NParam ] <- FUN2(ParamIn[, (NParam-1):NParam ], Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- if (inherits(CalibOptions, "hysteresis") & inherits(CalibOptions, "SD")) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 2:(NParam-4)] <- FUN1( ParamIn[, 2:(NParam-4)], Direction)
- ParamOut[, (NParam-3):NParam] <- FUN2( ParamIn[, (NParam-3):NParam], Direction)
- ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1 ]), Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
- if (!inherits(CalibOptions, "hysteresis") & inherits(CalibOptions, "SD")) {
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- if (!Bool) {
- ParamIn <- rbind(ParamIn)
- }
- ParamOut <- NA * ParamIn
- NParam <- ncol(ParamIn)
- ParamOut[, 2:(NParam-2)] <- FUN1( ParamIn[, 2:(NParam-2)], Direction)
- ParamOut[, (NParam-1):NParam] <- FUN2( ParamIn[, (NParam-1):NParam], Direction)
- ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1 ]), Direction)
- if (!Bool) {
- ParamOut <- ParamOut[1, ]
- }
- return(ParamOut)
- }
- }
-
- }
- if (is.null(FUN_TRANSFO)) {
- stop("'FUN_TRANSFO' was not found (in 'Calibration' function)")
- }
}
-
- ##_variables_initialisation
+
+ ##_variables_initialisation
ParamFinalR <- NULL
ParamFinalT <- NULL
CritFinal <- NULL
@@ -219,20 +88,20 @@ Calibration_Michel <- function(InputsModel,
CritOptim <- +1e100
##_temporary_change_of_Outputs_Sim
RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
-
-
-
+
+
+
##_____Parameter_Grid_Screening____________________________________________________________
-
-
+
+
##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
ProposeCandidatesGrid <- function(DistribParam) {
NewCandidates <- expand.grid(lapply(seq_len(ncol(DistribParamR)), function(x) DistribParam[, x]))
NewCandidates <- unique(NewCandidates) # to avoid duplicates when a parameter is set
Output <- list(NewCandidates = NewCandidates)
- }
-
-
+ }
+
+
##Creation_of_new_candidates_______________________________________________
OptimParam <- is.na(CalibOptions$FixedParam)
if (PrefilteringType == 1) {
@@ -253,7 +122,7 @@ Calibration_Michel <- function(InputsModel,
} else {
CandidatesParamR <- cbind(CandidatesParamR)
}
-
+
##Loop_to_test_the_various_candidates______________________________________
iNewOptim <- 0
Ncandidates <- nrow(CandidatesParamR)
@@ -272,12 +141,12 @@ Calibration_Michel <- function(InputsModel,
if (iNew == round(k / 10 * Ncandidates)) {
message(" ", 10 * k, "%", appendLF = FALSE)
}
- }
+ }
}
##Model_run
Param <- CandidatesParamR[iNew, ]
OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD)
-
+
##Calibration_criterion_computation
OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
if (!is.na(OutputsCrit$CritValue)) {
@@ -296,8 +165,8 @@ Calibration_Michel <- function(InputsModel,
if (verbose & Ncandidates > 1) {
message(" 100%)\n", appendLF = FALSE)
}
-
-
+
+
##End_of_first_step_Parameter_Screening____________________________________
ParamStartR <- CandidatesParamR[iNewOptim, ]
if (!is.matrix(ParamStartR)) {
@@ -320,13 +189,13 @@ Calibration_Michel <- function(InputsModel,
HistParamR[1, ] <- ParamStartR
HistParamT[1, ] <- ParamStartT
HistCrit[1, ] <- CritStart
-
-
-
-
+
+
+
+
##_____Steepest_Descent_Local_Search_______________________________________________________
-
-
+
+
##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace) {
##Format_checking
@@ -377,11 +246,11 @@ Calibration_Michel <- function(InputsModel,
Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
return(Output)
}
-
-
+
+
##Initialisation_of_variables
if (verbose) {
- message("Steepest-descent local search in progress")
+ message("Steepest-descent local search in progress")
}
Pace <- 0.64
PaceDiag <- rep(0, NParam)
@@ -393,18 +262,18 @@ Calibration_Michel <- function(InputsModel,
RangesT <- FUN_TRANSFO(RangesR, "RT")
NewParamOptimT <- ParamStartT
OldParamOptimT <- ParamStartT
-
-
+
+
##START_LOOP_ITER_________________________________________________________
for (ITER in 1:(100 * NParam)) {
-
-
+
+
##Exit_loop_when_Pace_becomes_too_small___________________________________
if (Pace < 0.01) {
break
}
-
-
+
+
##Creation_of_new_candidates______________________________________________
CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
@@ -418,8 +287,8 @@ Calibration_Michel <- function(InputsModel,
} else {
CandidatesParamR <- cbind(CandidatesParamR)
}
-
-
+
+
##Loop_to_test_the_various_candidates_____________________________________
iNewOptim <- 0
for (iNew in 1:nrow(CandidatesParamR)) {
@@ -427,7 +296,7 @@ Calibration_Michel <- function(InputsModel,
Param <- CandidatesParamR[iNew, ]
OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD)
##Calibration_criterion_computation
- OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
+ OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
if (!is.na(OutputsCrit$CritValue)) {
if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
@@ -436,8 +305,8 @@ Calibration_Michel <- function(InputsModel,
}
}
NRuns <- NRuns + nrow(CandidatesParamR)
-
-
+
+
##When_a_progress_has_been_achieved_______________________________________
if (iNewOptim != 0) {
##We_store_the_optimal_set
@@ -452,7 +321,7 @@ Calibration_Michel <- function(InputsModel,
##We_update_PaceDiag
VectPace <- NewParamOptimT-OldParamOptimT
for (iC in 1:NParam) {
- if (OptimParam[iC]) {
+ if (OptimParam[iC]) {
PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
}
}
@@ -461,8 +330,8 @@ Calibration_Michel <- function(InputsModel,
Pace <- Pace / 2
Compt <- 0
}
-
-
+
+
##Test_of_an_additional_candidate_using_diagonal_progress_________________
if (ITER > 4 * NParam) {
NRuns <- NRuns + 1
@@ -498,34 +367,34 @@ Calibration_Michel <- function(InputsModel,
OldParamOptimT <- NewParamOptimT
NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
}
-
+
}
-
-
+
+
##Results_archiving_______________________________________________________
NewParamOptimR <- FUN_TRANSFO(NewParamOptimT, "TR")
HistParamR[ITER+1, ] <- NewParamOptimR
HistParamT[ITER+1, ] <- NewParamOptimT
HistCrit[ITER+1, ] <- CritOptim
### if (verbose) { cat(paste("\t Iter ",formatC(ITER,format="d",width=3), " Crit ",formatC(CritOptim,format="f",digits=4), " Pace ",formatC(Pace,format="f",digits=4), "\n",sep=""))}
-
-
-
+
+
+
} ##END_LOOP_ITER_________________________________________________________
ITER <- ITER - 1
-
-
+
+
##Case_when_the_starting_parameter_set_remains_the_best_solution__________
- if (CritOptim == CritStart & verbose) {
+ if (CritOptim == CritStart & verbose) {
message("\t No progress achieved")
}
-
+
##End_of_Steepest_Descent_Local_Search____________________________________
ParamFinalR <- NewParamOptimR
ParamFinalT <- NewParamOptimT
CritFinal <- CritOptim
NIter <- 1 + ITER
- if (verbose) {
+ if (verbose) {
message(sprintf("\t Calibration completed (%s iterations, %s runs)", NIter, NRuns))
message("\t Param = ", paste(sprintf("%8.3f", ParamFinalR), collapse = ", "))
message(sprintf("\t Crit. %-12s = %.4f", CritName, CritFinal * Multiplier))
@@ -547,13 +416,13 @@ Calibration_Michel <- function(InputsModel,
colnames(HistParamT) <- paste0("Param", 1:NParam)
HistCrit <- cbind(HistCrit[1:NIter, ])
###colnames(HistCrit) <- paste("HistCrit")
-
+
BoolCrit_Actual <- InputsCrit$BoolCrit
BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
-
-
+
+
##_____Output______________________________________________________________________________
OutputsCalib <- list(ParamFinalR = as.double(ParamFinalR), CritFinal = CritFinal * Multiplier,
NIter = NIter, NRuns = NRuns,
@@ -562,7 +431,7 @@ Calibration_Michel <- function(InputsModel,
CritName = CritName, CritBestValue = CritBestValue)
class(OutputsCalib) <- c("OutputsCalib", "HBAN")
return(OutputsCalib)
-
-
-
+
+
+
}
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 5815beb7b86c6306c09b4e9b3959660c37866ec7..4d0edcc2a00705e3f4a02cb8165fa0e44bb59280 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -105,46 +105,46 @@ CreateCalibOptions <- function(FUN_MOD,
##_set_FUN1
if (identical(FUN_MOD, RunModel_GR4H) |
identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
- FUN1 <- TransfoParam_GR4H
+ FUN_GR <- TransfoParam_GR4H
}
if (identical(FUN_MOD, RunModel_GR5H) |
identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
- FUN1 <- TransfoParam_GR5H
+ FUN_GR <- TransfoParam_GR5H
}
if (identical(FUN_MOD, RunModel_GR4J) |
identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
- FUN1 <- TransfoParam_GR4J
+ FUN_GR <- TransfoParam_GR4J
}
if (identical(FUN_MOD, RunModel_GR5J) |
identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
- FUN1 <- TransfoParam_GR5J
+ FUN_GR <- TransfoParam_GR5J
}
if (identical(FUN_MOD, RunModel_GR6J) |
identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- FUN1 <- TransfoParam_GR6J
+ FUN_GR <- TransfoParam_GR6J
}
if (identical(FUN_MOD, RunModel_GR2M)) {
- FUN1 <- TransfoParam_GR2M
+ FUN_GR <- TransfoParam_GR2M
}
if (identical(FUN_MOD, RunModel_GR1A)) {
- FUN1 <- TransfoParam_GR1A
+ FUN_GR <- TransfoParam_GR1A
}
if (identical(FUN_MOD, RunModel_CemaNeige)) {
if (IsHyst) {
- FUN1 <- TransfoParam_CemaNeigeHyst
+ FUN_GR <- TransfoParam_CemaNeigeHyst
} else {
- FUN1 <- TransfoParam_CemaNeige
+ FUN_GR <- TransfoParam_CemaNeige
}
}
- if (is.null(FUN1)) {
- stop("'FUN1' was not found")
+ if (is.null(FUN_GR)) {
+ stop("'FUN_GR' was not found")
return(NULL)
}
##_set_FUN2
if (IsHyst) {
- FUN2 <- TransfoParam_CemaNeigeHyst
+ FUN_NEIGE <- TransfoParam_CemaNeigeHyst
} else {
- FUN2 <- TransfoParam_CemaNeige
+ FUN_NEIGE <- TransfoParam_CemaNeige
}
##_set_FUN_LAG
if (IsSD) {
@@ -153,7 +153,7 @@ CreateCalibOptions <- function(FUN_MOD,
##_set_FUN_TRANSFO
if (sum(ObjectClass %in% c("GR4H", "GR5H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige")) > 0) {
if (!IsSD) {
- FUN_TRANSFO <- FUN1
+ FUN_TRANSFO <- FUN_GR
} else {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
@@ -162,7 +162,7 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 2:NParam] <- FUN1(ParamIn[, 2:NParam], Direction)
+ ParamOut[, 2:NParam] <- FUN_GR(ParamIn[, 2:NParam], Direction)
ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
@@ -179,8 +179,8 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam - 4) ] <- FUN1(ParamIn[, 1:(NParam - 4) ], Direction)
- ParamOut[, (NParam - 3):NParam] <- FUN2(ParamIn[, (NParam - 3):NParam], Direction)
+ ParamOut[, 1:(NParam - 4) ] <- FUN_GR(ParamIn[, 1:(NParam - 4) ], Direction)
+ ParamOut[, (NParam - 3):NParam] <- FUN_NEIGE(ParamIn[, (NParam - 3):NParam], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -196,11 +196,11 @@ CreateCalibOptions <- function(FUN_MOD,
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
if (NParam <= 3) {
- ParamOut[, 1:(NParam - 2)] <- FUN1(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
+ ParamOut[, 1:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
} else {
- ParamOut[, 1:(NParam - 2)] <- FUN1( ParamIn[, 1:(NParam - 2)], Direction)
+ ParamOut[, 1:(NParam - 2)] <- FUN_GR( ParamIn[, 1:(NParam - 2)], Direction)
}
- ParamOut[, (NParam - 1):NParam] <- FUN2(ParamIn[, (NParam - 1):NParam], Direction)
+ ParamOut[, (NParam - 1):NParam] <- FUN_NEIGE(ParamIn[, (NParam - 1):NParam], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -215,8 +215,8 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 2:(NParam - 4) ] <- FUN1( ParamIn[, 2:(NParam - 4) ], Direction)
- ParamOut[, (NParam - 3):NParam] <- FUN2( ParamIn[, (NParam - 3):NParam], Direction)
+ ParamOut[, 2:(NParam - 4) ] <- FUN_GR( ParamIn[, 2:(NParam - 4) ], Direction)
+ ParamOut[, (NParam - 3):NParam] <- FUN_NEIGE( ParamIn[, (NParam - 3):NParam], Direction)
ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1 ]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
@@ -233,11 +233,11 @@ CreateCalibOptions <- function(FUN_MOD,
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
if (NParam <= 3) {
- ParamOut[, 2:(NParam - 2)] <- FUN1(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
+ ParamOut[, 2:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
} else {
- ParamOut[, 2:(NParam - 2)] <- FUN1( ParamIn[, 2:(NParam - 2)], Direction)
+ ParamOut[, 2:(NParam - 2)] <- FUN_GR( ParamIn[, 2:(NParam - 2)], Direction)
}
- ParamOut[, (NParam - 1):NParam] <- FUN2( ParamIn[, (NParam - 1):NParam], Direction)
+ ParamOut[, (NParam - 1):NParam] <- FUN_NEIGE( ParamIn[, (NParam - 1):NParam], Direction)
ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
@@ -472,7 +472,7 @@ CreateCalibOptions <- function(FUN_MOD,
##Create_CalibOptions
- CalibOptions <- list(FixedParam = FixedParam, SearchRanges = SearchRanges)
+ CalibOptions <- list(FixedParam = FixedParam, SearchRanges = SearchRanges, FUN_TRANSFO = FUN_TRANSFO)
if (!is.null(StartParamList)) {
CalibOptions <- c(CalibOptions, list(StartParamList = StartParamList))
diff --git a/R/CreateIniStates.R b/R/CreateIniStates.R
index d3a6c9481cf7bbe77cecb536fbe7d1345154980b..6b311d414f2e3f5fd486056b5ace76e0b605c724 100644
--- a/R/CreateIniStates.R
+++ b/R/CreateIniStates.R
@@ -3,17 +3,18 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
+ SD = NULL,
verbose = TRUE) {
-
-
+
+
ObjectClass <- NULL
-
+
UH1n <- 20L
UH2n <- UH1n * 2L
-
+
nameFUN_MOD <- as.character(substitute(FUN_MOD))
FUN_MOD <- match.fun(FUN_MOD)
-
+
## check FUN_MOD
BOOL <- FALSE
if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H)) {
@@ -56,7 +57,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & IsIntStore) {
stop("'IsIntStore' cannot be TRUE if GR5H is not used in 'FUN_MOD'")
}
-
+
## check InputsModel
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
@@ -68,13 +69,13 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
!inherits(InputsModel, "CemaNeige")) {
stop("'InputsModel' must be of class 'CemaNeige'")
}
-
-
+
+
## check states
if (any(eTGCemaNeigeLayers > 0)) {
stop("Positive values are not allowed for 'eTGCemaNeigeLayers'")
- }
-
+ }
+
if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
if (is.null(ExpStore)) {
stop("'RunModel_*GR6J' need an 'ExpStore' value")
@@ -85,7 +86,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
}
ExpStore <- Inf
}
-
+
if (identical(FUN_MOD, RunModel_GR2M)) {
if (!is.null(UH1)) {
if (verbose) {
@@ -100,20 +101,20 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
UH2 <- rep(Inf, UH2n)
}
}
-
+
if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !is.null(UH1)) {
if (verbose) {
warning(sprintf("'%s' does not require 'UH1'. Values set to NA", nameFUN_MOD))
}
UH1 <- rep(Inf, UH1n)
- }
+ }
if ((!identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & !is.null(IntStore)) {
if (verbose) {
warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", nameFUN_MOD))
}
IntStore <- Inf
}
-
+
if ("CemaNeige" %in% ObjectClass & ! "GR" %in% ObjectClass) {
if (!is.null(ProdStore)) {
if (verbose) {
@@ -170,7 +171,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", nameFUN_MOD))
}
GthrCemaNeigeLayers <- Inf
- GlocmaxCemaNeigeLayers <- Inf
+ GlocmaxCemaNeigeLayers <- Inf
}
if(!"CemaNeige" %in% ObjectClass &
(!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers) | !is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
@@ -180,24 +181,24 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
GCemaNeigeLayers <- Inf
eTGCemaNeigeLayers <- Inf
GthrCemaNeigeLayers <- Inf
- GlocmaxCemaNeigeLayers <- Inf
+ GlocmaxCemaNeigeLayers <- Inf
}
-
-
+
+
## set states
if("CemaNeige" %in% ObjectClass) {
NLayers <- length(InputsModel$LayerPrecip)
} else {
NLayers <- 1
}
-
-
+
+
## manage NULL values
if (is.null(ExpStore)) {
- ExpStore <- Inf
+ ExpStore <- Inf
}
if (is.null(IntStore)) {
- IntStore <- Inf
+ IntStore <- Inf
}
if (is.null(UH1)) {
if ("hourly" %in% ObjectClass) {
@@ -232,16 +233,16 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
}
if (any(is.infinite(GlocmaxCemaNeigeLayers))) {
GlocmaxCemaNeigeLayers <- rep(Inf, NLayers)
- }
-
+ }
+
# check negative values
if (any(ProdStore < 0) | any(RoutStore < 0) | any(IntStore < 0) |
any(UH1 < 0) | any(UH2 < 0) |
any(GCemaNeigeLayers < 0)) {
stop("Negative values are not allowed for any of 'ProdStore', 'RoutStore', 'IntStore', 'UH1', 'UH2', 'GCemaNeigeLayers'")
}
-
-
+
+
## check length
if (!is.numeric(ProdStore) || length(ProdStore) != 1L) {
stop("'ProdStore' must be numeric of length one")
@@ -281,7 +282,23 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
}
}
-
+
+ # SD model state handling
+ if(!is.null(SD)) {
+ if(!inherits(InputsModel, "SD")) {
+ stop("'SD' argument provided and 'InputsModel' is not of class 'SD'")
+ }
+ if(!is.list(SD)) {
+ stop("'SD' argument must be a list")
+ }
+ lapply(SD, function(x) {
+ if(!is.numeric(x)) stop("Each item of 'SD' list argument must be numeric")
+ })
+ if(length(SD) != length(InputsModel$LengthHydro)) {
+ stop("Number of items of 'SD' list argument must be the same as the number of upstream connections",
+ sprintf(" (%i required, found %i)", length(InputsModel$LengthHydro), length(SD)))
+ }
+ }
## format output
IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore, Int = IntStore),
@@ -291,7 +308,11 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
IniStatesNA <- unlist(IniStates)
IniStatesNA[is.infinite(IniStatesNA)] <- NA
IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
-
+
+ if(!is.null(SD)) {
+ IniStatesNA$SD <- SD
+ }
+
class(IniStatesNA) <- c("IniStates", ObjectClass)
if(IsHyst) {
class(IniStatesNA) <- c(class(IniStatesNA), "hysteresis")
@@ -299,8 +320,8 @@ CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = F
if(IsIntStore) {
class(IniStatesNA) <- c(class(IniStatesNA), "interception")
}
-
+
return(IniStatesNA)
-
-
+
+
}
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index d82a0241ec42abb8a607d291d60f588885f8defb..60c32a26898d8fa164b78d46d02cce8e7267419b 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -192,7 +192,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
## Output data preparation
## OutputsModel only
- if (!ExportDatesR== FALSE & !ExportStateEnd) {
+ if (!ExportDatesR & !ExportStateEnd) {
OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
list(CemaNeigeLayers))
names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod], NameCemaNeigeLayers)
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 29f7a475f7a08f2f812af9a64e1dfb4533d3a3d9..8bfd1f568e73a5a259545047926aaa151a52cfc9 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -1,14 +1,13 @@
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'")
}
@@ -19,23 +18,30 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
}
if (is.null(InputsModel$OutputsModel)) {
- stop("'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment")
+ stop(
+ "'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment"
+ )
}
if (is.null(InputsModel$OutputsModel$Qsim)) {
- stop("'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
+ stop(
+ "'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment"
+ )
}
if (sum(!is.na(InputsModel$OutputsModel$Qsim)) != length(RunOptions$IndPeriod_Run)) {
- stop("'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA")
+ stop(
+ "'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA"
+ )
}
-
+
OutputsModel <- InputsModel$OutputsModel
OutputsModel$QsimDown <- OutputsModel$Qsim
-
- if (inherits(InputsModel, "daily")) {
- TimeStep <- 60 * 60 * 24
- }
+
if (inherits(InputsModel, "hourly")) {
TimeStep <- 60 * 60
+ } else if (inherits(InputsModel, "daily")) {
+ TimeStep <- 60 * 60 * 24
+ } else {
+ stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
}
# propagation time from upstream meshes to outlet
@@ -44,19 +50,46 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
NbUpBasins <- length(InputsModel$LengthHydro)
LengthTs <- length(OutputsModel$QsimDown)
- OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+ OutputsModel$Qsim <-
+ OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+
+ IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
+ if (length(IniSD) > 0) {
+ if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
+ stop(
+ sprintf(
+ "SD initial states has a length of %i and a length of %i is required",
+ length(IniSD),
+ sum(floor(PT)) + NbUpBasins
+ )
+ )
+ }
+ 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]
+ Qupstream <-
+ InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
# Upstream flow with area needs to be converted to m3 by time step
- Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
+ 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 +100,12 @@ 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[x])):LengthTs]
+ })
+ }
+
return(OutputsModel)
-}
\ No newline at end of file
+}
diff --git a/man/CreateIniStates.Rd b/man/CreateIniStates.Rd
index 5e2fed8357751f995973ab695373951c5a90d669..76a5d9f3b760ac97e5e0fb729c756973c429baef 100644
--- a/man/CreateIniStates.Rd
+++ b/man/CreateIniStates.Rd
@@ -19,7 +19,7 @@ CreateIniStates(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = FALSE,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
- verbose = TRUE)
+ SD = NULL, verbose = TRUE)
}
@@ -52,6 +52,8 @@ CreateIniStates(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = FALSE,
\item{GlocmaxCemaNeigeLayers}{(optional) [numeric] local melt threshold for hysteresis [mm], possibly used to create the CemaNeige model initial state in case the Linear Hysteresis version is used}
+\item{SD}{(optional) [list] of [numeric] states of delayed upstream flows for semi-distributed models, the nature of the state and the unit depend on the model and the unit of the upstream flow}
+
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
}
diff --git a/man/RunModel_CemaNeige.Rd b/man/RunModel_CemaNeige.Rd
index f430f8136a1b0072d279a31b926503b238d953a7..0d0afd20b223563fb7399e31008e1ecf145bc93d 100644
--- a/man/RunModel_CemaNeige.Rd
+++ b/man/RunModel_CemaNeige.Rd
@@ -40,7 +40,7 @@ CemaNeige X4 \tab (optional) percentage (between 0 and 1) of annual snowfall def
\emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
\emph{$CemaNeigeLayers[[iLayer]]$Pliq } \tab [numeric] series of liquid precip. [mm/time step] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Psol } \tab [numeric] series of solid precip. [mm/time step] \cr
- \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack [mm] \cr
+ \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack (snow water equivalent) [mm] \cr
\emph{$CemaNeigeLayers[[iLayer]]$ThermalState} \tab [numeric] series of snow pack thermal state [°C] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Gratio } \tab [numeric] series of Gratio [0-1] \cr
\emph{$CemaNeigeLayers[[iLayer]]$PotMelt } \tab [numeric] series of potential snow melt [mm/time step] \cr
diff --git a/man/RunModel_CemaNeigeGR4H.Rd b/man/RunModel_CemaNeigeGR4H.Rd
index e97b5cd30742cad3d04b77180fafc502cdbf3423..0204920e38388cd0ca1048d9dc91d2516a8d9d92 100644
--- a/man/RunModel_CemaNeigeGR4H.Rd
+++ b/man/RunModel_CemaNeigeGR4H.Rd
@@ -40,29 +40,29 @@ CemaNeige X4 \tab (optional) percentage (between 0 and 1) of annual snowfall def
\value{
[list] list containing the function outputs organised as follows:
\tabular{ll}{
- \emph{$DatesR } \tab [POSIXlt] series of dates \cr
- \emph{$PotEvap } \tab [numeric] series of input potential evapotranspiration [mm/h] \cr
- \emph{$Precip } \tab [numeric] series of input total precipitation [mm/h] \cr
- \emph{$Prod } \tab [numeric] series of production store level [mm] \cr
- \emph{$Pn } \tab [numeric] series of net rainfall [mm/h] \cr
- \emph{$Ps } \tab [numeric] series of the part of Pn filling the production store [mm/h] \cr
- \emph{$AE } \tab [numeric] series of actual evapotranspiration [mm/h] \cr
- \emph{$Perc } \tab [numeric] series of percolation (PERC) [mm/h] \cr
- \emph{$PR } \tab [numeric] series of PR=Pn-Ps+Perc [mm/h] \cr
- \emph{$Q9 } \tab [numeric] series of UH1 outflow (Q9) [mm/h] \cr
- \emph{$Q1 } \tab [numeric] series of UH2 outflow (Q1) [mm/h] \cr
- \emph{$Rout } \tab [numeric] series of routing store level [mm] \cr
- \emph{$Exch } \tab [numeric] series of potential semi-exchange between catchments [mm/h] \cr
- \emph{$AExch1 } \tab [numeric] series of actual exchange between catchments for branch 1 [mm/h] \cr
- \emph{$AExch2 } \tab [numeric] series of actual exchange between catchments for branch 2 [mm/h] \cr
- \emph{$AExch } \tab [numeric] series of actual exchange between catchments (1+2) [mm/h] \cr
- \emph{$QR } \tab [numeric] series of routing store outflow (QR) [mm/h] \cr
- \emph{$QD } \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/h] \cr
- \emph{$Qsim } \tab [numeric] series of simulated discharge [mm/h] \cr
- \emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
+ \emph{$DatesR } \tab [POSIXlt] series of dates \cr
+ \emph{$PotEvap } \tab [numeric] series of input potential evapotranspiration [mm/h] \cr
+ \emph{$Precip } \tab [numeric] series of input total precipitation [mm/h] \cr
+ \emph{$Prod } \tab [numeric] series of production store level [mm] \cr
+ \emph{$Pn } \tab [numeric] series of net rainfall [mm/h] \cr
+ \emph{$Ps } \tab [numeric] series of the part of Pn filling the production store [mm/h] \cr
+ \emph{$AE } \tab [numeric] series of actual evapotranspiration [mm/h] \cr
+ \emph{$Perc } \tab [numeric] series of percolation (PERC) [mm/h] \cr
+ \emph{$PR } \tab [numeric] series of PR=Pn-Ps+Perc [mm/h] \cr
+ \emph{$Q9 } \tab [numeric] series of UH1 outflow (Q9) [mm/h] \cr
+ \emph{$Q1 } \tab [numeric] series of UH2 outflow (Q1) [mm/h] \cr
+ \emph{$Rout } \tab [numeric] series of routing store level [mm] \cr
+ \emph{$Exch } \tab [numeric] series of potential semi-exchange between catchments [mm/h] \cr
+ \emph{$AExch1 } \tab [numeric] series of actual exchange between catchments for branch 1 [mm/h] \cr
+ \emph{$AExch2 } \tab [numeric] series of actual exchange between catchments for branch 2 [mm/h] \cr
+ \emph{$AExch } \tab [numeric] series of actual exchange between catchments (1+2) [mm/h] \cr
+ \emph{$QR } \tab [numeric] series of routing store outflow (QR) [mm/h] \cr
+ \emph{$QD } \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/h] \cr
+ \emph{$Qsim } \tab [numeric] series of simulated discharge [mm/h] \cr
+ \emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
\emph{$CemaNeigeLayers[[iLayer]]$Pliq } \tab [numeric] series of liquid precip. [mm/h] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Psol } \tab [numeric] series of solid precip. [mm/h] \cr
- \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack [mm] \cr
+ \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack (snow water equivalent)[mm] \cr
\emph{$CemaNeigeLayers[[iLayer]]$ThermalState } \tab [numeric] series of snow pack thermal state [°C] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Gratio } \tab [numeric] series of Gratio [0-1] \cr
\emph{$CemaNeigeLayers[[iLayer]]$PotMelt } \tab [numeric] series of potential snow melt [mm/h] \cr
diff --git a/man/RunModel_CemaNeigeGR4J.Rd b/man/RunModel_CemaNeigeGR4J.Rd
index 99e09453d52bf998bd3530c674eeb0ae513446d6..7c2aef74403789f3da96d4b1ce174d9e1e2092f3 100644
--- a/man/RunModel_CemaNeigeGR4J.Rd
+++ b/man/RunModel_CemaNeigeGR4J.Rd
@@ -40,29 +40,29 @@ CemaNeige X4 \tab (optional) percentage (between 0 and 1) of annual snowfall def
\value{
[list] list containing the function outputs organised as follows:
\tabular{ll}{
- \emph{$DatesR } \tab [POSIXlt] series of dates \cr
- \emph{$PotEvap } \tab [numeric] series of input potential evapotranspiration [mm/d] \cr
- \emph{$Precip } \tab [numeric] series of input total precipitation [mm/d] \cr
- \emph{$Prod } \tab [numeric] series of production store level [mm] \cr
- \emph{$Pn } \tab [numeric] series of net rainfall [mm/d] \cr
- \emph{$Ps } \tab [numeric] series of the part of Pn filling the production store [mm/d] \cr
- \emph{$AE } \tab [numeric] series of actual evapotranspiration [mm/d] \cr
- \emph{$Perc } \tab [numeric] series of percolation (PERC) [mm/d] \cr
- \emph{$PR } \tab [numeric] series of PR=Pn-Ps+Perc [mm/d] \cr
- \emph{$Q9 } \tab [numeric] series of UH1 outflow (Q9) [mm/d] \cr
- \emph{$Q1 } \tab [numeric] series of UH2 outflow (Q1) [mm/d] \cr
- \emph{$Rout } \tab [numeric] series of routing store level [mm] \cr
- \emph{$Exch } \tab [numeric] series of potential semi-exchange between catchments [mm/d] \cr
- \emph{$AExch1 } \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d] \cr
- \emph{$AExch2 } \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d] \cr
- \emph{$AExch } \tab [numeric] series of actual exchange between catchments (1+2) [mm/d] \cr
- \emph{$QR } \tab [numeric] series of routing store outflow (QR) [mm/d] \cr
- \emph{$QD } \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d] \cr
- \emph{$Qsim } \tab [numeric] series of simulated discharge [mm/d] \cr
- \emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
+ \emph{$DatesR } \tab [POSIXlt] series of dates \cr
+ \emph{$PotEvap } \tab [numeric] series of input potential evapotranspiration [mm/d] \cr
+ \emph{$Precip } \tab [numeric] series of input total precipitation [mm/d] \cr
+ \emph{$Prod } \tab [numeric] series of production store level [mm] \cr
+ \emph{$Pn } \tab [numeric] series of net rainfall [mm/d] \cr
+ \emph{$Ps } \tab [numeric] series of the part of Pn filling the production store [mm/d] \cr
+ \emph{$AE } \tab [numeric] series of actual evapotranspiration [mm/d] \cr
+ \emph{$Perc } \tab [numeric] series of percolation (PERC) [mm/d] \cr
+ \emph{$PR } \tab [numeric] series of PR=Pn-Ps+Perc [mm/d] \cr
+ \emph{$Q9 } \tab [numeric] series of UH1 outflow (Q9) [mm/d] \cr
+ \emph{$Q1 } \tab [numeric] series of UH2 outflow (Q1) [mm/d] \cr
+ \emph{$Rout } \tab [numeric] series of routing store level [mm] \cr
+ \emph{$Exch } \tab [numeric] series of potential semi-exchange between catchments [mm/d] \cr
+ \emph{$AExch1 } \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d] \cr
+ \emph{$AExch2 } \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d] \cr
+ \emph{$AExch } \tab [numeric] series of actual exchange between catchments (1+2) [mm/d] \cr
+ \emph{$QR } \tab [numeric] series of routing store outflow (QR) [mm/d] \cr
+ \emph{$QD } \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d] \cr
+ \emph{$Qsim } \tab [numeric] series of simulated discharge [mm/d] \cr
+ \emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
\emph{$CemaNeigeLayers[[iLayer]]$Pliq } \tab [numeric] series of liquid precip. [mm/d] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Psol } \tab [numeric] series of solid precip. [mm/d] \cr
- \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack [mm] \cr
+ \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack (snow water equivalent) [mm] \cr
\emph{$CemaNeigeLayers[[iLayer]]$ThermalState } \tab [numeric] series of snow pack thermal state [°C] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Gratio } \tab [numeric] series of Gratio [0-1] \cr
\emph{$CemaNeigeLayers[[iLayer]]$PotMelt } \tab [numeric] series of potential snow melt [mm/d] \cr
diff --git a/man/RunModel_CemaNeigeGR5H.Rd b/man/RunModel_CemaNeigeGR5H.Rd
index f7059f4c36606f043be10219e9dcfd6035706eb3..3e9353a01989b05a138928d94b07e6947dad7df3 100644
--- a/man/RunModel_CemaNeigeGR5H.Rd
+++ b/man/RunModel_CemaNeigeGR5H.Rd
@@ -65,7 +65,7 @@ CemaNeige X4 \tab (optional) percentage (between 0 and 1) of annual snowfall def
\emph{$Qsim } \tab [numeric] series of simulated discharge [mm/h] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Pliq } \tab [numeric] series of liquid precip. [mm/h] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Psol } \tab [numeric] series of solid precip. [mm/h] \cr
- \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack [mm] \cr
+ \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack (snow water equivalent) [mm] \cr
\emph{$CemaNeigeLayers[[iLayer]]$ThermalState } \tab [numeric] series of snow pack thermal state [°C] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Gratio } \tab [numeric] series of Gratio [0-1] \cr
\emph{$CemaNeigeLayers[[iLayer]]$PotMelt } \tab [numeric] series of potential snow melt [mm/h] \cr
diff --git a/man/RunModel_CemaNeigeGR5J.Rd b/man/RunModel_CemaNeigeGR5J.Rd
index 6f5b2bf4b333f3fb07ddc56a0f5406a19e5d4fbe..3572376e5c04311ea5c70027e5888ef927105910 100644
--- a/man/RunModel_CemaNeigeGR5J.Rd
+++ b/man/RunModel_CemaNeigeGR5J.Rd
@@ -41,29 +41,29 @@ CemaNeige X4 \tab (optional) percentage (between 0 and 1) of annual snowfall def
\value{
[list] list containing the function outputs organised as follows:
\tabular{ll}{
- \emph{$DatesR } \tab [POSIXlt] series of dates \cr
- \emph{$PotEvap } \tab [numeric] series of input potential evapotranspiration [mm/d] \cr
- \emph{$Precip } \tab [numeric] series of input total precipitation [mm/d] \cr
- \emph{$Prod } \tab [numeric] series of production store level [mm] \cr
- \emph{$Pn } \tab [numeric] series of net rainfall [mm/d] \cr
- \emph{$Ps } \tab [numeric] series of the part of Pn filling the production store [mm/d] \cr
- \emph{$AE } \tab [numeric] series of actual evapotranspiration [mm/d] \cr
- \emph{$Perc } \tab [numeric] series of percolation (PERC) [mm/d] \cr
- \emph{$PR } \tab [numeric] series of PR=Pn-Ps+Perc [mm/d] \cr
- \emph{$Q9 } \tab [numeric] series of UH1 outflow (Q9) [mm/d] \cr
- \emph{$Q1 } \tab [numeric] series of UH2 outflow (Q1) [mm/d] \cr
- \emph{$Rout } \tab [numeric] series of routing store level [mm] \cr
- \emph{$Exch } \tab [numeric] series of potential semi-exchange between catchments [mm/d] \cr
- \emph{$AExch1 } \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d] \cr
- \emph{$AExch2 } \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d] \cr
- \emph{$AExch } \tab [numeric] series of actual exchange between catchments (1+2) [mm/d] \cr
- \emph{$QR } \tab [numeric] series of routing store outflow (QR) [mm/d] \cr
- \emph{$QD } \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d] \cr
- \emph{$Qsim } \tab [numeric] series of simulated discharge [mm/d] \cr
- \emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
+ \emph{$DatesR } \tab [POSIXlt] series of dates \cr
+ \emph{$PotEvap } \tab [numeric] series of input potential evapotranspiration [mm/d] \cr
+ \emph{$Precip } \tab [numeric] series of input total precipitation [mm/d] \cr
+ \emph{$Prod } \tab [numeric] series of production store level [mm] \cr
+ \emph{$Pn } \tab [numeric] series of net rainfall [mm/d] \cr
+ \emph{$Ps } \tab [numeric] series of the part of Pn filling the production store [mm/d] \cr
+ \emph{$AE } \tab [numeric] series of actual evapotranspiration [mm/d] \cr
+ \emph{$Perc } \tab [numeric] series of percolation (PERC) [mm/d] \cr
+ \emph{$PR } \tab [numeric] series of PR=Pn-Ps+Perc [mm/d] \cr
+ \emph{$Q9 } \tab [numeric] series of UH1 outflow (Q9) [mm/d] \cr
+ \emph{$Q1 } \tab [numeric] series of UH2 outflow (Q1) [mm/d] \cr
+ \emph{$Rout } \tab [numeric] series of routing store level [mm] \cr
+ \emph{$Exch } \tab [numeric] series of potential semi-exchange between catchments [mm/d] \cr
+ \emph{$AExch1 } \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d] \cr
+ \emph{$AExch2 } \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d] \cr
+ \emph{$AExch } \tab [numeric] series of actual exchange between catchments (1+2) [mm/d] \cr
+ \emph{$QR } \tab [numeric] series of routing store outflow (QR) [mm/d] \cr
+ \emph{$QD } \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d] \cr
+ \emph{$Qsim } \tab [numeric] series of simulated discharge [mm/d] \cr
+ \emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
\emph{$CemaNeigeLayers[[iLayer]]$Pliq } \tab [numeric] series of liquid precip. [mm/d] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Psol } \tab [numeric] series of solid precip. [mm/d] \cr
- \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack [mm] \cr
+ \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack (snow water equivalent) [mm] \cr
\emph{$CemaNeigeLayers[[iLayer]]$ThermalState } \tab [numeric] series of snow pack thermal state [°C] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Gratio } \tab [numeric] series of Gratio [0-1] \cr
\emph{$CemaNeigeLayers[[iLayer]]$PotMelt } \tab [numeric] series of potential snow melt [mm/d] \cr
diff --git a/man/RunModel_CemaNeigeGR6J.Rd b/man/RunModel_CemaNeigeGR6J.Rd
index f1d936c98e6d06fd5ea8ba9c1cba3484b4564df9..99b37e9ded7f78bbd65dd275ad13e23d147d2a2d 100644
--- a/man/RunModel_CemaNeigeGR6J.Rd
+++ b/man/RunModel_CemaNeigeGR6J.Rd
@@ -30,7 +30,7 @@ GR6J X2 \tab intercatchment exchange coefficient [mm/d]
GR6J X3 \tab routing store capacity [mm] \cr
GR6J X4 \tab unit hydrograph time constant [d] \cr
GR6J X5 \tab intercatchment exchange threshold [-] \cr
-GR6J X6 \tab coefficient for emptying exponential store [mm] \cr
+GR6J X6 \tab exponential store depletion coefficient [mm] \cr
CemaNeige X1 \tab weighting coefficient for snow pack thermal state [-] \cr
CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d] \cr
CemaNeige X3 \tab (optional) accumulation threshold [mm] (needed if \code{IsHyst = TRUE}) \cr
@@ -66,7 +66,7 @@ CemaNeige X4 \tab (optional) percentage (between 0 and 1) of annual snowfall def
\emph{$CemaNeigeLayers} \tab [list] list of CemaNeige outputs (1 list per layer) \cr
\emph{$CemaNeigeLayers[[iLayer]]$Pliq } \tab [numeric] series of liquid precip. [mm/d] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Psol } \tab [numeric] series of solid precip. [mm/d] \cr
- \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack [mm] \cr
+ \emph{$CemaNeigeLayers[[iLayer]]$SnowPack } \tab [numeric] series of snow pack (snow water equivalent) [mm] \cr
\emph{$CemaNeigeLayers[[iLayer]]$ThermalState } \tab [numeric] series of snow pack thermal state [°C] \cr
\emph{$CemaNeigeLayers[[iLayer]]$Gratio } \tab [numeric] series of Gratio [0-1] \cr
\emph{$CemaNeigeLayers[[iLayer]]$PotMelt } \tab [numeric] series of potential snow melt [mm/d] \cr
diff --git a/man/RunModel_GR6J.Rd b/man/RunModel_GR6J.Rd
index e9e7b57e6a3914fde4f864fffe4604c65f107cd5..4b7174e5cf841fdd54d93d0879842582d1d60c59 100644
--- a/man/RunModel_GR6J.Rd
+++ b/man/RunModel_GR6J.Rd
@@ -25,12 +25,12 @@ RunModel_GR6J(InputsModel, RunOptions, Param)
\item{Param}{[numeric] vector of 6 parameters
\tabular{ll}{
- GR6J X1 \tab production store capacity [mm] \cr
- GR6J X2 \tab intercatchment exchange coefficient [mm/d] \cr
- GR6J X3 \tab routing store capacity [mm] \cr
- GR6J X4 \tab unit hydrograph time constant [d] \cr
- GR6J X5 \tab intercatchment exchange threshold [-] \cr
- GR6J X6 \tab coefficient for emptying exponential store [mm] \cr
+ GR6J X1 \tab production store capacity [mm] \cr
+ GR6J X2 \tab intercatchment exchange coefficient [mm/d] \cr
+ GR6J X3 \tab routing store capacity [mm] \cr
+ GR6J X4 \tab unit hydrograph time constant [d] \cr
+ GR6J X5 \tab intercatchment exchange threshold [-] \cr
+ GR6J X6 \tab exponential store depletion coefficient [mm] \cr
}}
}
diff --git a/tests/testthat/test-CreateRunOptions.R b/tests/testthat/test-CreateRunOptions.R
new file mode 100644
index 0000000000000000000000000000000000000000..91723330598bf17fe8d5dca46a00074d62d92168
--- /dev/null
+++ b/tests/testthat/test-CreateRunOptions.R
@@ -0,0 +1,31 @@
+context("CreateRunOptions")
+
+test_that("Warm start of GR4J should give same result as warmed model", {
+ data(L0123001)
+ InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P, PotEvap = BasinObs$E)
+ Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+ 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_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions, Param = Param)
+ # 1990
+ RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel, IndPeriod_Run = Ind_Run1))
+ OutputsModel1 <- RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions1, Param = Param)
+ # Warm start 1991
+ RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel, IndPeriod_Run = Ind_Run2,
+ IndPeriod_WarmUp = 0L,
+ IniStates = OutputsModel1$StateEnd)
+ OutputsModel2 <- RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions2, Param = Param)
+ # Compare 1991 Qsim from warm started and from 1990-1991
+ expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730])
+})
diff --git a/tests/testthat/test-CreateiniStates.R b/tests/testthat/test-CreateiniStates.R
new file mode 100644
index 0000000000000000000000000000000000000000..58dbcf8a7ee9fd26967338907048774f82bd0ee5
--- /dev/null
+++ b/tests/testthat/test-CreateiniStates.R
@@ -0,0 +1,102 @@
+context("CreateIniStates on SD model")
+
+data(L0123001)
+
+test_that("Error: SD argument provided on non-SD 'InputsModel'", {
+ InputsModel <-
+ CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E
+ )
+ expect_error(
+ IniStates <-
+ CreateIniStates(
+ FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ ProdStore = 0,
+ RoutStore = 0,
+ ExpStore = NULL,
+ UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+ UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+ SD = list(rep(0, 10))
+ ),
+ regexp = "'SD' argument provided and"
+ )
+})
+
+BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
+
+# Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
+Qupstream <-
+ floor((sin((
+ seq_along(BasinObs$Qmm) / 365 * 2 * 3.14
+ )) + 1) * mean(BasinObs$Qmm, na.rm = TRUE))
+
+InputsModel <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(Qupstream, ncol = 1),
+ LengthHydro = 1000,
+ BasinAreas = BasinAreas
+)
+
+test_that("Error: Non-list 'SD' argument", {
+ expect_error(
+ IniStates <-
+ CreateIniStates(
+ FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ ProdStore = 0,
+ RoutStore = 0,
+ ExpStore = NULL,
+ UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+ UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+ SD = rep(0, 10)
+ ),
+ regexp = "'SD' argument must be a list"
+ )
+})
+
+test_that("Error: Non-numeric items in 'SD' list argument", {
+ lapply(list(list(list(rep(0, 10))), list(toto = NULL)),
+ function(x) {
+ expect_error(
+ IniStates <-
+ CreateIniStates(
+ FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ ProdStore = 0,
+ RoutStore = 0,
+ ExpStore = NULL,
+ UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+ UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+ SD = x
+ ),
+ regexp = "Each item of 'SD' list argument must be numeric"
+ )
+ })
+})
+
+test_that("Error: Number of items not equal to number of upstream connections", {
+ lapply(list(list(), list(rep(0, 10), rep(0, 10))),
+ function(x) {
+ expect_error(
+ IniStates <-
+ CreateIniStates(
+ FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ ProdStore = 0,
+ RoutStore = 0,
+ ExpStore = NULL,
+ UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+ UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+ SD = x
+ ),
+ regexp = "list argument must be the same as the number of upstream"
+ )
+ })
+})
diff --git a/tests/testthat/test-RunModel_LAG.R b/tests/testthat/test-RunModel_LAG.R
index 880f33c311e904065db8dac0d56a665a4f6ca683..00d0104afc576bf336f5730af9d11df00ceeff85 100644
--- a/tests/testthat/test-RunModel_LAG.R
+++ b/tests/testthat/test-RunModel_LAG.R
@@ -50,13 +50,13 @@ 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(
- RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+ RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
regexp = "'InputsModel' should contain an 'OutputsModel' key"
)
})
@@ -134,7 +134,6 @@ test_that("1 input with lag of 0.5 time step delay out gives an output delayed o
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,
@@ -168,13 +167,60 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
InputsModel$BasinAreas <- c(NA, BasinAreas[2L])
# Convert upstream flow to m3/day
InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinAreas[1L] * 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
Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
-
+
expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
+
+# *** IniStates tests ***
+IM <- InputsModel
+IM$BasinAreas <- rep(BasinInfo$BasinArea, 3)
+IM$Qupstream <- cbind(IM$Qupstream, IM$Qupstream)
+IM$LengthHydro <- c(1000, 1500)
+
+PSDini <- ParamSD
+PSDini[1] <- PSDini[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
+RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = IM, IndPeriod_Run = Ind_Run1))
+OutputsModel1 <- RunModel(InputsModel = IM,
+ RunOptions = RunOptions1, Param = PSDini, FUN_MOD = RunModel_GR4J)
+# 1990-1991
+RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = IM, IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
+OutputsModel <- RunModel(InputsModel = IM,
+ RunOptions = RunOptions, Param = PSDini, FUN_MOD = RunModel_GR4J)
+
+test_that("Warm start should give same result as warmed model", {
+ # Warm start 1991
+ RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = IM, IndPeriod_Run = Ind_Run2,
+ IndPeriod_WarmUp = 0L,
+ IniStates = OutputsModel1$StateEnd)
+ OutputsModel2 <- RunModel(InputsModel = IM,
+ RunOptions = RunOptions2, Param = PSDini, 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])
+})
+
+test_that("Error on Wrong length of iniState$SD", {
+ OutputsModel1$StateEnd$SD[[1]] <- c(1,1)
+ RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ 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)
+ )
+})