Commit 477a3a2e authored by Dorchies David's avatar Dorchies David
Browse files

test: add benchmarking tests on RunModel

Refs #136
Showing with 207 additions and 123 deletions
+207 -123
......@@ -23,6 +23,10 @@ default:
stage: check
script:
- R -e 'rcmdcheck::rcmdcheck(args = ifelse(as.logical(Sys.getenv("NOT_CRAN")), "", "--as-cran"), error_on = "warning")'
- R -e 'install.packages("airGR", repos = "http://cran.r-project.org")'
- Rscript tests/scheduled_tests/benchmarkRunModel.R
- R CMD INSTALL .
- Rscript tests/scheduled_tests/benchmarkRunModel.R
.check_not_cran:
variables:
......
library(airGR)
source("tests/testthat/helper_scheduled_Calibration.R")
BenchmarkRunModel <- function(model) {
e <- PrepareCalibration(model)
for(n in ls(e, all.names=TRUE)) assign(n, get(n, e))
# RunOptions calibration configuration
RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal
mbm <- microbenchmark::microbenchmark(RunModel =
RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = as.numeric(strsplit(model$ParamFinalR, ";")[[1]]), FUN_MOD = sModel)
)
mbm$time
}
Args <- commandArgs(trailingOnly = TRUE)
sModelNames <- paste0(dfModels$name,
ifelse(as.logical(dfModels$IsHyst), "_Hyst", ""))
dfBM <- as.data.frame(apply(dfModels, 1, BenchmarkRunModel))
colnames(dfBM) <- sModelNames
dfBM <- cbind(version = as.character(packageVersion('airGR')), dfBM)
file <- "tests/tmp/benchmark.txt"
write.table(dfBM, file = file,
row.names = FALSE, col.names = !file.exists(file), quote = FALSE,
sep = "\t", append = TRUE)
df <- read.table(file = file, sep = "\t", header = TRUE)
if(length(unique(df$version)) > 1) {
lV <- lapply(unique(df$version), function(version) {
apply(df[df$version == version, -1] / 1E6, 2, mean)
})
names(lV) <- unique(df$version)
dfMean <- cbind(model = sModelNames, as.data.frame(t(do.call(rbind, lV))))
dfMean$evolution <- (dfMean[,3] - dfMean[,2]) / dfMean[,2]
write.table(dfMean, "tests/tmp/mean_execution_time.txt", row.names = FALSE, quote = F, sep = "\t")
res <- testthat::test_file("tests/testthat/benchmark_RunModel.R")
dRes <- as.data.frame(res)
if (any(dRes[, "failed"] > 0) | any(dRes[, "error"])) {
quit(status = 1)
}
}
......@@ -34,13 +34,15 @@ StoreExampleResults <- function(package, path, run.dontrun = FALSE, run.donttest
# Get the list of documentation pages
rd <- unique(readRDS(system.file("help", "aliases.rds", package = package)))
dir.create(path, showWarnings = FALSE)
unlink(path, recursive = TRUE)
dir.create(path)
lapply(
rd,
StoreTopicResults,
package, path, run.dontrun = run.dontrun, run.donttest = run.donttest
)
}
StoreTopicResults <- function(topic, package, path, run.dontrun = TRUE, run.donttest = TRUE) {
......@@ -54,12 +56,21 @@ StoreTopicResults <- function(topic, package, path, run.dontrun = TRUE, run.dont
varBefore <- c()
varBefore <- ls(envir = globalenv())
start_time = Sys.time()
example(
topic, package = package, character.only = TRUE, echo = FALSE, ask = FALSE, local = FALSE, setRNG = TRUE,
run.dontrun = run.dontrun, run.donttest = run.donttest
)
end_time = Sys.time()
dev.off()
write.table(data.frame(topic = topic, time = end_time - start_time),
file.path(path, "timing.csv"),
row.names = FALSE, col.names = FALSE, quote = FALSE,
sep = "\t", append = TRUE)
varAfter <- ls(envir = globalenv())
varToSave <- setdiff(varAfter, varBefore)
......
df <- read.table("../tmp/mean_execution_time.txt", sep = "\t", header = T)
lapply(df$model, function(model) {
test_that(paste(model, ": RunModel should be as fast as CRAN version"), {
expect_lt(df$evolution[df$model == model], 1.1)
})
})
# Calibrations on v1.6.12
sModelCalibrations <- c(
"name IsHyst data aggreg ParamFinalR",
"GR1A FALSE L0123001 %Y 0.91125",
"GR2M FALSE L0123001 %Y%m 259.8228;0.9975",
"GR4J FALSE L0123001 NA 223.6315877;0.5781516;97.5143942;2.2177177",
"GR5J FALSE L0123001 NA 220.3863609;0.8944531;93.5640705;1.7628720;0.4846427",
"GR6J FALSE L0123001 NA 192.8761657;0.6933087;49.1783293;2.2145422;0.5088240;6.8146261",
"CemaNeigeGR4J FALSE L0123001 NA 2.043839e+02;5.781516e-01;1.025141e+02;2.217718e+00;1.501502e-03;1.432036e+01",
"CemaNeigeGR5J FALSE L0123001 NA 1.983434e+02;8.747758e-01;9.849443e+01;1.768769e+00;4.829830e-01;1.501502e-03;1.432036e+01",
"CemaNeigeGR6J FALSE L0123001 NA 184.9341841;0.5551637;59.7398917;2.2177177;0.4760000;6.0496475;0.0000000;14.4642868",
"CemaNeigeGR4J TRUE L0123001 NA 208.5127103;0.5781516;102.5140641;2.2274775;0.0000000;6.7644613;8.0000000;1.0000000",
"CemaNeigeGR5J TRUE L0123001 NA 202.350228;0.901525;98.494430;1.788288;0.483984;0.000000;7.401500;6.100000;1.000000",
"CemaNeigeGR6J TRUE L0123001 NA 188.67010241;0.56662930;60.34028760;2.22747748;0.47600000;5.98945247;0.03203203;7.93816892;10.80000000;1.00000000",
"GR4H FALSE L0123003 NA 711.676649;-1.158469;150.556095;4.686093",
"GR5H FALSE L0123003 NA 804.0021672;-0.1898488;137.7524699;3.0436628;0.1951163",
"CemaNeigeGR4H FALSE L0123003 NA 1.595685e+03;-8.183484e-01;2.320697e+02;5.000000e-01;5.005005e-04;9.342369e+01",
"CemaNeigeGR5H FALSE L0123003 NA 33.34921883;-4.98925432;332.00673122;1.58534106;0.20792716;0.02214393;4.28498513",
"CemaNeigeGR4H TRUE L0123003 NA 1.766316e+03;-6.920667e-01;2.192034e+02;3.451688e+00;5.005005e-04;4.869585e+01;1.111447e+01;5.064090e-01",
"CemaNeigeGR5H TRUE L0123003 NA 66.6863310;-1.4558128;138.3795123;2.6499450;0.2325000;0.0000000;0.3017014;48.4000000;0.9914915"
)
dfModels <- read.table(text = paste(sModelCalibrations, collapse = "\n"), header = TRUE)
PrepareCalibration <- function(model) {
model <- as.list(model)
sModel <- paste0("RunModel_", model$name)
sIM_FUN_MOD <- sModel
if (model$data == "L0123003") {
# hourly time step database
dates <- c("2004-01-01 00:00", "2004-12-31 23:00", "2005-01-01 00:00", "2008-12-31 23:00")
date_format = "%Y-%m-%d %H:%M"
TempMean <- fakeHourlyTemp()
} else {
# yearly, monthly, daily time step databases
dates <- c("1985-01-01", "1985-12-31", "1986-01-01", "2012-12-31")
date_format <- "%Y-%m-%d"
if (!is.na(model$aggreg)) {
# Aggregation on monthly and yearly databases
sIM_FUN_MOD <- "RunModel_GR4J" # CreateInputsModel with daily data
date_format <- model$aggreg
}
}
## loading catchment data
data(list = model$data)
if (model$data != "L0123003") TempMean <- BasinObs$T
# preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = sIM_FUN_MOD,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E,
TempMean = TempMean,
ZInputs = median(BasinInfo$HypsoData),
HypsoData = BasinInfo$HypsoData,
NLayers = 5)
if (!is.na(model$aggreg)) {
# conversion of InputsModel to target time step
InputsModel <- SeriesAggreg(InputsModel, Format = model$aggreg)
dfQobs <- SeriesAggreg(data.frame(DatesR = BasinObs$DatesR, Qmm = BasinObs$Qmm),
Format = model$aggreg, ConvertFun = "sum")
Obs <- dfQobs$Qmm
} else {
Obs <- BasinObs$Qmm
}
# calibration period selection
dates <- sapply(dates, function(x) format(as.Date(x), format = date_format))
Ind_WarmUp <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[1]),
which(format(InputsModel$DatesR, format = date_format)==dates[2])
)
Ind_Run <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[3]),
which(format(InputsModel$DatesR, format = date_format)==dates[4])
)
# preparation of the RunOptions object
suppressWarnings(
RunOptions <- CreateRunOptions(
FUN_MOD = sModel,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
IndPeriod_WarmUp = Ind_WarmUp,
IsHyst = as.logical(model$IsHyst)
)
)
# calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = Obs[Ind_Run])
# preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(sModel, IsHyst = as.logical(model$IsHyst))
return(environment())
}
ModelCalibration <- function(model) {
e <- PrepareCalibration(model)
for(n in ls(e, all.names=TRUE)) assign(n, get(n, e))
# calibration
suppressWarnings(OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = sModel))
OutputsCalib$ParamFinalR
}
#' Create Fake hourly temperature from daily temperatures in L0123001
#'
#' @param start_date [character] start date in format "%Y-%m-%d"
#' @param end_date [character] end date in format "%Y-%m-%d"
#' @return [numeric] hourly temperature time series between `start_date` and `end_date`
fakeHourlyTemp <- function(start_date = "2004-01-01", end_date = "2008-12-31") {
dates <- as.POSIXct(c(start_date, end_date), tz = "UTC")
data(L0123002)
indJ <- seq.int(which(BasinObs$DatesR == as.POSIXct(dates[1])),
which(BasinObs$DatesR == as.POSIXct(dates[2])))
TJ <- BasinObs$T[indJ]
TH <- approx((seq.int(length(TJ)) - 1) * 24,TJ,
seq.int(length(TJ) * 24 ) - 1,
rule = 2)$y
varT_1J <- -sin(0:23/24 * 2 * pi) # Temp min at 6 and max at 18
varT <- rep(varT_1J, length(TJ))
TH <- TH + varT * 5 # For a mean daily amplitude of 10°
TH
}
context("Calibration")
sModels <- c(
"name IsHyst data aggreg ParamFinalR",
"GR1A FALSE L0123001 %Y 0.91125",
"GR2M FALSE L0123001 %Y%m 259.8228;0.9975",
"GR4J FALSE L0123001 NA 223.6315877;0.5781516;97.5143942;2.2177177",
"GR5J FALSE L0123001 NA 220.3863609;0.8944531;93.5640705;1.7628720;0.4846427",
"GR6J FALSE L0123001 NA 192.8761657;0.6933087;49.1783293;2.2145422;0.5088240;6.8146261",
"CemaNeigeGR4J FALSE L0123001 NA 2.043839e+02;5.781516e-01;1.025141e+02;2.217718e+00;1.501502e-03;1.432036e+01",
"CemaNeigeGR5J FALSE L0123001 NA 1.983434e+02;8.747758e-01;9.849443e+01;1.768769e+00;4.829830e-01;1.501502e-03;1.432036e+01",
"CemaNeigeGR6J FALSE L0123001 NA 184.9341841;0.5551637;59.7398917;2.2177177;0.4760000;6.0496475;0.0000000;14.4642868",
"CemaNeigeGR4J TRUE L0123001 NA 208.5127103;0.5781516;102.5140641;2.2274775;0.0000000;6.7644613;8.0000000;1.0000000",
"CemaNeigeGR5J TRUE L0123001 NA 202.350228;0.901525;98.494430;1.788288;0.483984;0.000000;7.401500;6.100000;1.000000",
"CemaNeigeGR6J TRUE L0123001 NA 188.67010241;0.56662930;60.34028760;2.22747748;0.47600000;5.98945247;0.03203203;7.93816892;10.80000000;1.00000000",
"GR4H FALSE L0123003 NA 711.676649;-1.158469;150.556095;4.686093",
"GR5H FALSE L0123003 NA 804.0021672;-0.1898488;137.7524699;3.0436628;0.1951163",
"CemaNeigeGR4H FALSE L0123003 NA 1.595685e+03;-8.183484e-01;2.320697e+02;5.000000e-01;5.005005e-04;9.342369e+01",
"CemaNeigeGR5H FALSE L0123003 NA 33.34921883;-4.98925432;332.00673122;1.58534106;0.20792716;0.02214393;4.28498513",
"CemaNeigeGR4H TRUE L0123003 NA 1.766316e+03;-6.920667e-01;2.192034e+02;3.451688e+00;5.005005e-04;4.869585e+01;1.111447e+01;5.064090e-01",
"CemaNeigeGR5H TRUE L0123003 NA 66.6863310;-1.4558128;138.3795123;2.6499450;0.2325000;0.0000000;0.3017014;48.4000000;0.9914915"
)
dfModels <- read.table(text = paste(sModels, collapse = "\n"), header = TRUE)
# helpers functions and dfModels variable are in tests/testthat/helper_scheduled_Calibration.R
ModelCalibration <- function(model) {
sModel <- paste0("RunModel_", model$name)
sIM_FUN_MOD <- sModel
if (model$data == "L0123003") {
# hourly time step database
dates <- c("2004-01-01 00:00", "2004-12-31 23:00", "2005-01-01 00:00", "2008-12-31 23:00")
date_format = "%Y-%m-%d %H:%M"
TempMean <- fakeHourlyTemp()
} else {
# yearly, monthly, daily time step databases
dates <- c("1985-01-01", "1985-12-31", "1986-01-01", "2012-12-31")
date_format <- "%Y-%m-%d"
if (!is.na(model$aggreg)) {
# Aggregation on monthly and yearly databases
sIM_FUN_MOD <- "RunModel_GR4J" # CreateInputsModel with daily data
date_format <- model$aggreg
}
}
## loading catchment data
data(list = model$data)
if (model$data != "L0123003") TempMean <- BasinObs$T
# preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = sIM_FUN_MOD,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E,
TempMean = TempMean,
ZInputs = median(BasinInfo$HypsoData),
HypsoData = BasinInfo$HypsoData,
NLayers = 5)
if (!is.na(model$aggreg)) {
# conversion of InputsModel to target time step
InputsModel <- SeriesAggreg(InputsModel, Format = model$aggreg)
dfQobs <- SeriesAggreg(data.frame(DatesR = BasinObs$DatesR, Qmm = BasinObs$Qmm),
Format = model$aggreg, ConvertFun = "sum")
Obs <- dfQobs$Qmm
} else {
Obs <- BasinObs$Qmm
}
# calibration period selection
dates <- sapply(dates, function(x) format(as.Date(x), format = date_format))
Ind_WarmUp <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[1]),
which(format(InputsModel$DatesR, format = date_format)==dates[2])
)
Ind_Run <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[3]),
which(format(InputsModel$DatesR, format = date_format)==dates[4])
)
# preparation of the RunOptions object
suppressWarnings(
RunOptions <- CreateRunOptions(
FUN_MOD = sModel,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
IndPeriod_WarmUp = Ind_WarmUp,
IsHyst = as.logical(model$IsHyst)
)
)
# calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = Obs[Ind_Run])
# preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(sModel, IsHyst = as.logical(model$IsHyst))
# calibration
suppressWarnings(OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = sModel))
OutputsCalib$ParamFinalR
}
TestModelCalibration <- function(model) {
......@@ -116,26 +17,4 @@ TestModelCalibration <- function(model) {
})
}
#' Create Fake hourly temperature from daily temperatures in L0123001
#'
#' @param start_date [character] start date in format "%Y-%m-%d"
#' @param end_date [character] end date in format "%Y-%m-%d"
#' @return [numeric] hourly temperature time series between `start_date` and `end_date`
fakeHourlyTemp <- function(start_date = "2004-01-01", end_date = "2008-12-31") {
dates <- as.POSIXct(c(start_date, end_date), tz = "UTC")
data(L0123002)
indJ <- seq.int(which(BasinObs$DatesR == as.POSIXct(dates[1])),
which(BasinObs$DatesR == as.POSIXct(dates[2])))
TJ <- BasinObs$T[indJ]
TH <- approx((seq.int(length(TJ)) - 1) * 24,TJ,
seq.int(length(TJ) * 24 ) - 1,
rule = 2)$y
varT_1J <- -sin(0:23/24 * 2 * pi) # Temp min at 6 and max at 18
varT <- rep(varT_1J, length(TJ))
TH <- TH + varT * 5 # For a mean daily amplitude of 10°
TH
}
apply(dfModels, 1, TestModelCalibration)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment