Commit e31bbbfc authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Merge branch '101-fix-vignette-plugging-in-new-calibration-algorithms-in-airgr' into 'dev'

Resolve "Fix vignette "Plugging in new calibration algorithms in airGR""

Closes #101

See merge request !33
parents a9edee4b fc918ee3
Pipeline #21707 passed with stages
in 8 minutes and 37 seconds
...@@ -45,6 +45,13 @@ Please note that the calibration period is defined in the `CreateRunOptions()` f ...@@ -45,6 +45,13 @@ Please note that the calibration period is defined in the `CreateRunOptions()` f
example("Calibration_Michel") example("Calibration_Michel")
``` ```
In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
```{r, results='hide', eval=FALSE}
RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
Outputs_Sim = "Qsim")
```
Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used. Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used.
...@@ -58,7 +65,6 @@ Here we choose to minimize the root mean square error. ...@@ -58,7 +65,6 @@ Here we choose to minimize the root mean square error.
The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model. The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model.
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
OptimGR4J <- function(ParamOptim) { OptimGR4J <- function(ParamOptim) {
## Transformation of the parameter set to real space ## Transformation of the parameter set to real space
...@@ -78,6 +84,7 @@ OptimGR4J <- function(ParamOptim) { ...@@ -78,6 +84,7 @@ OptimGR4J <- function(ParamOptim) {
In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space: In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space:
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
lowerGR4J <- rep(-9.99, times = 4) lowerGR4J <- rep(-9.99, times = 4)
upperGR4J <- rep(+9.99, times = 4) upperGR4J <- rep(+9.99, times = 4)
...@@ -86,6 +93,7 @@ upperGR4J <- rep(+9.99, times = 4) ...@@ -86,6 +93,7 @@ upperGR4J <- rep(+9.99, times = 4)
# Local optimization # Local optimization
We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space: We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space:
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
startGR4J <- c(4.1, 3.9, -0.9, -8.7) startGR4J <- c(4.1, 3.9, -0.9, -8.7)
optPORT <- stats::nlminb(start = startGR4J, optPORT <- stats::nlminb(start = startGR4J,
...@@ -93,13 +101,17 @@ optPORT <- stats::nlminb(start = startGR4J, ...@@ -93,13 +101,17 @@ optPORT <- stats::nlminb(start = startGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
control = list(trace = 1)) control = list(trace = 1))
``` ```
The RMSE value reaches a local minimum value after 35 iterations. The RMSE value reaches a local minimum value after 35 iterations.
We can also try a multi-start approach to test the consistency of the local optimization. We can also try a multi-start approach to test the consistency of the local optimization.
Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function). Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function).
For each starting point, a local optimization is performed. For each starting point, a local optimization is performed.
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
startGR4J <- expand.grid(data.frame(CalibOptions$StartParamDistrib)) StarttGR4JDistrib <- TransfoParam_GR4J(ParamIn = CalibOptions$StartParamDistrib,
Direction = "RT")
startGR4J <- expand.grid(data.frame(StarttGR4JDistrib))
optPORT_ <- function(x) { optPORT_ <- function(x) {
opt <- stats::nlminb(start = x, opt <- stats::nlminb(start = x,
objective = OptimGR4J, objective = OptimGR4J,
...@@ -110,6 +122,7 @@ listOptPORT <- apply(startGR4J, MARGIN = 1, FUN = optPORT_) ...@@ -110,6 +122,7 @@ listOptPORT <- apply(startGR4J, MARGIN = 1, FUN = optPORT_)
``` ```
We can then extract the best parameter sets and the value of the performance criteria: We can then extract the best parameter sets and the value of the performance criteria:
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
parPORT <- t(sapply(listOptPORT, function(x) x$par)) parPORT <- t(sapply(listOptPORT, function(x) x$par))
objPORT <- sapply(listOptPORT, function(x) x$objective) objPORT <- sapply(listOptPORT, function(x) x$objective)
...@@ -117,6 +130,7 @@ resPORT <- data.frame(parPORT, RMSE = objPORT) ...@@ -117,6 +130,7 @@ resPORT <- data.frame(parPORT, RMSE = objPORT)
``` ```
As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets. As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets.
```{r, warning=FALSE} ```{r, warning=FALSE}
summary(resPORT) summary(resPORT)
``` ```
...@@ -134,6 +148,7 @@ Here we use the following R implementation of some popular strategies: ...@@ -134,6 +148,7 @@ Here we use the following R implementation of some popular strategies:
* [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains) * [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains)
## Differential Evolution ## Differential Evolution
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
optDE <- DEoptim::DEoptim(fn = OptimGR4J, optDE <- DEoptim::DEoptim(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
...@@ -142,6 +157,7 @@ optDE <- DEoptim::DEoptim(fn = OptimGR4J, ...@@ -142,6 +157,7 @@ optDE <- DEoptim::DEoptim(fn = OptimGR4J,
## Particle Swarm ## Particle Swarm
```{r, warning=FALSE, results='hide', message=FALSE, eval=FALSE} ```{r, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J, optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
...@@ -149,6 +165,7 @@ optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J, ...@@ -149,6 +165,7 @@ optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
``` ```
## MA-LS-Chains ## MA-LS-Chains
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
optMALS <- Rmalschains::malschains(fn = OptimGR4J, optMALS <- Rmalschains::malschains(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
...@@ -214,6 +231,7 @@ MOptimGR4J <- function(i) { ...@@ -214,6 +231,7 @@ MOptimGR4J <- function(i) {
``` ```
## caRamel ## caRamel
caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm. caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm.
```{r, warning=FALSE, results='hide', eval=FALSE} ```{r, warning=FALSE, results='hide', eval=FALSE}
......
...@@ -38,6 +38,14 @@ We use the **GR4J** model and we assume that the R global environment contains d ...@@ -38,6 +38,14 @@ We use the **GR4J** model and we assume that the R global environment contains d
example("Calibration_Michel") example("Calibration_Michel")
``` ```
In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
```{r, results='hide', eval=FALSE}
RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
Outputs_Sim = "Qsim")
```
Please refer to the [2.1 Plugging in new calibration](V02.1_param_optim.html) vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**. Please refer to the [2.1 Plugging in new calibration](V02.1_param_optim.html) vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**.
...@@ -49,6 +57,7 @@ We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modM ...@@ -49,6 +57,7 @@ We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modM
First, we need to define a function that returns twice the opposite of the log-likelihood for a given parameter set. First, we need to define a function that returns twice the opposite of the log-likelihood for a given parameter set.
Nota: in the `LogLikeGR4J()` function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation of the two following lines. Nota: in the `LogLikeGR4J()` function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation of the two following lines.
```{r, echo=TRUE, eval=FALSE, purl=FALSE} ```{r, echo=TRUE, eval=FALSE, purl=FALSE}
Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2) Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2)
LogLike <- -2 * log(Likelihood) LogLike <- -2 * log(Likelihood)
...@@ -96,8 +105,10 @@ Nota: in this example, there are relatively few iterations (2000), in order to l ...@@ -96,8 +105,10 @@ Nota: in this example, there are relatively few iterations (2000), in order to l
With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied. With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied.
```{r, results='hide', eval=FALSE} ```{r, results='hide', eval=FALSE}
iniParPORT <- data.frame(Chain1 = iniParPORT, Chain2 = iniParPORT, Chain3 = iniParPORT, iniParPORT <- data.frame(Chain1 = iniParPORT,
row.names = paste0("X", 1:4)) Chain2 = iniParPORT,
Chain3 = iniParPORT,
row.names = paste0("X", 1:4))
iniParPORT <- sweep(iniParPORT, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*") iniParPORT <- sweep(iniParPORT, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*")
iniParPORT[iniParPORT < -9.9] <- -9.9 iniParPORT[iniParPORT < -9.9] <- -9.9
iniParPORT[iniParPORT > +9.9] <- +9.9 iniParPORT[iniParPORT > +9.9] <- +9.9
...@@ -111,8 +122,8 @@ mcmcDRAM <- apply(iniParPORT, MARGIN = 2, FUN = function(iIniParPORT) { ...@@ -111,8 +122,8 @@ mcmcDRAM <- apply(iniParPORT, MARGIN = 2, FUN = function(iIniParPORT) {
jump = 0.01, jump = 0.01,
outputlength = 2000, outputlength = 2000,
burninlength = 0, burninlength = 0,
updatecov = 100, ## Adaptative Metropolis updatecov = 100, ## adaptative Metropolis (AM)
ntrydr = 2) ## Delayed Rejection ntrydr = 2) ## delayed rejection (RD)
}) })
``` ```
...@@ -143,7 +154,7 @@ In addition, graphical tools can be used, with for example the [ggmcmc](https:// ...@@ -143,7 +154,7 @@ In addition, graphical tools can be used, with for example the [ggmcmc](https://
First, the evolution of the Markov chains can be seen with a traceplot: First, the evolution of the Markov chains can be seen with a traceplot:
```{r, fig.width=6, fig.height=9, warning=FALSE} ```{r, fig.width=6, fig.height=9, warning=FALSE}
parDRAM <- ggmcmc::ggs(multDRAM) ## to convert objet for using by all ggs_* graphical functions parDRAM <- ggmcmc::ggs(multDRAM) ## to convert object for using by all ggs_* graphical functions
ggmcmc::ggs_traceplot(parDRAM) ggmcmc::ggs_traceplot(parDRAM)
``` ```
......
Markdown is supported
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