More than 1/8th of Earth’s land area is within some form of protected area; however, protected area (PA) networks often do a poor job of representing the biodiversity that they ostensibly protect. In large part, this is because new parks and PAs historically haven’t been placed in the optimal locations for conservation. Often they are placed opportunistically, where the land is cheap and unlikely to experience land use change anyway. Many of the most famous parks in North America (e.g. Yosemite or Banff) were chosen more for scenic and recreation value than conservation potential. Or, more subtly, PAs may not complement each other, with some taxa or habitats over-represented and others underrepresented.

Enter the field of systematic conservation planning, which offers a set of tools and techniques for systematically identifying optimal locations for new PAs in a way that maximizes conservation benefit and minimizes socio-economic cost. These tools have changed the way conservation decisions are made, and how and where conservation resources are allocated. Marxan is the most widely used tool for systematic reserve design. Conservation researchers and decision makers all over the world use Marxan in both the terrestrial and marine realms, for projects large and small.

A detailed understanding of how Marxan works isn’t required to use the software; however, I hate a black box so, as I started using Marxan, I began digging into the documentation and source code to see exactly what it was doing. This led me to emulating Marxan within R as a means of checking that I was understanding everything correctly. Of course, Marxan is highly optimized in C and very feature rich, and this R emulator could never be used in a real planning exercise. In contrast, I’ve made absolutely no effort to optimize my code, but I’m sharing my efforts here anyway in the hopes that it may be instructive.

Required packages

library(plyr)
library(dplyr)
library(tidyr)
library(assertthat)
library(sp)
library(raster)
library(rgeos)
library(rasterVis)
library(viridis)
library(gstat)
library(marxan)
library(stringr)
library(Rcpp)
library(ggplot2)
library(gridExtra)
library(scales)
set.seed(1)

Background

Resources

To write this post, I used the following resources:

Marxan

A Marxan exercise starts by dividing the study region into planning units (typically square or hexagonal cells) and, for each planning unit, assigning values that quantify socio-economic cost and conservation benefit for a set of conservation features. The cost can be the acquisition cost of the land, the costs of management, the opportunity cost of foregone commercial activities (e.g. from logging or oil palm development), or simply the area. The conservation features are typically species (e.g. Clouded Leopard) or habitats (e.g. mangroves or cloud forest). The benefit that each features derives from a planning unit can take a variety of forms, but is typically either occupancy (i.e. presence or absence) or area of occurrence within each planning unit. Finally, for each conservation feature, representation targets must be set, such as 20% of the historical extent of cloud forest or 10,000km2 of Clouded Leopard habitat.

The ultimate goal of a conservation planning exercise is to ensure the long term persistence of the conservation features. When planning for persistence, it’s important to consider the total area of protected habitat and the spatial configuration of that habitat. A reserve network that is highly fragmented, with poor connectivity between patches, will not be able to support species in the long term. To account for this, Marxan can be configured to look for more compact reserve networks by minimizing the total boundary length of the reserve.

Given these inputs, Marxan finds the set of planning units that meets the representation targets while minimizing cost and boundary length. This is a form of the minimum set cover problem and is formulated mathematically, for \( n \) planning units and \( m \) conservation features, as:

where \( x_i \) is a binary decision variable specifying whether a planning unit has been selected (1) or not (0), \( c_i \) is the cost of planning unit \( i \), \( r_{ij} \) is the representation level of feature \( j \) in planning unit \( i \), and \( T_j \) is the target for feature \( j \). \( \nu_{ij} \) is a matrix where off diagonal components are the lengths of the shared boundaries between planning units and diagonal components are external boundaries of planning units (i.e. those that are not shared). Finally, \( b \) is known as the Boundary Length Modifier (BLM) and determines how much emphasis should be placed on producing compact solutions relative to meeting targets and minimizing cost. Note that the second and third terms together are just the total length of the perimeter of the reserve network. In all the Marxan documentation I’ve seen the third term is absent; however, when Marxan is run, this term is included.

Finally, a status code can be assigned to each planning unit to specify whether it is locked in or locked out the of solution. The codes and their meanings are as follows:

Status Meaning
0 Not guaranteed in initial reserve system
1 Included in initial reserve system, but not guaranteed in final
2 Locked in to solution, i.e. guaranteed protected
3 Lockout out of solution, i.e. guaranteed unprotected

Marxan requires 5 input files in order to run:

  1. Planning Units (pu.dat): units to be considered for the reserve system, including cost and status.
  2. Conservation Features (spec.dat): features and corresponding targets.
  3. Planning Units vs. Conservation Features (puvspr.dat): representation level of each feature within each planning unit.
  4. Boundary Length (bound.dat): pairwise shared boundaries between planning units.
  5. Parameters (input.dat): parameters and settings for the Marxan run.

The marxan R package acts as an interface between R and Marxan. It automates the process of generating these input files from R objects and provides an R wrapper for the Marxan executable.

Simulated Annealing

The main method that Marxan uses for solving reserve design problems is simulated annealing, a stochastic heuristic for approximating global optima of complex functions with many local optima. Simulated annealing is a highly flexible optimization method that minimizes a given objective function over a set of decision variables by stochastically exploring the state space of the decision variables. Starting from a randomly generated initial state, simulated annealing proceeds iteratively by:

  1. Randomly picking a neighbouring state from a well-defined set of candidates, and calculating the value of the objective function at this state.
  2. If the candidate state reduces the objective function, accept it. If the candidate state increases the objective function, accept it according to an acceptance probability, which depends on the change in objective function value and a global temperature parameter.
  3. As the heuristic progresses, the temperature parameter gradually decreases according to an annealing schedule. As a result, changes that increase the objective function become less likely to be accepted. Initially permitting “bad” changes ensures that the heuristic is unlikely to get caught in a local minimum.

