## S4 method for signature 'ConservationProblem,ANY,ANY,matrix'
add_connectivity_penalties(x, penalty, zones, data)

## S4 method for signature 'ConservationProblem,ANY,ANY,Matrix'
add_connectivity_penalties(x, penalty, zones, data)

## S4 method for signature 'ConservationProblem,ANY,ANY,data.frame'
add_connectivity_penalties(x, penalty, zones, data)

## S4 method for signature 'ConservationProblem,ANY,ANY,dgCMatrix'
add_connectivity_penalties(x, penalty, zones, data)

## S4 method for signature 'ConservationProblem,ANY,ANY,array'
add_connectivity_penalties(x, penalty, zones, data)

# set seed for reproducibility
set.seed(600)

# load Matrix package for visualizing matrices
require(Matrix)

# load data
data(sim_pu_polygons, sim_pu_zones_stack, sim_features, sim_features_zones)

# define function to rescale values between zero and one so that we
# can compare solutions from different connectivity matrices
rescale <- function(x, to = c(0, 1), from = range(x, na.rm = TRUE)) {
  (x - from[1]) / diff(from) * diff(to) + to[1]
}

# create basic problem
p1 <- problem(sim_pu_polygons, sim_features, "cost") %>%
      add_min_set_objective() %>%
      add_relative_targets(0.2) %>%
      add_default_solver(verbose = FALSE)

# create a symmetric connectivity matrix where the connectivity between
# two planning units corresponds to their shared boundary length
b_matrix <- boundary_matrix(sim_pu_polygons)

# standardize matrix values to lay between zero and one
b_matrix[] <- rescale(b_matrix[])

# visualize connectivity matrix
## Not run: 
image(b_matrix)

## End(Not run)
# create a symmetric connectivity matrix where the connectivity between
# two planning units corresponds to their spatial proximity
# i.e. planning units that are further apart share less connectivity
centroids <- rgeos::gCentroid(sim_pu_polygons, byid = TRUE)
d_matrix <- (1 / (as(dist(centroids@coords), "Matrix") + 1))

# standardize matrix values to lay between zero and one
d_matrix[] <- rescale(d_matrix[])

# remove connections between planning units without connectivity to
# reduce run-time
d_matrix[d_matrix < 0.7] <- 0

# visualize connectivity matrix
## Not run: 
image(d_matrix)

## End(Not run)
# create a symmetric connectivity matrix where the connectivity
# between adjacent two planning units corresponds to their combined
# value in a field in the planning unit attribute data
# for example, this field could describe the extent of native vegetation in
# each planning unit and we could use connectivity penalties to identify
# solutions that cluster planning units together that both contain large
# amounts of native vegetation
c_matrix <- connectivity_matrix(sim_pu_polygons, "cost")

# standardize matrix values to lay between zero and one
c_matrix[] <- rescale(c_matrix[])

# visualize connectivity matrix
## Not run: 
image(c_matrix)

## End(Not run)
# create an asymmetric connectivity matrix. Here, connectivity occurs between
# adjacent planning units and, due to rivers flowing southwards
# through the study area, connectivity from northern planning units to
# southern planning units is ten times stronger than the reverse.
ac_matrix <- matrix(0, length(sim_pu_polygons), length(sim_pu_polygons))
ac_matrix <- as(ac_matrix, "Matrix")
adjacent_units <- rgeos::gIntersects(sim_pu_polygons, byid = TRUE)
for (i in seq_len(length(sim_pu_polygons))) {
  for (j in seq_len(length(sim_pu_polygons))) {
    # find if planning units are adjacent
    if (adjacent_units[i, j]) {
      # find if planning units lay north and south of each other
      # i.e. they have the same x-coordinate
      if (centroids@coords[i, 1] == centroids@coords[j, 1]) {
        if (centroids@coords[i, 2] > centroids@coords[j, 2]) {
          # if i is north of j add 10 units of connectivity
          ac_matrix[i, j] <- ac_matrix[i, j] + 10
        } else if (centroids@coords[i, 2] < centroids@coords[j, 2]) {
          # if i is south of j add 1 unit of connectivity
          ac_matrix[i, j] <- ac_matrix[i, j] + 1
        }
      }
    }
  }
}

# standardize matrix values to lay between zero and one
ac_matrix[] <- rescale(ac_matrix[])

# visualize asymmetric connectivity matrix
## Not run: 
image(ac_matrix)

## End(Not run)
# create penalties
penalties <- c(10, 25)

# create problems using the different connectivity matrices and penalties
p2 <- list(p1,
           p1 %>% add_connectivity_penalties(penalties[1], data = b_matrix),
           p1 %>% add_connectivity_penalties(penalties[2], data = b_matrix),
           p1 %>% add_connectivity_penalties(penalties[1], data = d_matrix),
           p1 %>% add_connectivity_penalties(penalties[2], data = d_matrix),
           p1 %>% add_connectivity_penalties(penalties[1], data = c_matrix),
           p1 %>% add_connectivity_penalties(penalties[2], data = c_matrix),
           p1 %>% add_connectivity_penalties(penalties[1], data = ac_matrix),
           p1 %>% add_connectivity_penalties(penalties[2], data = ac_matrix))

