---
title: "belg: Boltzmann Entropy of a Landscape Gradient"
author: "Jakub Nowosad"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{belg: Boltzmann Entropy of a Landscape Gradient}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(ggplot2)
library(rasterVis)
theme_clean = function (base_size = 12, base_family = ""){
  theme_grey(base_size = base_size, base_family = base_family) %+replace% 
    theme(axis.line = element_blank(), axis.text.x = element_blank(), 
      axis.text.y = element_blank(), axis.ticks = element_blank(), 
      axis.title.x = element_blank(), axis.title.y = element_blank(), 
      axis.ticks.length = grid::unit(0, "lines"), legend.position = "none", 
      panel.background = element_blank(), panel.border = element_blank(), 
      panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
      panel.spacing = grid::unit(0, "lines"), plot.background = element_blank(), 
      strip.background = element_rect(colour = "#ffffff", 
        fill = "#eeeeee"))
}
theme_set(theme_clean())
```

Boltzmann entropy (also called configurational entropy) has been recently adopted to analyze entropy of landscape gradients (Gao et al. (2017), Gao et al. (2018)).
The goal of **belg** is to provide an efficient C++ implementation of this method in R.
It also extend the original idea by allowing calculations on data with missing values.

# Basic example

```{r, message=FALSE}
library(raster)
library(belg)
```

Let's take two small rasters - `complex_land` representing a complex landscape and `simple_land` representing a simple landscape.

```{r, fig.height=4, echo=FALSE}
gplot(complex_land) + 
  geom_tile(aes(fill = value), color = "black", size = 2) +
  geom_text(aes(label = value), size = 4, color = "black") +
  coord_equal() +
  scale_fill_distiller(palette = "RdYlBu") +
  labs(title = "Complex landscape")
gplot(simple_land) + 
  geom_tile(aes(fill = value), color = "black", size = 2) +
  geom_text(aes(label = value), size = 4, color = "black") +
  coord_equal() +
  scale_fill_distiller(palette = "RdYlBu") +
  labs(title = "Simple landscape")
```

The main function in this package, `get_boltzmann()`, calculates the Boltzmann entropy of a landscape gradient:

```{r}
get_boltzmann(complex_land, method = "hierarchy")
get_boltzmann(simple_land, method = "hierarchy")
```

The results, unsurprisingly, showed that the complex landscape has a larger value of the Boltzmann entropy than the simple one.

The `get_boltzmann()` function accepts a `RasterLayer`, `RasterStack`, `RasterBrick`, `matrix`, or `array` object as an input.
As a default, it uses a logarithm of base 10 (`log10`), however `log` and `log2` are also available options for the `base` argument.

```{r}
get_boltzmann(complex_land, method = "hierarchy") # log10
get_boltzmann(complex_land, method = "hierarchy", base = "log")
get_boltzmann(complex_land, method = "hierarchy", base = "log2")
```

It also allows for calculation of the relative (the `relative` argument equal to `TRUE`) and absolute Boltzmann entropy of a landscape gradient.

# Relative Boltzmann entropy of a landscape gradient

The main idea behind the Boltzmann entropy of a landscape gradient is to calculate an entropy in a sliding window of 2 x 2 pixels.
The relative configurational entropy is a sum of entropies for all windows of the original data.

```{r}
get_boltzmann(complex_land, method = "hierarchy", relative = TRUE)
```

# Absolute Boltzmann entropy of a landscape gradient

It is possible to calculate an average value for each sliding window of 2 x 2 pixels and therefore create a resampled version of the original dataset:

```{r, fig.height=4, echo=FALSE}
complex_land_l1 = complex_land
raster_template = raster(ncols = 7, nrows = 5, xmn = 0, xmx = 7, ymn = 0, ymx = 5)
complex_land_l2 = raster(matrix(c(53, 32, 50, 53, 32, 
                                  81, 58, 69, 81, 58,
                                  76, 80, 79, 76, 80,
                                  44, 59, 62, 44, 59, 
                                  41, 47, 56, 41, 47, 
                                  57, 50, 69, 57, 50, 
                                  51, 49, 59, 51, 49), ncol = 7),
                         template = raster_template)
gplot(complex_land_l1) + 
  geom_tile(aes(fill = value), color = "black", size = 2) +
  geom_text(aes(label = value), size = 4, color = "black") +
  coord_equal() +
  scale_fill_distiller(palette = "RdYlBu", limits = c(12, 98)) +
  labs(title = "Original dataset (Level 1)")
gplot(complex_land_l2) + 
  geom_tile(aes(fill = value), color = "black", size = 2) +
  geom_text(aes(label = value), size = 4, color = "black") +
  coord_equal() +
  scale_fill_distiller(palette = "RdYlBu", limits = c(12, 98)) +
  labs(title = "Resampled dataset (Level 2)")
