
---
title: "Nelder-Mead Optimization"
---

```{r ch_opt, child="include/chunk_options.rmd", cache=FALSE}
```

```{r, echo = FALSE}
library(nlexperiment)
nl_netlogo_path("c:/Program Files/NetLogo 6.0.1/app") 
```

This example is using NetLogo Flocking model (Wilensky, 1998) to
demonstrate model fitting with Nelder-Mead optimization method.

<aside>
Wilensky, U. (1998). 
[NetLogo Flocking model](http://ccl.northwestern.edu/netlogo/models/Flocking). 
Center for Connected Learning and Computer-Based Modeling, 
Northwestern University, Evanston, IL.
</aside>

There are many R packages for solving optimization problems 
(see [CRAN Task View](https://cran.rstudio.com/web/views/Optimization.html)).
This example uses Nelder-Mead optimization algorithm for derivative-free optimization
implemented in `nmkb` function from __dfoptim__ package.


<aside>
Ravi Varadhan, Johns Hopkins University, Hans W. Borchers and ABB Corporate Research. (2011). dfoptim:
Derivative-free Optimization. R package version 2011.8-1. http://CRAN.R-project.org/package=dfoptim
</aside>


```{r flockNMModel, title = "Define the experiment"}
experiment <- nl_experiment( 
  model_file = "models/Sample Models/Biology/Flocking.nlogo",

  setup_commands = c("setup", "repeat 100 [go]"),
  iterations = 5,

  param_values = list(
    world_size = 50,
    population = 80,
    vision = 6,
    min_separation = seq(from = 0, to = 4, by = 0.5),
    max_align_turn = seq(from = 0, to = 20, by = 2.5)
  ),
  mapping = c(
    min_separation = "minimum-separation",
    max_align_turn = "max-align-turn"),

  step_measures = measures(
    converged = "1 - 
      (standard-deviation [dx] of turtles + 
       standard-deviation [dy] of turtles) / 2",
    mean_crowding = 
      "mean [count flockmates + 1] of turtles"
  ),
  eval_criteria = criteria(
    c_converged = mean(step$converged),
    c_mcrowding = mean(step$mean_crowding)
  ),

  repetitions = 10,                 # repeat simulations 10 times

  eval_aggregate_fun = mean,        # aggregate over repetitions

  eval_mutate = criteria(           # evaluation criterium
    eval_value = 
      sqrt((c_mcrowding - 8)^2 + 400*(c_converged - 1)^2)
  )
)
```
<aside>
Note that best-fit criteria function is defined in the experiment.
$$ \sqrt{(crowding - 8)^2 + 800 \times (converged - 1)^2} $$
</aside>

In this example the `nl_eval_run` function is used instead of `nl_run`. 

While `nl_run` expects predefined parameter sets, `nl_run_eval` accepts a single parameter set and returns a numeric value. That makes it compatible with optimization
functions from different R packages where next parameter combination depends on 
the result of previous simulation.

But note that `nl_eval_run` requires started NetLogo instance. 
User have to take care to 
initialize NetLogo and load the model before optimization begins and 
close NetLogo when it is no longer needed 
(see `nl_eval_init` and `nl_eval_close` in package documentation).

Use `nl_eval_run` _parallel_ option when optimizing stochastic models
with more than a few repetitions needed to evaluate one parameter set. 


```{r flockNMOptimRun, title = "Run optim function", cache=TRUE}
library(dfoptim)

cl <- nl_eval_init(experiment, parallel = TRUE)
trace <- nl_eval_tracer(verbose = FALSE)
param_range <- nl_get_param_range(experiment)   
set.seed(1) 

o_result <- nmkb(
  par = (param_range$upper + param_range$lower)/2,
  fn = nl_eval_run, 
    experiment = experiment, 
    criteria = "eval_value", 
    call_back = trace$add, 
    parallel = TRUE, cluster = cl,
    param_names = names(param_range$lower),
  lower = param_range$lower, 
  upper = param_range$upper, 
  control = list(maxfeval = 200)
)

nl_eval_close(parallel = TRUE, cl)
```
<aside>
Function `nl_eval_init` starts processes with NetLogo instance

Function `nl_get_param_range` returns parameters' lower and upper values
as they are defined in experiment object.
This range will be used as optimization constraints in L-BFGS-B method.
</aside>


```{r flockNMGetData, title = "Optim results"}
o_result
```


The `trace$add` function colected every iteration of `nmkb`. 

```{r pFlockNMPlot, title = "Plot tracing data", dev.args = list(bg = 'transparent')}
tr <- trace$get()
tr$current_best <- cummin(tr$result)
library(ggplot2)
ggplot(tr, aes(x = iter_id,  y = result)) +
  geom_point() +
  geom_step(aes(y = current_best), alpha  =0.3) +
  theme_minimal()
```


```{r pFlockNMPlot2, title = "Plot tracing data", fig.cap="Where did `nmkb` look for best value?", dev.args = list(bg = 'transparent')}
library(ggplot2)
ggplot(tr, aes(x = min_separation,  y = max_align_turn)) +
  geom_path(alpha = 0.1)+
  geom_point(shape = 4, size = 5, alpha = 0.3) +
  scale_x_continuous(limits = c(0, 4)) +
  scale_y_continuous(limits = c(0, 20)) +
  theme_minimal()

```

## Related

[Simulated Annealing](flocking_sa.html) demonstrates optimization
  with simulated annealing.

[Genetic Algorithm](flocking_ga.html) demonstrates optimization
  with genetic algorithm.

[Random Search](flocking_random_search.html) demonstrates searching for optimal parameters
  using random search.
  
<asideclose>
Get [rmarkdown source](`r knitr::current_input()`) of this page.
</asideclose>