# assign names to the problems
names(p2) <- c("basic problem",
               paste0("b_matrix (", penalties,")"),
               paste0("d_matrix (", penalties,")"),
               paste0("c_matrix (", penalties,")"),
               paste0("ac_matrix (", penalties,")"))
## Not run: 
# solve problems
s2 <- lapply(p2, solve)

# plot solutions
par(mfrow = c(3, 3))
for (i in seq_along(s2)) {
  plot(s2[[i]], main = names(p2)[i], cex = 1.5, col = "white")
  plot(s2[[i]][s2[[i]]$solution_1 == 1, ], col = "darkgreen", add = TRUE)
}

## End(Not run)

# create minimal multi-zone problem and limit solver to one minute
# to obtain solutions in a short period of time
p3 <- problem(sim_pu_zones_stack, sim_features_zones) %>%
      add_min_set_objective() %>%
      add_relative_targets(matrix(0.15, nrow = 5, ncol = 3)) %>%
      add_binary_decisions() %>%
      add_default_solver(time_limit = 60, verbose = FALSE)

# create matrix showing which planning units are adjacent to other units
a_matrix <- adjacency_matrix(sim_pu_zones_stack)

# visualize matrix
## Not run: 
image(a_matrix)

## End(Not run)
# create a zone matrix where connectivities are only present between
# planning units that are allocated to the same zone
zm1 <- as(diag(3), "Matrix")

# print zone matrix
print(zm1)

# create a zone matrix where connectivities are strongest between
# planning units allocated to different zones
zm2 <- matrix(1, ncol = 3, nrow = 3)
diag(zm2) <- 0
zm2 <- as(zm2, "Matrix")

# print zone matrix
print(zm2)

# create a zone matrix that indicates that connectivities between planning
# units assigned to the same zone are much higher than connectivities
# assigned to different zones
zm3 <- matrix(0.1, ncol = 3, nrow = 3)
diag(zm3) <- 1
zm3 <- as(zm3, "Matrix")

# print zone matrix
print(zm3)

# create a zone matrix that indicates that connectivities between planning
# units allocated to zone 1 are very high, connectivities between planning
# units allocated to zones 1 and 2 are moderately high, and connectivities
# planning units allocated to other zones are low
zm4 <- matrix(0.1, ncol = 3, nrow = 3)
zm4[1, 1] <- 1
zm4[1, 2] <- 0.5
zm4[2, 1] <- 0.5
zm4 <- as(zm4, "Matrix")

# print zone matrix
print(zm4)

# create a zone matrix with strong connectivities between planning units
# allocated to the same zone, moderate connectivities between planning
# unit allocated to zone 1 and zone 2, and negative connectivities between
# planning units allocated to zone 3 and the other two zones
zm5 <- matrix(-1, ncol = 3, nrow = 3)
zm5[1, 2] <- 0.5
zm5[2, 1] <- 0.5
diag(zm5) <- 1
zm5 <- as(zm5, "Matrix")

# print zone matrix
print(zm5)

# create vector of penalties to use creating problems
penalties2 <- c(5, 15)

# create multi-zone problems using the adjacent connectivity matrix and
# different zone matrices
p4 <- list(
  p3,
  p3 %>% add_connectivity_penalties(penalties2[1], zm1, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[2], zm1, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[1], zm2, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[2], zm2, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[1], zm3, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[2], zm3, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[1], zm4, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[2], zm4, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[1], zm5, a_matrix),
  p3 %>% add_connectivity_penalties(penalties2[2], zm5, a_matrix))

# assign names to the problems
names(p4) <- c("basic problem",
               paste0("zm", rep(seq_len(5), each = 2), " (",
                      rep(penalties2, 2), ")"))
## Not run: 
# solve problems
s4 <- lapply(p4, solve)
s4 <- lapply(s4, category_layer)
s4 <- stack(s4)

# plot solutions
plot(s4, main = names(p4), axes = FALSE, box = FALSE)

## End(Not run)

# create an array to manually specify the connectivities between
# each planning unit when they are allocated to each different zone
# for real-world problems, these connectivities would be generated using
# data - but here these connectivity values are assigned as random
# ones or zeros
c_array <- array(0, c(rep(ncell(sim_pu_zones_stack[[1]]), 2), 3, 3))
for (z1 in seq_len(3))
  for (z2 in seq_len(3))
    c_array[, , z1, z2] <- round(runif(ncell(sim_pu_zones_stack[[1]]) ^ 2,
                                       0, 0.505))

# create a problem with the manually specified connectivity array
# note that the zones argument is set to NULL because the connectivity
# data is an array
p5 <- list(p3,
           p3 %>% add_connectivity_penalties(15, zones = NULL, c_array))


# assign names to the problems
names(p5) <- c("basic problem", "connectivity array")
## Not run: 
# solve problems
s5 <- lapply(p5, solve)
s5 <- lapply(s5, category_layer)
s5 <- stack(s5)

# plot solutions
plot(s5, main = names(p5), axes = FALSE, box = FALSE)

## End(Not run)