```

The absolute configurational entropy is a sum of relative configurational entropies for all levels, starting from the original data to the resampled dataset with at least two rows or columns.

```{r}
get_boltzmann(complex_land, method = "hierarchy", relative = FALSE)
```

# Calculation of the configurational entropy in a sliding window 

Determining the number of microstates belonging to a defined macrostate in a crucial concept for calculation of the configurational entropy.
We explore this topic using five different cases of 2 x 2 windows:

```{r}
win_1 = raster(matrix(c(1, 3, 3, 4), ncol = 2))
win_2 = raster(matrix(c(1, 3, 3, NA), ncol = 2))
win_3 = raster(matrix(c(1, 3, NA, NA), ncol = 2))
win_4 = raster(matrix(c(1, NA, NA, NA), ncol = 2))
win_5 = raster(matrix(c(NA, NA, NA, NA), ncol = 2))
```

## Data without missing values

The configurational entropy for data without missing values is calculated using the analytical method by Gao et al. (2018).

```{r, echo = FALSE}
val_cols = c("1" = "#fc8d59", "3" = "#ffffbf", "4" = "#91bfdb")
gplot(win_1) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = value), size = 6) +
  scale_fill_manual(values = val_cols) +
  labs(title = "Window 1")
```

Twenty-four different microstate are possible in the above case. 
The common (base 10) logarithm of 24 is equal to `r round(log10(24), 6)`.
We can compare this result to the `get_boltzmann()` output:

```{r}
get_boltzmann(win_1, method = "hierarchy")
```

The generalized (resampled) version of this window has one value, `3`, which is a rounded average of the four original values.

## Data with missing values

The papers of Gao et al. (2017, 2018) only considered data without missing values.
However, the **belg** package provides a modification allowing for calculation also for data with missing values.
Cells with `NA` are not considered when calculating microstates.

```{r, echo = FALSE}
gplot(win_2) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = c(1, 3, 3, "NA")), size = 6) +
  scale_fill_manual(values = val_cols) +
  labs(title = "Window 2")
```

For example, three microstates are possible for the above case:

```{r, echo = FALSE}
win_2_1 = win_2
win_2_2 = raster(matrix(c(3, 1, 3, NA), ncol = 2))
win_2_3 = raster(matrix(c(3, 3, 1, NA), ncol = 2))

gplot(win_2_1) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = c(1, 3, 3, "NA")), size = 6) +
  scale_fill_manual(values = val_cols) +
  labs(title = "Window 2", subtitle = "Microstate I")
gplot(win_2_2) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = c(3, 3, 1, "NA")), size = 6) +
  scale_fill_manual(values = val_cols) +
  labs(title = "Window 2", subtitle = "Microstate II")
gplot(win_2_3) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = c(3, 1, 3, "NA")), size = 6) +
  scale_fill_manual(values = val_cols) +
  labs(title = "Window 2", subtitle = "Microstate III")
```

The common (base 10) logarithm of 3 is equal to `r round(log10(3), 6)`.

```{r}
get_boltzmann(win_2, method = "hierarchy")
```

The generalized (resampled) version of this window is 2.

```{r, echo = FALSE}
win_3_1 = win_3
win_3_2 = raster(matrix(c(3, 1, NA, NA), ncol = 2))
gplot(win_3_1) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = c(1, "NA", 3, "NA")), size = 6) +
  scale_fill_manual(values = c("#fc8d59", "#ffffbf", "#91bfdb")) +
  labs(title = "Window 3", subtitle = "Microstate I")
gplot(win_3_2) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = c(3, "NA", 1, "NA")), size = 6) +
  scale_fill_manual(values = c("#fc8d59", "#ffffbf", "#91bfdb")) +
  labs(title = "Window 3", subtitle = "Microstate II")
```

The third window has two combinations.
The common logarithm of 2 is equal to `r round(log10(2), 6)`.

```{r}
get_boltzmann(win_3, method = "hierarchy")
```

The generalized (resampled) version of this window is also 2.

```{r, echo = FALSE}
gplot(win_4) + 
  geom_tile(aes(fill = as.factor(value)), color = "black", size = 2) +
  geom_text(aes(label = c(1, "NA", "NA", "NA")), size = 6) +
  scale_fill_manual(values = c("#fc8d59", "#ffffbf", "#91bfdb")) +
  labs(title = "Window 4")
```

The fourth window has only one microstate, therefore its common logarithm equals to `r round(log10(1), 6)`.

```{r}
get_boltzmann(win_4, method = "hierarchy")
```

The generalized (resampled) version of this window is the same as only existing value - 1.

```{r, echo = FALSE}
gplot(win_5) + 
  geom_tile(fill = "white", color = "black", size = 2) +
  geom_text(aes(label = c("NA", "NA", "NA", "NA")), size = 6) +
  labs(title = "Window 5")
```

Finally, the last window consists of four missing values.
In these cases, the configurational entropy is zero.

```{r}
get_boltzmann(win_5, method = "hierarchy")
```

Importantly, the generalized version of this window is represented by NA.

## References

- Gao, Peichao, Hong Zhang, and Zhilin Li. "An efficient analytical method for computing the Boltzmann entropy of a landscape gradient." Transactions in GIS (2018).
- Gao, Peichao, Hong Zhang, and Zhilin Li. "A hierarchy-based solution to calculate the configurational entropy of landscape gradients." Landscape Ecology 32(6) (2017): 1133-1146.
