diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 599ec847b10e5b32ad7701752284926fa025ca7d..73d31f89ad7f752944b114e78ad61f76eab5b84a 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -9,6 +9,7 @@ default:
- echo "setwd(\"$(pwd)\")" > .Rprofile
- PATH=~/R/sources/R-${R_VERSION}/bin:$PATH
- rename "s/${R_VERSION}.airGR/airGR/" *.tar.gz
+ - R -e 'chooseCRANmirror(graphics = FALSE, ind = 1); pkg <- "caRamel"; pkgInst <- installed.packages()[, "Package"]; pkgMiss <- setdiff(pkg, pkgInst); if (length(pkgMiss) > 0) install.packages(pkgMiss)'
.update_packages:
stage: update_packages
diff --git a/DESCRIPTION b/DESCRIPTION
index 547f14b5c8914a2a5c42d5fb87d43e1b71c6ab52..202669ea1e813c56b219b1c3ed1e839553ce43fe 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.6.9.36
+Version: 1.6.10.0
Date: 2021-01-28
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
@@ -28,7 +28,8 @@ Imports:
utils
Suggests:
knitr, rmarkdown,
- coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, imputeTS, Rmalschains,
+ caRamel, coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, imputeTS, Rmalschains,
+ ggplot2, GGally,
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
diff --git a/NEWS.md b/NEWS.md
index 13d0bb6ee6644fa4c44acb8efd26b448790015b0..978dccf3659b6f9a6239165fb9047fa935bb9285 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,11 @@
-### 1.6.9.36 Release Notes (2021-01-28)
+### 1.6.10.0 Release Notes (2021-01-28)
+
+#### New features
+
+- Added a section 'param_optim' vignette to explain how to manage with multiobjective optimization using the 'CaRamel' package. ([#61](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/61))
#### Major user-visible changes
diff --git a/inst/vignettesData/vignetteParamOptimCaramel.rda b/inst/vignettesData/vignetteParamOptimCaramel.rda
new file mode 100644
index 0000000000000000000000000000000000000000..7d2e9fcc41220cab22fa5b1e57e7c666ef74803e
Binary files /dev/null and b/inst/vignettesData/vignetteParamOptimCaramel.rda differ
diff --git a/tests/testthat/test-vignettes.R b/tests/testthat/test-vignettes.R
index 497cb7d8d04c53ef56cebacaa6aefbe456e8f76d..88de5d118cbdd54a94c6bd298065df54db0a664a 100644
--- a/tests/testthat/test-vignettes.R
+++ b/tests/testthat/test-vignettes.R
@@ -11,11 +11,14 @@ test_that("V02.1_param_optim works", {
skip_on_cran()
rm(list = ls())
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+ load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
rda_resGLOB <- resGLOB
rda_resPORT <- resPORT
+ rda_optMO <- optMO
expect_true(RunVignetteChunks("V02.1_param_optim"))
expect_equal(summary(resGLOB), summary(rda_resGLOB), tolerance = 1e-7)
expect_equal(resGLOB[, -1], rda_resGLOB[, -1], tolerance = 1e-2) # High tolerance due to randomisation in optimisations
+ expect_equal(summary(optMO$parameters), summary(rda_optMO$parameters), tolerance = 1e-7)
})
test_that("V02.2_param_mcmc works", {
diff --git a/vignettes/V02.1_param_optim.Rmd b/vignettes/V02.1_param_optim.Rmd
index 0063d9151e5c4c90de00a1d61aaa0374ebac4e31..2dd121e92d89e6838e3c7db3af73322e74080400 100644
--- a/vignettes/V02.1_param_optim.Rmd
+++ b/vignettes/V02.1_param_optim.Rmd
@@ -1,6 +1,6 @@
---
title: "Plugging in new calibration algorithms in airGR"
-author: "François Bourgin"
+author: "François Bourgin, Guillaume Thirel"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
@@ -15,9 +15,13 @@ library(airGR)
library(DEoptim)
library(hydroPSO) # Needs R version >= 3.6 or latticeExtra <= 0.6-28 on R 3.5
library(Rmalschains)
+library(caRamel)
+library(ggplot2)
+library(GGally)
# source("airGR.R")
set.seed(321)
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
```
@@ -158,7 +162,7 @@ As it can be seen in the table below, the four additional optimization strategie
```{r, warning=FALSE, echo=FALSE, eval=FALSE}
resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
round(rbind(
- OutputsCalib$ParamFinalR ,
+ OutputsCalib$ParamFinalR,
airGR::TransfoParam_GR4J(ParamIn = optPORT$par , Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optPSO$par) , Direction = "TR"),
@@ -172,4 +176,100 @@ resGLOB
<!-- This is an expected result because the response surface for quadratic performance criteria of the **GR4J** model is generally sufficiently well defined in the transformed parameter space to allow using a local optimization strategy instead of a more time consuming global optimization strategy. -->
+# Multiobjective optimization
+
+Multiobjective optimization is used to explore possible trade-offs between different performances criteria.
+Here we use the following R implementation of an efficient strategy:
+* [caRamel: Automatic Calibration by Evolutionary Multi Objective Algorithm](https://cran.r-project.org/package=caRamel)
+
+Motivated by using the rainfall-runoff model for low flow simulation, we explore the trade-offs between the KGE values obtained without any data transformation and with the inverse transformation.
+
+First, the OptimGR4J function previously used is modified to return two values.
+
+```{r, warning=FALSE, results='hide', eval=FALSE}
+InputsCrit_inv <- InputsCrit
+InputsCrit_inv$transfo <- "inv"
+
+MOptimGR4J <- function(i) {
+ if (algo == "caRamel") {
+ ParamOptim <- x[i, ]
+ }
+ ## Transformation of the parameter set to real space
+ RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
+ Direction = "TR")
+ ## Simulation given a parameter set
+ OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Param = RawParamOptim)
+ ## Computation of the value of the performance criteria
+ OutputsCrit1 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit,
+ OutputsModel = OutputsModel,
+ verbose = FALSE)
+ ## Computation of the value of the performance criteria
+ OutputsCrit2 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit_inv,
+ OutputsModel = OutputsModel,
+ verbose = FALSE)
+ return(c(OutputsCrit1$CritValue, OutputsCrit2$CritValue))
+}
+```
+
+## caRamel
+caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm.
+
+```{r, warning=FALSE, results='hide', eval=FALSE}
+algo <- "caRamel"
+optMO <- caRamel::caRamel(nobj = 2,
+ nvar = 4,
+ minmax = rep(TRUE, 2),
+ bounds = matrix(c(lowerGR4J, upperGR4J), ncol = 2),
+ func = MOptimGR4J,
+ popsize = 100,
+ archsize = 100,
+ maxrun = 15000,
+ prec = rep(1.e-3, 2),
+ carallel = FALSE,
+ graph = FALSE)
+```
+
+The algorithm returns parameter sets that describe the pareto front, illustrating the trade-off between overall good performance and good performance for low flow.
+
+```{r, fig.width=6, fig.height=6, warning=FALSE}
+ggplot() +
+ geom_point(aes(optMO$objectives[, 1], optMO$objectives[, 2])) +
+ coord_equal(xlim = c(0.4, 0.9), ylim = c(0.4, 0.9)) +
+ xlab("KGE") + ylab("KGE [1/Q]") +
+ theme_bw()
+```
+
+The pameter sets can be viewed in the parameter space, illustrating different populations.
+
+```{r fig.height=6, fig.width=6, message=FALSE, warning=FALSE}
+param_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
+ airGR::TransfoParam_GR4J(x, Direction = "TR")
+ })
+GGally::ggpairs(data.frame(t(param_optMO)), diag = NULL) + theme_bw()
+```
+
+```{r fig.height=6, fig.width=12, message=FALSE, warning=FALSE}
+RunOptions$Outputs_Sim <- "Qsim"
+run_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
+ airGR::RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Param = x)
+ }$Qsim)
+run_optMO <- data.frame(run_optMO)
+
+ggplot() +
+ geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
+ y = run_optMO$X1)) +
+ geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
+ y = run_optMO$X54),
+ colour = "darkred") +
+ scale_x_datetime(limits = c(as.POSIXct("1998-01-01"), NA)) +
+ ylab("Discharge [mm/d]") + xlab("Date") +
+ scale_y_sqrt() +
+ theme_bw()
+```
+
+
diff --git a/vignettes/V02.2_param_mcmc.Rmd b/vignettes/V02.2_param_mcmc.Rmd
index 02e9c5520c37833183840c03c4cb226d062b8888..11bdad00f73c7b6f617794b1031929b2aecb6a24 100644
--- a/vignettes/V02.2_param_mcmc.Rmd
+++ b/vignettes/V02.2_param_mcmc.Rmd
@@ -1,6 +1,5 @@
---
title: "Parameter estimation within a Bayesian MCMC framework"
-author: "François Bourgin"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}