diff --git a/.github/workflows/MacOS.yaml b/.github/workflows/MacOS.yaml index eeb6db9..f44ef01 100644 --- a/.github/workflows/MacOS.yaml +++ b/.github/workflows/MacOS.yaml @@ -46,6 +46,7 @@ jobs: any::rcmdcheck any::remotes any::lpsymphony + any::slam needs: check - uses: r-lib/actions/check-r-package@v2 diff --git a/.github/workflows/Ubuntu.yaml b/.github/workflows/Ubuntu.yaml index 03878a2..93f329b 100644 --- a/.github/workflows/Ubuntu.yaml +++ b/.github/workflows/Ubuntu.yaml @@ -45,7 +45,7 @@ jobs: sudo apt-get -y install \ libcurl4-gnutls-dev coinor-libsymphony-dev \ libudunits2-dev libgdal-dev libgeos-dev libproj-dev libglpk-dev \ - libgmp3-dev libmpfr-dev chromium-browser + libgmp3-dev libmpfr-dev - uses: r-lib/actions/setup-r-dependencies@v2 with: @@ -53,6 +53,7 @@ jobs: any::rcmdcheck any::remotes any::lpsymphony + any::slam needs: check - uses: r-lib/actions/check-r-package@v2 diff --git a/.github/workflows/Windows.yaml b/.github/workflows/Windows.yaml index e46efd1..0c14c96 100644 --- a/.github/workflows/Windows.yaml +++ b/.github/workflows/Windows.yaml @@ -35,8 +35,6 @@ jobs: - uses: r-lib/actions/setup-pandoc@v2 - - uses: browser-actions/setup-chrome@v1 - - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} @@ -50,6 +48,7 @@ jobs: any::rcmdcheck any::remotes any::lpsymphony + any::slam needs: check - uses: r-lib/actions/check-r-package@v2 diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index fd20dfa..6df7cd0 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -30,7 +30,7 @@ jobs: sudo apt-get -y install \ libcurl4-gnutls-dev coinor-libsymphony-dev \ libudunits2-dev libgdal-dev libgeos-dev libproj-dev libglpk-dev \ - libgmp3-dev libmpfr-dev chromium-browser + libgmp3-dev libmpfr-dev - name: Install R package dependencies uses: r-lib/actions/setup-r-dependencies@v2 @@ -40,6 +40,7 @@ jobs: any::xml2 any::remotes any::lpsymphony + any::slam needs: coverage - name: Test coverage diff --git a/DESCRIPTION b/DESCRIPTION index 6f9060e..8e9f343 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,7 +41,9 @@ Imports: prioritizr (>= 7.0.0), rlang, stats, - magrittr + magrittr, + igraph, + Matrix Suggests: dplyr, ggplot2, @@ -52,7 +54,8 @@ Suggests: patchwork, rmarkdown, stars, - stringr + stringr, + parallelly VignetteBuilder: knitr URL: https://github.com/SpatialPlanning/minpatch BugReports: https://github.com/SpatialPlanning/minpatch/issues diff --git a/R/data_structures.R b/R/data_structures.R index ca46822..d622468 100644 --- a/R/data_structures.R +++ b/R/data_structures.R @@ -159,28 +159,47 @@ initialize_minpatch_data <- function(solution, planning_units, targets, costs, #' Create boundary matrix from planning units #' -#' Creates a matrix of shared boundary lengths between adjacent planning units +#' Creates a sparse matrix of shared boundary lengths between adjacent planning units. +#' Returns a Matrix::sparseMatrix for efficient storage and operations. +#' This optimized version supports parallel processing via the parallelly package. +#' When n_cores = 1, runs sequentially with no parallel overhead. #' #' @param planning_units sf object with planning unit geometries +#' @param verbose Logical, whether to print progress +#' @param n_cores Integer, number of cores to use. If NULL, uses availableCores(omit=1). +#' Set to 1 for sequential processing. #' -#' @return Named list where each element contains neighbors and shared boundary lengths +#' @return Matrix::dgCMatrix sparse matrix where [i,j] is the shared boundary length #' @keywords internal -create_boundary_matrix <- function(planning_units, verbose = TRUE) { +create_boundary_matrix <- function(planning_units, verbose = TRUE, n_cores = NULL) { n_units <- nrow(planning_units) - boundary_matrix <- vector("list", n_units) - names(boundary_matrix) <- as.character(seq_len(n_units)) - # Initialize empty lists for each planning unit - for (i in seq_len(n_units)) { - boundary_matrix[[i]] <- list() + # Determine number of cores + if (is.null(n_cores)) { + if (requireNamespace("parallelly", quietly = TRUE)) { + n_cores <- parallelly::availableCores(omit = 2) + } else { + n_cores <- 1 + } + } + # Only use parallel for larger datasets (overhead not worth it for small ones) + if (n_units < 500) { + n_cores <- 1 + } else { + n_cores <- min(n_cores, n_units) } - # Find adjacent planning units and calculate shared boundary lengths - # This is computationally intensive, so we'll use sf::st_touches for adjacency - # and sf::st_intersection for boundary lengths + # Final safety check: ensure n_cores is always between 1 and n_units + n_cores <- max(1, min(n_cores, n_units)) - if (verbose) cat("Calculating boundary matrix (this may take a while)...\n") + if (verbose) { + if (n_cores > 1) { + cat("Calculating boundary matrix using", n_cores, "cores...\n") + } else { + cat("Calculating boundary matrix (optimized version)...\n") + } + } # Check for invalid geometries and repair if needed if (any(!sf::st_is_valid(planning_units))) { @@ -188,49 +207,113 @@ create_boundary_matrix <- function(planning_units, verbose = TRUE) { planning_units <- sf::st_make_valid(planning_units) } - # Get adjacency matrix using a more robust method - # sf::st_touches() can be unreliable due to precision issues - # Use st_intersects() with boundaries instead + # Pre-compute all boundaries once (major optimization) boundaries <- sf::st_boundary(planning_units) - touches <- sf::st_intersects(boundaries, boundaries, sparse = FALSE) - - # Remove self-intersections (diagonal) - diag(touches) <- FALSE + # Pre-compute all perimeters once for diagonal + perimeters <- as.numeric(sf::st_length(boundaries)) + + # Get sparse adjacency list (much more efficient than dense matrix) + touches_sparse <- sf::st_intersects(boundaries, boundaries) + + # Split work into chunks - handle edge cases properly + if (n_cores == 1) { + # Single core: all units in one chunk + chunks <- list(seq_len(n_units)) + } else { + # Multiple cores: split evenly + # Ensure we don't try to create more chunks than units + actual_cores <- min(n_cores, n_units) + if (actual_cores >= n_units) { + # If cores >= units, each unit gets its own chunk + chunks <- as.list(seq_len(n_units)) + } else { + # Normal case: split into chunks + chunks <- split(seq_len(n_units), cut(seq_len(n_units), actual_cores, labels = FALSE)) + } + } - # Calculate shared boundaries - for (i in seq_len(n_units)) { - for (j in seq_len(n_units)) { - if (i != j && touches[i, j]) { - # Calculate shared boundary length (suppress sf warnings) - intersection <- suppressWarnings(sf::st_intersection( - sf::st_boundary(planning_units[i, ]), - sf::st_boundary(planning_units[j, ]) - )) - - if (nrow(intersection) > 0) { - shared_length <- sum(as.numeric(sf::st_length(intersection))) - # Use tolerance for very small shared lengths (floating-point precision issues) - if (shared_length > 1e-10) { - boundary_matrix[[i]][[as.character(j)]] <- shared_length - } else if (shared_length > 0) { - # For very small but non-zero lengths, use a minimal positive value - boundary_matrix[[i]][[as.character(j)]] <- 1e-6 + # Function to process a chunk of units + process_chunk <- function(unit_indices) { + local_i <- integer() + local_j <- integer() + local_lengths <- numeric() + + for (i in unit_indices) { + neighbors <- touches_sparse[[i]] + neighbors <- neighbors[neighbors != i] + + if (length(neighbors) > 0) { + for (j in neighbors) { + if (i < j) { # Only process each pair once + intersection <- suppressWarnings(sf::st_intersection( + boundaries[i, ], + boundaries[j, ] + )) + + if (nrow(intersection) > 0) { + shared_length <- sum(as.numeric(sf::st_length(intersection))) + if (shared_length > 1e-10) { + local_i <- c(local_i, i, j) + local_j <- c(local_j, j, i) + local_lengths <- c(local_lengths, shared_length, shared_length) + } else if (shared_length > 0) { + local_i <- c(local_i, i, j) + local_j <- c(local_j, j, i) + local_lengths <- c(local_lengths, 1e-6, 1e-6) + } + } } } } } - # Add self-boundary (external edge) - approximate as perimeter - perimeter <- as.numeric(sf::st_length(sf::st_boundary(planning_units[i, ]))) - boundary_matrix[[i]][[as.character(i)]] <- perimeter + list(i = local_i, j = local_j, x = local_lengths) + } - if (verbose && i %% 100 == 0) { - cat("Processed", i, "of", n_units, "planning units\n") - } + # Process chunks (parallel if n_cores > 1, sequential if n_cores = 1) + if (n_cores > 1 && requireNamespace("parallelly", quietly = TRUE)) { + # Parallel processing + cl <- parallelly::makeClusterPSOCK(n_cores, autoStop = TRUE, verbose = FALSE) + on.exit(parallel::stopCluster(cl), add = TRUE) + + parallel::clusterExport(cl, c("boundaries", "touches_sparse"), + envir = environment()) + parallel::clusterEvalQ(cl, library(sf)) + + if (verbose) cat("Processing chunks in parallel...\n") + results <- parallel::parLapply(cl, chunks, process_chunk) + } else { + # Sequential processing + results <- lapply(chunks, function(chunk) { + result <- process_chunk(chunk) + if (verbose && max(chunk) %% 100 == 0) { + cat("Processed", max(chunk), "of", n_units, "planning units\n") + } + result + }) } - return(boundary_matrix) + # Combine results + if (verbose && n_cores > 1) cat("Combining results...\n") + i_indices <- unlist(lapply(results, function(r) r$i)) + j_indices <- unlist(lapply(results, function(r) r$j)) + boundary_lengths <- unlist(lapply(results, function(r) r$x)) + + # Add perimeters on diagonal + i_indices <- c(i_indices, seq_len(n_units)) + j_indices <- c(j_indices, seq_len(n_units)) + boundary_lengths <- c(boundary_lengths, perimeters) + + # Create sparse matrix + Matrix::sparseMatrix( + i = i_indices, + j = j_indices, + x = boundary_lengths, + dims = c(n_units, n_units), + dimnames = list(as.character(seq_len(n_units)), + as.character(seq_len(n_units))) + ) } #' Create abundance matrix from planning units @@ -282,7 +365,8 @@ create_abundance_matrix <- function(planning_units, prioritizr_problem) { #' Create patch radius dictionary #' -#' For each planning unit, find all units within the specified patch radius +#' For each planning unit, find all units within the specified patch radius. +#' Optimized version computes full distance matrix once instead of n times. #' #' @param planning_units sf object with planning unit geometries #' @param patch_radius radius for patch creation @@ -299,16 +383,21 @@ create_patch_radius_dict <- function(planning_units, patch_radius, verbose = TRU centroids <- sf::st_centroid(planning_units %>% dplyr::select("geometry")) - if (verbose) cat("Creating patch radius dictionary...\n") + if (verbose) cat("Creating patch radius dictionary (optimized)...\n") - for (i in seq_len(n_units)) { - # Find all units within patch_radius of unit i - distances <- sf::st_distance(centroids[i, ], centroids) - # Convert both to numeric to avoid units mismatch - distances_numeric <- as.numeric(distances) - patch_radius_numeric <- as.numeric(patch_radius) - within_radius <- which(distances_numeric <= patch_radius_numeric & seq_len(n_units) != i) + # OPTIMIZATION: Compute full distance matrix ONCE instead of n times + # This changes from O(n²) distance calculations to O(n²/2) calculations + dist_matrix <- sf::st_distance(centroids, centroids) + dist_matrix_numeric <- as.numeric(dist_matrix) + patch_radius_numeric <- as.numeric(patch_radius) + # Create matrix of dimensions n x n + dist_mat <- matrix(dist_matrix_numeric, nrow = n_units, ncol = n_units) + + # For each unit, find neighbors within radius + for (i in seq_len(n_units)) { + # Use vectorized comparison on pre-computed distances + within_radius <- which(dist_mat[i, ] <= patch_radius_numeric & seq_len(n_units) != i) patch_radius_dict[[i]] <- as.character(within_radius) if (verbose && i %% 100 == 0) { diff --git a/R/minpatch.R b/R/minpatch.R index 0c73372..3a541cc 100644 --- a/R/minpatch.R +++ b/R/minpatch.R @@ -99,6 +99,8 @@ run_minpatch <- function(prioritizr_problem, solution_column = "solution_1", verbose = TRUE) { + # Stage 0: Checks and data preparation ----- + # Check if prioritizr is available if (!requireNamespace("prioritizr", quietly = TRUE)) { stop("prioritizr package is required for this function") @@ -206,7 +208,9 @@ run_minpatch <- function(prioritizr_problem, minpatch_data <- calculate_patch_stats(minpatch_data) initial_patch_stats <- minpatch_data$patch_stats - # Stage 1: Remove small patches (conditional) + + + # Stage 1: Remove small patches (conditional) ----- if (remove_small_patches) { if (verbose) cat("Stage 1: Removing small patches...\n") minpatch_data <- remove_small_patches_from_solution(minpatch_data) @@ -231,15 +235,13 @@ run_minpatch <- function(prioritizr_problem, minpatch_data$prioritizr_solution$minpatch <- create_solution_vector(minpatch_data$unit_dict) } - # Stage 2: Add new patches to meet targets (conditional) + # Stage 2: Add new patches to meet targets (conditional) ---- unmet_targets <- character(0) if (add_patches) { if (verbose) cat("Stage 2: Adding new patches...\n") - minpatch_data <- add_new_patches(minpatch_data, verbose) - # Check final unmet targets unmet_targets <- identify_unmet_targets(minpatch_data) @@ -253,7 +255,7 @@ run_minpatch <- function(prioritizr_problem, if (verbose) cat("Stage 2: Skipping addition of new patches...\n") } - # Stage 3: Simulated whittling (conditional) + # Stage 3: Simulated whittling (conditional) ---- if (whittle_patches) { if (verbose) cat("Stage 3: Removing unnecessary planning units...\n") minpatch_data <- simulated_whittling(minpatch_data, verbose) diff --git a/R/patch_functions.R b/R/patch_functions.R index 21ea8e3..4f6e050 100644 --- a/R/patch_functions.R +++ b/R/patch_functions.R @@ -1,6 +1,8 @@ #' Create patch dictionary from unit dictionary #' -#' Identifies connected components (patches) in the current solution +#' Identifies connected components (patches) in the current solution using igraph +#' and sparse matrix operations. This implementation follows the wheretowork approach +#' for efficient patch identification using matrix subsetting. #' #' @param minpatch_data List containing all MinPatch data structures #' @@ -11,53 +13,41 @@ make_patch_dict <- function(minpatch_data) { unit_dict <- minpatch_data$unit_dict boundary_matrix <- minpatch_data$boundary_matrix - # Get all selected planning units (status = 1 or 2) - selected_units <- names(unit_dict)[sapply(unit_dict, function(x) x$status %in% c(1, 2))] + # Get indices of selected planning units (status = 1 or 2) + selected_idx <- which(sapply(unit_dict, function(x) x$status %in% c(1, 2))) - - if (length(selected_units) == 0) { + if (length(selected_idx) == 0) { return(list()) } - patch_dict <- list() - patch_id <- 1 - remaining_units <- selected_units - - while (length(remaining_units) > 0) { - # Start a new patch with the first remaining unit - current_patch <- character(0) - units_to_check <- remaining_units[1] - - - # Find all connected units using breadth-first search - while (length(units_to_check) > 0) { - current_unit <- units_to_check[1] - units_to_check <- units_to_check[-1] + # Subset boundary matrix to only selected units (WTW approach) + adj_matrix <- boundary_matrix[selected_idx, selected_idx] + + # Convert to binary adjacency (1 if boundary exists, 0 otherwise) + adj_matrix@x <- rep(1, length(adj_matrix@x)) + + # Set diagonal to 1 (self-connections) + Matrix::diag(adj_matrix) <- 1 - if (!current_unit %in% current_patch) { - current_patch <- c(current_patch, current_unit) + # Create graph from adjacency matrix + g <- igraph::graph_from_adjacency_matrix(adj_matrix, mode = "undirected") - # Find neighbors of current unit that are also selected - neighbors <- names(boundary_matrix[[current_unit]]) - selected_neighbors <- intersect(neighbors, remaining_units) - selected_neighbors <- selected_neighbors[!selected_neighbors %in% current_patch] + # Identify components (patches) using igraph + clu <- igraph::components(g) + # Create patch dictionary from component membership + patch_dict <- list() + unit_names <- names(unit_dict)[selected_idx] - units_to_check <- unique(c(units_to_check, selected_neighbors)) - } - } - + for (patch_id in seq_len(clu$no)) { + # Get unit IDs in this component + unit_ids <- unit_names[clu$membership == patch_id] - # Store patch information patch_dict[[as.character(patch_id)]] <- list( area = 0, # Will be calculated later - unit_count = length(current_patch), - unit_ids = current_patch + unit_count = length(unit_ids), + unit_ids = unit_ids ) - - # Remove processed units from remaining units - remaining_units <- setdiff(remaining_units, current_patch) - patch_id <- patch_id + 1 } return(patch_dict) diff --git a/R/whittling_functions.R b/R/whittling_functions.R index e8a926f..7e68849 100644 --- a/R/whittling_functions.R +++ b/R/whittling_functions.R @@ -21,6 +21,9 @@ simulated_whittling <- function(minpatch_data, verbose = TRUE) { iteration <- 0 max_iterations <- 10000 # Prevent infinite loops keystone_pu_id_set <- character(0) # Initialize keystone set + + # OPTIMIZATION: Cache feature amounts and only update when unit is removed + feature_amounts <- calculate_feature_conservation(minpatch_data) while (iteration < max_iterations) { iteration <- iteration + 1 @@ -45,8 +48,8 @@ simulated_whittling <- function(minpatch_data, verbose = TRUE) { break # No more edge units to consider } - # Calculate whittling scores for edge units - whittle_scores_raw <- calculate_whittle_scores(edge_units, minpatch_data) + # Calculate whittling scores for edge units (pass cached feature amounts) + whittle_scores_raw <- calculate_whittle_scores(edge_units, minpatch_data, feature_amounts) # Separate keystone units (needed for targets) from scoreable units whittle_scores <- list() @@ -74,11 +77,19 @@ simulated_whittling <- function(minpatch_data, verbose = TRUE) { # Find the unit with the lowest whittling score (most suitable for removal) candidate_unit <- names(whittle_scores)[which.min(unlist(whittle_scores))] - # Check if removing this unit is acceptable - removal_result <- can_remove_unit(candidate_unit, minpatch_data) + # Check if removing this unit is acceptable (pass cached feature amounts) + removal_result <- can_remove_unit(candidate_unit, minpatch_data, feature_amounts) if (removal_result) { # Remove the unit unit_dict[[candidate_unit]]$status <- 0 + + # OPTIMIZATION: Update cached feature amounts incrementally instead of recalculating + unit_abundances <- abundance_matrix[[candidate_unit]] + for (feat_id in names(unit_abundances)) { + if (feat_id %in% names(feature_amounts)) { + feature_amounts[feat_id] <- feature_amounts[feat_id] - unit_abundances[[feat_id]] + } + } if (verbose && iteration <= 10) { cat(" Removed unit", candidate_unit, "at iteration", iteration, "\n") @@ -107,7 +118,7 @@ simulated_whittling <- function(minpatch_data, verbose = TRUE) { #' Find edge planning units #' #' Identifies planning units that are on the edge of selected areas -#' (have at least one unselected neighbor) +#' (have at least one unselected neighbor). Optimized with vector pre-allocation. #' #' @param minpatch_data List containing all MinPatch data structures #' @@ -117,55 +128,72 @@ find_edge_units <- function(minpatch_data) { unit_dict <- minpatch_data$unit_dict boundary_matrix <- minpatch_data$boundary_matrix - - edge_units <- character(0) + + n_units <- length(unit_dict) + + # Pre-allocate for worst case (all selected units are edge units) + edge_units <- character(n_units) + edge_count <- 0 for (unit_id in names(unit_dict)) { if (unit_dict[[unit_id]]$status %in% c(1, 2)) { # Consider both selected (1) and conserved (2) units - # Check if this unit has any unselected neighbors - neighbors <- names(boundary_matrix[[unit_id]]) + # Get neighbors from sparse matrix (units with non-zero boundary) + unit_idx <- as.integer(unit_id) + neighbor_indices <- which(boundary_matrix[unit_idx, ] > 0) + neighbor_ids <- colnames(boundary_matrix)[neighbor_indices] - for (neighbor_id in neighbors) { + is_edge <- FALSE + for (neighbor_id in neighbor_ids) { if (neighbor_id != unit_id && neighbor_id %in% names(unit_dict)) { neighbor_status <- unit_dict[[neighbor_id]]$status # If neighbor is available (0) or excluded (3), this is an edge unit if (neighbor_status %in% c(0, 3)) { - edge_units <- c(edge_units, unit_id) + is_edge <- TRUE break } } else if (neighbor_id == unit_id) { - # Self-reference indicates external boundary - edge_units <- c(edge_units, unit_id) + # Self-reference (diagonal) indicates external boundary + is_edge <- TRUE break } } + + if (is_edge) { + edge_count <- edge_count + 1 + edge_units[edge_count] <- unit_id + } } } - - return(unique(edge_units)) + + # Trim to actual size and return unique values + return(unique(edge_units[seq_len(edge_count)])) } #' Calculate whittling scores for edge units #' #' Calculates the "Low Relevance" score for each edge unit based on -#' feature importance (Equation A2 from the original paper) +#' feature importance (Equation A2 from the original paper). +#' Optimized to accept pre-computed feature amounts to avoid redundant calculations. #' #' @param edge_units Character vector of edge unit IDs #' @param minpatch_data List containing all MinPatch data structures +#' @param feature_amounts Optional pre-computed feature conservation amounts #' #' @return Named vector of whittling scores #' @keywords internal -calculate_whittle_scores <- function(edge_units, minpatch_data) { +calculate_whittle_scores <- function(edge_units, minpatch_data, feature_amounts = NULL) { unit_dict <- minpatch_data$unit_dict abundance_matrix <- minpatch_data$abundance_matrix target_dict <- minpatch_data$target_dict - # Calculate current feature conservation amounts - feature_amounts <- calculate_feature_conservation(minpatch_data) + # Calculate current feature conservation amounts only if not provided + if (is.null(feature_amounts)) { + feature_amounts <- calculate_feature_conservation(minpatch_data) + } scores <- list() @@ -230,7 +258,7 @@ calculate_whittle_scores <- function(edge_units, minpatch_data) { #' #' @return Logical indicating if unit can be removed #' @keywords internal -can_remove_unit <- function(unit_id, minpatch_data) { +can_remove_unit <- function(unit_id, minpatch_data, feature_amounts = NULL) { unit_dict <- minpatch_data$unit_dict min_patch_size <- minpatch_data$min_patch_size @@ -240,8 +268,8 @@ can_remove_unit <- function(unit_id, minpatch_data) { return(FALSE) } - # Check if removal would violate conservation targets - if (removal_violates_targets(unit_id, minpatch_data)) { + # Check if removal would violate conservation targets (pass cached feature amounts) + if (removal_violates_targets(unit_id, minpatch_data, feature_amounts)) { return(FALSE) } @@ -265,19 +293,24 @@ can_remove_unit <- function(unit_id, minpatch_data) { #' Check if removing unit would violate conservation targets #' +#' Optimized to accept pre-computed feature amounts to avoid redundant calculations. +#' #' @param unit_id ID of unit to potentially remove #' @param minpatch_data List containing all MinPatch data structures +#' @param feature_amounts Optional pre-computed feature conservation amounts #' #' @return Logical indicating if removal would violate targets #' @keywords internal -removal_violates_targets <- function(unit_id, minpatch_data) { +removal_violates_targets <- function(unit_id, minpatch_data, feature_amounts = NULL) { unit_dict <- minpatch_data$unit_dict abundance_matrix <- minpatch_data$abundance_matrix target_dict <- minpatch_data$target_dict - # Calculate current feature amounts - feature_amounts <- calculate_feature_conservation(minpatch_data) + # Calculate current feature amounts only if not provided + if (is.null(feature_amounts)) { + feature_amounts <- calculate_feature_conservation(minpatch_data) + } # Get features in this unit unit_abundances <- abundance_matrix[[unit_id]] @@ -354,14 +387,19 @@ removal_increases_cost <- function(unit_id, minpatch_data) { unit_cost <- unit_dict[[unit_id]]$cost + # Get neighbors from sparse matrix + unit_idx <- as.integer(unit_id) + neighbor_indices <- which(boundary_matrix[unit_idx, ] > 0) + neighbor_ids <- colnames(boundary_matrix)[neighbor_indices] + # Calculate change in boundary cost - neighbors <- names(boundary_matrix[[unit_id]]) boundary_change <- 0 - for (neighbor_id in neighbors) { + for (neighbor_id in neighbor_ids) { if (neighbor_id %in% names(unit_dict)) { neighbor_status <- unit_dict[[neighbor_id]]$status - boundary_length <- boundary_matrix[[unit_id]][[neighbor_id]] + neighbor_idx <- as.integer(neighbor_id) + boundary_length <- boundary_matrix[unit_idx, neighbor_idx] if (neighbor_status %in% c(1, 2)) { # Removing unit increases boundary (neighbor becomes edge) diff --git a/README.Rmd b/README.Rmd index 9d36e6f..d1555cc 100644 --- a/README.Rmd +++ b/README.Rmd @@ -66,7 +66,6 @@ pak::pak("SpatialPlanning/minpatch") - **Flexible Parameters**: Control minimum patch sizes, patch radius, and boundary penalties - **Comprehensive Reporting**: Detailed statistics and comparisons - **Visualization Support**: Plot results with ggplot2 (optional) -- **Well Documented**: Extensive documentation and examples ## Algorithm Details diff --git a/README.md b/README.md index 8611a54..281fe26 100644 --- a/README.md +++ b/README.md @@ -76,7 +76,6 @@ pak::pak("SpatialPlanning/minpatch") and boundary penalties - **Comprehensive Reporting**: Detailed statistics and comparisons - **Visualization Support**: Plot results with ggplot2 (optional) -- **Well Documented**: Extensive documentation and examples ## Algorithm Details diff --git a/Untitled.Rprofvis b/Untitled.Rprofvis new file mode 100644 index 0000000..f06a1f6 --- /dev/null +++ b/Untitled.Rprofvis @@ -0,0 +1,4637 @@ + + + + +profvis + + + + + + + + + + + + + +
+
+
+ + + + diff --git a/docs/articles/minpatch.html b/docs/articles/minpatch.html index 442c63f..de067dd 100644 --- a/docs/articles/minpatch.html +++ b/docs/articles/minpatch.html @@ -57,10 +57,7 @@

