#' @include internal.R ConservationProblem-class.R OptimizationProblem-class.R compile.R problem.R solve.R presolve_check.R
NULL
#' Evaluate solution importance using replacement cost scores
#'
#' Calculate importance scores for planning units selected in a solution
#' based on the replacement cost method (Cabeza and Moilanen 2006).
#'
#' @inheritParams eval_cost_summary
#'
#' @param rescale `logical` flag indicating if replacement cost
#' values -- excepting infinite (`Inf`) and zero values -- should be
#' rescaled to range between 0.01 and 1. Defaults to `TRUE`.
#'
#' @param run_checks `logical` flag indicating whether presolve checks
#' should be run prior solving the problem. These checks are performed using
#' the [presolve_check()] function. Defaults to `TRUE`.
#' Skipping these checks may reduce run time for large problems.
#'
#' @param force `logical` flag indicating if an attempt should be
#' made to solve the problem even if potential issues were detected during
#' the presolve checks. Defaults to `FALSE`.
#'
#' @param threads `integer` number of threads to use for the
#' optimization algorithm. Defaults to 1 such that only a single
#' thread is used.
#'
#' @param ... not used.
#'
#' @details
#' This function implements a modified version of the
#' replacement cost method (Cabeza and Moilanen 2006).
#' Specifically, the score for each planning unit is calculated
#' as the difference in the objective value of a solution when each planning
#' unit is locked out and the optimization processes rerun with all other
#' selected planning units locked in. In other words, the replacement cost
#' metric corresponds to change in solution quality incurred if a given
#' planning unit cannot be acquired when implementing the solution and the
#' next best planning unit (or set of planning units) will need to be
#' considered instead. Thus planning units with a higher score are more
#' important (and irreplaceable).
#' For example, when using the minimum set objective function
#' ([add_min_set_objective()]), the replacement cost scores
#' correspond to the additional costs needed to meet targets when each
#' planning unit is locked out. When using the maximum utility
#' objective function ([add_max_utility_objective()], the
#' replacement cost scores correspond to the reduction in the utility when
#' each planning unit is locked out. Infinite values mean that no feasible
#' solution exists when planning units are locked out---they are
#' absolutely essential for obtaining a solution (e.g., they contain rare
#' species that are not found in any other planning units or were locked in).
#' Zeros values mean that planning units can be swapped with other planning
#' units and this will have no effect on the performance of the solution at all
#' (e.g., because they were only selected due to spatial fragmentation
#' penalties).
#'
#' These calculations can take a long time to complete for large
#' or complex conservation planning problems. As such, we recommend using this
#' method for small or moderate-sized conservation planning problems
#' (e.g., < 30,000 planning units). To reduce run time, we
#' recommend calculating these scores without additional penalties (e.g.,
#' [add_boundary_penalties()]) or spatial constraints (e.g.,
#' [add_contiguity_constraints()]). To further reduce run time,
#' we recommend using proportion-type decisions when calculating the scores
#' (see below for an example).
#'
#' @inheritSection eval_cost_summary Solution format
#'
#' @return A `numeric`, `matrix`, `data.frame`,
#' [terra::rast()], or [sf::sf()] object
#' containing the importance scores for each planning
#' unit in the solution. Specifically, the returned object is in the
#' same format as the planning unit data in the argument to `x`.
#'
#' @seealso
#' See [importance] for an overview of all functions for evaluating
#' the importance of planning units selected in a solution.
#.
#' @family importances
#'
#' @examples
#' \dontrun{
#' # seed seed for reproducibility
#' set.seed(600)
#'
#' # load data
#' sim_pu_raster <- get_sim_pu_raster()
#' sim_pu_polygons <- get_sim_pu_polygons()
#' sim_features <- get_sim_features()
#' sim_zones_pu_raster <- get_sim_zones_pu_raster()
#' sim_zones_features <- get_sim_zones_features()
#'
#' # create minimal problem with binary decisions
#' p1 <-
#' problem(sim_pu_raster, sim_features) %>%
#' add_min_set_objective() %>%
#' add_relative_targets(0.1) %>%
#' add_binary_decisions() %>%
#' add_default_solver(gap = 0, verbose = FALSE)
#'
#' # solve problem
#' s1 <- solve(p1)
#'
#' # print solution
#' print(s1)
#'
#' # plot solution
#' plot(s1, main = "solution", axes = FALSE)
#'
#' # calculate importance scores
#' rc1 <- eval_replacement_importance(p1, s1)
#'
#' # print importance scores
#' print(rc1)
#'
#' # plot importance scores
#' plot(rc1, main = "replacement cost", axes = FALSE)
#'
#' # since replacement cost scores can take a long time to calculate with
#' # binary decisions, we can calculate them using proportion-type
#' # decision variables. Note we are still calculating the scores for our
#' # previous solution (s1), we are just using a different optimization
#' # problem when calculating the scores.
#' p2 <-
#' problem(sim_pu_raster, sim_features) %>%
#' add_min_set_objective() %>%
#' add_relative_targets(0.1) %>%
#' add_proportion_decisions() %>%
#' add_default_solver(gap = 0, verbose = FALSE)
#'
#' # calculate importance scores using proportion type decisions
#' rc2 <- eval_replacement_importance(p2, s1)
#'
#' # print importance scores based on proportion type decisions
#' print(rc2)
#'
#' # plot importance scores based on proportion type decisions
#' # we can see that the importance values in rc1 and rc2 are similar,
#' # and this confirms that the proportion type decisions are a good
#' # approximation
#' plot(rc2, main = "replacement cost", axes = FALSE)
#'
#' # create minimal problem with polygon planning units
#' p3 <-
#' problem(sim_pu_polygons, sim_features, cost_column = "cost") %>%
#' add_min_set_objective() %>%
#' add_relative_targets(0.05) %>%
#' add_binary_decisions() %>%
#' add_default_solver(gap = 0, verbose = FALSE)
#'
#' # solve problem
#' s3 <- solve(p3)
#'
#' # print solution
#' print(s3)
#'
#' # plot solution
#' plot(s3[, "solution_1"], main = "solution")
#'
#' # calculate importance scores
#' rc3 <- eval_rare_richness_importance(p3, s3[, "solution_1"])
#'
#' # plot importance scores
#' plot(rc3, main = "replacement cost")
#'
#' # build multi-zone conservation problem with raster data
#' p4 <-
#' problem(sim_zones_pu_raster, sim_zones_features) %>%
#' add_min_set_objective() %>%
#' add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>%
#' add_binary_decisions() %>%
#' add_default_solver(gap = 0, verbose = FALSE)
#'
#' # solve the problem
#' s4 <- solve(p4)
#' names(s4) <- paste0("zone ", seq_len(terra::nlyr(s4)))
#'
#' # print solution
#' print(s4)
#'
#' # plot solution
#' # each panel corresponds to a different zone, and data show the
#' # status of each planning unit in a given zone
#' plot(s4, axes = FALSE)
#'
#' # calculate importance scores
#' rc4 <- eval_replacement_importance(p4, s4)
#' names(rc4) <- paste0("zone ", seq_len(terra::nlyr(s4)))
#'
#' # plot importance
#' # each panel corresponds to a different zone, and data show the
#' # importance of each planning unit in a given zone
#' plot(rc4, axes = FALSE)
#' }
#'
#' @references
#' Cabeza M and Moilanen A (2006) Replacement cost: A practical measure of site
#' value for cost-effective reserve planning. *Biological Conservation*,
#' 132: 336--342.
#'
#' @aliases eval_replacement_importance,ConservationProblem,numeric-method eval_replacement_importance,ConservationProblem,matrix-method eval_replacement_importance,ConservationProblem,data.frame-method eval_replacement_importance,ConservationProblem,Spatial-method eval_replacement_importance,ConservationProblem,sf-method eval_replacement_importance,ConservationProblem,Raster-method eval_replacement_importance,ConservationProblem,SpatRaster-method
#'
#' @name eval_replacement_importance
#'
#' @rdname eval_replacement_importance
#'
#' @exportMethod eval_replacement_importance
methods::setGeneric("eval_replacement_importance",
function(x, solution, ...) {
assert_required(x)
assert_required(solution)
assert(
is_conservation_problem(x),
is_inherits(
solution,
c(
"numeric", "data.frame", "matrix", "sf", "SpatRaster",
"Spatial", "Raster"
)
)
)
standardGeneric("eval_replacement_importance")
}
)
#' @name eval_replacement_importance
#' @usage \S4method{eval_replacement_importance}{ConservationProblem,numeric}(x, solution, rescale, run_checks, force, threads, ...)
#' @rdname eval_replacement_importance
methods::setMethod("eval_replacement_importance",
methods::signature("ConservationProblem", "numeric"),
function(x, solution, rescale = TRUE, run_checks = TRUE, force = FALSE,
threads = 1L, ...) {
# assert valid arguments
assert(is.numeric(solution))
assert_dots_empty()
# extract planning unit solution status
status <- planning_unit_solution_status(x, solution)
# subset planning units with finite cost values
idx <- x$planning_unit_indices()
pos <- which(status > 1e-10)
# calculate replacement costs
v <- internal_eval_replacement_importance(
x, pos, rescale, run_checks, force, threads
)
# return replacement costs
out <- rep(NA_real_, x$number_of_total_units())
out[idx] <- 0
out[idx[pos]] <- c(v)
out
}
)
#' @name eval_replacement_importance
#' @usage \S4method{eval_replacement_importance}{ConservationProblem,matrix}(x, solution, rescale, run_checks, force, threads, ...)
#' @rdname eval_replacement_importance
methods::setMethod("eval_replacement_importance",
methods::signature("ConservationProblem", "matrix"),
function(x, solution, rescale = TRUE, run_checks = TRUE, force = FALSE,
threads = 1L, ...) {
# assert valid arguments
assert(
is.matrix(solution),
is.numeric(solution)
)
assert_dots_empty()
# extract data
status <- planning_unit_solution_status(x, solution)
# extract planning units in solution
pos <- which(status > 1e-10)
# calculate replacement costs
v <- internal_eval_replacement_importance(
x, pos, rescale, run_checks, force, threads
)
# initialize matrix
m_total <- matrix(
NA_real_,
nrow = x$number_of_total_units(),
ncol = x$number_of_zones()
)
m_pu <- matrix(
0,
nrow = x$number_of_planning_units(),
ncol = x$number_of_zones()
)
m_pu[pos] <- c(v)
m_pu[is.na(x$planning_unit_costs())] <- NA_real_
m_total[x$planning_unit_indices(), ] <- m_pu
# add column names to matrix
if (x$number_of_zones() > 1) {
colnames(m_total) <- paste0("rc_", x$zone_names())
} else {
colnames(m_total) <- "rc"
}
# return result
m_total
}
)
#' @name eval_replacement_importance
#' @usage \S4method{eval_replacement_importance}{ConservationProblem,data.frame}(x, solution, rescale, run_checks, force, threads, ...)
#' @rdname eval_replacement_importance
methods::setMethod("eval_replacement_importance",
methods::signature("ConservationProblem", "data.frame"),
function(x, solution, rescale = TRUE, run_checks = TRUE, force = FALSE,
threads = 1L, ...) {
# assert valid arguments
assert(is.data.frame(solution))
assert_dots_empty()
# extract data
status <- planning_unit_solution_status(x, solution)
# extract planning units in solution
pos <- which(status > 1e-10)
# calculate replacement costs
v <- internal_eval_replacement_importance(
x, pos, rescale, run_checks, force, threads
)
# initialize matrix
m_total <- matrix(
NA_real_,
nrow = x$number_of_total_units(),
ncol = x$number_of_zones()
)
m_pu <- matrix(
0,
nrow = x$number_of_planning_units(),
ncol = x$number_of_zones()
)
m_pu[pos] <- c(v)
m_pu[is.na(x$planning_unit_costs())] <- NA_real_
m_total[x$planning_unit_indices(), ] <- m_pu
# add column names to matrix
if (x$number_of_zones() > 1) {
colnames(m_total) <- paste0("rc_", x$zone_names())
} else {
colnames(m_total) <- "rc"
}
# return result
tibble::as_tibble(m_total)
}
)
#' @name eval_replacement_importance
#' @usage \S4method{eval_replacement_importance}{ConservationProblem,Spatial}(x, solution, rescale, run_checks, force, threads, ...)
#' @rdname eval_replacement_importance
methods::setMethod("eval_replacement_importance",
methods::signature("ConservationProblem", "Spatial"),
function(x, solution, rescale = TRUE, run_checks = TRUE, force = FALSE,
threads = 1L, ...) {
# assert valid arguments
assert(
is_inherits(
solution,
c(
"SpatialPointsDataFrame", "SpatialLinesDataFrame",
"SpatialPolygonsDataFrame"
)
)
)
assert_dots_empty()
# extract data
status <- planning_unit_solution_status(x, solution)
# extract planning units in solution
pos <- which(status > 1e-10)
# calculate replacement costs
v <- internal_eval_replacement_importance(
x, pos, rescale, run_checks, force, threads
)
# initialize matrix
m_total <- matrix(
NA_real_,
nrow = x$number_of_total_units(),
ncol = x$number_of_zones()
)
m_pu <- matrix(
0,
nrow = x$number_of_planning_units(),
ncol = x$number_of_zones()
)
m_pu[pos] <- c(v)
m_pu[is.na(x$planning_unit_costs())] <- NA_real_
m_total[x$planning_unit_indices(), ] <- m_pu
# add column names to matrix
if (x$number_of_zones() > 1) {
colnames(m_total) <- paste0("rc_", x$zone_names())
} else {
colnames(m_total) <- "rc"
}
# return result
out <- as.data.frame(m_total)
rownames(out) <- rownames(solution@data)
solution@data <- out
solution
}
)
#' @name eval_replacement_importance
#' @usage \S4method{eval_replacement_importance}{ConservationProblem,sf}(x, solution, rescale, run_checks, force, threads, ...)
#' @rdname eval_replacement_importance
methods::setMethod("eval_replacement_importance",
methods::signature("ConservationProblem", "sf"),
function(x, solution, rescale = TRUE, run_checks = TRUE, force = FALSE,
threads = 1L, ...) {
# assert valid arguments
assert(inherits(solution, "sf"))
assert_dots_empty()
# extract data
status <- planning_unit_solution_status(x, solution)
# extract planning units in solution
pos <- which(status > 1e-10)
# calculate replacement costs
v <- internal_eval_replacement_importance(
x, pos, rescale, run_checks, force, threads
)
# initialize matrix
m_total <- matrix(
NA_real_,
nrow = x$number_of_total_units(),
ncol = x$number_of_zones()
)
m_pu <- matrix(
0,
nrow = x$number_of_planning_units(),
ncol = x$number_of_zones()
)
m_pu[pos] <- c(v)
m_pu[is.na(x$planning_unit_costs())] <- NA_real_
m_total[x$planning_unit_indices(), ] <- m_pu
# add column names to matrix
if (x$number_of_zones() > 1) {
colnames(m_total) <- paste0("rc_", x$zone_names())
} else {
colnames(m_total) <- "rc"
}
# return result
out <- tibble::as_tibble(as.data.frame(m_total))
out$geometry <- sf::st_geometry(x$data$cost)
sf::st_sf(out, crs = sf::st_crs(x$data$cost))
}
)
#' @name eval_replacement_importance
#' @usage \S4method{eval_replacement_importance}{ConservationProblem,Raster}(x, solution, rescale, run_checks, force, threads, ...)
#' @rdname eval_replacement_importance
methods::setMethod("eval_replacement_importance",
methods::signature("ConservationProblem", "Raster"),
function(x, solution, rescale = TRUE, run_checks = TRUE, force = FALSE,
threads = 1L, ...) {
assert(inherits(solution, "Raster"))
assert_dots_empty()
# extract data
status <- planning_unit_solution_status(x, solution)
# extract planning units in solution
pos <- which(status > 1e-10)
# calculate replacement costs
v <- internal_eval_replacement_importance(
x, pos, rescale, run_checks, force, threads
)
# initialize matrix
m_pu <- matrix(
0,
nrow = x$number_of_planning_units(),
ncol = x$number_of_zones()
)
m_pu[pos] <- c(v)
m_pu[is.na(x$planning_unit_costs())] <- NA_real_
# add column names to matrix
if (x$number_of_zones() > 1) {
colnames(m_pu) <- paste0("rc_", x$zone_names())
} else {
colnames(m_pu) <- "rc"
}
# return result
out <- raster::as.list(solution)
for (i in seq_along(out)) {
out[[i]][x$planning_unit_indices()] <- m_pu[, i]
out[[i]][raster::Which(is.na(solution[[i]]), cells = TRUE)] <- NA_real_
}
if (length(out) > 1) {
out <- raster::stack(out)
} else {
out <- out[[1]]
}
names(out) <- colnames(m_pu)
out
}
)
#' @name eval_replacement_importance
#' @usage \S4method{eval_replacement_importance}{ConservationProblem,SpatRaster}(x, solution, rescale, run_checks, force, threads, ...)
#' @rdname eval_replacement_importance
methods::setMethod("eval_replacement_importance",
methods::signature("ConservationProblem", "SpatRaster"),
function(x, solution, rescale = TRUE, run_checks = TRUE, force = FALSE,
threads = 1L, ...) {
assert(inherits(solution, "SpatRaster"))
assert_dots_empty()
# extract data
status <- planning_unit_solution_status(x, solution)
# extract planning units in solution
pos <- which(status > 1e-10)
# calculate replacement costs
v <- internal_eval_replacement_importance(
x, pos, rescale, run_checks, force, threads
)
# initialize matrix
m_pu <- matrix(
0,
nrow = x$number_of_planning_units(),
ncol = x$number_of_zones()
)
m_pu[pos] <- c(v)
m_pu[is.na(x$planning_unit_costs())] <- NA_real_
# add column names to matrix
if (x$number_of_zones() > 1) {
colnames(m_pu) <- paste0("rc_", x$zone_names())
} else {
colnames(m_pu) <- "rc"
}
# return result
out <- terra::as.list(solution)
for (i in seq_along(out)) {
out[[i]][x$planning_unit_indices()] <- m_pu[, i]
out[[i]][is.na(solution[[i]])] <- NA_real_
}
out <- terra::rast(out)
names(out) <- colnames(m_pu)
out
}
)
internal_eval_replacement_importance <- function(x, indices, rescale,
run_checks, force,
threads = 1,
call = fn_caller_env()) {
# assert valid arguments
assert(
is_conservation_problem(x),
is.integer(indices),
length(indices) > 0,
assertthat::is.flag(rescale),
assertthat::is.flag(run_checks),
assertthat::is.flag(force),
is_thread_count(threads),
call = call
)
# assign default solver
if (inherits(x$solver, "Waiver"))
x <- add_default_solver(x)
# overwrite portfolio
x <- add_shuffle_portfolio(x, 1)
# compile problem
opt <- compile.ConservationProblem(x)
# run presolve check to try to identify potential problems
if (run_checks) {
## run checks
presolve_res <- internal_presolve_check(opt)
## prepare message
msg <- presolve_res$msg
if (!isTRUE(force)) {
msg <- c(
msg,
"i" = paste(
"To ignore checks and attempt optimization anyway,",
"use {.code solve(force = TRUE)}."
)
)
}
## determine if error or warning should be thrown
if (!isTRUE(force)) {
f <- assert
} else {
f <- verify
}
## throw error or warning if checks failed
f(isTRUE(presolve_res$pass), call = parent.frame(), msg = msg)
}
# solve problem
x$solver$calculate(opt)
modelsense <- opt$modelsense()
rm(opt)
# lock in decisions in the solution
x$solver$set_variable_lb(indices, rep.int(1, length(indices)))
x$solver$set_variable_ub(indices, rep.int(1, length(indices)))
# generate objective value for solution if unknown
solution_obj <- x$solver$run()
if (is.null(solution_obj) || is.null(solution_obj$x))
cli::cli_abort(
"{.arg solution} is not feasible solution for the problem {.arg x}",
call = call
)
solution_obj <- solution_obj[[2]]
# define function for processing
calculate_alt_solution_obj <- function(y) {
vapply(
indices[y],
FUN.VALUE = numeric(1),
function(i) {
# lock out i'th selected planning unit in solution
x$solver$set_variable_lb(i, 0)
x$solver$set_variable_ub(i, 0)
# solve problem
sol <- x$solver$run()
# check that solution is valid
if (is.null(sol) || is.null(sol$x)) {
out <- Inf
} else {
out <- sol[[2]]
}
# reset upper bound
x$solver$set_variable_lb(i, 1)
x$solver$set_variable_ub(i, 1)
# return result
out
}
)
}
# iterate over decision variables in solution and store new objective values
if (isTRUE(threads > 1L)) {
## initialize cluster
cl <- parallel::makeCluster(threads, "PSOCK")
## prepare cluster clean up
on.exit(try(cl <- parallel::stopCluster(cl), silent = TRUE))
## move data to workers
parallel::clusterExport(
cl,
c("indices", "x", "calculate_alt_solution_obj"),
envir = environment()
)
## main processing
alt_solution_obj <- parallel::parLapply(
cl,
parallel::splitIndices(length(indices), threads),
calculate_alt_solution_obj
)
} else {
## main processing
alt_solution_obj <- lapply(
parallel::splitIndices(length(indices), threads),
calculate_alt_solution_obj
)
}
# reformat output
alt_solution_obj <- unlist(
alt_solution_obj, recursive = TRUE, use.names = FALSE
)
# calculate replacement costs
out <- alt_solution_obj - solution_obj
if (identical(modelsense, "max")) # rescale values if maximization problem
out <- out * -1
# rescale values if specified
if (rescale) {
rescale_ind <- is.finite(out) & (abs(out) > 1e-10)
out[rescale_ind] <- rescale(out[rescale_ind], to = c(0.01, 1))
}
# return result
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.