---
title: "3. Calculate the lowest elevation path and barrier height between stable states"
author: "Jingmeng Cui"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{3. Calculate the lowest elevation path and barrier height between stable states}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(simlandr)
```

An important property of the states in a landscape is their (kinetic) stability, characterized by the barrier height between these states and other adjacent states. `simlandr` also provides tools to calculate the barrier heights from landscapes.

You can use the general function `calculate_barrier()` to calculate the barrier for most landscapes. There are also specific `calculate_barrier_*()` functions available. The output of these functions is a `barrier object.

The `barrier` objects contain the potential function and the position of both states and the saddle point. For 3D landscapes, the minimum energy path (MEP) is also provided. `barrier`s can also be calculated for landscapes from multiple simulations. In this case, remember to set `individual_landscape = TRUE` in landscape construction functions.

The local minimums are searched in a square space around a given point. The point with the lowest potential value in the given region is set as the position of the stable state. If all the potential values in the region are equal to `Umax` (which represents ~`Inf`), the barrier calculation functions will expand the searching area automatically. Use `expand = FALSE` to disable this feature.

For landscapes from multiple simulations, the searching regions for their starting and ending points can be different. `simlandr` provides `make_barrier_grid_2d()` and `make_barrier_grid_3d()` functions to help you put these settings into a data frame with the correct format.

The `barrier` objects also provide a `ggplot` geom object that can be added to the landscape plots to show the starting (white), end (white), and saddle (red) points, as well as the MEP (white line, only for 3d landscapes). Use `autolayer(b)` to access those geoms.

Below are examples of different barrier calculations. See the help documents of those functions for further details.

------

Prepare data sets and landscapes (see `vignette("landscape")`)
```{r}
single_test <- sim_fun_test(
  arg1 = list(ele1 = 1),
  arg2 = list(ele2 = 1, ele3 = 0)
)

l_single_2d <- make_2d_static(single_test, x = "out1")
l_single_3d <- make_3d_static(single_test, x = "out1", y = "out2")

batch_test <- new_arg_set()
batch_test <- batch_test %>%
  add_arg_ele("arg2", "ele3", 0.2, 0.5, 0.1)
batch_test_grid <- make_arg_grid(batch_test)
batch_test_result <- batch_simulation(batch_test_grid, sim_fun_test,
  default_list = list(
    arg1 = list(ele1 = 0),
    arg2 = list(ele2 = 0, ele3 = 0)
  ),
  bigmemory = FALSE
)
batch_test2 <- new_arg_set()
batch_test2 <- batch_test2 %>%
  add_arg_ele("arg1", "ele1", 0.2, 0.6, 0.2) %>%
  add_arg_ele("arg2", "ele2", 0.2, 0.6, 0.2)
batch_test_grid2 <- make_arg_grid(batch_test2)
batch_test_result2 <- batch_simulation(batch_test_grid2, sim_fun_test,
  default_list = list(
    arg1 = list(ele1 = 0),
    arg2 = list(ele2 = 0, ele3 = 0)
  ),
  bigmemory = FALSE
)

l_batch_3d_m1 <- make_3d_matrix(batch_test_result, x = "out1", y = "out2", cols = "ele3")
l_batch_2d_m2 <- make_2d_matrix(batch_test_result2, x = "out1", rows = "ele1", cols = "ele2", individual_landscape = TRUE)
l_batch_3d_m2 <- make_3d_matrix(batch_test_result2, x = "out1", y = "out2", rows = "ele1", cols = "ele2", Umax = 10, individual_landscape = TRUE)
```

------

------

Frequently used parameters for the family of barrier functions:

`start_location_value`,`end_location_value`: the initial position (in value) for searching the start/end point;
`start_r`,`end_r`: the searching (L1) radius for searching the start/end point.

------


# Barrier calculation for 2d single landscape
```{r}
b_single_2d <- calculate_barrier(l_single_2d, start_location_value = -2, end_location_value = 2, start_r = 1, end_r = 1)

b_single_2d$local_min_start
b_single_2d$local_min_end
b_single_2d$saddle_point

get_barrier_height(b_single_2d)

plot(l_single_2d) + autolayer(b_single_2d)
```


# Barrier calculation for 3d single landscape
```{r}
b_single_3d <- calculate_barrier(l_single_3d, start_location_value = c(-2.5, -2), end_location_value = c(2.5, 0), start_r = 0.3, end_r = 0.3)
plot(l_single_3d, 2) + autolayer(b_single_3d)
```

# Barrier calculation for 2d batch landscape
```{r}
b_batch_2d_m2 <- calculate_barrier(l_batch_2d_m2, start_location_value = -1, end_location_value = 1, start_r = 0.99, end_r = 0.99)
plot(l_batch_2d_m2) + autolayer(b_batch_2d_m2)
```

# Barrier calculation for 3d batch landscape
```{r}
b_batch_3d_m2 <- calculate_barrier(l_batch_3d_m2, start_location_value = c(-1, -1), end_location_value = c(1, 1), start_r = 0.9, end_r = 0.9)
plot(l_batch_3d_m2) + autolayer(b_batch_3d_m2)
```

# Set the starting and end values for each landscape

```{r}
b_batch_3d_m1 <- calculate_barrier(l_batch_3d_m1, start_location_value = c(0, 0), end_location_value = c(2, 1), start_r = 0.3, end_r = 0.6)
plot(l_batch_3d_m1) + autolayer(b_batch_3d_m1)

## This barrier calculation doesn't find proper local minimums for several landscapes. Specify the searching parameters per landscape manually.
## First, print a template of the data format.

make_barrier_grid_3d(batch_test_grid, start_location_value = c(0, 0), end_location_value = c(2, 1), start_r = 0.3, end_r = 0.6, print_template = TRUE)

## Then, modify the parameters as you want, and send this `barrier_grid` to the barrier calculation function.
b_batch_3d_m1 <- calculate_barrier(
  l_batch_3d_m1,
  make_barrier_grid_3d(batch_test_grid,
    df =
      structure(list(start_location_value = list(
        c(0, 0), c(0, 0), c(0, 0), c(0, 0)
      ), start_r = list(c(0.2, 0.2), c(0.3, 0.3), c(0.3, 0.3), c(0.3, 0.3)), end_location_value = list(c(1, 0.5), c(1.8, 0.8), c(2, 1), c(2, 1)), end_r = c(
        0.6, 0.6, 0.6, 0.6
      )), row.names = c(NA, -4L), class = c(
        "arg_grid",
        "data.frame"
      ))
  )
)
plot(l_batch_3d_m1) + autolayer(b_batch_3d_m1)

## Now it works well.
```