Simulated annealing does not find the true minimum of the objective function. Rather it finds a near optimal solution and, since it is stochastic, each time simulated annealing is run, a different near optimal solution is returned.

In the context of Marxan, the decision variable is a binary vector indicating which planning units are included in the reserve, and the states are different reserve configurations. Choosing a neighbouring state involves randomly selecting a planning unit and switching its status from selected to not selected or vice versa. In a typical Marxan run, simulated annealing is performed many times to generate a suite of possible reserve networks. Thus, multiple options can be presented to decision makers, and the importance of particular planning units can be measured by the selection frequency (i.e. number of solutions containing that planning unit).

Example data

First I create some very simple example data: 4 species distributions, a cost layer, and a grid of 25km2 planning units, all defined on a 100km x 100km study area.

Conservation features

I create 4 artificial 1km resolution raster layers that I’ll use to represent species distributions. To work with data that are at least somewhat realistic I introduce some auto-correlation into the layers by generating spatially autocorrelated Gaussian random fields. I also give each of the four species different levels of rarity.

# raster template
utm10 <- crs('+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=km +no_defs')
r <- extent(c(0, 100, 0, 100)) %>% 
    raster(nrows = 100, ncols = 100, crs = utm10, vals = 1)

gaussian_field <- function(r, range, pct, beta = c(1, 0, 0)) {
  gsim <- gstat(formula = (z ~ x + y), dummy = TRUE, beta = beta, nmax = 20,
               model = vgm(psill = 1, range = range, model = 'Exp'))
  vals <- rasterToPoints(r, spatial = TRUE) %>% 
    geometry %>% 
    predict(gsim, newdata = ., nsim = 1) %>% 
    {scale(.$sim1)}
  if (!missing(pct)) {
    vals <- as.integer(vals < quantile(vals, pct))
  }
  r[] <- vals
  r
}

# conservation features
species <- mapply(function(x, y, r) gaussian_field(r = r, range = x, pct = y, 
                                                   beta = c(1, 0.01, 0.01)),
                  15, c(0.01, 0.05, 0.1, 0.15),
                  MoreArgs = list(r = r)) %>% 
  stack %>% 
  setNames(letters[1:4])
levelplot(species, main='Feature Distribution', layout = c(2, 2),
          scales=list(draw=FALSE),
          col.regions = c("grey20", "#fd9900"), colorkey = FALSE)

plot of chunk features

Cost

I use a similar approach to generating a spatially auto-correlated cost layer.

cost <- gaussian_field(r, 10, beta = c(1, 0, 0)) %>% 
  stretch(10, 1000) %>% 
  setNames("cost")
levelplot(cost, main = "Cost", margin = FALSE,
          col.regions = viridis, at = seq(0, 1000, len = 256))

plot of chunk cost

Planning units

In a previous post, I discussed the merits of hexagonal grids and defined a function to create such grids. I now use that function to generate a grid (see previous post for definition of make_grid). Since I’ll be performing the simulated annealing within R, I want to keep the parameter space small so the problem is feasible. Therefore, I divide the study area up in to just 378 25km2 planning units.

For completeness I randomly assign a patch of 7 planning units to be locked in and a patch of 7 to be locked out. Locked in planning units are typically those already within protected areas, while locked out planning units could be developed for another land use (e.g. urban or plantation) and therefore be unavailable for protection.

sa_border <- as(extent(r), "SpatialPolygons")
projection(sa_border) <- projection(r)
pu <- make_grid(sa_border, "hexagonal", cell_area = 25, clip = FALSE)
pu <- pu[gContains(sa_border, pu, byid = TRUE), ]
pu$id <- 1:length(pu)
row.names(pu) <- as.character(pu$id)
pu$status <- 0L
# randomly create patches with status = 2 and 3
inner_cells <- gContainsProperly(gUnaryUnion(pu), pu, byid = TRUE) %>% 
  which %>% 
  sample(2)
locked_in <- gIntersects(pu[inner_cells[1],], pu, byid = TRUE)
pu$status[locked_in] <- 2L
locked_out <- gIntersects(pu[inner_cells[2],], pu, byid = TRUE)
pu$status[locked_out] <- 3L

Now I average the cost raster layer over the hexagonal planning units.

pu <- raster::extract(cost, pu, fun = mean, na.rm = TRUE, sp = TRUE)
spplot(pu, "cost", main = "Planning Unit Cost", col.regions = viridis(256),
       at = seq(floor(min(pu$cost)), ceiling(max(pu$cost)), len = 256))

plot of chunk pu-cost

Implementation in R

Now that I’ve given an overview of the problem, and generated some data, it’s time to explore the details of Marxan and reconstruct them in R. I break down each component in turn, defining functions for each task.

Pre-processing

Some pre-processing needs to be done to get the spatial data into a form useful for simulated annealing.

Feature representation

The spatial distribution of features over the landscape typically comes in the form of a series of raster layers. Here I create a function that summarizes these layers to a data frame of species representation within each planning unit, which is essentially the \( r_{ij} \) term in the above objective function and the puvspr.dat Marxan input file.

summarize_features <- function(features, pu) {
  raster::extract(features, pu, fun = sum, na.rm = TRUE, sp = FALSE) %>% 
    data.frame(pu = pu$id, .) %>% 
    gather(feature, amount, -pu) %>% 
    filter(amount > 0)
}
puvspr <- summarize_features(species, pu)
head(puvspr, 5) %>% 
  kable
