From f406b70adc9aaeb93aae5487aebf5abc4c22c507 Mon Sep 17 00:00:00 2001 From: Jason Everett Date: Mon, 24 Nov 2025 21:58:19 +1100 Subject: [PATCH 1/8] Fix errors with checks --- .github/workflows/MacOS.yaml | 1 + .github/workflows/Ubuntu.yaml | 1 + .github/workflows/Windows.yaml | 1 + 3 files changed, 3 insertions(+) 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..562a505 100644 --- a/.github/workflows/Ubuntu.yaml +++ b/.github/workflows/Ubuntu.yaml @@ -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..76154cf 100644 --- a/.github/workflows/Windows.yaml +++ b/.github/workflows/Windows.yaml @@ -50,6 +50,7 @@ jobs: any::rcmdcheck any::remotes any::lpsymphony + any::slam needs: check - uses: r-lib/actions/check-r-package@v2 From 48a6bd797ff2036ddf434733bed4a70483f49085 Mon Sep 17 00:00:00 2001 From: Jason Everett Date: Tue, 25 Nov 2025 08:03:53 +1100 Subject: [PATCH 2/8] Website updates --- R/minpatch.R | 12 +++++++----- README.Rmd | 1 - README.md | 1 - docs/articles/minpatch.html | 12 +++--------- docs/articles/minpatch.md | 3 --- docs/articles/minpatchTasmania.html | 9 +++------ docs/articles/minpatchTasmania.md | 3 --- docs/pkgdown.yml | 2 +- man/minpatch-package.Rd | 4 ++-- vignettes/minpatch.Rmd | 6 ------ vignettes/minpatchTasmania.Rmd | 8 -------- 11 files changed, 16 insertions(+), 45 deletions(-) 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/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/docs/articles/minpatch.html b/docs/articles/minpatch.html index 442c63f..5b2685c 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,
@@ -301,8 +296,7 @@ 
Step 4: Analyze the ResultsFeature Representation Analysis
 

 
-
-# Create solution data for prioritizr analysis
+# Create solution data for prioritizr analysis
 minpatch_solution_data <- result$solution[c("minpatch")]
 
 # Use prioritizr functions for accurate feature representation analysis
diff --git a/docs/articles/minpatch.md b/docs/articles/minpatch.md
index cd46381..e723364 100644
--- a/docs/articles/minpatch.md
+++ b/docs/articles/minpatch.md
@@ -22,7 +22,6 @@ library(patchwork)
 
 ``` r
 
-
 dat <- c(get_sim_pu_raster(), get_sim_features()) %>% 
   as.polygons(dissolve = FALSE, values = TRUE) %>%
   sf::st_as_sf() %>% 
@@ -92,7 +91,6 @@ Run MinPatch with automatic data extraction from prioritizr objects
 
 ``` r
 