MinPatch with prioritizr

-

MinPatch -Package Authors

-

2025-11-24

Source: vignettes/minpatch.Rmd
minpatch.Rmd
@@ -87,8 +84,7 @@

IntroductionStep 1: Load and Examine the Data 

-
-dat <- c(get_sim_pu_raster(), get_sim_features()) %>% 
+dat <- c(get_sim_pu_raster(), get_sim_features()) %>% 
   as.polygons(dissolve = FALSE, values = TRUE) %>%
   sf::st_as_sf() %>% 
   dplyr::rename(cost = layer)
@@ -146,8 +142,7 @@ 
Step 3: Run MinPatch
-
-result <- run_minpatch(
+result <- run_minpatch(
   prioritizr_problem = p,
   prioritizr_solution = s,
   min_patch_size = min_patch_size,
@@ -160,8 +155,8 @@ 
Step 3: Run MinPatch)
 #> Validating inputs...
 #> Initializing data structures...
-#> Calculating boundary matrix (this may take a while)...
-#> Creating patch radius dictionary...
+#> Calculating boundary matrix (optimized version)...
+#> Creating patch radius dictionary (optimized)...
 #> Calculating initial patch statistics...
 #> Stage 1: Removing small patches...
 #> Stage 2: Adding new patches...
@@ -301,8 +296,7 @@ 
Step 4: Analyze the ResultsFeature Representation Analysis
 

 
+#>  [1] gtable_0.3.6         xfun_0.54            bslib_0.9.0         
+#>  [4] raster_3.6-32        htmlwidgets_1.6.4    lattice_0.22-7      
+#>  [7] vctrs_0.6.5          tools_4.5.0          generics_0.1.4      
+#> [10] parallel_4.5.0       tibble_3.3.0         proxy_0.4-27        
+#> [13] pkgconfig_2.0.3      Matrix_1.7-4         KernSmooth_2.23-26  
+#> [16] RColorBrewer_1.1-3   S7_0.2.0             desc_1.4.3          
+#> [19] assertthat_0.2.1     lifecycle_1.0.4      compiler_4.5.0      
+#> [22] farver_2.1.2         stringr_1.6.0        textshaping_1.0.4   
+#> [25] codetools_0.2-20     htmltools_0.5.8.1    class_7.3-23        
+#> [28] sass_0.4.10          yaml_2.3.10          pillar_1.11.1       
+#> [31] pkgdown_2.2.0        exactextractr_0.10.0 jquerylib_0.1.4     
+#> [34] rcbc_0.1.0.9003      classInt_0.4-11      cachem_1.1.0        
+#> [37] nlme_3.1-168         parallelly_1.45.1    tidyselect_1.2.1    
+#> [40] digest_0.6.38        stringi_1.8.7        fastmap_1.2.0       
+#> [43] grid_4.5.0           cli_3.6.5            magrittr_2.0.4      
+#> [46] dichromat_2.0-0.1    e1071_1.7-16         ape_5.8-1           
+#> [49] withr_3.0.2          scales_1.4.0         sp_2.2-0            
+#> [52] rmarkdown_2.30       igraph_2.2.1         ragg_1.5.0          
+#> [55] evaluate_1.0.5       knitr_1.50           rlang_1.1.6         
+#> [58] Rcpp_1.1.0           glue_1.8.0           DBI_1.2.3           
+#> [61] rstudioapi_0.17.1    jsonlite_2.0.0       R6_2.6.1            
+#> [64] systemfonts_1.3.1    fs_1.6.6             units_1.0-0