pu feature amount
7 a 1
9 a 6
21 a 1
119 a 3
132 a 2

Setting targets

Prior to a Marxan run, absolute representation targets for each conservation feature must be set. This may occur through expert elicitation, modeling exercise to determine the amount of habitat necessary to ensure persistence, or more informal methods. Regardless, it’s very common to set proportional targets, such as protecting 20% of the range of a given species or habitat type. I define a function that converts proportional targets to absolute representation targets.

In addition, this function sets the Species Penalty Factor (SPF) for each feature. I will explain these in detail below but, for now, they can be thought of as measuring the relative importance of meeting the representation target for each feature.

The data frame output by this function corresponds to the spec.dat Marxan input file.

prepare_features <- function(rij, target_prop, spf = 1) {
  assert_that(length(target_prop) == 1 || is.data.frame(target_prop),
              length(spf) == 1 || is.data.frame(spf))
  
  rij_sum <- group_by(rij, feature) %>% 
    summarize(amount = sum(amount))
  
  # set targets
  if (is.data.frame(target_prop)) {
    assert_that("feature" %in% names(target_prop),
                all(rij_sum$feature %in% target_prop$feature))
    features <- inner_join(rij_sum, target_prop, by = "feature") %>% 
      mutate(target = prop * amount)
  } else {
    features <- mutate(rij_sum, prop = target_prop, target = prop * amount)
  }
  
  # set spf
  if (is.data.frame(spf)) {
    assert_that("feature" %in% names(spf),
                all(target_prop$feature %in% spf$feature))
    features <- inner_join(features, spf, by = "feature")
  } else {
    features <- mutate(features, spf = spf)
  }
  dplyr::select(features, feature, prop, total_amount = amount, target, spf)
}
# 20% target for all species 
prepare_features(puvspr, 0.2, 10) %>% 
  kable
feature prop total_amount target spf
a 0.2 100 20.0 10
b 0.2 489 97.8 10
c 0.2 929 185.8 10
d 0.2 1407 281.4 10
# species specific targets
ex_targets <- data_frame(feature = letters[1:4], prop = c(0.4, 0.3, 0.2, 0.1))
ex_spf <- data_frame(feature = letters[1:4], spf = rep(1, 4))
spec <- prepare_features(puvspr, ex_targets, ex_spf)
kable(spec)
feature prop total_amount target spf
a 0.4 100 40.0 1
b 0.3 489 146.7 1
c 0.2 929 185.8 1
d 0.1 1407 140.7 1

Calculate shared boundaries

In addition to minimizing cost, Marxan can also be instructed to minimize the length of the external boundary of the reserve network. This external boundary can be calculated explicitly at each simulated annealing iteration; however, it’s more efficient to calculate a dataframe of all the pairwise shared boundaries once, then use this to calculate the overall boundary length at each iteration.

calculate_boundary <- function(pu) {
  row.names(pu) <- as.character(pu$id)
  # calculate all boundaries
  boundary <- geometry(pu) %>% 
    as("SpatialLines") %>% 
    gIntersection(., ., byid = TRUE) %>% 
    gLength(byid = TRUE) %>% 
    data_frame(pid = names(.), boundary = .) %>% 
    separate(pid, c("id1", "id2"), sep = " ") %>% 
    mutate(id1 = as(id1, class(pu$id)), id2 = as(id2, class(pu$id))) %>% 
    # don't double count
    filter(id1 <= id2)
  # used to avoid rounding errors
  min_boundary <- 1e-6 * min(boundary$boundary)
  # internal boundaries
  internal_boundary <- filter(boundary, id1 != id2) %>% 
    rbind(., mutate(., id1 = id2)) %>% 
    group_by(id1) %>% 
    summarize(internal = sum(boundary))
  # subtract internal boundaries from boundary of pu with itself
  # to get external boundary
  external_boundary <- filter(boundary, id1 == id2) %>% 
    left_join(internal_boundary, "id1") %>% 
    mutate(boundary = (boundary - ifelse(is.na(internal), 0, internal))) %>% 
    dplyr::select(-internal) %>% 
    mutate(boundary = ifelse(boundary < min_boundary, 0, boundary))
  boundary <- filter(boundary, id1 != id2) %>% 
    rbind(external_boundary) %>% 
    arrange(id1, id2) %>% 
    filter(boundary > 0)
  return(boundary)
}
boundary <- calculate_boundary(pu)
head(boundary) %>% 
  kable(digits = 3)
id1 id2 boundary
1 1 9.306
1 2 3.102
1 19 3.102
1 20 3.102
2 2 6.204
2 3 3.102

Simulated annealing

Next, I’ll delve into simulated annealing, which I introduced above. The function optim in the stats package provides implementations for a variety of optimization methods, including simulated annealing. For efficiency, each of the optimization methods is actually implemented in C and called via R. The C source code is available on GitHub.

Annealing schedule

Recall that the key to simulated annealing is that “bad” changes are accepted with some probability that declines as the heuristic progresses. This acceptance probability is given by:

where \( \vec{x}_n \) is the binary decision variable at annealing iteration \( n \) and \( T_m \) is the temperature parameter after \( m \) decreases. The two different indices, \( n \) and \( m \), arise because, at each temperature, there are typically multiple annealing iterations performed. The annealing schedule determines the rate at which temperature decreases, and alternate implementations of simulated annealing use different functional forms for the schedule. In particular, Marxan uses

and optim uses

where \( T_0 \) is the starting temperature and \( \alpha \) is a cooling factor, both of which are user-defined parameters.

Vanilla R implementation

