if (FALSE) {
# load package
require(ape)
# load data
sim_pu_raster <- get_sim_pu_raster()
sim_features <- get_sim_features()
sim_phylogeny <- get_sim_phylogeny()
sim_zones_pu_raster <- get_sim_zones_pu_raster()
sim_zones_features <- get_sim_zones_features()
# create minimal problem that aims to maximize the number of features
# adequately conserved given a total budget of 3800. Here, each feature
# needs 20% of its habitat for it to be considered adequately conserved
p1 <-
problem(sim_pu_raster, sim_features) %>%
add_max_features_objective(budget = 3800) %>%
add_relative_targets(0.2) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# create weights that assign higher importance to features with less
# suitable habitat in the study area
w2 <- exp((1 / terra::global(sim_features, "sum", na.rm = TRUE)[[1]]) * 200)
# create problem using rarity weights
p2 <- p1 %>% add_feature_weights(w2)
# create manually specified weights that assign higher importance to
# certain features. These weights could be based on a pre-calculated index
# (e.g., an index measuring extinction risk where higher values
# denote higher extinction risk)
w3 <- c(0, 0, 0, 100, 200)
p3 <- p1 %>% add_feature_weights(w3)
# solve problems
s1 <- c(solve(p1), solve(p2), solve(p3))
names(s1) <- c("equal weights", "rarity weights", "manual weights")
# plot solutions
plot(s1, axes = FALSE)
# plot the example phylogeny
par(mfrow = c(1, 1))
plot(sim_phylogeny, main = "simulated phylogeny")
# create problem with a maximum phylogenetic diversity objective,
# where each feature needs 10% of its distribution to be secured for
# it to be adequately conserved and a total budget of 1900
p4 <-
problem(sim_pu_raster, sim_features) %>%
add_max_phylo_div_objective(1900, sim_phylogeny) %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# solve problem
s4 <- solve(p4)
# plot solution
plot(s4, main = "solution", axes = FALSE)
# find out which features have their targets met
r4 <- eval_target_coverage_summary(p4, s4)
print(r4, width = Inf)
# plot the example phylogeny and color the represented features in red
plot(
sim_phylogeny, main = "represented features",
tip.color = replace(
rep("black", terra::nlyr(sim_features)), which(r4$met), "red"
)
)
# we can see here that the third feature ("layer.3", i.e.,
# sim_features[[3]]) is not represented in the solution. Let us pretend
# that it is absolutely critical this feature is adequately conserved
# in the solution. For example, this feature could represent a species
# that plays important role in the ecosystem, or a species that is
# important commercial activities (e.g., eco-tourism). So, to generate
# a solution that conserves the third feature whilst also aiming to
# maximize phylogenetic diversity, we will create a set of weights that
# assign a particularly high weighting to the third feature
w5 <- c(0, 0, 10000, 0, 0)
# we can see that this weighting (i.e., w5[3]) has a much higher value than
# the branch lengths in the phylogeny so solutions that represent this
# feature be much closer to optimality
print(sim_phylogeny$edge.length)
# create problem with high weighting for the third feature and solve it
s5 <- p4 %>% add_feature_weights(w5) %>% solve()
# plot solution
plot(s5, main = "solution", axes = FALSE)
# find which features have their targets met
r5 <- eval_target_coverage_summary(p4, s5)
print(r5, width = Inf)
# plot the example phylogeny and color the represented features in red
# here we can see that this solution only adequately conserves the
# third feature. This means that, given the budget, we are faced with the
# trade-off of conserving either the third feature, or a phylogenetically
# diverse set of three different features.
plot(
sim_phylogeny, main = "represented features",
tip.color = replace(
rep("black", terra::nlyr(sim_features)), which(r5$met), "red"
)
)
# create multi-zone problem with maximum features objective,
# with 10% representation targets for each feature, and set
# a budget such that the total maximum expenditure in all zones
# cannot exceed 3000
p6 <-
problem(sim_zones_pu_raster, sim_zones_features) %>%
add_max_features_objective(3000) %>%
add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# create weights that assign equal weighting for the representation
# of each feature in each zone except that it does not matter if
# feature 1 is represented in zone 1 and it really important
# that feature 3 is really in zone 1
w7 <- matrix(1, ncol = 3, nrow = 5)
w7[1, 1] <- 0
w7[3, 1] <- 100
# create problem with weights
p7 <- p6 %>% add_feature_weights(w7)
# solve problems
s6 <- solve(p6)
s7 <- solve(p7)
# convert solutions to category layers
c6 <- category_layer(s6)
c7 <- category_layer(s7)
# plot solutions
plot(c(c6, c7), main = c("equal weights", "manual weights"), axes = FALSE)
# create minimal problem to show the correct method for setting
# weights for problems with manual targets
p8 <-
problem(sim_pu_raster, sim_features) %>%
add_max_features_objective(budget = 3000) %>%
add_manual_targets(
data.frame(
feature = c("feature_1", "feature_4"),
type = "relative",
target = 0.1)
) %>%
add_feature_weights(matrix(c(1, 200), ncol = 1)) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# solve problem
s8 <- solve(p8)
# plot solution
plot(s8, main = "solution", axes = FALSE)
}
Run the code above in your browser using DataLab