---
title: "Getting started with the `calibrar` package"
author: "Ricardo Oliveros-Ramos"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with the `calibrar` package}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: inline
---

```{r setup, include = FALSE}
library(calibrar)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Introduction

This package allows the parameter estimation (i.e. calibration) of complex models, including stochastic ones. It implements generic functions that can be used for fitting any type of models, especially those with non-differentiable objective functions, with the same syntax as `base::optim`. It supports multiple phases estimation (sequential parameter masking), constrained optimization (bounding box restrictions) and automatic parallel computation of numerical gradients. Some common maximum likelihood estimation methods and automated construction of the objective function from simulated model outputs is provided.

## Basic usage

This vignette covers the basic usage of the package, introducing the functions `optim2()`, `optimh()` and `calibrate()`.

### optim2()

As the name sugests, `optim2()` is intended to extend the functionality of `stats::optim()` and it uses the same arguments (with some additions):

```{r, echo=TRUE, eval=FALSE, results='markup'}
optim2(
  par,
  fn,
  gr = NULL,
  ...,
  method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent", "nlm", "nlminb",
    "Rcgmin", "Rvmmin", "hjn", "spg", "LBFGSB3", "AHR-ES"),
  lower = -Inf,
  upper = +Inf,
  active = NULL,
  control = list(),
  hessian = FALSE,
  parallel = FALSE
)
```

The first difference is the possible values for the `method` argument. In addition to the first six methods, also available in `optim()`, `optim2()` gives access to `stats::nlm()` and `stats::nlminb()` but with the same syntax as `optim()` to make them easy to use. In addition, three methods from the `optimr` package (`Rcgmin`, `Rvmmin`, `hjn`), the L-BFGS-B v3 implemented in the `bfgsb3c` package and the `AHR-ES` (Adaptative Hierarchical Recombination Evolutionary Strategy) implemented in this package.

In the next example, we compare the outputs of `optim()` and `optim2()`:

```{r, echo=TRUE, results='markup'}
library(calibrar)
optim(par=rep(1, 5), fn=function(x) sum(x^2))
```

```{r, echo=TRUE, results='markup'}
optim2(par=rep(1, 5), fn=function(x) sum(x^2))
```

The results are identical, as here `optim2()` acts just as a wrapper for `optim()`. Now, we can run the same example with two other methods:

```{r, echo=TRUE, results='markup'}
optim2(par=rep(1, 5), fn=function(x) sum(x^2), method="nlm")
```

```{r, echo=TRUE, results='markup'}
set.seed(880820) # for reproducibility
optim2(par=rep(1, 5), fn=function(x) sum(x^2), method="AHR-ES")
```

The second difference is the new argument `active`, which is a vector indicating if a parameters will be optimized (i.e. active) or fixed to a constant value during the optimization process. In the next example, we will fix the third and fourth parameters to its initial values:

```{r, echo=TRUE, results='markup'}
optim2(par=rep(1, 5), fn=function(x) sum(x^2), 
       active=c(TRUE, TRUE, FALSE, FALSE, TRUE))
```

As we can see, in the final solution, the `par` value keep the values at 1, the initial value provided. All the numerical gradients computed internally have also 'masked' this parameters and the derivatives are not computed for them to speed up computation time.

Finally, the third difference is the new argument `parallel`, that active the parallel computation of the numerical gradient, when `gr` is not supplied:

```{r, echo=TRUE, results='markup'}
optim2(par=rep(1, 5), fn=function(x) sum(x^2), parallel=TRUE)
```

This last option will increase performance when the computation time of `fn` in considerable. We will explain in detail this feature in a following section.

Additionally, the method for the computation of the numerical gradient can be chosen within the `control` list:

```{r, echo=TRUE, results='hide', eval=FALSE}
optim2(par=rep(0.5, 5), fn=function(x) sum(2*x^(3.1*x)), control=list(gr.method="richardson"))
optim2(par=rep(0.5, 5), fn=function(x) sum(2*x^(3.1*x)), control=list(gr.method="central"))
optim2(par=rep(0.5, 5), fn=function(x) sum(2*x^(3.1*x)), control=list(gr.method="forward"))