-
 result <- run_minpatch(
   prioritizr_problem = p,
   prioritizr_solution = s,
@@ -251,7 +249,6 @@ print(comparison$summary)
 
 ``` r
 
-
 # Create solution data for prioritizr analysis
 minpatch_solution_data <- result$solution[c("minpatch")]
 
diff --git a/docs/articles/minpatchTasmania.html b/docs/articles/minpatchTasmania.html
index ed1b1bb..4955424 100644
--- a/docs/articles/minpatchTasmania.html
+++ b/docs/articles/minpatchTasmania.html
@@ -157,8 +157,7 @@ 
Choose a patch sizeRun minpatch
 

 
-
-result <- run_minpatch(
+result <- run_minpatch(
   prioritizr_problem = p,
   prioritizr_solution = s,
   min_patch_size = min_patch_size,
@@ -456,8 +455,7 @@ 
Run different patch sizes
 
-
-# Calculate reasonable parameters based on planning unit characteristics
+# Calculate reasonable parameters based on planning unit characteristics
 median_area <- median(st_area(tas))
 min_patch_size <- median_area * 10
 patch_radius <- sqrt(median_area * 10)
@@ -523,8 +521,7 @@ 
Run different Boundary Penalties
 
-
-# Calculate reasonable parameters based on planning unit characteristics
+# Calculate reasonable parameters based on planning unit characteristics
 median_area <- median(st_area(tas))
 min_patch_size <- median_area * 5
 patch_radius <- sqrt(median_area * 10)
diff --git a/docs/articles/minpatchTasmania.md b/docs/articles/minpatchTasmania.md
index 3bfcc93..23c7238 100644
--- a/docs/articles/minpatchTasmania.md
+++ b/docs/articles/minpatchTasmania.md
@@ -97,7 +97,6 @@ cat("- This means patches must be at least", round(min_patch_size/median_area, 3
 
 ``` r
 
-
 result <- run_minpatch(
   prioritizr_problem = p,
   prioritizr_solution = s,
@@ -400,7 +399,6 @@ areas.
 
 ``` r
 
-
 # Calculate reasonable parameters based on planning unit characteristics
 median_area <- median(st_area(tas))
 min_patch_size <- median_area * 10
@@ -470,7 +468,6 @@ needed to see meaningful differences in such cases.
 
 ``` r
 
-
 # Calculate reasonable parameters based on planning unit characteristics
 median_area <- median(st_area(tas))
 min_patch_size <- median_area * 5
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
index 7db4512..777aa0f 100644
--- a/docs/pkgdown.yml
+++ b/docs/pkgdown.yml
@@ -4,4 +4,4 @@ pkgdown_sha: ~
 articles:
   minpatch: minpatch.html
   minpatchTasmania: minpatchTasmania.html
-last_built: 2025-11-24T09:41Z
+last_built: 2025-11-24T11:05Z
diff --git a/man/minpatch-package.Rd b/man/minpatch-package.Rd
index 017bcb3..527e320 100644
--- a/man/minpatch-package.Rd
+++ b/man/minpatch-package.Rd
@@ -13,8 +13,8 @@ thresholds, following the methodology described in Smith et al. (2010).
 \seealso{
 Useful links:
 \itemize{
-  \item \url{https://github.com/yourusername/minpatch}
-  \item Report bugs at \url{https://github.com/yourusername/minpatch/issues}
+  \item \url{https://github.com/SpatialPlanning/minpatch}
+  \item Report bugs at \url{https://github.com/SpatialPlanning/minpatch/issues}
 }
 
 }
diff --git a/vignettes/minpatch.Rmd b/vignettes/minpatch.Rmd
index 7848e95..955529f 100644
--- a/vignettes/minpatch.Rmd
+++ b/vignettes/minpatch.Rmd
@@ -1,7 +1,5 @@
 ---
 title: "MinPatch with prioritizr"
-author: "MinPatch Package Authors"
-date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteIndexEntry{MinPatch with prioritizr}
@@ -36,7 +34,6 @@ library(patchwork)
 ## Step 1: Load and Examine the Data
 
 ```{r load-data}
-
 dat <- c(get_sim_pu_raster(), get_sim_features()) %>% 
   as.polygons(dissolve = FALSE, values = TRUE) %>%
   sf::st_as_sf() %>% 
@@ -95,7 +92,6 @@ cat("- This means patches must be at least", round(min_patch_size/median_area, 3
 Run MinPatch with automatic data extraction from prioritizr objects
 
 ```{r run-minpatch}
-
 result <- run_minpatch(
   prioritizr_problem = p,
   prioritizr_solution = s,
@@ -143,7 +139,6 @@ print(comparison$summary)
 ### Feature Representation Analysis
 
 ```{r feature-analysis}
-
 # Create solution data for prioritizr analysis
 minpatch_solution_data <- result$solution[c("minpatch")]
 
@@ -196,7 +191,6 @@ Let's create maps to visualize the changes MinPatch made:
 
 ```{r visualization}
 plot_minpatch(result, title = "MinPatch Results")
-
 ```
 
 ## Understanding the Results
diff --git a/vignettes/minpatchTasmania.Rmd b/vignettes/minpatchTasmania.Rmd
index 06e1642..b539f69 100644
--- a/vignettes/minpatchTasmania.Rmd
+++ b/vignettes/minpatchTasmania.Rmd
@@ -56,7 +56,6 @@ tas <- tas %>%
 # Run prioritizr analysis
 
 ```{r}
-
 p <- problem(tas, features = features, cost_column = "cost") %>%
   add_min_set_objective() %>%
   add_relative_targets(0.30) %>%  # 30% of each feature
@@ -64,14 +63,11 @@ p <- problem(tas, features = features, cost_column = "cost") %>%
   add_default_solver(verbose = FALSE)
 
 s <- solve(p)
-
 ```
 
 ## Plot prioritizr solution
 ```{r}
-
 plot_prioritizr(s)
-
 ```
 
 
@@ -80,7 +76,6 @@ plot_prioritizr(s)
 ## Choose a patch size
 
 ```{r}
-
 # Calculate reasonable parameters based on planning unit characteristics
 median_area <- median(st_area(tas))
 
@@ -102,7 +97,6 @@ cat("- This means patches must be at least", round(min_patch_size/median_area, 3
 ## Run minpatch
 
 ```{r}
-
 result <- run_minpatch(
   prioritizr_problem = p,
   prioritizr_solution = s,
@@ -156,7 +150,6 @@ The minimum patch size parameter is the core constraint that drives MinPatch beh
 Larger minimum patch sizes result in fewer, bigger patches with potentially higher total area, as MinPatch must ensure every patch meets the size requirement. Smaller minimum patch sizes allow more flexibility, potentially resulting in more patches that are closer to the original *prioritizr* solution. The choice of minimum patch size should reflect ecological or management considerations - for example, larger patches may be needed to support viable populations or reduce edge effects, while smaller patches may be acceptable in highly connected landscapes or for features that don't require large contiguous areas.
 
 ```{r}
-
 # Calculate reasonable parameters based on planning unit characteristics
 median_area <- median(st_area(tas))
 min_patch_size <- median_area * 10
@@ -208,7 +201,6 @@ patchwork::wrap_plots(plot_minpatch(result, title = "Patch Size x5"),
 The boundary penalty controls how much MinPatch prioritizes spatial compactness during the "simulated whittling" stage. During whittling, MinPatch considers removing planning units from patch edges, but only if doing so doesn't increase the total *priortizr* cost. The boundary penalty affects this decision by penalizing fragmented solutions - higher penalties favour more compact patches by making it costly to create additional "edge" between selected and unselected areas. When MinPatch evaluates whether to remove a unit, it calculates the change in boundary length (units with selected neighbours increase boundary when removed, while units with unselected neighbours decrease boundary) and multiplies this by the boundary penalty. If the resulting boundary cost change exceeds the unit's cost, the unit cannot be removed. In datasets like Tasmania with long planning unit boundaries relative to unit costs, even small boundary penalties can be highly influential, potentially preventing most unit removals and resulting in similar solutions across different penalty values. Very small penalties (e.g., 1e-10) may be needed to see meaningful differences in such cases.
 
 ```{r}
-
 # Calculate reasonable parameters based on planning unit characteristics
 median_area <- median(st_area(tas))
 min_patch_size <- median_area * 5

From 8b08e3dbdbc3b61337338dfa3a58e332f38d44a2 Mon Sep 17 00:00:00 2001
From: Jason Everett 
Date: Tue, 25 Nov 2025 08:10:24 +1100
Subject: [PATCH 3/8] Add slam to checks

---
 .github/workflows/test-coverage.yaml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index fd20dfa..d7f4474 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -40,6 +40,7 @@ jobs:
             any::xml2
             any::remotes
             any::lpsymphony
+            any::slam
           needs: coverage
 
       - name: Test coverage

From ebb51750e6367eb7348e332e77c0afeb6fb7f22e Mon Sep 17 00:00:00 2001
From: Jason Everett 
Date: Tue, 25 Nov 2025 08:14:35 +1100
Subject: [PATCH 4/8] Updates to include igraph

---
 DESCRIPTION         |  3 ++-
 R/patch_functions.R | 66 +++++++++++++++++++++++----------------------
 2 files changed, 36 insertions(+), 33 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 6f9060e..cdd2777 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -41,7 +41,8 @@ Imports:
     prioritizr (>= 7.0.0),
     rlang,
     stats,
-    magrittr
+    magrittr,
+    igraph
 Suggests:
     dplyr,
     ggplot2,
diff --git a/R/patch_functions.R b/R/patch_functions.R
index 21ea8e3..6aedb30 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.
+#' This implementation uses graph theory algorithms (via igraph) for faster
+#' performance, especially with large planning problems.
 #'
 #' @param minpatch_data List containing all MinPatch data structures
 #'
@@ -14,50 +16,50 @@ make_patch_dict <- function(minpatch_data) {
   # 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))]
 
-
   if (length(selected_units) == 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]
+  # Create adjacency matrix for selected units
+  n_selected <- length(selected_units)
+  adj_matrix <- matrix(0, nrow = n_selected, ncol = n_selected)
+  rownames(adj_matrix) <- selected_units
+  colnames(adj_matrix) <- selected_units
+
+  # Fill adjacency matrix based on boundary_matrix
+  for (i in seq_along(selected_units)) {
+    unit_i <- selected_units[i]
+    neighbors <- names(boundary_matrix[[unit_i]])
+
+    for (j in seq_along(selected_units)) {
+      unit_j <- selected_units[j]
+      if (unit_j %in% neighbors) {
+        adj_matrix[i, j] <- 1
+      }
+    }
+  }
 
-      if (!current_unit %in% current_patch) {
-        current_patch <- c(current_patch, current_unit)
+  # Set diagonal to 1 (self-connections)
+  diag(adj_matrix) <- 1
 
-        # 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]
+  # Create graph from adjacency matrix
+  g <- igraph::graph_from_adjacency_matrix(adj_matrix, mode = "undirected")
 
+  # Identify components (patches) using igraph
+  clu <- igraph::components(g)
 
-        units_to_check <- unique(c(units_to_check, selected_neighbors))
-      }
-    }
+  # Create patch dictionary from component membership
+  patch_dict <- list()
 
+  for (patch_id in seq_len(clu$no)) {
+    # Get unit IDs in this component
+    unit_ids <- selected_units[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)

From 7887afbbe52df3d5fff1bae50c451412126e2e66 Mon Sep 17 00:00:00 2001
From: Jason Everett 
Date: Tue, 25 Nov 2025 11:58:11 +1100
Subject: [PATCH 5/8] Remove bottlenecks and add parallel

---
 DESCRIPTION             |    6 +-
 R/data_structures.R     |  195 +-
 R/patch_functions.R     |   42 +-
 R/whittling_functions.R |   94 +-
 Untitled.Rprofvis       | 4637 +++++++++++++++++++++++++++++++++++++++
 5 files changed, 4864 insertions(+), 110 deletions(-)
 create mode 100644 Untitled.Rprofvis

diff --git a/DESCRIPTION b/DESCRIPTION
index cdd2777..8e9f343 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -42,7 +42,8 @@ Imports:
     rlang,
     stats,
     magrittr,
-    igraph
+    igraph,
+    Matrix
 Suggests:
     dplyr,
     ggplot2,
@@ -53,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/patch_functions.R b/R/patch_functions.R
index 6aedb30..4f6e050 100644
--- a/R/patch_functions.R
+++ b/R/patch_functions.R
@@ -1,8 +1,8 @@
 #' Create patch dictionary from unit dictionary
 #'
-#' Identifies connected components (patches) in the current solution using igraph.
-#' This implementation uses graph theory algorithms (via igraph) for faster
-#' performance, especially with large planning problems.
+#' 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
 #'
@@ -13,34 +13,21 @@ 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())
   }
 
-  # Create adjacency matrix for selected units
-  n_selected <- length(selected_units)
-  adj_matrix <- matrix(0, nrow = n_selected, ncol = n_selected)
-  rownames(adj_matrix) <- selected_units
-  colnames(adj_matrix) <- selected_units
-
-  # Fill adjacency matrix based on boundary_matrix
-  for (i in seq_along(selected_units)) {
-    unit_i <- selected_units[i]
-    neighbors <- names(boundary_matrix[[unit_i]])
-
-    for (j in seq_along(selected_units)) {
-      unit_j <- selected_units[j]
-      if (unit_j %in% neighbors) {
-        adj_matrix[i, j] <- 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)
-  diag(adj_matrix) <- 1
+  Matrix::diag(adj_matrix) <- 1
 
   # Create graph from adjacency matrix
   g <- igraph::graph_from_adjacency_matrix(adj_matrix, mode = "undirected")
@@ -50,10 +37,11 @@ make_patch_dict <- function(minpatch_data) {
 
   # Create patch dictionary from component membership
   patch_dict <- list()
+  unit_names <- names(unit_dict)[selected_idx]
 
   for (patch_id in seq_len(clu$no)) {
     # Get unit IDs in this component
-    unit_ids <- selected_units[clu$membership == patch_id]
+    unit_ids <- unit_names[clu$membership == patch_id]
 
     patch_dict[[as.character(patch_id)]] <- list(
       area = 0,  # Will be calculated later
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/Untitled.Rprofvis b/Untitled.Rprofvis
new file mode 100644
index 0000000..f06a1f6
--- /dev/null
+++ b/Untitled.Rprofvis
@@ -0,0 +1,4637 @@
+
+
+
+
+profvis
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+  

+

+
+
+
+

From 081e52b3fe007a2cd844dd6e293b761093b62136 Mon Sep 17 00:00:00 2001
From: Jason Everett 
Date: Tue, 25 Nov 2025 12:33:13 +1100
Subject: [PATCH 6/8] Update website

---
 docs/articles/minpatch.html                | 48 +++++++++++-----------
 docs/articles/minpatch.md                  | 48 +++++++++++-----------
 docs/articles/minpatchTasmania.html        | 26 ++++--------
 docs/articles/minpatchTasmania.md          | 20 ++-------
 docs/index.html                            |  2 -
 docs/index.md                              |  1 -
 docs/llms.txt                              |  1 -
 docs/pkgdown.yml                           |  2 +-
 docs/reference/minpatch-package.html       |  4 +-
 docs/reference/minpatch-package.md         |  4 +-
 docs/reference/print_minpatch_summary.html |  4 +-
 docs/reference/print_minpatch_summary.md   |  4 +-
 docs/reference/run_minpatch.html           |  8 ++--
 docs/reference/run_minpatch.md             |  8 ++--
 docs/search.json                           |  2 +-
 15 files changed, 77 insertions(+), 105 deletions(-)

diff --git a/docs/articles/minpatch.html b/docs/articles/minpatch.html
index 5b2685c..de067dd 100644
--- a/docs/articles/minpatch.html
+++ b/docs/articles/minpatch.html
@@ -155,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...
@@ -617,28 +617,28 @@ 
Session Information#> [5] sf_1.0-22        prioritizr_8.1.0 minpatch_0.1.0  
 #> 
 #> loaded via a namespace (and not attached):
-#>  [1] generics_0.1.4       sass_0.4.10          class_7.3-23        
-#>  [4] KernSmooth_2.23-26   stringi_1.8.7        lattice_0.22-7      
-#>  [7] digest_0.6.38        magrittr_2.0.4       RColorBrewer_1.1-3  
-#> [10] evaluate_1.0.5       grid_4.5.0           fastmap_1.2.0       
-#> [13] jsonlite_2.0.0       Matrix_1.7-4         ape_5.8-1           
-#> [16] e1071_1.7-16         DBI_1.2.3            scales_1.4.0        
-#> [19] codetools_0.2-20     textshaping_1.0.4    jquerylib_0.1.4     
-#> [22] cli_3.6.5            rlang_1.1.6          units_1.0-0         
-#> [25] withr_3.0.2          cachem_1.1.0         yaml_2.3.10         
-#> [28] tools_4.5.0          raster_3.6-32        parallel_4.5.0      
-#> [31] rcbc_0.1.0.9003      assertthat_0.2.1     exactextractr_0.10.0
-#> [34] vctrs_0.6.5          R6_2.6.1             proxy_0.4-27        
-#> [37] lifecycle_1.0.4      classInt_0.4-11      stringr_1.6.0       
-#> [40] fs_1.6.6             htmlwidgets_1.6.4    ragg_1.5.0          
-#> [43] pkgconfig_2.0.3      desc_1.4.3           gtable_0.3.6        
-#> [46] pkgdown_2.2.0        bslib_0.9.0          pillar_1.11.1       
-#> [49] glue_1.8.0           Rcpp_1.1.0           systemfonts_1.3.1   
-#> [52] tidyselect_1.2.1     xfun_0.54            tibble_3.3.0        
-#> [55] dichromat_2.0-0.1    rstudioapi_0.17.1    knitr_1.50          
-#> [58] farver_2.1.2         htmltools_0.5.8.1    nlme_3.1-168        
-#> [61] rmarkdown_2.30       compiler_4.5.0       S7_0.2.0            
-#> [64] sp_2.2-0

+#>  [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