Since optim uses a different annealing schedule, and the inner workings are obscured because it uses C, I implement my own Marxan-like simulated annealing in plain R.

anneal <- function(x_init, objective, neighbour, n, ntemp, t_init, alpha) {
  assert_that(is.count(n), is.count(ntemp),
              ntemp > 1, ntemp <= n, 
              n %% ntemp == 0)
  
  n_per_temp <- floor(n / ntemp)
  t <- t_init
  x <- x_init
  for (i in 1:n) {
    if (i %% n_per_temp == 0) {
      # decrease temperature
      t <- alpha * t
    }
    x_try <- neighbour(x)
    delta <- objective(x_try) - objective(x)
    if (delta <= 0) {
      # "good" change
      x <- x_try
    } else if (runif(1) < exp(-delta / t)) {
      # "bad" change
      x <- x_try
    }
  }
  return(x)
}

In this implementation, objective and neighbour are the objective function and neighbour selection function respectively. x_init is a vector of initial values for the decision variable, n is the number of iterations, ntemp is the number of temperature decreases, t_init is the starting temperature (\( T_0 \)), and alpha is the cooling factor (\( \alpha \)).

To keep this simulated annealing implementation as flexible as possible, the only argument the objective and neighbour selection functions take is the vector of decision variables. In general, other data or parameters will be required to evaluate these functions, but this additional information must be stored within the enclosing environment of the function. This can be accomplished using closures, functions created by other functions, which have access to all variables in the parent function. This will become more clear when I define these closure below. The Advanced R book by Hadley Wickham has great coverage of this topic, particularly the sections on function environments and closures

Initial reserve system

Simulated annealing needs to start somewhere. Marxan randomly picks a user-defined proportion of planning units to be included in the initial reserve system. I define a function that, given a set of planning units and a proportion, returns a logical vector indicating which units are selected for the initial reserve. In addition, the planning unit status described above is taken into account when choosing which units to include in the initial reserve.

init_res <- function(pu, prop = 0) {
  n <- round(prop * nrow(pu))
  # account for locked in planning units
  x <- pu$status %in% 1:2
  if (sum(x) >= n) {
    return(x)
  }
  # add more to meet specified proportion
  n <- n - sum(x)
  add <- which(pu$status == 0) %>% 
    sample(n)
  x[add] <- TRUE
  return(x)
}
sum(init_res(pu, 0.5)) / length(pu)
#> [1] 0.5

Neighbour selection function

At each iteration, simulated annealing must have some means of moving from the current candidate state to a neighbouring candidate state. In the case of Marxan, this is as simple as switching the binary decision variable for a single, randomly chosen planning unit. Again, planning unit status must be taken into account since some units are not permitted to be included.

As noted above, the neighbour function must be a closure taking just one argument: the vector of decision variables. The list of planning units that can be switched (i.e. those with status 0 or 1) is stored in the enclosing environment. To make this all happen, I define a generator function for Marxan-like neighbour selection closures.

generate_nsf <- function(pu) {
  x_valid <- (1:nrow(pu))[pu$status %in% 0:1]
  rm(pu)
  function(x) {
    i <- sample(x_valid, 1)
    x[i] <- !x[i]
    return(x)
  }
}
ex_nsf <- generate_nsf(pu)
ex_res <- init_res(pu, 0)
neigh <- ex_nsf(ex_res)
sum(ex_res - neigh)
#> [1] -1

Objective function

Recall that the problem Marxan seeks to solve can be formulated mathematically in terms of the minimization of an objective function subject to the constraint that all representation targets are met:

The fist terms in the objective function correspond to minimizing the cost, while second and third correspond to minimizing the perimeter of the reserve network. Rather than strictly enforcing the representation constraint, Marxan incorporates the constraint into the objective function as a target shortfall penalty. This approach is more efficient because the constraint does not have to be enforced in each simulated annealing iteration. Furthermore, it allows more flexibility within simulated annealing, by temporarily allowing solutions that don’t meet targets.

The final term is the shortfall penalty. For each feature for which the target isn’t met, this penalty is proportional to the relative shortfall to the target. The parameter \( s_j \) is the Species Penalty Factor (SPF), which determines the magnitude of the penalty to apply if the target isn’t met for a given feature. A larger value of the SPF will put more emphasis on meeting targets for that feature relative to minimizing cost and boundary length. Finally, \( \gamma_j \) is the base shortfall penalty, which essentially acts to convert the shortfall penalty into units of cost. Marxan calculates these values as the cost of meeting a given feature’s target in the cheapest way possible, ignoring all other features. This results in the shortfall penalty being an approximation of the cost required to raise a feature’s level of representation to the target level.

Base shortfall penalty

Marxan uses a simple greedy algorithm for estimating \( \gamma_j \). For each feature, planning units are ranked in descending order of efficiency, i.e. the ratio of representation level for that feature and cost (economic cost plus boundary length). Planning units are then added in order until the target is met, and the resulting cost is \( \gamma_j \). This doesn’t guarantee the cheapest way to meet a feature’s target, but it give a good approximation.

The approach described in the previous paragraph is what the Marxan User’s Manual says Marxan is doing; however, looking at the source code, there’s an additional wrinkle. Rather that just adding planning units in order of efficiency, Marxan will also look to see if, at any point, the addition of a single planning unit will result in the target being met. If so, and the cost of this planning unit is lower than the cost of the most efficient unit, than this planning unit is chosen instead of the most efficient one.

