R/solve.R

In prioritizr/prioritizr: Systematic Conservation Prioritization in R

#' @include internal.R ConservationProblem-class.R OptimizationProblem-class.R compile.R presolve_check.R
NULL

#' Solve
#'
#' Solve a conservation planning problem.
#'
#' @param a [problem()] object.
#'
#' @param b missing.
#'
#' @param ... arguments passed to [compile()].
#'
#' @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 to should be
#'   made to solve the problem even if potential issues were detected during
#'   the presolve checks. Defaults to `FALSE`.
#'
#' @details
#' After formulating a conservation planning [problem()],
#' it can be solved using an exact algorithm solver (see [solvers]
#' for available solvers). If no solver has been explicitly specified,
#' then the best available exact algorithm solver will be used by default
#' (see [add_default_solver()]). Although these exact algorithm
#' solvers will often display a lot of information that isn't really that
#' helpful (e.g., nodes, cutting planes), they do display information
#' about the progress they are making on solving the problem (e.g., the
#' performance of the best solution found at a given point in time). If
#' potential issues were detected during the
#' presolve checks (see [presolve_check()])
#' and the problem is being forcibly solved (i.e., with `force = TRUE`),
#' then it is also worth checking for any warnings displayed by the solver
#' to see if these potential issues are actually causing issues
#' (e.g., *Gurobi* can display warnings that include
#' `"Warning: Model contains large matrix coefficient range"` and
#' `"Warning: Model contains large rhs"`).
#'
#' @section Output format:
#' This function will output solutions in a similar format to the
#' planning units associated with `a`. Specifically, it will return
#' solutions based on the following types of planning units.
#'
#'   \describe{
#'
#'   \item{`a` has `numeric` planning units}{The solution will be
#'    returned as a `numeric` vector. Here, each element in the vector
#'     corresponds to a different planning unit.
#'     Note that if a portfolio is used to generate multiple solutions,
#'     then a `list` of such `numeric` vectors will be returned.}
#'
#'   \item{`a` has `matrix` planning units}{The solution will be
#'     returned as a `matrix` object.
#'     Here, rows correspond to different planning units,
#'     and columns correspond to different  management zones.
#'     Note that if a portfolio is used to generate multiple solutions,
#'     then a `list` of such `matrix` objects will be returned.}
#'
#'   \item{`a` has [terra::rast()] planning units}{The solution
#'     will be returned as a [terra::rast()] object.
#'     If the argument to `x` contains multiple zones, then the object
#'     will have a different layer for each management zone.
#'     Note that if a portfolio is used to generate multiple solutions,
#'     then a `list` of [terra::rast()] objects will be returned.}
#'
#'   \item{`a` has [sf::sf()], or `data.frame` planning units}{
#'     The solution will be returned in the same data format as the planning
#'     units.
#'     Here, each row corresponds to a different planning unit,
#'     and columns contain solutions.
#'     If the argument to `a` contains a single zone, then the solution object
#'     will contain columns named by solution.
#'     Specifically, the column names containing the solution values
#'     be will named as `"solution_XXX"` where `"XXX"` corresponds to a solution
#'     identifier (e.g., `"solution_1"`).
#'     If the argument to `a` contains multiple zones, then the columns
#'     containing solutions will be named as `"solution_XXX_YYY"` where
#'     `"XXX"` corresponds to the solution identifier and `"YYY"` is the name
#'     of the management zone (e.g., `"solution_1_zone1"`).}
#'
#'   }
#'
#' @return
#' A `numeric`, `matrix`, `data.frame`, [sf::st_sf()], or
#' [terra::rast()] object containing the solution to the problem.
#' Additionally, the returned object will have the following additional
#' attributes: `"objective"` containing the solution's objective,
#' `"runtime"` denoting the number of seconds that elapsed while solving
#' the problem, and `"status"` describing the status of the solution
#' (e.g., `"OPTIMAL"` indicates that the optimal solution was found).
#'
#' @seealso
#' See [problem()] to create conservation planning problems, and
#' [presolve_check()] to check problems for potential issues.
#' Also, see the [category_layer()] and [category_vector()] function to
#' reformat solutions that contain multiple zones.
#'
#' @examples
#' \dontrun{
#' # set seed for reproducibility
#' set.seed(500)
#'
#' # 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_pu_polygons <- get_sim_zones_pu_polygons()
#' sim_zones_features <- get_sim_zones_features()
#'
#' # build minimal conservation problem with raster data
#' p1 <-
#'   problem(sim_pu_raster, sim_features) %>%
#'   add_min_set_objective() %>%
#'   add_relative_targets(0.1) %>%
#'   add_binary_decisions() %>%
#'   add_default_solver(verbose = FALSE)
#'
#' # solve the problem
#' s1 <- solve(p1)
#'
#' # print solution
#' print(s1)
#'
#' # print attributes describing the optimization process and the solution
#' print(attr(s1, "objective"))
#' print(attr(s1, "runtime"))
#' print(attr(s1, "status"))
#'
#' # calculate feature representation in the solution
#' r1 <- eval_feature_representation_summary(p1, s1)
#' print(r1)
#'
#' # plot solution
#' plot(s1, main = "solution", axes = FALSE)
#'
#' # build minimal conservation problem with polygon data
#' p2 <-
#'   problem(sim_pu_polygons, sim_features, cost_column = "cost") %>%
#'   add_min_set_objective() %>%
#'   add_relative_targets(0.1) %>%
#'   add_binary_decisions() %>%
#'   add_default_solver(verbose = FALSE)
#'
#' # solve the problem
#' s2 <- solve(p2)
#'
#' # print solution
#' print(s2)
#'
#' # calculate feature representation in the solution
#' r2 <- eval_feature_representation_summary(p2, s2[, "solution_1"])
#' print(r2)
#'
#' # plot solution
#' plot(s2[, "solution_1"], main = "solution", axes = FALSE)
#'
#' # build multi-zone conservation problem with raster data
#' p3 <-
#'   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(verbose = FALSE)
#'
#' # solve the problem
#' s3 <- solve(p3)
#'
#' # print solution
#' print(s3)
#'
#' # calculate feature representation in the solution
#' r3 <- eval_feature_representation_summary(p3, s3)
#' print(r3)
#'
#' # plot solution
#' plot(category_layer(s3), main = "solution", axes = FALSE)
#'
#' # build multi-zone conservation problem with polygon data
#' p4 <-
#'   problem(
#'     sim_zones_pu_polygons, sim_zones_features,
#'     cost_column = c("cost_1", "cost_2", "cost_3")
#'   ) %>%
#'   add_min_set_objective() %>%
#'   add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>%
#'   add_binary_decisions() %>%
#'   add_default_solver(verbose = FALSE)
#'
#' # solve the problem
#' s4 <- solve(p4)
#'
#' # print solution
#' print(s4)
#'
#' # calculate feature representation in the solution
#' r4 <- eval_feature_representation_summary(
#'   p4, s4[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")]
#' )
#' print(r4)
#'
#' # create new column representing the zone id that each planning unit
#' # was allocated to in the solution
#' s4$solution <- category_vector(
#'   s4[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")]
#' )
#' s4$solution <- factor(s4$solution)
#'
#' # plot solution
#' plot(s4[, "solution"])
#' }
#' @name solve
NULL