```

### optimh()

The function `optimh()` has a similar functionality as `optim2()` but acts as a wrapper for several heuristic optimization algorithms implemented in several packages: `dfoptim`, `optimr`, `minqa`, `cmaes`, `genSA`, `DEoptim`, `soma`, `rgenoud` and `psoptim`. The `optimh()` function standardizes the inputs and outputs to those of `optim()`, providing a more convenient user interface. All specific arguments of this methods can be passed to the original function using the `control` argument.

```{r, echo=TRUE, eval=FALSE, results='markup'}
optimh(
  par,
  fn,
  gr = NULL,
  ...,
  method = c("AHR-ES", "Nelder-Mead", "SANN", "hjn", "CMA-ES", "genSA", "DE", "soma",
    "genoud", "PSO", "hybridPSO", "mads", "hjk", "hjkb", "nmk", "nmkb"),
  lower = -Inf,
  upper = +Inf,
  active = NULL,
  control = list(),
  hessian = FALSE,
  parallel = FALSE
)
```

```{r, echo=TRUE, results='markup'}
# Covariance Matrix Adaptation Evolutionary Strategy
set.seed(880820) # for reproducibility
optimh(par=rep(1, 5), fn=function(x) sum(x^2), method="CMA-ES",
       control=list(maxit=200))
```

```{r, echo=TRUE, results='markup'}
# Generalized Simulated Anneling
set.seed(880820) # for reproducibility
optimh(par=rep(1, 5), fn=function(x) sum(x^2), method="genSA", 
       lower=rep(-100, 5), upper=rep(100, 5),
       control=list(maxit=200, temperature=6000))
```

```{r, echo=TRUE, results='markup'}
# Self-Organising Migrating Algorithm
set.seed(880820) # for reproducibility
optimh(par=rep(1, 5), fn=function(x) sum(x^2), method="soma",
       lower=rep(-100, 5), upper=rep(100, 5),
       control=list(maxit=200))
```
The `maxit` control argument has been standardized to work with all methods and to represent the maximum number of iterations of the algorithm. However, the specific number of function evaluations per iteration may vary between methods. You can refer to the help pages from each package to have details about every specific method and its control arguments.

### Running in parallel

Most algorithms implemented in `optim2()` and some in `optimh` can benefit of parallel computation. For the methods that uses the numerical computation of the gradient, this will be calculated in parallel. In order to support any type of parallel implementation, the parallel setup is NOT automatic, and must be done by the user previous to executed the optimization, as described in the following example:

```{r, echo=FALSE, results='hide', eval=TRUE, message=FALSE}
library(parallel)
```

```{r, echo=TRUE, results='markup', eval=FALSE}
library(parallel)
ncores = detectCores() - 1 # number of cores to be used
cl = makeCluster(ncores)
# this is slower than sequential for very fast models (like this one)
optim2(par=rep(0.5, 5), fn=function(x) sum(x^2), 
               control=list(ncores=ncores), parallel=TRUE)
stopCluster(cl) # close the parallel connections
```

## Sequential parameter estimation

### calibrate()

The `calibrate()` function implements an automatic sequential parameter estimation, meaning parameters can be un-masked (set active) progressively during sequential phases of the calibration process.

```{r, echo=TRUE, eval=FALSE, results='markup'}
calibrate(
  par,
  fn,
  gr = NULL,
  ...,
  method = NULL,
  lower = NULL,
  upper = NULL,
  phases = NULL,
  control = list(),
  hessian = FALSE,
  replicates = 1,
  parallel = FALSE
)
```

With the basic syntax, the `calibrate()` function will work similarly to `optim2()` and `optimh()`, performing a simple optimization:

```{r, echo=TRUE, results='markup'}
calibrate(par=c(1,2,3,NA,NA), fn=function(x) sum(x^2))
```

If `upper` and `lower` bounds are provided, the calibrate function can take `NA` as starting values for the optimization, to some or all of the parameters:

```{r, echo=TRUE, results='markup'}
calibrate(par=c(1,2,3,NA,5), fn=function(x) sum(x^2),
          lower=rep(-100, 5), upper=rep(100, 5))