base_penalty <- function(pu, rij, features, blm) {
  # calculate boundary of each planning unit
  pu$boundary <- gLength(pu, byid = TRUE)
  # calculate efficiency = representation / cost for each planning unit
  remaining <- as.data.frame(pu) %>% 
    mutate(total_cost = (cost + blm * boundary)) %>% 
    inner_join(rij, ., by = c("pu" = "id")) %>% 
    filter(amount > 0, status != 3) %>% 
    mutate(efficiency = amount / total_cost)
  # start by including free planning units and those with status = 2
  p <- filter(remaining, status == 2 | total_cost == 0) %>% 
    group_by(feature) %>% 
    summarize(cumulative_amount = sum(amount),
              penalty = sum(total_cost)) %>% 
    right_join(dplyr::select(features, feature, target), by = "feature") %>% 
    mutate(shortfall = (target - cumulative_amount),
           shortfall = ifelse(is.na(shortfall), target, shortfall),
           penalty = ifelse(is.na(penalty), 0, penalty)) %>% 
    dplyr::select(feature, penalty, shortfall)
  remaining <- filter(remaining, status != 2, total_cost != 0)
  for (f in p$feature) {
    # find and add next best planning unit
    r <- filter(remaining, feature == f) %>% 
      arrange(desc(efficiency))
    while(p$shortfall[p$feature == f] > 0 && nrow(r) > 0) {
      # most efficient planning unit
      best <- head(r, 1)
      # cheapest planning unit that meets target
      meets_target <- r %>% 
        filter(amount >= p$shortfall[p$feature == f],
               total_cost < best$total_cost) %>% 
        arrange(total_cost) %>% 
        head(1)
      if (nrow(meets_target) > 0) {
        best <- meets_target
      }
      p$shortfall[p$feature == f] <- p$shortfall[p$feature == f] - best$amount
      p$penalty[p$feature == f] <- p$penalty[p$feature == f] + best$total_cost
      r <- filter(r, pu != best$pu)
    }
  }
  dplyr::select(p, feature, penalty)
}
base_penalty(pu, puvspr, spec, blm = 10) %>% 
  kable(digits = 1)
feature penalty
a 2171.8
b 3820.4
c 6527.3
d 6846.2

Getting the base penalty exactly right was fairly challenging, since what Marxan actually does is more complicated than what is outlined in the documentation. I modified the Marxan source code to output the set of planning units that it chooses as the cheapest way to meet each single-species target. Then I tweaked and debugged my function until it found the same set of reserves. Along the way I found what I believe is a minor bug in Marxan’s source code, which causes the final planning unit added for each species to be incorrect. This results in a slightly higher base penalty than intended. Although I believe this has no impact on Marxan’s ability to find optimal reserve configurations, I’ve corrected the bug in the above function, which means there will be slight differences in the objective function values between my R emulator and Marxan.

Also, Marxan has two different functions for calculating the base penalty: CalcPenalties and CalcPenaltiesOptimise. The former is used by default and the latter is used only if the file puvspr_sporder.dat is provided as an input. This additional file contains the same information as puvspr.dat (i.e. representation level of each species in each planning unit), but it is ordered by species ID as opposed to planning unit ID. Apparently, CalcPenaltiesOptimise is much faster for large problems, but the catch is they don’t’ give the same penalties except in the case where no reserves are locked in. Since the R package marxan always provides a puvspr_sporder.dat file and therefore Marxan always uses CalcPenaltiesOptimise when being run through R, I’ve written the above function to match CalcPenaltiesOptimise, except for the issue with the final planning unit noted above.

Objective function implementation

As with the neighbour selection function defined earlier, I implement the objective function as a closure that takes the decision variable as an argument and store all other required data and parameters in the enclosing environment. Again, to make this happen, I define a generator function for Marxan-like objective function closures.

generate_objective <- function(pu, rij, features, boundary, blm = 0) {
  # calculate base shortfall penalties
  features <- base_penalty(pu = pu, rij = rij, features = features, blm = blm) %>% 
    inner_join(features, by = "feature")
  function(x, breakdown = FALSE) {
    if (all(!x)) {
      cost <- 0
      boundary <- 0
      shortfall <- sum(targets$spf * features$penalty)
    } else {
      pu_selected <- pu[x,]
      # soci-economic cost
      cost <- sum(pu_selected$cost)
      # boundary length
      if (blm == 0) {
        boundary_cost <- 0
      } else {
        boundary_cost <- filter(boundary, 
            (id1 %in% pu_selected$id & !id2 %in% pu_selected$id) |
            (id2 %in% pu_selected$id & !id1 %in% pu_selected$id) |
            (id1 %in% pu_selected$id & id1 == id2)) %>% 
          {blm * sum(.$boundary)}
      }
      # shortfall penalty
      shortfall <- data_frame(pu = pu_selected$id) %>% 
        inner_join(rij, by = "pu") %>% 
        group_by(feature) %>% 
        summarize(amount = sum(amount)) %>% 
        right_join(features, by = "feature") %>% 
        mutate(amount = ifelse(is.na(amount), 0, amount),
               shortfall = spf * penalty * pmax(1 - amount / target, 0)) %>% 
               {sum(.$shortfall)}
    }
    score <- cost + boundary_cost + shortfall
    if (breakdown) {
      return(data.frame(score = score, cost = cost,
                        boundary = boundary_cost, 
                        shortfall = shortfall))
    } else {
      return(score)
    }
  }
}
ex_objective <- generate_objective(pu, puvspr, spec, boundary, blm = 1)
ex_objective(init_res(pu, 0.5))
#> [1] 96087.66
ex_objective(init_res(pu, 0.5), breakdown = TRUE)
#>      score     cost boundary shortfall
#> 1 99428.08 96922.43 1867.414  638.2427