#' @rdname solve
#' @method solve ConservationProblem
#' @export solve.ConservationProblem
#' @export
solve.ConservationProblem <- function(a, b, ...,
                                      run_checks = TRUE, force = FALSE) {
  # assert arguments are valid
  assert_required(a)
  assert(
    assertthat::is.flag(run_checks),
    assertthat::noNA(run_checks),
    assertthat::is.flag(force),
    assertthat::noNA(force)
  )
  if (!rlang::is_missing(b)) {
    cli::cli_abort("{.arg b} must not be specified.")
  }
  # compile optimization problem
  opt <- compile.ConservationProblem(a, ...)
  # 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), msg = msg)
  }
  # solve problem
  sol <- a$portfolio$run(opt, a$solver)
  # prepare message
  msg <- c(
    "Can't find a solution!",
    "i" = paste(
      "This is because it is impossible to meet the",
      "targets, budgets, or constraints."
    )
  )
  if (isTRUE(a$solver$data$time_limit < 1e5)) {
    msg <- c(
      msg,
      "i" = paste(
        "It could also be because the {.arg time_limit}",
        "is too low."
      )
    )
  }
  # check that solution is valid
  assert(
    !is.null(sol) && !is.null(sol[[1]]$x),
    msg = msg
  )
  # check that desired number of solutions were found
  portfolio_number_solutions <- a$portfolio$get_data("number_solutions")
  if (!is.Waiver(portfolio_number_solutions)) {
    if (length(sol) != portfolio_number_solutions) {
      cli_warning(
        paste(
          "Portfolio could only find",
          "{.val {length(sol)}} out of",
          "{.val {portfolio_number_solutions}}",
          "solution{?s}."
        )
      )
    }
  }
  ## format solutions
  # format solutions into planning unit by zones matrix
  na_pos <- which(is.na(a$planning_unit_costs()), arr.ind = TRUE)
  sol_status <- lapply(sol, function(x) {
    m <- matrix(
      x[[1]][seq_len(a$number_of_planning_units() * a$number_of_zones())],
      nrow = a$number_of_planning_units(),
      ncol = a$number_of_zones()
    )
    m[na_pos] <- NA_real_
    m
  })
  # create solution data
  pu <- a$data$cost
  if (inherits(pu, c("SpatRaster", "Raster"))) {
    # SpatRaster or Raster planning units
    pos <- a$planning_unit_indices()
    pu <- terra::rast(pu)
    pu <- suppressWarnings(terra::setValues(pu[[1]], NA))
    ret <- lapply(sol_status, function(s) {
      ret <- lapply(seq_len(ncol(s)), function(z) {
        pu[pos] <- s[, z]
        pu
      })
      ret <- terra::rast(ret)
      names(ret) <- a$zone_names()
      ret
    })
    # convert to RasterStack or RasterLayer if needed
    if (inherits(a$data$cost, c("RasterStack", "RasterBrick"))) {
      ret <- lapply(ret, raster::stack)
    } else if (inherits(a$data$cost, "RasterLayer")) {
      ret <- lapply(ret, raster::raster)
    }
    names(ret) <- paste0("solution_", seq_along(sol))
  } else if (inherits(pu, c("data.frame", "Spatial", "sf"))) {
    # Spatial* or data.frame planning units
    sol_status <- do.call(cbind, sol_status)
    if (a$number_of_zones() == 1) {
      colnames(sol_status) <- paste0("solution_", seq_along(sol))
    } else {
      colnames(sol_status) <- paste0(
        "solution_",
        rep(seq_along(sol), each = a$number_of_zones()),
        "_",
        rep(a$zone_names(), length(sol))
      )
    }
    # add in NA values for planning units that contained NA values in
    # all zones that were discarded from the mathematical formulation
    # to reduce overheads
    pos <- a$planning_unit_indices()
    if (!identical(pos, seq_len(a$number_of_total_units()))) {
      sol_status2 <- matrix(
        NA_real_,
        nrow = a$number_of_total_units(),
        ncol = ncol(sol_status)
      )
      sol_status2[pos, ] <- sol_status
      dimnames(sol_status2) <- dimnames(sol_status)
    } else {
      sol_status2 <- sol_status
    }
    # cbind solutions to planning unit data
    sol_status2 <- as.data.frame(sol_status2)
    if (inherits(pu, "Spatial")) {
      ret <- pu
      ret@data <- cbind(ret@data, sol_status2)
    } else if (inherits(pu, "sf")) {
      ret <- tibble::as_tibble(data.frame(pu, sol_status2))
      sf_col <- attr(pu, "sf_column")
      ret <- ret[, c(setdiff(names(ret), sf_col), sf_col), drop = FALSE]
      ret <- ret <- sf::st_sf(ret)
    } else {
      ret <- cbind(pu, sol_status2)
      if (inherits(pu, "tbl_df")) {
        pu <- tibble::as_tibble(pu)
      }
    }
  } else if (is.matrix(pu)) {
    # matrix planning units
    # add in NA values for planning units that contained NA values in
    # all zones that were discarded from the mathematical formulation
    # to reduce overheads
    pos <- a$planning_unit_indices()
    pu[] <- NA
    colnames(pu) <- a$zone_names()
    ret <- lapply(sol_status, function(s) {
      pu[pos, ] <- s
      pu
    })
    names(ret) <- paste0("solution_", seq_along(sol))
  } else {
    # nocov start
    cli::cli_abort(
      "Planning unit data is of an unrecognized class.",
      .internal = TRUE
    )
    # nocov end
  }
  # if ret is a list of matrices with a single column then convert to numeric
  if (is.matrix(ret[[1]]) && ncol(ret[[1]]) == 1) {
    ret <- lapply(ret, as.numeric)
  }
  # if ret is a list with a single element then extract the element
  if (length(ret) == 1) {
    ret <- ret[[1]]
  }
  # add attributes
  attr(ret, "objective") <- stats::setNames(
    vapply(sol, `[[`, numeric(1), 2), paste0("solution_", seq_along(sol))
  )
  attr(ret, "status") <- stats::setNames(
    vapply(sol, `[[`, character(1), 3), paste0("solution_", seq_along(sol))
  )
  attr(ret, "runtime") <- stats::setNames(
    vapply(sol, `[[`, numeric(1), 4), paste0("solution_", seq_along(sol))
  )
  # return object
  ret
}