```

### Setting up a parameter estimation with multiple phases.

Multiple `phases` can be set up by selecting a different one for each parameter. The starting value of the optimization for each phase is updated with the best parameters found in the previous phase:

```{r, echo=TRUE, results='markup'}
calibrate(par=c(1,2,3,NA,5), fn=function(x) sum(x^2),
          lower=rep(-100, 5), upper=rep(100, 5),
          phases=c(1,2,3,2,1))
```

When a phase is set to a negative number, the parameter is fixed at its initial value during all the calibration and it is never optimized:

```{r, echo=TRUE, results='markup'}
calibrate(par=c(1,2,3,NA,5), fn=function(x) sum(x^2),
          lower=rep(-100, 5), upper=rep(100, 5),
          phases=c(1,2,-1,2,1))
```

### Dealing with stochastic functions

When dealing with stochastic functions, the argument `replicates` can be helpful, as it allows to evaluate the objective function several times, taking the average value of them as the actual function value (as an approximation of the expected value of the function). When `replicates=1`, the algorithm used by default is "LBFGSB3", but when replicates is greater than 1, "AHR-ES" is used. The next examples use the function `sphereN()`, which computes the Euclidean distance from a point `x` to the origin of coordinates after a random displacement of its position (see `?sphereN` for details). We will set the maximum number of iterations to 1000 to speed up the execution of the vignette.

```{r, echo=TRUE, results='markup'}
calibrate(par=c(1,2,3,NA,5), fn=sphereN,
          lower=rep(-100, 5), upper=rep(100, 5),
          phases=c(1,2,3,2,1), replicates=3, control=list(maxit=1000))
```

The number of replicates can be one single value or a vector with the length equal to the number of phases:

```{r, echo=TRUE, results='markup'}
calibrate(par=c(1,2,3,NA,5), fn=sphereN,
          lower=rep(-100, 5), upper=rep(100, 5),
          phases=c(1,2,3,2,1), replicates=c(1,1,5), control=list(maxit=1000))
```

### Parameters as lists

```{r, echo=TRUE, results='markup'}
calibrate(par=list(par1=c(1,2,3), par2=NA, par3=5), fn=sphereN,
          lower=rep(-100, 5), upper=rep(100, 5),
          phases=c(1,2,-3,2,1), replicates=c(1,5), control=list(maxit=1000))
```

Note that the function `fn` must be able to take a list as parameter set, and the user must ensure this works beforehand.

### Running in parallel

Most algorithms implemented in the `calibrate` function can benefit of parallel computation. For the methods that uses the numerical computation of the gradient, this will be calculated in parallel as in `optim2()`. In order to support any type of parallel implementation, the parallel setup is NOT automatic, and must be done by the user previous to executed the optimization, as described in the following example:

```{r, echo=TRUE, results='hide', eval=FALSE}
library(parallel)
ncores = detectCores() - 1 # number of cores to be used
cl = makeCluster(ncores)
# this is slower than sequential for very fast models (like this one)
calib = calibrate(par=rep(0.5, 5), fn=sphereN,
                  replicates=3, 
                  lower=rep(-5, 5), 
                  upper=rep(+5, 5), 
                  phases=c(1,1,1,2,3), 
                  control=list(parallel=TRUE, ncores=ncores))
stopCluster(cl) # close the parallel connections
```

Please, refer to `vignette(package="calibrar")` for additional vignettes or to the [calibrar website](https://roliveros-ramos.github.io/calibrar/) for more details.