Adaptive annealing

In the simulated annealing formulation Marxan uses, there are two parameters that define the annealing schedule: \( T_0 \) is the initial temperature and \( \alpha \) is the cooling factor. The optimal values for these parameters depend on the specific objective function being minimized. Given this, Marxan offers the option of setting these parameters adaptively based on the objective function. In the language of Marxan, this is referred to adaptive annealing in contrast to fixed annealing in which the user sets these parameters directly.

Marxan chooses sensible values for these parameters by iteratively flipping the status of one planning unit at a time, calculating the resulting change in objective function, and keeping track of the maximum and minimum size of “bad” changes (i.e. those that increased the objective function). The initial temperature is set to the maximum bad change and the final temperature is set to 10% of the minimum bad change. The cooling factor is then inferred from the initial and final temperatures.

adaptive_schedule <- function(n, objective, nsf, x_init) {
  epsilon <- 1e-10
  scores <- numeric(n)
  x <- x_init
  scores[1] <- objective(x)
  for (i in 2:n) {
    x <- nsf(x)
    scores[i] <- objective(x)
  }
  delta <- diff(scores)
  delta <- delta[delta > epsilon]
  c(t_initial = max(delta), t_final = 0.1 * min(delta))
}
adaptive_schedule(1000, ex_objective, ex_nsf, init_res(pu, 0.5))
#>  t_initial    t_final 
#> 864.751510   8.428172

Marxan!

All the pieces are now in place and it’s time to build the Marxan emulator!

marxan_emulator <- function(pu, feature_stack, target_prop,
                            spf = 1, blm = 0, init_prop = 0,
                            nrep = 1, n = 10000, ntemp = n / 10) {
  # representation level of features in planning units
  rij <- summarize_features(feature_stack, pu)
  # set targets and spf
  features <- prepare_features(rij, target_prop, spf)
  # calculate boundaries
  boundary <- calculate_boundary(pu)
  # neighbour function
  nsf <- generate_nsf(pu)
  # objective function
  objective <- generate_objective(pu, rij, features, boundary, blm)
  # adaptive annealing, find annealing parameters
  temp <- adaptive_schedule(max(1000, n / 100), objective, nsf,
                                x_init = init_res(pu, prop = init_prop))
  alpha <- unname(exp(log(temp["t_final"] / temp["t_initial"]) / ntemp))
  s <- data.frame(run = integer(), pu = integer(), x = logical())
  # run simulated annealing for multiple repetitions to get a suite of
  # candidate reserve networks to analyze
  for (i in 1:nrep) {
    #print(paste0("Annealing repetition ", i))
    x_result <- anneal(x_init = init_res(pu, prop = init_prop),
           objective = objective, neighbour = nsf, n = n, 
           ntemp = ntemp, t_init = temp["t_initial"], alpha = alpha)
    s <- rbind(s, data.frame(run = i, pu = pu$id, x = x_result))
  }
  selection_summary <- s %>% 
    arrange(run, pu) %>% 
    group_by(run) %>% 
    do(objective(.$x, breakdown = TRUE))
  list(pu = pu, features = features, boundary = boundary,
       objective = objective, nsf = nsf,
       params = c(temp, alpha = alpha, nrep = nrep, n = n, ntemp = ntemp),
       blm = blm,
       selections = s,
       summary = selection_summary)
}

Reserve selection

Now that I’ve created an R function that emulates Marxan, let’s to put it to the test. I’ll compare the results of reserve selection exercises using Marxan and R with identical parameters. Since there is inherent stochasticity in simulated annealing, the results won’t be identical, but they should be quite similar if I’ve implemented everything correctly.

Preparations

In generating the species distributions at the start of this tutorial, I created four levels of rarity: species cover 1%, 5%, 10%, or 15% of the study area. Working under the assumption that rare species need more of their range protected, I set percent representation targets for species in the different rarity categories at 40%, 30%, 20%, and 10%, respectively.

target_prop <- data_frame(feature = letters[1:4], prop = c(0.4, 0.3, 0.2, 0.1))
kable(target_prop)
feature prop
a 0.4
b 0.3
c 0.2
d 0.1

The Marxan Good Practices Handbook outlines an iterative approach for determining the optimal Species Penalty Factors. Essentially, one runs Marxan many times for different values of the SPF to find the minimum value for each species below which targets stop being met. For this tutorial, I won’t delve into calibration of the SPF, rather I’ve done this outside the tutorial and enter some sensible values here.

spfs <- data_frame(feature = letters[1:4], spf = c(10, 25, 15, 15))
kable(spfs)
feature spf
a 10
b 25
c 15
d 15

The optimal Boundary Length Modifier can be found in a similar iterative fashion. Again, this tutorial isn’t about Marxan best practices, so I won’t delve into calibrating the BLM.

best_blm <- 10L

Emulator

I run the Marxan emulator for 100 simulated annealing repetitions with 10,000 iterations each. Typically more iterations are used, but since this is just for demonstration purposes, and the R emulator is super inefficient, I keep the number of iterations low.

res_emulator <- marxan_emulator(pu, species, target_prop, spfs,
                                blm = best_blm, nrep = 100, n = 10000)

Now I map the best reserve, i.e. the one with the lowest objective function:

map_selection <- function(x, pu, title = NULL) {
  status_names <- c("Excluded", "Selected", "Locked in", "Locked out")
  status_pal <- c("grey40", "#4DAF4A", "#377EB8", "black")
  pu$selected <- ifelse(pu$status %in% 2:3,
                        pu$status,
                        as.integer(x)) %>% 
    factor(levels = 0:3, labels = status_names)
  spplot(pu, "selected", main = title, col.regions = status_pal, 
         colorkey = list(space = "bottom", height = 1)) %>% 
    print
}
best <- which.max(res_emulator$summary$cost)
filter(res_emulator$selections, run == best) %>% 
  {.$x} %>% 
  map_selection(res_emulator$pu, title = "Best Reserve (Emulator)")

