diff --git a/tests/testthat/test-RunModel_LLR.R b/tests/testthat/test-RunModel_LLR.R new file mode 100644 index 0000000000000000000000000000000000000000..7ef9a1028077b0550ef4184f31501f1fc2c5f664 --- /dev/null +++ b/tests/testthat/test-RunModel_LLR.R @@ -0,0 +1,350 @@ +context("RunModel_LLR") + +data(L0123001) + +test_that("'BasinAreas' must have one more element than 'LengthHydro'", { + expect_error( + InputsModel <- CreateInputsModel( + FUN_MOD = RunModel_GR4J, + DatesR = BasinObs$DatesR, + Precip = BasinObs$P, + PotEvap = BasinObs$E, + Qupstream = matrix(BasinObs$Qmm, ncol = 1), + LengthHydro = 1, + BasinAreas = 1 + ), + regexp = "'BasinAreas' must have one more element than 'LengthHydro'" + ) +}) + +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 = 1, + BasinAreas = BasinAreas +) + +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 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J, + InputsModel = InputsModel, + IndPeriod_Run = Ind_Run, FUN_SD = RunModel_LLR)) + +test_that("QcontribDown parameter should be a numeric vector or an OutputModel object", { + regexp = "'QcontribDown' must be a numeric vector or a 'OutputsModel' object" + expect_error( + RunModel_LLR(InputsModel = InputsModel, RunOptions = RunOptions, Param = c(1,0), QcontribDown = "A"), + regexp = regexp + ) + expect_error( + RunModel_LLR(InputsModel = InputsModel, RunOptions = RunOptions, Param = c(1,0), QcontribDown = NULL), + regexp = regexp + ) + expect_error( + RunModel_LLR(InputsModel = InputsModel, RunOptions = RunOptions, Param = c(1,0), QcontribDown = matrix(1, ncol = 1)), + regexp = regexp + ) +}) + +ParamGR4J <- c(257.237556, 1.012237, 88.234673, 2.207958) # From vignettes/V01_get_started + +RunOptionsGR4J <- RunOptions +RunOptionsGR4J$FeatFUN_MOD$NbParam <- 4 +OutputsGR4JOnly <- RunModel_GR4J(InputsModel = InputsModel, + RunOptions = RunOptionsGR4J, + Param = ParamGR4J) + +test_that("QcontribDown should contain a Qsim key", { + QcontribDown <- OutputsGR4JOnly + QcontribDown$Qsim <- NULL + expect_error( + RunModel_LLR(InputsModel = InputsModel, RunOptions = RunOptions, Param = c(1,0), QcontribDown = QcontribDown), + regexp = "should contain a key 'Qsim'" + ) +}) + +test_that("'QcontribDown$Qim' should have the same length as 'RunOptions$IndPeriod_Run'", { + QcontribDown <- OutputsGR4JOnly + QcontribDown$Qsim <- c(QcontribDown$Qsim, 0) + expect_error( + RunModel_LLR(InputsModel = InputsModel, RunOptions = RunOptions, Param = c(1,0), QcontribDown = QcontribDown), + regexp = "should have the same length as" + ) +}) +test_that("OutputsModel must have a item 'QsimDown' equal to GR4J Qsim contribution", { + QsimGR4J <- OutputsGR4JOnly$Qsim + QsimDownLLR <- RunModel_LLR(InputsModel = InputsModel, + RunOptions = RunOptions, + Param = c(1,0), + QcontribDown = OutputsGR4JOnly)$QsimDown + expect_equal(QsimGR4J, QsimDownLLR) +}) + +test_that("RunModel(FUN=RunModel_LLR) should give same result as RunModel_LLR", { + QcontribDown <- OutputsGR4JOnly + Output_RunModel_LLR <- RunModel_LLR(InputsModel = InputsModel, + RunOptions = RunOptions, + Param = c(1,0), + QcontribDown = QcontribDown) + Output_RunModel <- RunModel(InputsModel = InputsModel, + RunOptions = RunOptions, + Param = c(1,0), + FUN_MOD = RunModel_LLR, + QcontribDown = QcontribDown, FUN_SD =RunModel_LLR) + expect_equal(Output_RunModel, Output_RunModel_LLR) +}) + +test_that("'Qupstream' contain NA values", { + expect_warning( + InputsModel <- CreateInputsModel( + FUN_MOD = RunModel_GR4J, + DatesR = BasinObs$DatesR, + Precip = BasinObs$P, + PotEvap = BasinObs$E, + Qupstream = matrix(BasinObs$Qmm, ncol = 1), + LengthHydro = 1, + BasinAreas = BasinAreas + ), + regexp = "'Qupstream' contains NA values: model outputs will contain NAs" + ) + + RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J, + InputsModel = InputsModel, + IndPeriod_Run = Ind_Run, FUN_SD = RunModel_LLR)) + QcontribDown <- OutputsGR4JOnly + # Warning with RunModel_LLR + expect_warning( + RunModel_LLR(InputsModel = InputsModel, RunOptions = RunOptions, Param = c(1,0), QcontribDown = QcontribDown), + regexp = "time steps with NA values" + ) + # No warning during calibration + RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal + expect_warning( + RunModel_LLR(InputsModel = InputsModel, RunOptions = RunOptions, Param = c(1,0), QcontribDown = QcontribDown), + regexp = NA + ) +}) + +test_that("Upstream basin with nul area should return same Qdown as GR4J alone", { + InputsModel <- CreateInputsModel( + FUN_MOD = RunModel_GR4J, + DatesR = BasinObs$DatesR, + Precip = BasinObs$P, + PotEvap = BasinObs$E, + Qupstream = matrix(Qupstream, ncol = 1), + LengthHydro = 1, + BasinAreas = c(0,BasinAreas[2]) + ) + OutputsSD <- RunModel(InputsModel, + RunOptions, + Param = c(1, 0, ParamGR4J), + FUN_MOD = RunModel_GR4J, FUN_SD = RunModel_LLR) + expect_equal(OutputsGR4JOnly$Qsim, OutputsSD$Qsim) +}) + +## Ce test la ne marche pas +test_that("Downstream basin with nul 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, + Qupstream = matrix(Qupstream, ncol = 1), + LengthHydro = 0, + BasinAreas = c(BasinInfo$BasinArea, 0) + ) + OutputsSD <- RunModel(InputsModel, + RunOptions, + Param = c(1, 1, ParamGR4J), + FUN_MOD = RunModel_GR4J, + FUN_SD = RunModel_LLR) + expect_equal(OutputsSD$Qsim, Qupstream[Ind_Run]) +}) + +ParamSD <- c(InputsModel$LengthHydro * 1e3 / (24 * 60 * 60), 0, ParamGR4J) # Speed corresponding to one time step delay + +Qm3GR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e3 + +## différence de 0.000289 +test_that('RunModel_Lag is identical to RunModel_LLR with k = 0', { + ParamLLR <-c( 19.990, 0, 43.380, 1.602, 31.817, 2.647) + ParamLAG <- c( 19.990, 43.380, 1.602, 31.817, 2.647) + + RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, IndPeriod_Run = Ind_Run, + FUN_SD = RunModel_Lag)) + OutputLag <- suppressWarnings(RunModel(InputsModel = InputsModel, FUN_MOD = RunModel_GR4J, RunOptions = RunOptions, + Param = ParamLAG, FUN_SD = RunModel_Lag)) + RunOptions <- suppressWarnings((CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, IndPeriod_Run = Ind_Run, + FUN_SD = RunModel_LLR))) + OutputLLR <- RunModel(InputsModel = InputsModel, FUN_MOD = RunModel_GR4J, + RunOptions = RunOptions, Param = ParamLLR, FUN_SD = RunModel_LLR) + + + identical(OutputLLR$Qsim, OutputLag$Qsim) + expect_equal(OutputLLR$Qsim, OutputLag$Qsim, tolerance = 1e-3) +}) + +##relancer InputsModel RunOptions et ParamSD si ça passe pas + +test_that("Qupstream in different units should return the same result", { + OutputsSD_mm <- suppressWarnings(RunModel(InputsModel, RunOptions, + Param = ParamSD, + FUN_MOD = RunModel_GR4J, + FUN_SD = RunModel_LLR)) + InputsModel_m3 <- CreateInputsModel( + FUN_MOD = RunModel_GR4J, + DatesR = BasinObs$DatesR, + Precip = BasinObs$P, + PotEvap = BasinObs$E, + Qupstream = matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1e3, + LengthHydro = 1, + BasinAreas = BasinAreas, + QupstrUnit = "m3" + ) + OutputsSD_m3 <- suppressWarnings(RunModel(InputsModel_m3, RunOptions, + Param = ParamSD, + FUN_MOD = RunModel_GR4J, + FUN_SD = RunModel_LLR)) + expect_equal(OutputsSD_mm$Qsim, OutputsSD_m3$Qsim) + + InputsModel_m3s <- CreateInputsModel( + FUN_MOD = RunModel_GR4J, + DatesR = BasinObs$DatesR, + Precip = BasinObs$P, + PotEvap = BasinObs$E, + Qupstream = matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1e3 / 86400, + LengthHydro = 1, + BasinAreas = BasinAreas, + QupstrUnit = "m3/s" + ) + OutputsSD_m3s <- suppressWarnings(RunModel(InputsModel_m3s, RunOptions, + Param = ParamSD, + FUN_MOD = RunModel_GR4J, + FUN_SD = RunModel_LLR)) + expect_equal(OutputsSD_mm$Qsim, OutputsSD_m3s$Qsim) + + InputsModel_ls <- CreateInputsModel( + FUN_MOD = RunModel_GR4J, + DatesR = BasinObs$DatesR, + Precip = BasinObs$P, + PotEvap = BasinObs$E, + Qupstream = matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1e6 / 86400, + LengthHydro = 1, + BasinAreas = BasinAreas, + QupstrUnit = "L/s" + ) + OutputsSD_ls <- suppressWarnings(RunModel(InputsModel_ls, RunOptions, + Param = ParamSD, + FUN_MOD = RunModel_GR4J, + FUN_SD = RunModel_LLR)) + expect_equal(OutputsSD_mm$Qsim, OutputsSD_ls$Qsim) +}) + +InputsCrit <- CreateInputsCrit( + FUN_CRIT = ErrorCrit_NSE, + InputsModel = InputsModel, + RunOptions = RunOptions, + VarObs = "Q", + Obs = (Qupstream[Ind_Run - 1] * BasinAreas[1L] + + BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas) +) +ParamSD <- c(11.250, 1, 257.238, 1.012, 88.235, 2.208) +## ce test ne marche vraiment pas. Verifier la valeur de ParamSD en faisant retourner un modèle +test_that("Params from calibration with simulated data should be similar to initial params", { + CalibOptions <- CreateCalibOptions( + FUN_MOD = RunModel_GR4J, + FUN_CALIB = Calibration_Michel, + IsSD = TRUE, FUN_SD = RunModel_LLR + ) + OutputsCalib <- Calibration_Michel( + InputsModel = InputsModel, + RunOptions = RunOptions, + InputsCrit = InputsCrit, + CalibOptions = CalibOptions, + FUN_MOD = RunModel_GR4J + ) + expect_equal(OutputsCalib$ParamFinalR[3:6], ParamSD[3:6], tolerance = 1e-2) + expect_equal(OutputsCalib$ParamFinalR[1:2], ParamSD[1:2], tolerance = 2e-3) +}) + +## rajouter le test RunModel_Lag == RunModel_LLR avec k = 0 + +# *** IniStates tests *** +IM <- InputsModel +IM$BasinAreas <- rep(BasinInfo$BasinArea, 3) +IM$Qupstream <- cbind(IM$Qupstream, IM$Qupstream) +IM$LengthHydro <- c(1, 1.5) + +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, FUN_SD = RunModel_LLR)) +OutputsModel1 <- RunModel(InputsModel = IM, + RunOptions = RunOptions1, Param = PSDini, + FUN_MOD = RunModel_GR4J, FUN_SD = RunModel_LLR) +# 1990-1991 +RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J, + InputsModel = IM, + IndPeriod_Run = c(Ind_Run1, Ind_Run2), FUN_SD = RunModel_LLR)) +OutputsModel <- RunModel(InputsModel = IM, + RunOptions = RunOptions, + Param = PSDini, + FUN_MOD = RunModel_GR4J, FUN_SD = RunModel_LLR) + +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, FUN_SD = RunModel_LLR) + OutputsModel2 <- suppressWarnings(RunModel(InputsModel = IM, + RunOptions = RunOptions2, + Param = PSDini, + FUN_MOD = RunModel_GR4J, FUN_SD = RunModel_LLR)) + # Compare 1991 Qsim from warm started and from 1990-1991 + names(OutputsModel2$Qsim) <- NULL + expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730]) +}) + +## cette erreur marche pas non plus +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, FUN_SD = RunModel_LLR) + expect_error( + RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J, + FUN_SD = RunModel_LLR) + ) +}) +## ce test ne marche pas +test_that("First Qupstream time steps must be repeated if warm-up period is too short", { + IM2 <- IM[2558:3557] + IM2$BasinAreas[3] <- 0 + IM2$Qupstream <- matrix(rep(1:1000, 2), ncol = 2) + RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J, + InputsModel = IM2, IndPeriod_Run = seq_len(1000), + IndPeriod_WarmUp = 0L, FUN_SD = RunModel_LLR) + OM2 <- RunModel(InputsModel = IM2, + RunOptions = RunOptions2, + Param = PSDini, + FUN_MOD = RunModel_GR4J, FUN_SD = RunModel_LLR) + expect_equal(OM2$Qsim_m3[1:3], rep(2,3)) +})