plot of chunk map-emulator

And the selection frequency, which is the frequency with which each planning unit is selected, and gives a sense for the importance of each planning unit.

map_frequency <- function(selections, pu, title = "Selection Frequency") {
  pu <- group_by(selections, pu) %>% 
    summarize(frequency = sum(x) / n()) %>% 
    rename(id = pu) %>% 
    merge(pu, ., by = "id")
  spplot(pu, "frequency", main = title, col.regions = viridis(100),
         at = seq(0, 1.01, len = 100),
         colorkey = list(space = "bottom", height = 1))
}
map_frequency(res_emulator$selections, res_emulator$pu,
              title = "Selection Frequency (Emulator)")

plot of chunk map-emulator-sf

Marxan

Now I perform a normal Marxan run, with identical parameters, which I will compare to my R emulation. The R package marxan makes running Marxan from within R a breeze and has some great vignettes: browseVignettes("marxan").

res_marxan <- marxan(pu, species,
                     targets = res_emulator$features$target,
                     spf = res_emulator$features$spf,
                     BLM = 10, PROP = 0, VERBOSITY = 3L,
                     NUMREPS = 100L, NUMITNS = 10000L, NUMTEMP = 1000L)

Again, I map the best reserve and the selection frequency.

res_marxan@results@selections[res_marxan@results@best, ] %>% 
  map_selection(res_emulator$pu, title = "Best Reserve (Marxan)")

plot of chunk marxan-best

The best reserve is quite different between methods, which is to be expected given that simulated annealing is inherently stochastic and I’m using a small number of iterations. The selection frequency on the other hand should be fairly similar since it essentially averages over 100 solutions.

res_marxan@results@selections %>%
  data.frame(run = 1:nrow(.), .) %>%
  gather(pu, x, -run) %>% 
  mutate(pu = as.integer(sub("P", "", pu))) %>% 
  map_frequency(res_emulator$pu, title = "Selection Frequency (Marxan)")

plot of chunk marxan-freq

And, as expected, these results are very similar those from the emulator.

Comparison

Finally, I compare the results to demonstrate that the emulator I’ve created functions identically to Marxan. I’ll compare the two methods with respect.

Objective function

Getting the objective function right is the critical component to the whole reserve selection process. Recall that there are three components to the objective function: the socio-economic cost, the boundary length penalty, and the target shortfall penalty. I want to check that all three are being calculated correctly and, in addition to the total objective function value, Marxan does return the values of the individual components for the final reserve. However, in a typical Marxan run, the targets will all be met in the final reserve, so the shortfall term will be zero. To get around this, I re-run Marxan with a small number of simulated annealing iterations to intentionally produce a “bad” solution that doesn’t meet targets.

In the MarxanResults object that the marxan package returns, the individual components of the objective function appear in the summary table:

  • Score: full objective function
  • Cost: cost term
  • Connectivity: boundary length penalty prior to multiplication by BLM
  • Penalty: shortfall penalty
res_bad <- update(res_marxan, ~opt(NUMITNS = 500L, NUMTEMP = 100L))
marx_obj_components <- res_bad@results@summary %>% 
  dplyr::select(run = Run_Number, Score, Cost, Connectivity, Penalty)

Now I use my objective function to calculate these components and compare them to those from Marxan. Of the 100 runs, I display those for which there is a mismatch between Marxan and the emulator.

compare_components <- res_bad@results@selections %>% 
  adply(1, res_emulator$objective, breakdown = TRUE, .id = "run") %>%
  mutate(run = as.integer(run),
         # divide by blm to match marxan term
         boundary = boundary / best_blm) %>% 
  inner_join(marx_obj_components, ., by = "run") %>%
  # calcualte differences between components
  mutate(score_diff = (score - Score),
         cost_diff = (cost - Cost),
         boundary_diff = (boundary - Connectivity),
         shortfall_diff = (shortfall - Penalty)) %>% 
  dplyr::select(run, score_diff, cost_diff, boundary_diff, shortfall_diff) 
# identify mismatches between emulator and marxan proper
filter(compare_components, score_diff > 1e-6 | cost_diff > 1e-6 |
       boundary_diff > 1e-6 | shortfall_diff > 1e-6) %>% 
  kable
run score_diff cost_diff boundary_diff shortfall_diff
10 29.2609006 -3e-07 -3e-07 29.2609
16 29.2609009 -4e-07 1e-07 29.2609
37 0.0000011 -1e-07 -2e-07 0.0000
44 29.2609008 4e-07 1e-07 29.2609
58 29.2609005 -4e-07 -2e-07 29.2609
62 0.0000011 4e-07 2e-07 0.0000
71 0.0000011 4e-07 5e-07 0.0000
84 29.2609001 1e-07 -1e-07 29.2609

Apart from rounding errors, the only problems are with the target shortfall penalty. As noted above, Marxan appears to have a small bug in calculating the base shortfall penalty, which I’ve chosen to fix. This is the cause of the small difference in the shortfall penalty.

Adaptive annealing

Marxan adaptively sets the starting temperature and cooling factor, which determine the cooling schedule for simulated annealing. These parameters are not conveniently captured by the marxan R package; however, they are output to the screen as Marxan runs. So, I redirect this output to a file and extract the cooling parameters with regular expressions. Since there is stochasticity in the calculation of these parameters I use 100 Marxan runs to get a range of values.

# function to redirect C stdout
cppFunction('void redir(){FILE* F=freopen("/tmp/capture.txt","w+",stdout);}')
redir()
solve(res_marxan, force_reset = TRUE)
# read in marxan output
marx_stdout <- readLines("/tmp/capture.txt")
# extract lines with temperature parameters
marxan_temp <- marx_stdout[grepl("Tinit", marx_stdout)] %>% 
  str_match("^Run ([0-9]{1,3}) .* Tinit = ([.0-9]+) Tcool = ([.0-9]+)") %>% 
  data.frame(stringsAsFactors = FALSE) %>% 
  setNames(c("whole", "run", "t_initial", "alpha")) %>% 
  mutate(method = "marxan",
         t_initial = as.numeric(t_initial),
         alpha = as.numeric(alpha)) %>% 
  dplyr::select(method, run, t_initial, alpha)
unlink("/tmp/capture.txt")

I now get a sample of 100 cooling schedule parameters from the emulator for comparison, then plot a comparison of the results.

emulator_temp <- data.frame(run = 1:100) %>% 
  ddply("run", function(x) {
    adaptive_schedule(100, res_emulator$objective, res_emulator$nsf,
                      x_init = init_res(res_emulator$pu, 0))
  }) %>% 
  mutate(method = "emulator",
         alpha = exp(log(t_final / t_initial) / 1000)) %>%
  dplyr::select(method, run, t_initial, alpha)
  compare_temp <- rbind(marxan_temp, emulator_temp)
#plot distributions
p1 <- compare_temp %>%
  ggplot(aes(t_initial, ..scaled.., fill = method)) +
    geom_density(alpha = 0.75) +
    scale_fill_brewer(palette = "Set1") +
    scale_x_continuous(labels = comma) +
    labs(x = expression(Initial~Temperature~(T[0])), y = NULL) +
    theme(legend.position=c(.8, .8))
p2 <- compare_temp %>%
  ggplot(aes(alpha, ..scaled.., fill = method)) +
    geom_density(alpha = 0.75) +
    scale_fill_brewer(NULL, palette = "Set1", guide = FALSE) +
    labs(x = expression(Cooling~Factor~(alpha)), y = NULL)
grid.arrange(p1, p2, ncol = 2)

plot of chunk adaptive-check-emulator

Everything looks good, the emulator seems to match Marxan quite closely for both temperature parameters.

Selection frequency

The last thing I’ll check is the selection frequency. The stochastic nature of simulated annealing means each solution will be different, but the selection frequency should be similar.

# selection frequency - emulator
freq_compare <- group_by(res_emulator$selections, pu) %>% 
  summarize(freq_emulator = sum(x) / n())
# selection frequency - marxan
freq_compare <- res_marxan@results@selections %>%
  data.frame(run = 1:nrow(.), .) %>%
  gather(pu, x, -run) %>% 
  mutate(pu = as.integer(sub("P", "", pu))) %>% 
  group_by(pu) %>% 
  summarize(freq_marxan = sum(x) / n()) %>% 
  inner_join(freq_compare, by = "pu")
ggplot(freq_compare, aes(freq_marxan, freq_emulator)) +
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  geom_point() +
  geom_text(x = 0.75, y = 0.9, label = "y = x", color = "red", size = 6,
            family = "Helvetica Neue Light") +
  labs(x = "Marxan", y = "Emulator", title = "Selection Frequency") +
  theme(text=element_text(family = "Helvetica Neue Light"))

plot of chunk sf-compare

rename(freq_compare, id = pu, Marxan = freq_marxan, Emulator = freq_emulator) %>% 
  merge(pu, ., by = "id") %>% 
  spplot(c("Marxan", "Emulator"),
         main = "Selection Frequency", col.regions = viridis(100),
         at = seq(0, 1.01, len = 100),
         colorkey = list(space = "bottom", height = 1),
         par.strip.text = list(col = "white"),
         par.settings = list(
           strip.background = list(col = "grey40"))
         )

plot of chunk side-by-side

Looks good, we’ve got a pretty close match here.

Performance

Given that Marxan is written in C and presumably efforts were made to optimize the code, asking how this R emulator stacks up in terms of performance is almost pointless. However, I am curious just how slow the R code is, so let’s take a look. I’ll run 5 repetitions of the same problem we’ve already looked at: 378 planning units, 4 species, and 10,000 simulated annealing iterations. Note that by Marxan standards this is an extremely small problem, usually thousands of planning units, dozens of conservation features, and millions of iterations are used.

# emulator
time_emulator <- system.time(
  marxan_emulator(pu, species, target_prop, spfs,
                  blm = best_blm, nrep = 5, n = 10000)
)
# marxan
time_marxan <- system.time(
  marxan(pu, species,
         targets = res_emulator$features$target,
         spf = res_emulator$features$spf,
         BLM = 10, PROP = 0, NUMREPS = 5L, 
         NUMITNS = 10000L, NUMTEMP = 1000L)
)
time_emulator / time_marxan
#>     user   system  elapsed 
#> 140.4495 180.9157 158.8270

So, as expected, the performance is terrible. Marxan is 138 times faster!

Conclusions

In summary, I’ve built a series of pure R functions that solve reserve design problems using simulated annealing in a way that emulates the basic functionality of Marxan. This R emulator is of no practical value because it’s extremely inefficient. However, this isn’t the point of the emulator. Rather I wrote it for fun and as a means of learning what Marxan is doing under the hood. And, on those points, I think I’ve succeeded!