From eec7edb28dbb379f108f915f1cb4bb97f9064e57 Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Mon, 11 Aug 2025 09:49:28 +0930
Subject: [PATCH 1/9] Working on using factors for spatial columns

---
 R/plotting.R                   |  20 +++-
 tests/testthat/test-plotting.R | 207 +++++++++++++++++++++++++++++++++
 2 files changed, 223 insertions(+), 4 deletions(-)
 create mode 100644 tests/testthat/test-plotting.R

diff --git a/R/plotting.R b/R/plotting.R
index 40ef553..a580028 100644
--- a/R/plotting.R
+++ b/R/plotting.R
@@ -112,6 +112,15 @@ autoplot.design <- function(object, rotation = 0, size = 4, margin = FALSE, pale
     block_expr <- rlang::quo_name(block_expr)
     trt_expr <- rlang::quo_name(trt_expr)
 
+    # Handle factor columns for spatial coordinates
+    # Convert factor row/column to numeric while preserving order
+    if(is.factor(object[[row_expr]])) {
+        object[[row_expr]] <- as.numeric(as.character(object[[row_expr]]))
+    }
+    if(is.factor(object[[column_expr]])) {
+        object[[column_expr]] <- as.numeric(as.character(object[[column_expr]]))
+    }
+
     object[[trt_expr]] <- factor(as.character(object[[trt_expr]]), levels = unique(stringi::stri_sort(as.character(object[[trt_expr]]), numeric = TRUE)))
     ntrt <- nlevels(object[[trt_expr]])
 
@@ -177,10 +186,13 @@ autoplot.design <- function(object, rotation = 0, size = 4, margin = FALSE, pale
         }
         for (i in 1:nrow(blkdf)) {
             tmp <- object[object[[block_expr]] == blkdf$block[i], ]
-            blkdf[i, "ymin"] <- (min(tmp$row) - 0.5)
-            blkdf[i, "ymax"] <- (max(tmp$row) + 0.5)
-            blkdf[i, "xmin"] <- (min(tmp$col) - 0.5)
-            blkdf[i, "xmax"] <- (max(tmp$col) + 0.5)
+            # Ensure row and col are numeric for min/max calculations
+            tmp_row <- if(is.factor(tmp$row)) as.numeric(as.character(tmp$row)) else tmp$row
+            tmp_col <- if(is.factor(tmp$col)) as.numeric(as.character(tmp$col)) else tmp$col
+            blkdf[i, "ymin"] <- (min(tmp_row) - 0.5)
+            blkdf[i, "ymax"] <- (max(tmp_row) + 0.5)
+            blkdf[i, "xmin"] <- (min(tmp_col) - 0.5)
+            blkdf[i, "xmax"] <- (max(tmp_col) + 0.5)
         }
 
         plt <- ggplot2::ggplot(...) +
diff --git a/tests/testthat/test-plotting.R b/tests/testthat/test-plotting.R
new file mode 100644
index 0000000..047b384
--- /dev/null
+++ b/tests/testthat/test-plotting.R
@@ -0,0 +1,207 @@
+test_that("autoplot handles factor row and column variables", {
+  # Create a design with factor row and column variables
+  df <- data.frame(
+    row = factor(rep(1:3, each = 3), levels = 1:3),
+    col = factor(rep(1:3, times = 3), levels = 1:3),
+    treatment = rep(LETTERS[1:3], 3)
+  )
+
+  # Verify columns are factors
+  expect_true(is.factor(df$row))
+  expect_true(is.factor(df$col))
+
+  # Set seed for reproducibility
+  set.seed(42)
+
+  # Run speed optimization
+  result <- speed(df, "treatment", iterations = 100)
+
+  # Test that autoplot works with factor columns
+  expect_no_error(plot <- autoplot(result))
+  expect_s3_class(plot, "ggplot")
+
+  # Check that the plot data has numeric coordinates after conversion
+  plot_data <- plot$data
+  expect_true(is.numeric(plot_data$row))
+  expect_true(is.numeric(plot_data$col))
+})
+
+test_that("autoplot handles factor columns with non-sequential levels", {
+  # Create a design with factor columns that have non-sequential levels
+  df <- data.frame(
+    row = factor(rep(c(2, 4, 6), each = 3), levels = c(2, 4, 6)),
+    col = factor(rep(c(1, 3, 5), times = 3), levels = c(1, 3, 5)),
+    treatment = rep(LETTERS[1:3], 3)
+  )
+
+  # Verify columns are factors with non-sequential levels
+  expect_true(is.factor(df$row))
+  expect_true(is.factor(df$col))
+  expect_equal(levels(df$row), c("2", "4", "6"))
+  expect_equal(levels(df$col), c("1", "3", "5"))
+
+  # Set seed for reproducibility
+  set.seed(42)
+
+  # Run speed optimization
+  result <- speed(df, "treatment", iterations = 100)
+
+  # Test that autoplot works and converts factors appropriately
+  expect_no_error(plot <- autoplot(result))
+  expect_s3_class(plot, "ggplot")
+
+  # Check that the plot uses the correct coordinate ranges
+  plot_data <- plot$data
+  expect_true(all(plot_data$row %in% 1:3))  # Should be converted to 1,2,3
+  expect_true(all(plot_data$col %in% 1:3))  # Should be converted to 1,2,3
+})
+
+test_that("autoplot handles mixed factor and numeric columns", {
+  # Create a design with mixed factor and numeric spatial columns
+  df <- data.frame(
+    row = factor(rep(1:4, each = 2), levels = 1:4),  # Factor
+    col = rep(1:2, times = 4),                       # Numeric
+    treatment = rep(LETTERS[1:4], 2)
+  )
+
+  # Verify column types
+  expect_true(is.factor(df$row))
+  expect_true(is.numeric(df$col))
+
+  # Set seed for reproducibility
+  set.seed(42)
+
+  # Run speed optimization
+  result <- speed(df, "treatment", iterations = 100)
+
+  # Test that autoplot works with mixed column types
+  expect_no_error(plot <- autoplot(result))
+  expect_s3_class(plot, "ggplot")
+
+  # Check that both columns are numeric in the plot data
+  plot_data <- plot$data
+  expect_true(is.numeric(plot_data$row))
+  expect_true(is.numeric(plot_data$col))
+})
+
+test_that("autoplot handles factor columns with blocks", {
+  # Create a design with factor columns and blocks
+  df <- data.frame(
+    row = factor(rep(1:6, each = 4), levels = 1:6),
+    col = factor(rep(1:4, times = 6), levels = 1:4),
+    treatment = rep(LETTERS[1:8], 3),
+    block = rep(1:3, each = 8)
+  )
+
+  # Verify columns are factors
+  expect_true(is.factor(df$row))
+  expect_true(is.factor(df$col))
+
+  # Set seed for reproducibility
+  set.seed(42)
+
+  # Run speed optimization with blocks
+  result <- speed(df, "treatment", swap_within = "block", iterations = 100)
+
+  # Test that autoplot works with factor columns and blocks
+  expect_no_error(plot <- autoplot(result))
+  expect_s3_class(plot, "ggplot")
+
+  # Check that the plot data has numeric coordinates
+  plot_data <- plot$data
+  expect_true(is.numeric(plot_data$row))
+  expect_true(is.numeric(plot_data$col))
+})
+
+test_that("autoplot handles character columns converted to factors", {
+  # Create a design with character spatial columns that will be converted to factors
+  df <- data.frame(
+    row = rep(c("A", "B", "C"), each = 3),
+    col = rep(c("X", "Y", "Z"), times = 3),
+    treatment = rep(LETTERS[1:3], 3),
+    stringsAsFactors = FALSE
+  )
+
+  # Convert to factors manually to simulate user input
+  df$row <- factor(df$row, levels = c("A", "B", "C"))
+  df$col <- factor(df$col, levels = c("X", "Y", "Z"))
+
+  # Verify columns are factors
+  expect_true(is.factor(df$row))
+  expect_true(is.factor(df$col))
+
+  # Set seed for reproducibility
+  set.seed(42)
+
+  # Run speed optimization
+  result <- speed(df, "treatment", iterations = 100)
+
+  # Test that autoplot works with character-based factor levels
+  expect_no_error(plot <- autoplot(result))
+  expect_s3_class(plot, "ggplot")
+
+  # Check that the plot data has numeric coordinates
+  plot_data <- plot$data
+  expect_true(is.numeric(plot_data$row))
+  expect_true(is.numeric(plot_data$col))
+})
+
+test_that("autoplot preserves factor order when converting to numeric", {
+  # Create a design with factors in non-alphabetical order
+  df <- data.frame(
+    row = factor(rep(c("Third", "First", "Second"), each = 2),
+                 levels = c("First", "Second", "Third")),
+    col = factor(rep(c("B", "A"), times = 3),
+                 levels = c("A", "B")),
+    treatment = rep(LETTERS[1:3], 2)
+  )
+
+  # Set seed for reproducibility
+  set.seed(42)
+
+  # Run speed optimization
+  result <- speed(df, "treatment", iterations = 100)
+
+  # Test that autoplot works and preserves the factor level ordering
+  expect_no_error(plot <- autoplot(result))
+  expect_s3_class(plot, "ggplot")
+
+  # Check conversion preserves order
+  plot_data <- plot$data
+  expect_true(is.numeric(plot_data$row))
+  expect_true(is.numeric(plot_data$col))
+
+  # Verify the numeric values correspond to factor level positions
+  # "First" should be 1, "Second" should be 2, "Third" should be 3
+  # "A" should be 1, "B" should be 2
+  unique_rows <- unique(plot_data$row)
+  unique_cols <- unique(plot_data$col)
+  expect_true(all(unique_rows %in% 1:3))
+  expect_true(all(unique_cols %in% 1:2))
+})
+
+test_that("autoplot correctly converts factor labels to numeric values", {
+  # Create a design with factor columns that have numeric labels
+  # This tests that as.numeric(as.character()) is used instead of as.numeric()
+  df <- data.frame(
+    row = factor(rep(c("10", "20", "30"), each = 3), levels = c("10", "20", "30")),
+    col = factor(rep(c("5", "15", "25"), times = 3), levels = c("5", "15", "25")),
+    treatment = rep(LETTERS[1:3], 3)
+  )
+
+  # Create a design object directly (skip speed optimization for this test)
+  result <- structure(df, class = "design")
+
+  # Test that autoplot works and converts factor labels correctly
+  expect_no_error(plot <- autoplot(result))
+  expect_s3_class(plot, "ggplot")
+
+  # The plot data should have numeric values from the factor labels, not levels
+  plot_data <- plot$data
+  expect_true(all(plot_data$row %in% c(10, 20, 30)))
+  expect_true(all(plot_data$col %in% c(5, 15, 25)))
+
+  # Should NOT be the factor levels (1, 2, 3)
+  expect_false(any(plot_data$row %in% c(1, 2, 3)))
+  expect_false(any(plot_data$col %in% c(1, 2, 3)))
+})

From d746ca6f8bcc39f659785cc917277797ce4f86a7 Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Mon, 11 Aug 2025 09:55:31 +0930
Subject: [PATCH 2/9] First go at an incidence function

---
 R/incidence.R | 126 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 126 insertions(+)
 create mode 100644 R/incidence.R

diff --git a/R/incidence.R b/R/incidence.R
new file mode 100644
index 0000000..bb07e16
--- /dev/null
+++ b/R/incidence.R
@@ -0,0 +1,126 @@
+calculate_incidence <- function(df, row_col = "row", col_col = "col", var_col = "treatment") {
+
+  # Input validation
+  if (!all(c(row_col, col_col, var_col) %in% names(df))) {
+    stop("Specified columns not found in dataframe")
+  }
+
+  # Create contingency tables
+  row_incidence <- table(df[[row_col]], df[[var_col]])
+  col_incidence <- table(df[[col_col]], df[[var_col]])
+
+  # Convert to data frames for easier manipulation
+  row_df <- as.data.frame.matrix(row_incidence)
+  col_df <- as.data.frame.matrix(col_incidence)
+
+  # Add row/col identifiers
+  row_df$position <- as.numeric(rownames(row_df))
+  col_df$position <- as.numeric(rownames(col_df))
+
+  # Reshape to long format for easier analysis
+  row_long <- reshape2::melt(row_df, id.vars = "position",
+                            variable.name = "variable",
+                            value.name = "count")
+  row_long$type <- "row"
+
+  col_long <- reshape2::melt(col_df, id.vars = "position",
+                            variable.name = "variable",
+                            value.name = "count")
+  col_long$type <- "col"
+
+  # Combine results
+  result <- rbind(row_long, col_long)
+
+  # Return list with multiple formats
+  return(list(
+    row_incidence = row_incidence,
+    col_incidence = col_incidence,
+    row_matrix = row_df,
+    col_matrix = col_df,
+    combined_long = result,
+    summary = list(
+      total_by_row = rowSums(row_incidence),
+      total_by_col = rowSums(col_incidence),
+      total_by_variable = colSums(row_incidence)
+    )
+  ))
+}
+
+# Alternative function for proportions instead of counts
+calculate_incidence_prop <- function(df, row_col = "row", col_col = "col", var_col = "treatment") {
+  result <- calculate_incidence(df, row_col, col_col, var_col)
+
+  # Convert counts to proportions
+  row_prop <- prop.table(result$row_incidence, margin = 1)
+  col_prop <- prop.table(result$col_incidence, margin = 1)
+
+  return(list(
+    row_proportions = row_prop,
+    col_proportions = col_prop,
+    row_counts = result$row_incidence,
+    col_counts = result$col_incidence
+  ))
+}
+
+# # Example usage with your data
+# df <- data.frame(
+#   row = rep(1:10, times = 6),
+#   col = rep(1:6, each = 10),
+#   treatment = rep(LETTERS[1:10], 6)
+# )
+#
+# # Calculate incidence
+# incidence_result <- calculate_incidence(df)
+#
+# # Display results
+# print("Row incidence (treatments by row):")
+# print(incidence_result$row_incidence)
+#
+# print("\nColumn incidence (treatments by column):")
+# print(incidence_result$col_incidence)
+#
+# print("\nSummary statistics:")
+# print(incidence_result$summary)
+#
+# # Calculate proportions
+# prop_result <- calculate_incidence_prop(df)
+# print("\nRow proportions:")
+# print(prop_result$row_proportions)
+
+# Visualization function (requires ggplot2)
+plot_incidence <- function(incidence_result, type = "both") {
+  if (!require(ggplot2, quietly = TRUE)) {
+    stop("ggplot2 package required for plotting")
+  }
+
+  library(ggplot2)
+
+  if (type == "row" || type == "both") {
+    row_plot <- ggplot(incidence_result$combined_long[incidence_result$combined_long$type == "row",],
+                       aes(x = position, y = count, fill = variable)) +
+      geom_bar(stat = "identity", position = "stack") +
+      labs(title = "Variable Incidence by Row", x = "Row", y = "Count") +
+      theme_minimal()
+
+    if (type == "row") return(row_plot)
+  }
+
+  if (type == "col" || type == "both") {
+    col_plot <- ggplot(incidence_result$combined_long[incidence_result$combined_long$type == "col",],
+                       aes(x = position, y = count, fill = variable)) +
+      geom_bar(stat = "identity", position = "stack") +
+      labs(title = "Variable Incidence by Column", x = "Column", y = "Count") +
+      theme_minimal()
+
+    if (type == "col") return(col_plot)
+  }
+
+  if (type == "both") {
+    if (require(gridExtra, quietly = TRUE)) {
+      return(gridExtra::grid.arrange(row_plot, col_plot, ncol = 1))
+    } else {
+      print(row_plot)
+      print(col_plot)
+    }
+  }
+}

From c757566218a2add8145321b2cc91d5bb473f6f6f Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Wed, 10 Jun 2026 11:13:57 +0930
Subject: [PATCH 3/9] Updated instructions for NEWS entries.

---
 CLAUDE.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CLAUDE.md b/CLAUDE.md
index 4d256d4..c955963 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -56,7 +56,7 @@ quarto::quarto_render("vignettes/speed.qmd")  # render a vignette (Quarto, not k
 Code formatting is handled by [Air](https://posit-dev.github.io/air/) - see `air.toml` (80-col, 2-space indent).
 Per `CONTRIBUTING.md`, do **not** restyle code that is unrelated to your PR.
 
-User-facing changes should add a bullet to the top of `NEWS.md`.
+User-facing changes should add a bullet to the top of `NEWS.md`. Use headings Major Changes, Minor Changes and Bug Fixes to document changes appropriately. Link to an issue number or PR such as (#999) where relevant.
 
 ## Architecture
 

From 813766eab4b6394951f511d3a39f0f0fb6812d3d Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Wed, 10 Jun 2026 11:14:18 +0930
Subject: [PATCH 4/9] Ignoring more files

---
 .gitignore | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/.gitignore b/.gitignore
index 6a2325c..5b8c1ea 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,8 @@
 .Rproj.user
 .Rhistory
 .RData
+*\.RData
+*\.rds
 .Ruserdata
 lib/
 docs
@@ -12,3 +14,8 @@ inst/doc
 
 Rplots.pdf
 *.local.json
+.DS_Store
+Thumbs.db
+*\.txt
+*\.xlsx
+*\.csv

From a4eff499054356c1638050e5c1cf99c8e956123d Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Wed, 10 Jun 2026 11:15:35 +0930
Subject: [PATCH 5/9] Adding incidence functions

---
 DESCRIPTION                          |   2 +-
 NAMESPACE                            |   2 +
 NEWS.md                              |   7 +
 R/calculate_adjacency_score.R        |   8 +-
 R/design_utils.R                     |  28 +++
 R/incidence.R                        | 274 ++++++++++++++++-----------
 R/metrics.R                          |   7 +-
 man/build_design_matrix.Rd           |  28 +++
 man/calculate_pair_incidence.Rd      |  66 +++++++
 man/calculate_position_incidence.Rd  |  54 ++++++
 tests/testthat/test-pair-incidence.R | 114 +++++++++++
 tests/testthat/test-pos-incidence.R  | 103 ++++++++++
 12 files changed, 572 insertions(+), 121 deletions(-)
 create mode 100644 man/build_design_matrix.Rd
 create mode 100644 man/calculate_pair_incidence.Rd
 create mode 100644 man/calculate_position_incidence.Rd
 create mode 100644 tests/testthat/test-pair-incidence.R
 create mode 100644 tests/testthat/test-pos-incidence.R

diff --git a/DESCRIPTION b/DESCRIPTION
index 4365616..6d631fa 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: speed
 Title: Generate Spatially Efficient Experimental Designs
-Version: 0.0.8
+Version: 0.0.9
 Authors@R: 
     c(person(given = "Sam",
              family = "Rogers",
diff --git a/NAMESPACE b/NAMESPACE
index 923ae6e..d8a936c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,6 +9,8 @@ export(calculate_balance_score)
 export(calculate_ed)
 export(calculate_efficiency_factor)
 export(calculate_nb)
+export(calculate_pair_incidence)
+export(calculate_position_incidence)
 export(create_pair_mapping)
 export(get_edges)
 export(get_vertices)
diff --git a/NEWS.md b/NEWS.md
index b62a504..7fe54e8 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,10 @@
+# speed 0.0.9
+
+## Major Changes
+
+- Added `calculate_pair_incidence()` to return a symmetric treatment × treatment matrix of neighbour-pair counts.
+- Added `calculate_position_incidence()` to return treatment × row/column position count matrices.
+
 # speed 0.0.8
 
 ## Major Changes
diff --git a/R/calculate_adjacency_score.R b/R/calculate_adjacency_score.R
index c6bdb28..9bc3975 100755
--- a/R/calculate_adjacency_score.R
+++ b/R/calculate_adjacency_score.R
@@ -247,11 +247,9 @@ calculate_adjacency_score <- function(layout_df,
                                       relationship = NULL) {
   ring_type <- match.arg(ring_type)
 
-  design_matrix <- matrix(
-    layout_df[[swap]],
-    nrow = max(as_numeric_factor(layout_df[[row_column]]), na.rm = TRUE),
-    ncol = max(as_numeric_factor(layout_df[[col_column]]), na.rm = TRUE),
-    byrow = TRUE
+  design_matrix <- build_design_matrix(
+    layout_df, swap,
+    row_col = row_column, col_col = col_column
   )
 
   per_cell <- adjacency_score_vec(
diff --git a/R/design_utils.R b/R/design_utils.R
index 8e712f7..98ab9e0 100644
--- a/R/design_utils.R
+++ b/R/design_utils.R
@@ -623,3 +623,31 @@ random_initialise <- function(design, optimise, seed = NULL, ...) {
 #' @rdname initialise_design_df
 #' @export
 initialize_design_df <- initialise_design_df
+
+#' Build a Spatial Design Matrix from a Data Frame
+#'
+#' @description
+#' Places each treatment value at the grid position given by its `row_col` and
+#' `col_col` coordinates, returning a character matrix of dimensions
+#' `max(row) × max(col)`. Cells with no corresponding row in `df` are `NA`.
+#' Unlike filling via `matrix(..., byrow)`, this is robust to any row ordering
+#' of `df`.
+#'
+#' @param df A data frame with columns named by `swap`, `row_col`, `col_col`.
+#' @param swap Column name of the treatment variable.
+#' @param row_col Column name of the row position variable (default `"row"`).
+#' @param col_col Column name of the column position variable (default `"col"`).
+#'
+#' @return A character matrix of dimensions `max(row) × max(col)`.
+#'
+#' @keywords internal
+build_design_matrix <- function(df, swap, row_col = "row", col_col = "col") {
+  nr <- max(as_numeric_factor(df[[row_col]]), na.rm = TRUE)
+  nc <- max(as_numeric_factor(df[[col_col]]), na.rm = TRUE)
+  m <- matrix(NA_character_, nrow = nr, ncol = nc)
+  m[cbind(
+    as_numeric_factor(df[[row_col]]),
+    as_numeric_factor(df[[col_col]])
+  )] <- as.character(df[[swap]])
+  m
+}
diff --git a/R/incidence.R b/R/incidence.R
index bb07e16..e48c256 100644
--- a/R/incidence.R
+++ b/R/incidence.R
@@ -1,126 +1,178 @@
-calculate_incidence <- function(df, row_col = "row", col_col = "col", var_col = "treatment") {
-
-  # Input validation
-  if (!all(c(row_col, col_col, var_col) %in% names(df))) {
-    stop("Specified columns not found in dataframe")
+#' Calculate Treatment Pair Incidence
+#'
+#' @description
+#' Returns a symmetric matrix counting how many times each pair of treatments
+#' appears as neighbours in the design. Two plots are neighbours if they share
+#' an edge on the `row × col` grid (orthogonal / rook adjacency; diagonals
+#' excluded). Horizontal and vertical adjacencies are pooled.
+#'
+#' The diagonal holds self-adjacency counts — the quantity that
+#' [calculate_adjacency_score()] penalises. A well-optimised design will have
+#' zeros (or near-zeros) on the diagonal. Off-diagonal entry \eqn{(i, j)}
+#' is the total number of edges where treatment \eqn{i} and treatment \eqn{j}
+#' are immediate neighbours.
+#'
+#' @param design A `design` object returned by [speed()], or a plain data
+#'   frame containing at least the columns named by `swap`, `row_col`, and
+#'   `col_col`.
+#' @param swap Column name of the treatment variable (default `"treatment"`).
+#' @param row_col Column name of the row position variable (default `"row"`).
+#' @param col_col Column name of the column position variable (default `"col"`).
+#' @param as_list If `TRUE`, returns a named list instead of a matrix: one
+#'   entry per treatment, each a named integer vector of neighbour counts with
+#'   every other treatment (including itself on the diagonal). This is the
+#'   matrix split by row — identical data, different container.
+#'
+#' @return When `as_list = FALSE` (default): a symmetric integer matrix of
+#'   dimension \eqn{t \times t} where \eqn{t} is the number of unique
+#'   treatments, with treatment names as row and column names.
+#'
+#'   When `as_list = TRUE`: a named list of length \eqn{t}, where each element
+#'   is a named integer vector of neighbour counts for that treatment.
+#'
+#' @examples
+#' # 3x3 Latin square — all pairs equally adjacent, zero self-adjacency
+#' df <- initialise_design_df(
+#'   items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
+#'   nrows = 3,
+#'   ncols = 3
+#' )
+#' calculate_pair_incidence(df)
+#'
+#' # as named list
+#' calculate_pair_incidence(df, as_list = TRUE)
+#'
+#' @seealso [calculate_adjacency_score()], [calculate_position_incidence()]
+#'
+#' @export
+calculate_pair_incidence <- function(design,
+                                     swap = "treatment",
+                                     row_col = "row",
+                                     col_col = "col",
+                                     as_list = FALSE) {
+  df <- if (inherits(design, "design")) design$design_df else design
+  missing_cols <- setdiff(c(swap, row_col, col_col), names(df))
+  if (length(missing_cols)) {
+    stop(
+      "Column(s) not found in data: ",
+      paste(missing_cols, collapse = ", "),
+      call. = FALSE
+    )
   }
 
-  # Create contingency tables
-  row_incidence <- table(df[[row_col]], df[[var_col]])
-  col_incidence <- table(df[[col_col]], df[[var_col]])
-
-  # Convert to data frames for easier manipulation
-  row_df <- as.data.frame.matrix(row_incidence)
-  col_df <- as.data.frame.matrix(col_incidence)
-
-  # Add row/col identifiers
-  row_df$position <- as.numeric(rownames(row_df))
-  col_df$position <- as.numeric(rownames(col_df))
-
-  # Reshape to long format for easier analysis
-  row_long <- reshape2::melt(row_df, id.vars = "position",
-                            variable.name = "variable",
-                            value.name = "count")
-  row_long$type <- "row"
+  m <- build_design_matrix(df, swap, row_col, col_col)
+  nr <- nrow(m)
+  nc <- ncol(m)
 
-  col_long <- reshape2::melt(col_df, id.vars = "position",
-                            variable.name = "variable",
-                            value.name = "count")
-  col_long$type <- "col"
+  lvls <- stringi::stri_sort(
+    unique(as.character(df[[swap]][!is.na(df[[swap]])])),
+    numeric = TRUE
+  )
 
-  # Combine results
-  result <- rbind(row_long, col_long)
-
-  # Return list with multiple formats
-  return(list(
-    row_incidence = row_incidence,
-    col_incidence = col_incidence,
-    row_matrix = row_df,
-    col_matrix = col_df,
-    combined_long = result,
-    summary = list(
-      total_by_row = rowSums(row_incidence),
-      total_by_col = rowSums(col_incidence),
-      total_by_variable = colSums(row_incidence)
+  horiz <- if (nc >= 2) {
+    cbind(
+      as.vector(m[, seq_len(nc - 1L)]),
+      as.vector(m[, seq.int(2L, nc)])
     )
-  ))
-}
-
-# Alternative function for proportions instead of counts
-calculate_incidence_prop <- function(df, row_col = "row", col_col = "col", var_col = "treatment") {
-  result <- calculate_incidence(df, row_col, col_col, var_col)
-
-  # Convert counts to proportions
-  row_prop <- prop.table(result$row_incidence, margin = 1)
-  col_prop <- prop.table(result$col_incidence, margin = 1)
-
-  return(list(
-    row_proportions = row_prop,
-    col_proportions = col_prop,
-    row_counts = result$row_incidence,
-    col_counts = result$col_incidence
-  ))
-}
-
-# # Example usage with your data
-# df <- data.frame(
-#   row = rep(1:10, times = 6),
-#   col = rep(1:6, each = 10),
-#   treatment = rep(LETTERS[1:10], 6)
-# )
-#
-# # Calculate incidence
-# incidence_result <- calculate_incidence(df)
-#
-# # Display results
-# print("Row incidence (treatments by row):")
-# print(incidence_result$row_incidence)
-#
-# print("\nColumn incidence (treatments by column):")
-# print(incidence_result$col_incidence)
-#
-# print("\nSummary statistics:")
-# print(incidence_result$summary)
-#
-# # Calculate proportions
-# prop_result <- calculate_incidence_prop(df)
-# print("\nRow proportions:")
-# print(prop_result$row_proportions)
-
-# Visualization function (requires ggplot2)
-plot_incidence <- function(incidence_result, type = "both") {
-  if (!require(ggplot2, quietly = TRUE)) {
-    stop("ggplot2 package required for plotting")
+  } else {
+    matrix(character(0L), nrow = 0L, ncol = 2L)
+  }
+  vert <- if (nr >= 2) {
+    cbind(
+      as.vector(m[seq_len(nr - 1L), ]),
+      as.vector(m[seq.int(2L, nr), ])
+    )
+  } else {
+    matrix(character(0L), nrow = 0L, ncol = 2L)
   }
 
-  library(ggplot2)
+  edges <- rbind(horiz, vert)
+  keep <- !is.na(edges[, 1L]) & !is.na(edges[, 2L])
+  edges <- edges[keep, , drop = FALSE]
 
-  if (type == "row" || type == "both") {
-    row_plot <- ggplot(incidence_result$combined_long[incidence_result$combined_long$type == "row",],
-                       aes(x = position, y = count, fill = variable)) +
-      geom_bar(stat = "identity", position = "stack") +
-      labs(title = "Variable Incidence by Row", x = "Row", y = "Count") +
-      theme_minimal()
+  a_fac <- factor(edges[, 1L], levels = lvls)
+  b_fac <- factor(edges[, 2L], levels = lvls)
+  raw <- matrix(
+    as.integer(table(a_fac, b_fac)),
+    nrow = length(lvls),
+    dimnames = list(lvls, lvls)
+  )
 
-    if (type == "row") return(row_plot)
+  M <- raw + t(raw)
+  diag(M) <- diag(raw)
+
+  if (as_list) {
+    return(lapply(setNames(lvls, lvls), function(trt) M[trt, ]))
   }
+  M
+}
 
-  if (type == "col" || type == "both") {
-    col_plot <- ggplot(incidence_result$combined_long[incidence_result$combined_long$type == "col",],
-                       aes(x = position, y = count, fill = variable)) +
-      geom_bar(stat = "identity", position = "stack") +
-      labs(title = "Variable Incidence by Column", x = "Column", y = "Count") +
-      theme_minimal()
 
-    if (type == "col") return(col_plot)
+#' Calculate Positional Incidence of Treatments
+#'
+#' @description
+#' Returns how many times each treatment appears in each row position and each
+#' column position of the design. This is the human-readable decomposition of
+#' what [calculate_balance_score()] collapses to a scalar: uneven row or column
+#' counts reveal where spatial imbalance lies.
+#'
+#' @inheritParams calculate_pair_incidence
+#'
+#' @return A named list with two integer matrices, both with treatments as rows
+#'   and positions as columns:
+#'   \describe{
+#'     \item{`row`}{Treatment × row-position count matrix.}
+#'     \item{`col`}{Treatment × column-position count matrix.}
+#'   }
+#'   Indexing by treatment name is natural: `result$row["A", ]` gives the row
+#'   distribution of treatment A.
+#'
+#' @examples
+#' df <- initialise_design_df(
+#'   items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
+#'   nrows = 3,
+#'   ncols = 3
+#' )
+#' res <- calculate_position_incidence(df)
+#' res$row  # treatment × row counts
+#' res$col  # treatment × column counts
+#'
+#' @seealso [calculate_balance_score()], [calculate_pair_incidence()]
+#'
+#' @export
+calculate_position_incidence <- function(design,
+                                         swap = "treatment",
+                                         row_col = "row",
+                                         col_col = "col") {
+  df <- if (inherits(design, "design")) design$design_df else design
+  missing_cols <- setdiff(c(swap, row_col, col_col), names(df))
+  if (length(missing_cols)) {
+    stop(
+      "Column(s) not found in data: ",
+      paste(missing_cols, collapse = ", "),
+      call. = FALSE
+    )
   }
 
-  if (type == "both") {
-    if (require(gridExtra, quietly = TRUE)) {
-      return(gridExtra::grid.arrange(row_plot, col_plot, ncol = 1))
-    } else {
-      print(row_plot)
-      print(col_plot)
-    }
-  }
+  lvls <- stringi::stri_sort(
+    unique(as.character(df[[swap]][!is.na(df[[swap]])])),
+    numeric = TRUE
+  )
+  trt_fac <- factor(df[[swap]], levels = lvls)
+
+  row_tbl <- table(trt_fac, df[[row_col]])
+  col_tbl <- table(trt_fac, df[[col_col]])
+
+  row_mat <- matrix(
+    as.integer(row_tbl),
+    nrow = length(lvls),
+    dimnames = list(lvls, colnames(row_tbl))
+  )
+  col_mat <- matrix(
+    as.integer(col_tbl),
+    nrow = length(lvls),
+    dimnames = list(lvls, colnames(col_tbl))
+  )
+
+  list(row = row_mat, col = col_mat)
 }
diff --git a/R/metrics.R b/R/metrics.R
index c356786..47a8acf 100644
--- a/R/metrics.R
+++ b/R/metrics.R
@@ -219,10 +219,9 @@ objective_function_piepho <- function(design,
                                       row_column = "row",
                                       col_column = "col",
                                       ...) {
-  design_matrix <- matrix(
-    design[[swap]],
-    nrow = max(as_numeric_factor(design[[row_column]]), na.rm = TRUE),
-    ncol = max(as_numeric_factor(design[[col_column]]), na.rm = TRUE)
+  design_matrix <- build_design_matrix(
+    design, swap,
+    row_col = row_column, col_col = col_column
   )
 
   ed <- calculate_ed(design_matrix, current_score_obj$ed, swapped_items)
diff --git a/man/build_design_matrix.Rd b/man/build_design_matrix.Rd
new file mode 100644
index 0000000..469d041
--- /dev/null
+++ b/man/build_design_matrix.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/design_utils.R
+\name{build_design_matrix}
+\alias{build_design_matrix}
+\title{Build a Spatial Design Matrix from a Data Frame}
+\usage{
+build_design_matrix(df, swap, row_col = "row", col_col = "col")
+}
+\arguments{
+\item{df}{A data frame with columns named by \code{swap}, \code{row_col}, \code{col_col}.}
+
+\item{swap}{Column name of the treatment variable.}
+
+\item{row_col}{Column name of the row position variable (default \code{"row"}).}
+
+\item{col_col}{Column name of the column position variable (default \code{"col"}).}
+}
+\value{
+A character matrix of dimensions \verb{max(row) × max(col)}.
+}
+\description{
+Places each treatment value at the grid position given by its \code{row_col} and
+\code{col_col} coordinates, returning a character matrix of dimensions
+\verb{max(row) × max(col)}. Cells with no corresponding row in \code{df} are \code{NA}.
+Unlike filling via \code{matrix(..., byrow)}, this is robust to any row ordering
+of \code{df}.
+}
+\keyword{internal}
diff --git a/man/calculate_pair_incidence.Rd b/man/calculate_pair_incidence.Rd
new file mode 100644
index 0000000..86ff959
--- /dev/null
+++ b/man/calculate_pair_incidence.Rd
@@ -0,0 +1,66 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/incidence.R
+\name{calculate_pair_incidence}
+\alias{calculate_pair_incidence}
+\title{Calculate Treatment Pair Incidence}
+\usage{
+calculate_pair_incidence(
+  design,
+  swap = "treatment",
+  row_col = "row",
+  col_col = "col",
+  as_list = FALSE
+)
+}
+\arguments{
+\item{design}{A \code{design} object returned by \code{\link[=speed]{speed()}}, or a plain data
+frame containing at least the columns named by \code{swap}, \code{row_col}, and
+\code{col_col}.}
+
+\item{swap}{Column name of the treatment variable (default \code{"treatment"}).}
+
+\item{row_col}{Column name of the row position variable (default \code{"row"}).}
+
+\item{col_col}{Column name of the column position variable (default \code{"col"}).}
+
+\item{as_list}{If \code{TRUE}, returns a named list instead of a matrix: one
+entry per treatment, each a named integer vector of neighbour counts with
+every other treatment (including itself on the diagonal). This is the
+matrix split by row — identical data, different container.}
+}
+\value{
+When \code{as_list = FALSE} (default): a symmetric integer matrix of
+dimension \eqn{t \times t} where \eqn{t} is the number of unique
+treatments, with treatment names as row and column names.
+
+When \code{as_list = TRUE}: a named list of length \eqn{t}, where each element
+is a named integer vector of neighbour counts for that treatment.
+}
+\description{
+Returns a symmetric matrix counting how many times each pair of treatments
+appears as neighbours in the design. Two plots are neighbours if they share
+an edge on the \verb{row × col} grid (orthogonal / rook adjacency; diagonals
+excluded). Horizontal and vertical adjacencies are pooled.
+
+The diagonal holds self-adjacency counts — the quantity that
+\code{\link[=calculate_adjacency_score]{calculate_adjacency_score()}} penalises. A well-optimised design will have
+zeros (or near-zeros) on the diagonal. Off-diagonal entry \eqn{(i, j)}
+is the total number of edges where treatment \eqn{i} and treatment \eqn{j}
+are immediate neighbours.
+}
+\examples{
+# 3x3 Latin square — all pairs equally adjacent, zero self-adjacency
+df <- initialise_design_df(
+  items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
+  nrows = 3,
+  ncols = 3
+)
+calculate_pair_incidence(df)
+
+# as named list
+calculate_pair_incidence(df, as_list = TRUE)
+
+}
+\seealso{
+\code{\link[=calculate_adjacency_score]{calculate_adjacency_score()}}, \code{\link[=calculate_position_incidence]{calculate_position_incidence()}}
+}
diff --git a/man/calculate_position_incidence.Rd b/man/calculate_position_incidence.Rd
new file mode 100644
index 0000000..e19714b
--- /dev/null
+++ b/man/calculate_position_incidence.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/incidence.R
+\name{calculate_position_incidence}
+\alias{calculate_position_incidence}
+\title{Calculate Positional Incidence of Treatments}
+\usage{
+calculate_position_incidence(
+  design,
+  swap = "treatment",
+  row_col = "row",
+  col_col = "col"
+)
+}
+\arguments{
+\item{design}{A \code{design} object returned by \code{\link[=speed]{speed()}}, or a plain data
+frame containing at least the columns named by \code{swap}, \code{row_col}, and
+\code{col_col}.}
+
+\item{swap}{Column name of the treatment variable (default \code{"treatment"}).}
+
+\item{row_col}{Column name of the row position variable (default \code{"row"}).}
+
+\item{col_col}{Column name of the column position variable (default \code{"col"}).}
+}
+\value{
+A named list with two integer matrices, both with treatments as rows
+and positions as columns:
+\describe{
+\item{\code{row}}{Treatment × row-position count matrix.}
+\item{\code{col}}{Treatment × column-position count matrix.}
+}
+Indexing by treatment name is natural: \code{result$row["A", ]} gives the row
+distribution of treatment A.
+}
+\description{
+Returns how many times each treatment appears in each row position and each
+column position of the design. This is the human-readable decomposition of
+what \code{\link[=calculate_balance_score]{calculate_balance_score()}} collapses to a scalar: uneven row or column
+counts reveal where spatial imbalance lies.
+}
+\examples{
+df <- initialise_design_df(
+  items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
+  nrows = 3,
+  ncols = 3
+)
+res <- calculate_position_incidence(df)
+res$row  # treatment × row counts
+res$col  # treatment × column counts
+
+}
+\seealso{
+\code{\link[=calculate_balance_score]{calculate_balance_score()}}, \code{\link[=calculate_pair_incidence]{calculate_pair_incidence()}}
+}
diff --git a/tests/testthat/test-pair-incidence.R b/tests/testthat/test-pair-incidence.R
new file mode 100644
index 0000000..81d2c72
--- /dev/null
+++ b/tests/testthat/test-pair-incidence.R
@@ -0,0 +1,114 @@
+test_that("3x3 Latin square: symmetric matrix, zero diagonal, all pairs = 4", {
+  # Grid (column-major from initialise_design_df):
+  #   col 1  col 2  col 3
+  # row 1: A      B      C
+  # row 2: B      C      A
+  # row 3: C      A      B
+  df <- data.frame(
+    row = c(1, 2, 3, 1, 2, 3, 1, 2, 3),
+    col = c(1, 1, 1, 2, 2, 2, 3, 3, 3),
+    trt = c("A", "B", "C", "B", "C", "A", "C", "A", "B")
+  )
+  M <- calculate_pair_incidence(df, swap = "trt")
+
+  expect_true(isSymmetric(M))
+  expect_equal(diag(M), c(A = 0L, B = 0L, C = 0L))
+  expect_true(all(M[upper.tri(M)] == 4L))
+})
+
+test_that("as_list rows match matrix rows", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "B", "A")
+  )
+  M <- calculate_pair_incidence(df, swap = "trt")
+  lst <- calculate_pair_incidence(df, swap = "trt", as_list = TRUE)
+
+  expect_equal(lst[["A"]], M["A", ])
+  expect_equal(lst[["B"]], M["B", ])
+})
+
+test_that("self-adjacency recorded on diagonal", {
+  # 1x2 grid with two identical treatments: one self-adjacent edge
+  df <- data.frame(row = c(1, 1), col = c(1, 2), trt = c("A", "A"))
+  M <- calculate_pair_incidence(df, swap = "trt")
+
+  expect_equal(M["A", "A"], 1L)
+})
+
+test_that("checkerboard 2x2: zero diagonal, four A-B edges", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "B", "A")
+  )
+  M <- calculate_pair_incidence(df, swap = "trt")
+
+  expect_equal(diag(M), c(A = 0L, B = 0L))
+  expect_equal(M["A", "B"], 4L)
+})
+
+test_that("single row: only horizontal edges counted", {
+  df <- data.frame(row = c(1, 1, 1), col = c(1, 2, 3), trt = c("A", "B", "A"))
+  M <- calculate_pair_incidence(df, swap = "trt")
+
+  # Horizontal edges: (A,B), (B,A) -> A-B = 2; no A-A edges
+  expect_equal(M["A", "B"], 2L)
+  expect_equal(M["A", "A"], 0L)
+})
+
+test_that("single column: only vertical edges counted", {
+  df <- data.frame(row = c(1, 2, 3), col = c(1, 1, 1), trt = c("A", "B", "A"))
+  M <- calculate_pair_incidence(df, swap = "trt")
+
+  # Vertical edges: (A,B), (B,A) -> A-B = 2
+  expect_equal(M["A", "B"], 2L)
+  expect_equal(M["A", "A"], 0L)
+})
+
+test_that("NA cells are skipped", {
+  # 2x2 grid with one NA cell
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", NA, "B", "A")
+  )
+  M <- calculate_pair_incidence(df, swap = "trt")
+
+  # Edges involving the NA cell are dropped.
+  # Remaining edges: horiz row2 (B,A); vert col1 (A,B); vert col2 (NA,A) dropped
+  # So A-B = 2
+  expect_equal(M["A", "B"], 2L)
+  expect_true(!any(is.na(M)))
+})
+
+test_that("design object accepted as input", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "B", "A")
+  )
+  fake_design <- structure(list(design_df = df), class = "design")
+  M_df <- calculate_pair_incidence(df, swap = "trt")
+  M_obj <- calculate_pair_incidence(fake_design, swap = "trt")
+
+  expect_equal(M_df, M_obj)
+})
+
+test_that("treatments never adjacent still appear as zero row/col", {
+  # 1x3 grid: A at ends, B in middle — A never adjacent to itself
+  df <- data.frame(row = c(1, 1, 1), col = c(1, 2, 3), trt = c("A", "B", "A"))
+  M <- calculate_pair_incidence(df, swap = "trt")
+
+  expect_true(all(c("A", "B") %in% rownames(M)))
+  expect_equal(M["A", "A"], 0L)
+})
+
+test_that("missing column gives informative error", {
+  df <- data.frame(row = 1:3, trt = letters[1:3])
+  expect_error(
+    calculate_pair_incidence(df, swap = "trt"),
+    "col"
+  )
+})
diff --git a/tests/testthat/test-pos-incidence.R b/tests/testthat/test-pos-incidence.R
new file mode 100644
index 0000000..c3cde43
--- /dev/null
+++ b/tests/testthat/test-pos-incidence.R
@@ -0,0 +1,103 @@
+test_that("position incidence returns correct row and col counts", {
+  # 2x3 grid:
+  #   col 1  col 2  col 3
+  # row 1: A      A      B
+  # row 2: A      B      B
+  df <- data.frame(
+    row = c(1, 1, 1, 2, 2, 2),
+    col = c(1, 2, 3, 1, 2, 3),
+    trt = c("A", "A", "B", "A", "B", "B")
+  )
+  res <- calculate_position_incidence(df, swap = "trt")
+
+  # row counts
+  expect_equal(res$row["A", "1"], 2L)
+  expect_equal(res$row["A", "2"], 1L)
+  expect_equal(res$row["B", "1"], 1L)
+  expect_equal(res$row["B", "2"], 2L)
+
+  # col counts
+  expect_equal(res$col["A", "1"], 2L)
+  expect_equal(res$col["A", "2"], 1L)
+  expect_equal(res$col["A", "3"], 0L)
+  expect_equal(res$col["B", "1"], 0L)
+  expect_equal(res$col["B", "2"], 1L)
+  expect_equal(res$col["B", "3"], 2L)
+})
+
+test_that("returns list with row and col matrices", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "A", "B")
+  )
+  res <- calculate_position_incidence(df, swap = "trt")
+
+  expect_type(res, "list")
+  expect_named(res, c("row", "col"))
+  expect_true(is.matrix(res$row))
+  expect_true(is.matrix(res$col))
+})
+
+test_that("treatments as rows, positions as columns", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "A", "B")
+  )
+  res <- calculate_position_incidence(df, swap = "trt")
+
+  expect_equal(rownames(res$row), c("A", "B"))
+  expect_equal(colnames(res$row), c("1", "2"))
+  expect_equal(rownames(res$col), c("A", "B"))
+  expect_equal(colnames(res$col), c("1", "2"))
+})
+
+test_that("balanced design has equal counts everywhere", {
+  # Each treatment appears once per row and once per column (Latin square)
+  df <- data.frame(
+    row = c(1, 2, 3, 1, 2, 3, 1, 2, 3),
+    col = c(1, 1, 1, 2, 2, 2, 3, 3, 3),
+    trt = c("A", "B", "C", "B", "C", "A", "C", "A", "B")
+  )
+  res <- calculate_position_incidence(df, swap = "trt")
+
+  expect_true(all(res$row == 1L))
+  expect_true(all(res$col == 1L))
+})
+
+test_that("NA cells excluded consistently", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", NA, "B", "A")
+  )
+  res <- calculate_position_incidence(df, swap = "trt")
+
+  # NA treatment is dropped; row/col counts reflect only non-NA cells
+  expect_equal(res$row["A", "1"], 1L)
+  expect_equal(res$row["A", "2"], 1L)
+  expect_equal(res$row["B", "1"], 0L)
+  expect_equal(res$row["B", "2"], 1L)
+})
+
+test_that("design object accepted as input", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "A", "B")
+  )
+  fake_design <- structure(list(design_df = df), class = "design")
+  res_df <- calculate_position_incidence(df, swap = "trt")
+  res_obj <- calculate_position_incidence(fake_design, swap = "trt")
+
+  expect_equal(res_df, res_obj)
+})
+
+test_that("missing column gives informative error", {
+  df <- data.frame(row = 1:3, trt = letters[1:3])
+  expect_error(
+    calculate_position_incidence(df, swap = "trt"),
+    "col"
+  )
+})

From b7777003c343df75002ad94524722ce6a6ee9e0a Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Wed, 10 Jun 2026 11:40:26 +0930
Subject: [PATCH 6/9] Adding examples from a speed call

---
 R/incidence.R                       | 56 ++++++++++++++++-------------
 man/calculate_pair_incidence.Rd     | 30 +++++++++-------
 man/calculate_position_incidence.Rd | 24 +++++++------
 3 files changed, 63 insertions(+), 47 deletions(-)

diff --git a/R/incidence.R b/R/incidence.R
index e48c256..e2e22eb 100644
--- a/R/incidence.R
+++ b/R/incidence.R
@@ -13,11 +13,11 @@
 #' are immediate neighbours.
 #'
 #' @param design A `design` object returned by [speed()], or a plain data
-#'   frame containing at least the columns named by `swap`, `row_col`, and
-#'   `col_col`.
+#'   frame containing at least the columns named by `swap`, `row_column`, and
+#'   `col_column`.
 #' @param swap Column name of the treatment variable (default `"treatment"`).
-#' @param row_col Column name of the row position variable (default `"row"`).
-#' @param col_col Column name of the column position variable (default `"col"`).
+#' @param row_column Column name of the row position variable (default `"row"`).
+#' @param col_column Column name of the column position variable (default `"col"`).
 #' @param as_list If `TRUE`, returns a named list instead of a matrix: one
 #'   entry per treatment, each a named integer vector of neighbour counts with
 #'   every other treatment (including itself on the diagonal). This is the
@@ -31,27 +31,33 @@
 #'   is a named integer vector of neighbour counts for that treatment.
 #'
 #' @examples
-#' # 3x3 Latin square — all pairs equally adjacent, zero self-adjacency
-#' df <- initialise_design_df(
+#' # From a speed() result
+#' df <- data.frame(
+#'   row = rep(1:4, times = 3),
+#'   col = rep(1:3, each = 4),
+#'   treatment = rep(c("A", "B", "C"), 4)
+#' )
+#' result <- speed(df, swap = "treatment", seed = 42, iterations = 200)
+#' calculate_pair_incidence(result)
+#'
+#' # From a plain data frame; as named list
+#' df2 <- initialise_design_df(
 #'   items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
 #'   nrows = 3,
 #'   ncols = 3
 #' )
-#' calculate_pair_incidence(df)
-#'
-#' # as named list
-#' calculate_pair_incidence(df, as_list = TRUE)
+#' calculate_pair_incidence(df2, as_list = TRUE)
 #'
 #' @seealso [calculate_adjacency_score()], [calculate_position_incidence()]
 #'
 #' @export
 calculate_pair_incidence <- function(design,
                                      swap = "treatment",
-                                     row_col = "row",
-                                     col_col = "col",
+                                     row_column = "row",
+                                     col_column = "col",
                                      as_list = FALSE) {
   df <- if (inherits(design, "design")) design$design_df else design
-  missing_cols <- setdiff(c(swap, row_col, col_col), names(df))
+  missing_cols <- setdiff(c(swap, row_column, col_column), names(df))
   if (length(missing_cols)) {
     stop(
       "Column(s) not found in data: ",
@@ -60,7 +66,7 @@ calculate_pair_incidence <- function(design,
     )
   }
 
-  m <- build_design_matrix(df, swap, row_col, col_col)
+  m <- build_design_matrix(df, swap, row_column, col_column)
   nr <- nrow(m)
   nc <- ncol(m)
 
@@ -128,12 +134,14 @@ calculate_pair_incidence <- function(design,
 #'   distribution of treatment A.
 #'
 #' @examples
-#' df <- initialise_design_df(
-#'   items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
-#'   nrows = 3,
-#'   ncols = 3
+#' # From a speed() result
+#' df <- data.frame(
+#'   row = rep(1:4, times = 3),
+#'   col = rep(1:3, each = 4),
+#'   treatment = rep(c("A", "B", "C"), 4)
 #' )
-#' res <- calculate_position_incidence(df)
+#' result <- speed(df, swap = "treatment", seed = 42, iterations = 200)
+#' res <- calculate_position_incidence(result)
 #' res$row  # treatment × row counts
 #' res$col  # treatment × column counts
 #'
@@ -142,10 +150,10 @@ calculate_pair_incidence <- function(design,
 #' @export
 calculate_position_incidence <- function(design,
                                          swap = "treatment",
-                                         row_col = "row",
-                                         col_col = "col") {
+                                         row_column = "row",
+                                         col_column = "col") {
   df <- if (inherits(design, "design")) design$design_df else design
-  missing_cols <- setdiff(c(swap, row_col, col_col), names(df))
+  missing_cols <- setdiff(c(swap, row_column, col_column), names(df))
   if (length(missing_cols)) {
     stop(
       "Column(s) not found in data: ",
@@ -160,8 +168,8 @@ calculate_position_incidence <- function(design,
   )
   trt_fac <- factor(df[[swap]], levels = lvls)
 
-  row_tbl <- table(trt_fac, df[[row_col]])
-  col_tbl <- table(trt_fac, df[[col_col]])
+  row_tbl <- table(trt_fac, df[[row_column]])
+  col_tbl <- table(trt_fac, df[[col_column]])
 
   row_mat <- matrix(
     as.integer(row_tbl),
diff --git a/man/calculate_pair_incidence.Rd b/man/calculate_pair_incidence.Rd
index 86ff959..49f28a3 100644
--- a/man/calculate_pair_incidence.Rd
+++ b/man/calculate_pair_incidence.Rd
@@ -7,21 +7,21 @@
 calculate_pair_incidence(
   design,
   swap = "treatment",
-  row_col = "row",
-  col_col = "col",
+  row_column = "row",
+  col_column = "col",
   as_list = FALSE
 )
 }
 \arguments{
 \item{design}{A \code{design} object returned by \code{\link[=speed]{speed()}}, or a plain data
-frame containing at least the columns named by \code{swap}, \code{row_col}, and
-\code{col_col}.}
+frame containing at least the columns named by \code{swap}, \code{row_column}, and
+\code{col_column}.}
 
 \item{swap}{Column name of the treatment variable (default \code{"treatment"}).}
 
-\item{row_col}{Column name of the row position variable (default \code{"row"}).}
+\item{row_column}{Column name of the row position variable (default \code{"row"}).}
 
-\item{col_col}{Column name of the column position variable (default \code{"col"}).}
+\item{col_column}{Column name of the column position variable (default \code{"col"}).}
 
 \item{as_list}{If \code{TRUE}, returns a named list instead of a matrix: one
 entry per treatment, each a named integer vector of neighbour counts with
@@ -49,16 +49,22 @@ is the total number of edges where treatment \eqn{i} and treatment \eqn{j}
 are immediate neighbours.
 }
 \examples{
-# 3x3 Latin square — all pairs equally adjacent, zero self-adjacency
-df <- initialise_design_df(
+# From a speed() result
+df <- data.frame(
+  row = rep(1:4, times = 3),
+  col = rep(1:3, each = 4),
+  treatment = rep(c("A", "B", "C"), 4)
+)
+result <- speed(df, swap = "treatment", seed = 42)
+calculate_pair_incidence(result)
+
+# From a plain data frame; as named list
+df2 <- initialise_design_df(
   items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
   nrows = 3,
   ncols = 3
 )
-calculate_pair_incidence(df)
-
-# as named list
-calculate_pair_incidence(df, as_list = TRUE)
+calculate_pair_incidence(df2, as_list = TRUE)
 
 }
 \seealso{
diff --git a/man/calculate_position_incidence.Rd b/man/calculate_position_incidence.Rd
index e19714b..1eeefd9 100644
--- a/man/calculate_position_incidence.Rd
+++ b/man/calculate_position_incidence.Rd
@@ -7,20 +7,20 @@
 calculate_position_incidence(
   design,
   swap = "treatment",
-  row_col = "row",
-  col_col = "col"
+  row_column = "row",
+  col_column = "col"
 )
 }
 \arguments{
 \item{design}{A \code{design} object returned by \code{\link[=speed]{speed()}}, or a plain data
-frame containing at least the columns named by \code{swap}, \code{row_col}, and
-\code{col_col}.}
+frame containing at least the columns named by \code{swap}, \code{row_column}, and
+\code{col_column}.}
 
 \item{swap}{Column name of the treatment variable (default \code{"treatment"}).}
 
-\item{row_col}{Column name of the row position variable (default \code{"row"}).}
+\item{row_column}{Column name of the row position variable (default \code{"row"}).}
 
-\item{col_col}{Column name of the column position variable (default \code{"col"}).}
+\item{col_column}{Column name of the column position variable (default \code{"col"}).}
 }
 \value{
 A named list with two integer matrices, both with treatments as rows
@@ -39,12 +39,14 @@ what \code{\link[=calculate_balance_score]{calculate_balance_score()}} collapses
 counts reveal where spatial imbalance lies.
 }
 \examples{
-df <- initialise_design_df(
-  items = c("A", "B", "C", "B", "C", "A", "C", "A", "B"),
-  nrows = 3,
-  ncols = 3
+# From a speed() result
+df <- data.frame(
+  row = rep(1:4, times = 3),
+  col = rep(1:3, each = 4),
+  treatment = rep(c("A", "B", "C"), 4)
 )
-res <- calculate_position_incidence(df)
+result <- speed(df, swap = "treatment", seed = 42)
+res <- calculate_position_incidence(result)
 res$row  # treatment × row counts
 res$col  # treatment × column counts
 

From 5c5fcee3ea5b94d2f62b1b7797d4246cff6c372b Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Wed, 10 Jun 2026 11:57:23 +0930
Subject: [PATCH 7/9] Fix argument names

---
 R/calculate_adjacency_score.R |  2 +-
 R/design_utils.R              | 22 ++++++++++++----------
 R/metrics.R                   |  2 +-
 3 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/R/calculate_adjacency_score.R b/R/calculate_adjacency_score.R
index 9bc3975..bcea070 100755
--- a/R/calculate_adjacency_score.R
+++ b/R/calculate_adjacency_score.R
@@ -249,7 +249,7 @@ calculate_adjacency_score <- function(layout_df,
 
   design_matrix <- build_design_matrix(
     layout_df, swap,
-    row_col = row_column, col_col = col_column
+    row_column = row_column, col_column = col_column
   )
 
   per_cell <- adjacency_score_vec(
diff --git a/R/design_utils.R b/R/design_utils.R
index 98ab9e0..5bdf469 100644
--- a/R/design_utils.R
+++ b/R/design_utils.R
@@ -627,27 +627,29 @@ initialize_design_df <- initialise_design_df
 #' Build a Spatial Design Matrix from a Data Frame
 #'
 #' @description
-#' Places each treatment value at the grid position given by its `row_col` and
-#' `col_col` coordinates, returning a character matrix of dimensions
+#' Places each treatment value at the grid position given by its `row_column`
+#' and `col_column` coordinates, returning a character matrix of dimensions
 #' `max(row) × max(col)`. Cells with no corresponding row in `df` are `NA`.
 #' Unlike filling via `matrix(..., byrow)`, this is robust to any row ordering
 #' of `df`.
 #'
-#' @param df A data frame with columns named by `swap`, `row_col`, `col_col`.
+#' @param df A data frame with columns named by `swap`, `row_column`,
+#'   `col_column`.
 #' @param swap Column name of the treatment variable.
-#' @param row_col Column name of the row position variable (default `"row"`).
-#' @param col_col Column name of the column position variable (default `"col"`).
+#' @param row_column Column name of the row position variable (default `"row"`).
+#' @param col_column Column name of the column position variable
+#'   (default `"col"`).
 #'
 #' @return A character matrix of dimensions `max(row) × max(col)`.
 #'
 #' @keywords internal
-build_design_matrix <- function(df, swap, row_col = "row", col_col = "col") {
-  nr <- max(as_numeric_factor(df[[row_col]]), na.rm = TRUE)
-  nc <- max(as_numeric_factor(df[[col_col]]), na.rm = TRUE)
+build_design_matrix <- function(df, swap, row_column = "row", col_column = "col") {
+  nr <- max(as_numeric_factor(df[[row_column]]), na.rm = TRUE)
+  nc <- max(as_numeric_factor(df[[col_column]]), na.rm = TRUE)
   m <- matrix(NA_character_, nrow = nr, ncol = nc)
   m[cbind(
-    as_numeric_factor(df[[row_col]]),
-    as_numeric_factor(df[[col_col]])
+    as_numeric_factor(df[[row_column]]),
+    as_numeric_factor(df[[col_column]])
   )] <- as.character(df[[swap]])
   m
 }
diff --git a/R/metrics.R b/R/metrics.R
index 47a8acf..7604410 100644
--- a/R/metrics.R
+++ b/R/metrics.R
@@ -221,7 +221,7 @@ objective_function_piepho <- function(design,
                                       ...) {
   design_matrix <- build_design_matrix(
     design, swap,
-    row_col = row_column, col_col = col_column
+    row_column = row_column, col_column = col_column
   )
 
   ed <- calculate_ed(design_matrix, current_score_obj$ed, swapped_items)

From 2d23eb49cd3da7accdcb936ae74c05f69bc7f98e Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Thu, 11 Jun 2026 20:32:00 +0930
Subject: [PATCH 8/9] Updating docs

---
 man/build_design_matrix.Rd          | 14 ++++++++------
 man/calculate_pair_incidence.Rd     |  2 +-
 man/calculate_position_incidence.Rd |  2 +-
 3 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/man/build_design_matrix.Rd b/man/build_design_matrix.Rd
index 469d041..f104b19 100644
--- a/man/build_design_matrix.Rd
+++ b/man/build_design_matrix.Rd
@@ -4,23 +4,25 @@
 \alias{build_design_matrix}
 \title{Build a Spatial Design Matrix from a Data Frame}
 \usage{
-build_design_matrix(df, swap, row_col = "row", col_col = "col")
+build_design_matrix(df, swap, row_column = "row", col_column = "col")
 }
 \arguments{
-\item{df}{A data frame with columns named by \code{swap}, \code{row_col}, \code{col_col}.}
+\item{df}{A data frame with columns named by \code{swap}, \code{row_column},
+\code{col_column}.}
 
 \item{swap}{Column name of the treatment variable.}
 
-\item{row_col}{Column name of the row position variable (default \code{"row"}).}
+\item{row_column}{Column name of the row position variable (default \code{"row"}).}
 
-\item{col_col}{Column name of the column position variable (default \code{"col"}).}
+\item{col_column}{Column name of the column position variable
+(default \code{"col"}).}
 }
 \value{
 A character matrix of dimensions \verb{max(row) × max(col)}.
 }
 \description{
-Places each treatment value at the grid position given by its \code{row_col} and
-\code{col_col} coordinates, returning a character matrix of dimensions
+Places each treatment value at the grid position given by its \code{row_column}
+and \code{col_column} coordinates, returning a character matrix of dimensions
 \verb{max(row) × max(col)}. Cells with no corresponding row in \code{df} are \code{NA}.
 Unlike filling via \code{matrix(..., byrow)}, this is robust to any row ordering
 of \code{df}.
diff --git a/man/calculate_pair_incidence.Rd b/man/calculate_pair_incidence.Rd
index 49f28a3..9b50b86 100644
--- a/man/calculate_pair_incidence.Rd
+++ b/man/calculate_pair_incidence.Rd
@@ -55,7 +55,7 @@ df <- data.frame(
   col = rep(1:3, each = 4),
   treatment = rep(c("A", "B", "C"), 4)
 )
-result <- speed(df, swap = "treatment", seed = 42)
+result <- speed(df, swap = "treatment", seed = 42, iterations = 200)
 calculate_pair_incidence(result)
 
 # From a plain data frame; as named list
diff --git a/man/calculate_position_incidence.Rd b/man/calculate_position_incidence.Rd
index 1eeefd9..fa17173 100644
--- a/man/calculate_position_incidence.Rd
+++ b/man/calculate_position_incidence.Rd
@@ -45,7 +45,7 @@ df <- data.frame(
   col = rep(1:3, each = 4),
   treatment = rep(c("A", "B", "C"), 4)
 )
-result <- speed(df, swap = "treatment", seed = 42)
+result <- speed(df, swap = "treatment", seed = 42, iterations = 200)
 res <- calculate_position_incidence(result)
 res$row  # treatment × row counts
 res$col  # treatment × column counts

From 98653502ec8468e11fc6d8274b51186d0e6267ca Mon Sep 17 00:00:00 2001
From: Sam Rogers <7007561+rogerssam@users.noreply.github.com>
Date: Sun, 26 Jul 2026 07:43:15 +0930
Subject: [PATCH 9/9] Adding plans

---
 PLAN-adjacency-orientation-fix.md  | 516 +++++++++++++++++++++++++
 PR-91-alignment-and-integration.md | 303 +++++++++++++++
 PR-91-info-objective-review.md     | 193 ++++++++++
 PR-incidence-review.md             | 579 +++++++++++++++++++++++++++++
 4 files changed, 1591 insertions(+)
 create mode 100644 PLAN-adjacency-orientation-fix.md
 create mode 100644 PR-91-alignment-and-integration.md
 create mode 100644 PR-91-info-objective-review.md
 create mode 100644 PR-incidence-review.md

diff --git a/PLAN-adjacency-orientation-fix.md b/PLAN-adjacency-orientation-fix.md
new file mode 100644
index 0000000..bb004c0
--- /dev/null
+++ b/PLAN-adjacency-orientation-fix.md
@@ -0,0 +1,516 @@
+# Plan: make grid construction coordinate-based (fixes `objective_function_piepho()`)
+
+**Suggested branch:** `bugfix/grid-orientation` (off `main`) — matches the existing
+`bugfix/adjacency-sorting`, `bugfix/ed` naming.
+
+**Status of evidence:** every claim was measured at the console against R 4.6.1 with
+`pkgload::load_all()`. Test constants in §6 are hand-derivable, so the tests assert *correct* values
+rather than pinning current behaviour.
+
+> **This plan was revised after review feedback.** The first version claimed
+> `calculate_adjacency_score()` was broken for all non-square designs produced by `speed()`. That
+> was **wrong** — `speed()` sorts its input to row-major, which is what `byrow = TRUE` expects. The
+> corrected analysis is in §1; the genuinely broken function is `objective_function_piepho()`, and
+> it is broken for *all* grid shapes including square ones.
+
+---
+
+## 1. The bug
+
+### 1.1 What `speed()` already does
+
+`speed()` sorts its input before optimising ([R/speed.R:188-192](R/speed.R#L188-L192)):
+
+```r
+if (inferred$inferred) {
+  # Sort the data frame to start with to ensure consistency in calculating the adjacency later
+  data <- data[do.call(order, data[c(row_column, col_column)]), ]
+  rownames(data) <- seq_len(nrow(data))
+}
+```
+
+`order(row, col)` produces **row-major** order — row varies slowest, col fastest. `byrow = TRUE`
+fills row-major. So the sort and `calculate_adjacency_score()` are correctly paired, and the comment
+says so explicitly. Verified: on sorted data, `matrix(..., byrow = TRUE)` is `identical()` to
+coordinate placement for 3×3, 4×3, 3×4, 6×4, 4×4 and 12×5.
+
+**Designs produced by `speed()` with the default objective and numeric row/col columns were
+correct.** No regeneration needed.
+
+### 1.2 The actual bug: two functions, opposite assumptions
+
+Neither function reads coordinates; each hardcodes an ordering assumption, and they disagree:
+
+| Function | Fill | Assumes | Inside `speed()` (row-major) | Raw `initialise_design_df()` (column-major) |
+|---|---|---|---|---|
+| `calculate_adjacency_score()` | `matrix(..., byrow = TRUE)` | row-major | ✅ correct | ❌ wrong |
+| `objective_function_piepho()` | `matrix(...)` (column-major) | column-major | ❌ **wrong** | ✅ correct |
+
+Each is wrong exactly where the other is right. `initialise_design_df()` produces
+`expand.grid(row = 1:nrows, col = 1:ncols)` ([R/design_utils.R:304](R/design_utils.R#L304)), which is
+column-major — the opposite of what `speed()` feeds the objective.
+
+### 1.3 Impact 1 — `objective_function_piepho()` (all shapes) 🔴
+
+Piepho receives `speed()`'s row-major data and fills column-major. Verified **not identical** to
+coordinate placement for 3×3, 4×3, 6×4, 4×4 and 12×5 — square grids included, so the
+transpose-invariance that saved adjacency does not apply here. On a 4×3 design:
+
+```
+piepho old grid          correct grid
+  T1  T2  T3               T1  T1  T1
+  T1  T2  T4               T2  T2  T2
+  T1  T3  T4               T3  T3  T3
+  T2  T3  T4               T4  T4  T4
+```
+
+The multiset of treatments is preserved, so replication counts look plausible, but the spatial
+arrangement — the entire input to `calculate_ed()` — is scrambled. **Anyone who used
+`obj_function = objective_function_piepho` has designs optimised against a layout that isn't
+theirs.** These need regenerating.
+
+### 1.4 Impact 2 — `calculate_adjacency_score()` doesn't compose 🟠
+
+It is exported, and wrong when handed anything that isn't already row-major — including the output of
+the package's own `initialise_design_df()`:
+
+```r
+df <- initialise_design_df(items = rep(paste0("T", 1:4), length.out = 12),
+                           nrows = 4, ncols = 3)
+calculate_adjacency_score(df, "treatment")   # main: 0    correct: 8
+```
+
+Its documented examples use hand-written row-major data, so they pass and the inconsistency is
+invisible. Two exported functions in the same package that silently disagree about layout.
+
+### 1.5 Impact 3 — the sort is applied to factors 🟠
+
+`to_factor()` runs at [R/speed.R:185](R/speed.R#L185), *before* the sort, converting every column
+including `row`. `order()` on a factor follows **level order**. For an integer row column
+`as.factor(1:12)` sorts numerically and all is well — but a **character** row column with ≥10 rows
+gets levels `1, 10, 11, 12, 2, …`, and the sort follows that, so `byrow = TRUE` fills a permuted
+grid. Demonstrated on an 11×1 design with a genuine like-pair at true rows 1-2:
+
+```
+lexical row order : 1 10 11 2 3 4 5 6 7 8 9
+old grid    : A J K A C D E F G H I     score 0
+correct grid: A A C D E F G H I J K     score 1
+```
+
+So adjacency *inside* `speed()` is also wrong for this input shape — the sort silently does not
+deliver the invariant that `byrow = TRUE` depends on.
+
+---
+
+## 2. Root cause and fix
+
+`byrow` can only encode an assumption about data ordering; it cannot read coordinates. The sort at
+[R/speed.R:190](R/speed.R#L190) exists to manufacture that assumption, but it (a) only helps the one
+function whose assumption it matches, (b) is defeated by lexical factor levels, and (c) does nothing
+for direct calls to exported functions.
+
+Placing values by their actual `(row, col)` coordinates makes all call paths correct and makes the
+sort unnecessary for correctness.
+
+---
+
+## 3. Scope
+
+**In scope**
+
+1. Add internal `build_design_matrix()` to `R/design_utils.R`.
+2. Point `objective_function_piepho()` at it — **the bug fix**.
+3. Point `calculate_adjacency_score()` at it — robustness + composability.
+4. New `tests/testthat/test-build_design_matrix.R`.
+5. Regression tests: piepho grid orientation; adjacency direct-call and character-label cases.
+6. `NEWS.md` **Bug Fixes** entry; `DESCRIPTION` version bump.
+
+**Out of scope** (recorded so they aren't lost)
+
+- **Removing the sort at [R/speed.R:190](R/speed.R#L190).** Once grids are coordinate-based the sort
+  is no longer needed *for correctness*, but other code may rely on row order (`generate_neighbour`,
+  `random_initialise`, `print.design`, `autoplot`). Leave it; note it as a possible later
+  simplification. Do **not** bundle its removal with a bug fix.
+- **Duplicate `(row, col)` coordinates** — multi-site/MET designs reuse coordinates and
+  `build_design_matrix()` keeps only the last written. Separate pre-existing bug; see §8.
+- Performance work beyond keeping this cost-neutral (§5.1).
+- The new incidence functions — they stay on `feature/incidence`.
+
+---
+
+## 4. Branch setup
+
+`build_design_matrix()` was introduced on `feature/incidence` in commit `a4eff49`, mixed in with the
+incidence functions, so there is no clean cherry-pick. Start from `main`; the function is ~20 lines.
+
+```sh
+git checkout main
+git pull
+git checkout -b bugfix/grid-orientation
+```
+
+---
+
+## 5. Implementation
+
+### 5.1 Add `build_design_matrix()` to `R/design_utils.R`
+
+Append after `initialize_design_df`. Two differences from the `feature/incidence` version: the
+coordinate vectors are computed **once** into locals, and there is a clear error for non-coercible
+coordinates.
+
+```r
+#' Build a Spatial Design Matrix from a Data Frame
+#'
+#' @description
+#' Places each treatment value at the grid position given by its `row_column`
+#' and `col_column` coordinates, returning a character matrix of dimensions
+#' `max(row)` by `max(col)`. Cells with no corresponding row in `df` are `NA`.
+#' Unlike filling via `matrix(..., byrow)`, this is robust to any row ordering
+#' of `df`, and to factor columns whose level order is not numeric.
+#'
+#' @param df A data frame with columns named by `swap`, `row_column`,
+#'   `col_column`.
+#' @param swap Column name of the treatment variable.
+#' @param row_column Column name of the row position variable (default `"row"`).
+#' @param col_column Column name of the column position variable
+#'   (default `"col"`).
+#'
+#' @return A character matrix of dimensions `max(row)` by `max(col)`.
+#'
+#' @keywords internal
+build_design_matrix <- function(df,
+                                swap,
+                                row_column = "row",
+                                col_column = "col") {
+  rows <- as_numeric_factor(df[[row_column]])
+  cols <- as_numeric_factor(df[[col_column]])
+
+  if (all(is.na(rows)) || all(is.na(cols))) {
+    stop(
+      "Cannot place the design on a grid: `", row_column, "` and `",
+      col_column, "` must be numeric, or coercible to numeric.",
+      call. = FALSE
+    )
+  }
+
+  design_matrix <- matrix(
+    NA_character_,
+    nrow = max(rows, na.rm = TRUE),
+    ncol = max(cols, na.rm = TRUE)
+  )
+  design_matrix[cbind(rows, cols)] <- as.character(df[[swap]])
+  design_matrix
+}
+```
+
+**On the two locals** — this matters. The obvious implementation calls `as_numeric_factor()` four
+times (twice for `max()`, twice for the index). It is `as.numeric(as.character(x))`
+([R/utils.R:301](R/utils.R#L301)) — two allocating passes each — and this runs once per SA iteration
+(10,000 by default, per level). Hoisting keeps it at two passes, matching `main`. A bug fix should
+not be a performance regression.
+
+**On the error branch** — non-numeric coordinate labels currently give `NAs introduced by coercion`
+then `invalid 'nrow' value (too large or NA)`. Both pre-existing; a clear message is cheap here.
+
+### 5.2 `R/metrics.R` — `objective_function_piepho()` (the fix)
+
+```r
+  design_matrix <- build_design_matrix(
+    design, swap,
+    row_column = row_column, col_column = col_column
+  )
+```
+
+Confirmed safe downstream: `matrix()` on a factor already returns a **character** matrix in R 4.6.1,
+so `calculate_ed()`'s type-sensitive
+`design_matrix[!(design_matrix %in% swapped_items)] <- NA`
+([R/metrics.R:388](R/metrics.R#L388)) is unchanged, and `get_vertices()` calls
+`as.character(design_matrix)` anyway ([R/metrics.R:470](R/metrics.R#L470)).
+
+### 5.3 `R/calculate_adjacency_score.R`
+
+Same replacement. This is **not** a behaviour change for `speed()` with numeric row/col columns
+(verified identical on sorted data across six grid shapes). It fixes direct calls (§1.4) and the
+character-label case (§1.5).
+
+### 5.4 Formatting
+
+```sh
+air format R/
+```
+
+Note the single-line signature used on `feature/incidence`
+(`build_design_matrix <- function(df, swap, row_column = "row", col_column = "col") {`) is 83
+characters and violates `air.toml`'s 80-column setting — hence the wrapped signature above.
+
+---
+
+## 6. Tests
+
+### 6.1 New file: `tests/testthat/test-build_design_matrix.R`
+
+```r
+test_that("treatments are placed by coordinate, not by data-frame order", {
+  # Column-major input, as initialise_design_df() produces:
+  #   A A A
+  #   B B B
+  #   C C C
+  #   D D D
+  df <- data.frame(
+    row = rep(1:4, times = 3),
+    col = rep(1:3, each = 4),
+    swap = rep(c("A", "B", "C", "D"), times = 3)
+  )
+  expected <- matrix(
+    rep(c("A", "B", "C", "D"), times = 3),
+    nrow = 4, ncol = 3
+  )
+  expect_equal(build_design_matrix(df, "swap"), expected)
+})
+
+test_that("result is invariant to row ordering of the data frame", {
+  df <- data.frame(
+    row = rep(1:4, times = 3),
+    col = rep(1:3, each = 4),
+    swap = rep(c("A", "B", "C", "D"), times = 3)
+  )
+  withr::with_seed(9, shuffled <- df[sample(nrow(df)), ])
+  expect_equal(
+    build_design_matrix(shuffled, "swap"),
+    build_design_matrix(df, "swap")
+  )
+})
+
+test_that("factor columns with non-numeric level order are placed correctly", {
+  # as.factor() on character rows gives levels 1, 10, 11, 2, ... Coordinates
+  # must be read numerically, not by level order.
+  df <- data.frame(
+    row = factor(as.character(1:11)),
+    col = 1,
+    swap = c("A", "A", "C", "D", "E", "F", "G", "H", "I", "J", "K")
+  )
+  m <- build_design_matrix(df, "swap")
+  expect_equal(dim(m), c(11L, 1L))
+  expect_equal(as.vector(m), c("A", "A", "C", "D", "E", "F", "G", "H", "I", "J", "K"))
+})
+
+test_that("missing grid positions are NA", {
+  # rows 1, 2, 4, 5 present; row 3 absent
+  df <- data.frame(
+    row = c(1, 2, 4, 5),
+    col = c(1, 1, 1, 1),
+    swap = c("A", "A", "B", "B")
+  )
+  m <- build_design_matrix(df, "swap")
+  expect_equal(dim(m), c(5L, 1L))
+  expect_true(is.na(m[3, 1]))
+  expect_equal(as.vector(m), c("A", "A", NA, "B", "B"))
+})
+
+test_that("non-numeric coordinates give an informative error", {
+  df <- data.frame(row = c("R1", "R2"), col = c(1, 1), swap = c("A", "B"))
+  expect_error(
+    suppressWarnings(build_design_matrix(df, "swap")),
+    "must be numeric"
+  )
+})
+```
+
+### 6.2 `tests/testthat/test-objective_functions.R` — the piepho regression
+
+This is the test that would have caught the real bug. Requires `igraph` (a `Suggests`).
+
+```r
+test_that("objective_function_piepho scores the design's actual layout", {
+  skip_if_not_installed("igraph")
+
+  # Row-major data, exactly as speed() supplies after its internal sort.
+  #   A A A
+  #   B B B
+  #   C C C
+  # Piepho's pre-fix column-major fill produced a scrambled grid for this input.
+  df <- data.frame(
+    row = rep(1:3, each = 3),
+    col = rep(1:3, times = 3),
+    treatment = rep(c("A", "B", "C"), each = 3)
+  )
+  expect_equal(
+    build_design_matrix(df, "treatment"),
+    matrix(rep(c("A", "B", "C"), each = 3), nrow = 3, byrow = TRUE)
+  )
+
+  # Score must not depend on the row ordering of the input.
+  withr::with_seed(4, shuffled <- df[sample(nrow(df)), ])
+  expect_equal(
+    objective_function_piepho(shuffled, "treatment", c("row", "col"))$score,
+    objective_function_piepho(df, "treatment", c("row", "col"))$score
+  )
+})
+```
+
+### 6.3 `tests/testthat/test-calculate_adjacency_score.R` — composability
+
+```r
+test_that("adjacency score is correct for column-major input", {
+  # initialise_design_df() is column-major; the pre-fix byrow = TRUE fill
+  # scored this 0 instead of 8.
+  #   A A A     horizontal like-pairs: 2 per row x 4 rows = 8
+  #   B B B     vertical like-pairs:   0
+  #   C C C
+  #   D D D
+  design <- data.frame(
+    row = rep(1:4, times = 3),
+    col = rep(1:3, each = 4),
+    swap = rep(c("A", "B", "C", "D"), times = 3)
+  )
+  expect_equal(calculate_adjacency_score(design, "swap"), 8)
+})
+
+test_that("adjacency score composes with initialise_design_df()", {
+  # 12 x 5, 4 treatments cycling down each column of 12: because 12 is a
+  # multiple of 4 every column is identical, so every row is constant.
+  # horizontal like-pairs: 4 per row x 12 rows = 48; vertical: 0.
+  df <- initialise_design_df(
+    items = rep(paste0("T", 1:4), length.out = 60),
+    nrows = 12, ncols = 5
+  )
+  expect_equal(calculate_adjacency_score(df, "treatment"), 48)
+})
+
+test_that("adjacency score does not depend on data-frame row order", {
+  design <- data.frame(
+    row = rep(1:4, times = 3),
+    col = rep(1:3, each = 4),
+    swap = rep(c("A", "B", "C", "D"), times = 3)
+  )
+  withr::with_seed(9, shuffled <- design[sample(nrow(design)), ])
+  expect_equal(
+    calculate_adjacency_score(shuffled, "swap"),
+    calculate_adjacency_score(design, "swap")
+  )
+})
+
+test_that("adjacency score reads factor coordinates numerically", {
+  # Lexical levels 1, 10, 11, 2, ... must not permute the grid. True rows 1-2
+  # hold the only like-pair, so the score is 1; the pre-fix fill gave 0.
+  design <- data.frame(
+    row = factor(as.character(1:11)),
+    col = 1,
+    swap = c("A", "A", "C", "D", "E", "F", "G", "H", "I", "J", "K")
+  )
+  expect_equal(calculate_adjacency_score(design, "swap"), 1)
+})
+```
+
+### 6.4 Existing tests need no changes
+
+Confirmed: the four expectations in
+[test-calculate_adjacency_score.R:1-28](tests/testthat/test-calculate_adjacency_score.R#L1-L28)
+(6, 2, 7, 0) still hold — they use 3×3 designs with row-major data. The equivalence assertion at
+[test-objective_functions.R:68](tests/testthat/test-objective_functions.R#L68) also still holds (both
+sides give 6). **No test churn**, which keeps the diff honest and easy to review.
+
+Worth noting for the PR description: the suite's blind spot was that every adjacency case used
+row-major 3×3 data — the one shape and ordering for which both old implementations agreed with the
+truth.
+
+---
+
+## 7. NEWS and version
+
+`main` is at `0.0.8`. Bump `DESCRIPTION` to `0.0.9` and add to the top of `NEWS.md`:
+
+```markdown
+# speed 0.0.9
+
+## Bug Fixes
+
+- `objective_function_piepho()` built its design grid by filling column-wise from the data frame,
+  but `speed()` supplies rows in row-major order. The grid was therefore scrambled for all design
+  shapes, so designs optimised with `obj_function = objective_function_piepho` were scored against
+  a layout other than their own. **Designs produced with this objective should be regenerated.**
+  Designs using the default `objective_function()` were not affected. (#NNN)
+- `calculate_adjacency_score()` assumed its input was already sorted in row-major order, so calling
+  it directly on a data frame in any other order — including the column-major output of
+  `initialise_design_df()` — returned an incorrect score. It now places treatments by their `row`
+  and `col` coordinates and is independent of row ordering. (#NNN)
+- Grid construction now reads `row` and `col` numerically rather than by factor level order, fixing
+  designs with character row or column labels and ten or more rows, where lexical level ordering
+  (`1, 10, 11, 2, ...`) permuted the grid. (#NNN)
+```
+
+Replace `#NNN` with the issue/PR number per `CLAUDE.md`.
+
+**Version conflict:** `feature/incidence` already bumped to `0.0.9`. Whichever merges second needs
+`0.0.10` — expect a trivial `DESCRIPTION`/`NEWS.md` conflict.
+
+---
+
+## 8. Follow-up: duplicate coordinates (separate issue)
+
+Worth filing before this branch closes.
+
+`build_design_matrix()` silently keeps only the last value written to a repeated `(row, col)`.
+Multi-site designs reuse coordinates — the MET example in `?speed`
+([R/speed.R:118-122](R/speed.R#L118-L122)) `rbind`s five copies of a 28×5 site grid. Verified on a
+two-site reduction: 12 data rows collapse into 6 cells and only the last site survives. `main` is
+also wrong here, just differently (it truncates to the first `nrow × ncol` values).
+
+That example sets `adj_weight = 0`, and `objective_function()`'s `ifelse` leaves the score
+unevaluated, so the documented example doesn't hit it. Any multi-site design with adjacency enabled
+does.
+
+**Where the check belongs:** *not* in `build_design_matrix()` — anything that warns or errors there
+fires once per SA iteration. It belongs in `speed()`'s up-front validation as a `.verify_*` helper in
+`R/verify_utils.R`, where it runs once and can name the offending column and suggest `grid_factors`
+or a site-specific coordinate column.
+
+---
+
+## 9. Acceptance criteria
+
+```sh
+devcontainer-exec-here R -e "devtools::document()"
+devcontainer-exec-here R -e "devtools::test()"
+devcontainer-exec-here R -e "devtools::check()"
+```
+
+- [ ] `build_design_matrix()` added to `R/design_utils.R`; `man/` regenerated with
+      `devtools::document()` (never hand-edit `man/*.Rd`).
+- [ ] `objective_function_piepho()` and `calculate_adjacency_score()` both use it.
+- [ ] New `test-build_design_matrix.R` passes; new piepho and adjacency cases pass.
+- [ ] All pre-existing tests pass **unmodified**.
+- [ ] `devtools::check()` clean.
+- [ ] `air format R/` produces no further changes.
+- [ ] `NEWS.md` Bug Fixes entries with issue/PR links; `DESCRIPTION` bumped.
+- [ ] Confirm a `speed()` run with the **default** objective, numeric row/col, gives an unchanged
+      result for a fixed seed — this fix should be a no-op on that path.
+- [ ] Benchmark a representative `speed()` run against `main` to confirm no slowdown (§5.1).
+- [ ] Follow-up issue filed for duplicate coordinates (§8).
+
+`igraph` is a `Suggests` dependency needed by `calculate_ed()` for 4+ replicates; install it before
+running the piepho tests or they will error (or they will skip, per `skip_if_not_installed`).
+
+---
+
+## 10. Reviewer summary
+
+Suggested PR description:
+
+> `speed()` sorts its input to row-major order before optimising, and `calculate_adjacency_score()`
+> fills its grid with `byrow = TRUE` to match. But `objective_function_piepho()` fills column-major,
+> so it scored a scrambled grid for every design shape — including square ones. Designs optimised
+> with that objective were scored against a layout other than their own.
+>
+> The underlying problem is that both functions encode an assumption about data ordering instead of
+> reading coordinates. This replaces both with coordinate-based placement, which also fixes direct
+> calls to the exported `calculate_adjacency_score()` on column-major data (such as
+> `initialise_design_df()` output, previously scoring a 4×3 design 0 instead of 8), and designs whose
+> row labels are characters with ten or more rows, where lexical factor levels defeated the sort.
+>
+> No existing test changed — the suite's adjacency cases were all row-major 3×3, the one shape and
+> ordering where both old implementations happened to be right. Adds a test file for the new helper
+> plus regression cases for each affected path.
+>
+> The default `objective_function()` path is unaffected; piepho users should regenerate designs.
diff --git a/PR-91-alignment-and-integration.md b/PR-91-alignment-and-integration.md
new file mode 100644
index 0000000..139c89a
--- /dev/null
+++ b/PR-91-alignment-and-integration.md
@@ -0,0 +1,303 @@
+# PR #91 alignment + how `feature/incidence` fits on top
+
+Two questions:
+
+1. Beyond `calc_` → `calculate_`, what else should PR #91 change to match the rest of the codebase?
+2. Assuming #91 merges first, how should `feature/incidence` fit in?
+
+All claims verified at the console against R 4.6.1 (`pkgload::load_all()`), with #91's functions
+sourced from the `info-objective` branch.
+
+---
+
+## Part 1 — Aligning PR #91 with the codebase
+
+### 1.1 The rename is load-bearing, not cosmetic 🔴
+
+`calculate_` isn't just a house style — `_pkgdown.yml` builds its reference index from prefix
+patterns:
+
+```yaml
+  - subtitle: Calculation functions
+  - contents:
+      - starts_with("calculate")
+  - subtitle: Objective functions
+  - contents:
+      - starts_with("objective_function")
+```
+
+Which of #91's 8 exports match a pattern:
+
+| Export | Matches index? |
+|---|---|
+| `objective_function_info` | ✅ `objective_function*` |
+| `calculate_efficiency_factors` | ✅ `calculate*` |
+| `calc_info_matrix` | ❌ |
+| `calc_incidence_matrix` | ❌ |
+| `calc_concurrence_matrix` | ❌ |
+| `compute_L_projection` | ❌ |
+| `cor_ar1` | ❌ |
+| `cor_ar1_ar1` | ❌ |
+
+pkgdown **errors** when exported topics are absent from an explicit `reference:` index, so as it
+stands #91 breaks the docs site build. The `calc_` → `calculate_` rename fixes three of the five
+automatically. The remaining two groups still need attention:
+
+- `compute_L_projection` — either rename to `calculate_L_projection` (picked up automatically), or
+  add explicitly to the index. I'd rename: the package has no `compute_*` prefix anywhere, and this
+  function is doing the same kind of job as the `calculate_*` family.
+- `cor_ar1` / `cor_ar1_ar1` — these are genuinely a new category (covariance-structure
+  constructors). Add a new `_pkgdown.yml` section rather than contorting the names:
+
+  ```yaml
+    - title: Covariance structures
+      desc: Constructors for spatial correlation structures.
+    - contents:
+        - starts_with("cor_")
+  ```
+
+### 1.2 `treatment_column` → `swap` 🟠
+
+#91 is inconsistent *with itself*: `objective_function_info(layout_df, swap, spatial_cols, ...)`
+uses `swap`, but all five analysis helpers use `treatment_column`. The package uses `swap`
+universally — `calculate_adjacency_score(layout_df, swap, ...)`,
+`calculate_balance_score(layout_df, swap, spatial_cols)`, `objective_function(layout_df, swap, ...)`.
+
+Rename `treatment_column` → `swap` throughout. Note the one genuine exception already in the
+codebase: `calculate_efficiency_factor(design_df, item)` uses `item` *and* `substitute()` for
+unquoted input — that's an outlier, not a precedent to follow.
+
+While there: decide whether the treatment argument has a default. Currently there are three
+conventions in play — `calculate_adjacency_score()` has none, #91 defaults to `"treatment"`, and
+`feature/incidence` defaults to `"treatment"`. Pick one (I'd default to `"treatment"`, since it's
+what `initialise_design_df()` produces) and apply it consistently.
+
+### 1.3 `levels(factor(...))` → `stri_sort(numeric = TRUE)` 🟠
+
+`.compute_info()` and `calc_incidence_matrix()` both order treatments with `levels(factor(x))`,
+i.e. lexically. `initialise_design_df(items = 12)` produces `T1`…`T12`, so:
+
+```
+PR #91 rownames : T1 T10 T11 T12 T2 T3 T4 T5 T6 T7 T8 T9
+this branch     : T1 T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 T12
+```
+
+The package already imports `stringi` and uses `stri_sort(..., numeric = TRUE)` in `speed()` for
+exactly this reason. This affects #91's *output* (incidence/concurrence row order, and the row/col
+order of the information matrix), so it's user-visible, not internal. Worth a small shared internal
+helper — `.treatment_levels(x)` — used by #91, this branch, and `speed()`, so it can never drift
+again.
+
+### 1.4 Roxygen: markdown, not Rd macros 🟠
+
+`DESCRIPTION` sets `Roxygen: list(markdown = TRUE)` and `CLAUDE.md` calls markdown roxygen the
+source of truth. #91 uses `\code{...}` throughout and refers to other functions as
+`\code{compute_L_projection}` — which renders as code but is **not a link**. Package style is
+backticks and `[fn()]` auto-links:
+
+```r
+#' Use [calculate_L_projection()] to generate it.   # links
+#' Use \code{compute_L_projection} to generate it.   # doesn't
+```
+
+Keep `\eqn{}` / `\deqn{}` — there's no markdown equivalent and the maths is good. Just convert
+`\code{}` → backticks and cross-references → `[fn()]`.
+
+### 1.5 Missing `@return` on an exported function 🟠
+
+`calc_info_matrix()` has `@param` tags and `@export` but **no `@return`**
+([objectives.R:249-262](R/objectives.R)). `R CMD check` flags missing `\value` in Rd for exported
+objects, and `check()` must pass before merge. It returns a 5-element list
+(`info_matrix`, `eigenvalues`, `rank`, `A_value`, `D_value`) — worth an `\describe{}` block, as
+`calculate_position_incidence()` does on the other branch.
+
+Also `objective_function_info()` places `@return` *after* `@examples`; every other file puts it
+before `@export`. Harmless but inconsistent.
+
+### 1.6 The `swapped_items` / `current_score_obj` contract 🔴
+
+This is the biggest convention gap, and `CLAUDE.md` is explicit about it:
+
+> The objective function … returns a `list(score = ..., ...)`. The returned list is fed back as
+> `current_score_obj` on the next call so objective functions can incrementally update internal
+> state from `swapped_items` instead of recomputing from scratch — **any custom objective function
+> must honour this contract**.
+
+`objective_function_info()` accepts neither `current_score_obj` nor `swapped_items` (they vanish
+into `...`), and although it *returns* `info_matrix` and `eigenvalues`, it never consumes them. So
+every iteration rebuilds `X1` with an R `for` loop, forms `t(X1) %*% L %*% X1` against a dense
+n×n `L`, and runs `eigen()` — 10,000 times by default. A swap touches two plots, so the information
+matrix update is rank ≤ 4 and the whole thing should be a downdate/update.
+
+Compare `objective_function_piepho()`, which threads `current_score_obj$ed` and `swapped_items`
+into `calculate_ed()` precisely as the contract requires. #91 should follow that pattern, or at
+minimum document loudly that it is O(n²v + v³) per iteration.
+
+### 1.7 Vectorise the indicator-matrix loops 🟡
+
+Three places build indicator matrices with `for` loops:
+
+```r
+# .build_treatment_matrix
+for (i in seq_len(n)) X1[i, trt_idx[i]] <- 1
+
+# .build_L_from_df and compute_L_projection (duplicated)
+for (j in seq_len(b)) X2[blocks == levels(blocks)[j], j] <- 1
+```
+
+Both are one-liners with the matrix-index idiom that `feature/incidence` just introduced in
+`build_design_matrix()`:
+
+```r
+X1[cbind(seq_len(n), trt_idx)] <- 1
+X2[cbind(seq_along(blocks), as.integer(blocks))] <- 1
+```
+
+The `X2` construction is also duplicated verbatim between `.build_L_from_df()` and
+`compute_L_projection()` — factor into one internal helper. (`calculate_efficiency_factor()` on
+`main` has the same loop-built `X`; not #91's job to fix, but worth knowing the idiom is already
+inconsistent.)
+
+### 1.8 Reuse `pseudo_inverse()` 🟡
+
+`R/utils.R` already provides `pseudo_inverse()` (SVD-based Moore-Penrose with a rank guard).
+`calc_info_matrix()` and `objective_function_info()` reimplement the same rank-deficiency handling
+via `eigen()` with a `1e-10` relative tolerance. Reuse the helper, or explain in a comment why
+`eigen` is preferred here (it is symmetric, so it's defensible — but then say so).
+
+### 1.9 Smaller items 🟡
+
+- **`=` instead of `<-`** at `A_val = sum(1 / pos_eig)` / `D_val = -sum(log(pos_eig))`
+  ([objectives.R:288-289](R/objectives.R)). Package uses `<-` throughout; Air won't fix this.
+- **`Sigma` is capitalised.** Every other argument in the package is lower snake_case. Consider
+  `sigma` or `covariance`; `L_matrix` is borderline but at least reads as a proper noun.
+- **`stop()` without `call. = FALSE`.** `compute_L_projection()`'s six validation `stop()`s omit
+  it; `R/buffers.R` and `R/verify_utils.R` use `call. = FALSE`. Also consider moving them into
+  `R/verify_utils.R` as a `.verify_*` helper, which is where all other input validation lives.
+- **Bare `setNames()`** at [objectives.R:398](R/objectives.R) — resolves only because other files
+  carry `@importFrom stats setNames`. `R/buffers.R` uses explicit `stats::setNames()`; match that.
+- **`class(N) <- "matrix"`** leaves `names(dimnames(N))` as `c("treatments", "blocks")`, so the
+  matrix prints with stray header rows:
+
+  ```
+            blocks
+  treatments 1 2 3 4
+         T1  1 0 0 0
+  ```
+
+  Add `names(dimnames(N)) <- NULL`, or build with `matrix(as.integer(tbl), ...)` as
+  `feature/incidence` does. Pick one idiom for both branches.
+- **`spatial_cols` is accepted and never used.** `objective_function_info()` takes it to satisfy the
+  signature but ignores it entirely — block structure comes from `block_column` instead. A user who
+  sets `spatial_factors = ~ row + col + block` and switches to this objective gets silence. Document
+  it, or warn when `spatial_cols` is supplied and non-trivial.
+- **`calculate_efficiency_factors` vs existing `calculate_efficiency_factor`** — singular/plural
+  differing by one character, computing different things (canonical efficiency factors from the
+  information matrix vs. Piepho's single efficiency value). Rename to
+  `calculate_canonical_efficiency_factors()`.
+- **File placement** — every other `objective_function_*` lives in `R/metrics.R`. See §2.3 for a
+  suggested layout that resolves this alongside `feature/incidence`.
+- **`NEWS.md`** — 8 new exports and no entry.
+- **Unrelated `DESCRIPTION` whitespace churn** and the dead `.Rbuildignore` entry (both already
+  noted in `PR-91-info-objective-review.md`).
+
+---
+
+## Part 2 — How `feature/incidence` fits on top of #91
+
+### 2.1 The key finding: `calculate_position_incidence()` becomes redundant 🔴
+
+These are the **same computation**. Verified on a 12-treatment, 4×3 blocked design:
+
+```r
+calc_incidence_matrix(df, block_column = "row")   # PR #91
+calculate_position_incidence(df, swap)$row        # feature/incidence
+#> identical values; differ only in row ordering (§1.3)
+```
+
+Both reduce to `table(treatment, spatial_factor)`. #91's version takes one factor; this branch's
+takes exactly two (hardcoded `row` and `col`) and returns a list. Neither can do what users
+actually need, which is *arbitrary* factors — blocks, `site_col`, `site_block`.
+
+**So don't port `calculate_position_incidence()` on top of #91.** Instead generalise #91's function
+to accept multiple factors, mirroring `speed()`'s own interface:
+
+```r
+calculate_incidence(design, swap = "treatment", spatial_factors = ~ row + col)
+# -> named list of treatment x level matrices, one per factor
+```
+
+That single change:
+
+- removes the duplicate implementation,
+- fixes the hardcoded-row/col gap in this branch (review §2.4),
+- fixes the block gap in this branch (review §2.4),
+- and makes the argument name match `speed()`, so users write the same formula in both places.
+
+`calculate_position_incidence()` then disappears from the diff entirely. Nothing is lost — it was
+`calculate_incidence(design, spatial_factors = ~ row + col)`.
+
+### 2.2 What `feature/incidence` uniquely contributes
+
+Stripping out the redundancy, three things in this branch are genuinely additive and should survive:
+
+1. **`build_design_matrix()` and the adjacency orientation fix.** The most valuable thing in the
+   branch — `main` scores a 12×5 design as 0 adjacencies when it has 48 (see
+   `PR-incidence-review.md` §3.3). Independent of #91. **Split this into its own PR and merge it
+   first, ahead of both features.**
+2. **Adjacency-based pair counting.** #91's concurrence is block-based; nothing in #91 counts
+   *spatially adjacent* pairs. This is the real new capability.
+3. **The print/summary layer.** Neither branch has one, and both produce v×v or v×b matrices that
+   are unreadable at realistic sizes. Best built once, over a shared representation.
+
+### 2.3 Proposed combined surface
+
+Assuming #91 merges first with §1's renames applied:
+
+| Function | Origin | What it does |
+|---|---|---|
+| `calculate_incidence(design, swap, spatial_factors)` | #91, generalised (§2.1) | treatment × factor-level counts, one matrix per factor |
+| `calculate_concurrence(design, swap, spatial_factors)` | #91, generalised | pairs co-occurring in the same factor level (`N Nᵀ`) |
+| `calculate_adjacency_concurrence(design, swap, ...)` | this branch, renamed | pairs that are spatially adjacent |
+| `calculate_info_matrix()`, `calculate_L_projection()`, `calculate_canonical_efficiency_factors()`, `cor_ar1*()` | #91 | information-theoretic tooling |
+
+This resolves the terminology collision in `PR-incidence-review.md` §3.1: "incidence" keeps its
+established design-theory meaning (treatment × block), "concurrence" means co-occurrence, and the
+adjacency variant is explicitly labelled as such. `calculate_pair_incidence` — which is neither an
+incidence matrix nor block concurrence — goes away.
+
+Suggested file layout, resolving §1.9's placement point:
+
+- `R/metrics.R` — `objective_function_info()` joins the other `objective_function_*`.
+- `R/information.R` — #91's information-matrix machinery and covariance constructors.
+- `R/incidence.R` — `calculate_incidence()`, `calculate_concurrence()`,
+  `calculate_adjacency_concurrence()`, and their print methods.
+- `R/design_utils.R` — `build_design_matrix()` (already there).
+
+### 2.4 One thing the shared docs must explain
+
+The two diagonals mean completely different things, and users will compare them:
+
+- `diag(calculate_concurrence(...))` is `sum_j n_ij²` — for a binary design, the **replication
+  number** of treatment *i*.
+- `diag(calculate_adjacency_concurrence(...))` is the **self-adjacency count** — the thing
+  `calculate_adjacency_score()` penalises, and ideally 0.
+
+One is "how many plots does this treatment have", the other is "how many times did it end up next to
+itself". Same shape, same-looking output, opposite interpretation — a shared `@details` section or a
+short vignette comparing them is worth more than any amount of per-function documentation.
+
+### 2.5 Suggested merge order
+
+1. **`build_design_matrix()` + adjacency orientation fix + NEWS "Bug Fixes"** — split out of
+   `feature/incidence`, no dependency on #91, fixes a live correctness bug. Merge first.
+2. **PR #91** with §1 applied — the renames (§1.1, §1.2), sorting helper (§1.3), roxygen (§1.4),
+   `@return` (§1.5), and a decision on the incremental-update contract (§1.6).
+3. **`calculate_incidence()` / `calculate_concurrence()` generalised to `spatial_factors`** (§2.1) —
+   small follow-up to #91.
+4. **`feature/incidence` reduced** to `calculate_adjacency_concurrence()` plus the print layer, with
+   `calculate_position_incidence()` dropped and the blockers in `PR-incidence-review.md` §2.1–2.3
+   addressed.
+
+Step 3 is where the two branches actually meet, and it's small — which is the argument for
+sequencing it this way rather than merging both features and reconciling afterwards.
diff --git a/PR-91-info-objective-review.md b/PR-91-info-objective-review.md
new file mode 100644
index 0000000..92bc931
--- /dev/null
+++ b/PR-91-info-objective-review.md
@@ -0,0 +1,193 @@
+# Code Review — PR #91 "Info objective"
+
+- **Repository:** biometryhub/speed
+- **PR:** [#91 "Info objective"](https://github.com/biometryhub/speed/pull/91)
+- **Author:** kaibagley (fork `kaibagley/speed`, branch `info-objective`)
+- **Base:** `main` · **Footprint:** 11 commits, 12 files, +1079/-10
+- **Reviewed locally on branch:** `info-objective`
+- **Reviewer / date:** Sam Rogers, 2026-06-12
+
+---
+
+## 1. Context
+
+`speed` rearranges treatments via simulated annealing to minimise a weighted sum of an
+**adjacency score** and a **balance score**. Objective functions plug into the SA loop
+(`speed_hierarchical` in `R/speed.R`) and must return `list(score = ..., ...)`; the returned
+list is fed back as `current_score_obj` on the next call so objectives can *incrementally*
+update from `swapped_items` rather than recompute from scratch.
+
+This PR adds a fundamentally different objective: instead of the adjacency proxy, it scores a
+design by the **Fisher information for treatment contrasts** under a user-supplied spatial
+covariance — directly targeting the Cramér–Rao lower bound of the intended analysis model.
+
+### What it adds
+
+New file [`R/objectives.R`](R/objectives.R) (434 lines):
+
+- `objective_function_info()` — computes the treatment information matrix
+  $I = X_1^\top L X_1$ after projecting out nuisance fixed effects, scored under
+  **A**-optimality ($\mathrm{tr}(I^-) = \sum 1/\lambda_i$) or contrast-space
+  **D**-optimality ($-\log|I| = -\sum \log\lambda_i$) over the $v-1$ positive eigenvalues.
+- `compute_L_projection()` — builds the projection
+  $L = \Sigma^{-1} - \Sigma^{-1}X_2(X_2^\top\Sigma^{-1}X_2)^{-1}X_2^\top\Sigma^{-1}$
+  from a user covariance $\Sigma$ and block design $X_2$.
+- Analysis helpers: `calc_info_matrix()`, `calc_incidence_matrix()`,
+  `calc_concurrence_matrix()`, `calculate_efficiency_factors()`.
+- Correlation builders: `cor_ar1()`, `cor_ar1_ar1()`.
+- Internal helpers `.compute_info()`, `.build_L_from_df()`, `.build_treatment_matrix()`.
+
+Plus [`tests/testthat/test-objective_function_info.R`](tests/testthat/test-objective_function_info.R)
+(378 lines), 8 new `NAMESPACE` exports, `man/*` docs, and `DESCRIPTION` author additions
+(Kai Bagley / Curtin University).
+
+---
+
+## 2. Overall assessment
+
+**Recommendation: request changes.**
+
+The concept and statistical core are sound and worth merging eventually. However, there is a
+**silent-correctness defect** (§4.1) that currently affects the package's own example and
+integration tests and must be resolved before merge, plus several secondary issues around
+performance, conventions, and test rigour.
+
+---
+
+## 3. Strengths
+
+- **Sound statistical core.** Rank-deficiency is handled correctly: $I$ is genuinely rank
+  $v-1$ (treatment columns sum to the intercept, which $L$ removes), so filtering to positive
+  eigenvalues and summing over them is the correct contrast-space formulation. Relative
+  eigenvalue tolerance (`max_eig * 1e-10`) and the disconnected-design penalty are the right
+  instincts.
+- **Identity-Σ consistency** holds by construction and is tested:
+  `compute_L_projection(df, diag(n))` reduces to `.build_L_from_df`.
+- **No new hard dependencies** — base `solve`/`eigen` only.
+- **Valuable direction** — optimising against the actual analysis covariance is a real
+  improvement over the adjacency proxy.
+- Helper-level tests are thorough and readable.
+
+---
+
+## 4. Correctness concerns
+
+### 4.1 🔴 `L`/`Σ` ordering is not tied to design coordinates — and `speed()` reorders rows
+
+**Severity: high (blocking).**
+
+`compute_L_projection` and `.compute_info` use `L_matrix` purely **positionally** by the row
+order of `layout_df` ([objectives.R:188](R/objectives.R#L188)); the `spatial_cols` argument is
+accepted but never used to align anything. Meanwhile `speed()` re-sorts the data by
+`(row, col)` before optimising:
+
+```r
+# R/speed.R:190 — runs whenever row/col are present
+#   (infer_row_col returns inferred = TRUE, design_utils.R:167-173)
+data <- data[do.call(order, data[c(row_column, col_column)]), ]
+```
+
+`initialise_design_df` emits **column-major** order (`expand.grid(row, col)` varies `row`
+fastest), but `speed()` re-sorts to **row-major**. For any grid with `nrows > 1` and
+`ncols > 1` these orderings differ, so the SA loop optimises against a **scrambled** $\Sigma$ —
+silently, with no error.
+
+This is not hypothetical: the headline example ([objectives.R:57-76](R/objectives.R#L57)) and
+the 4×6 integration test ([test:324](tests/testthat/test-objective_function_info.R#L324)) both
+hit it. The test passes only because it asserts `rank == 5` and `is.finite(A_value)`, both of
+which are insensitive to the permutation.
+
+**Fix direction:** thread `spatial_cols` through so `L`/`Σ` are keyed to coordinates (reorder
+`L` to match the design, or sort `Σ`'s implied order), or have `speed()` pass the row
+permutation to the objective. At minimum, validate and document the required ordering loudly.
+
+### 4.2 🟠 `spatial_cols` is dead in the new objective
+
+`objective_function_info` requires `spatial_cols` (it is how `speed()` calls every objective)
+but the body ignores it ([objectives.R:85](R/objectives.R#L85)). This is the root cause of
+§4.1 and is misleading on its own: the signature implies spatial awareness, but the spatial
+information actually lives in a separately hand-built `L`.
+
+### 4.3 🟡 Documented `L` formula is missing a term
+
+The roxygen ([objectives.R:32-33](R/objectives.R#L32)) drops the trailing
+$X_2^\top\Sigma^{-1}$ term. The PR description and the code
+([objectives.R:246](R/objectives.R#L246)) are correct; only the man page is wrong.
+
+---
+
+## 5. Integration with `speed()`
+
+- **Call contract: OK.** The loop calls
+  `obj_function(design, swap, spatial_cols, adj_weight=, bal_weight=, current_score_obj=, swapped_items=, ...)`
+  ([speed.R:289](R/speed.R#L289)). User-supplied `criterion` / `L_matrix` / `block_column`
+  arrive by name through `speed(...)`'s dots and match correctly; `adj_weight` / `bal_weight` /
+  `current_score_obj` / `swapped_items` fall harmlessly into `...`. `random_initialise` also
+  forwards `...` ([design_utils.R:401](R/design_utils.R#L401)). `rlang::check_dots_used()` is
+  satisfied.
+- **🟠 Ignores the incremental-state contract → performance.** Unlike
+  `objective_function_piepho`, this recomputes from scratch each iteration: rebuilds `X1` in an
+  R `for`-loop, forms `t(X1) %*% L %*% X1` ($O(n^2v)$ with a **dense** $n \times n$ `L`), and
+  runs `eigen` ($O(v^3)$) — ×10 000 iterations by default. A swap touches only two rows; the
+  info-matrix update is rank-≤4. For field-trial sizes this will be far slower than the
+  adjacency objective. Not incorrect, but a real scalability weakness.
+
+---
+
+## 6. API, packaging & conventions
+
+- **8 new exports in one PR.** Several (`calc_incidence_matrix`, `calc_concurrence_matrix`,
+  `calculate_efficiency_factors`) are not used by the objective — they are analysis
+  conveniences. Worth confirming scope with the author; each is now a maintenance and
+  `R CMD check` surface.
+- **New `R/objectives.R` vs. existing `R/metrics.R`** — the other `objective_function_*` all
+  live in `metrics.R`. Splitting is defensible but should be deliberate.
+- **🟡 `NEWS.md` not updated** — CLAUDE.md requires a bullet for user-facing changes; 8 new
+  exports qualify.
+- **🟡 Unrelated whitespace churn in `DESCRIPTION`** (trailing spaces stripped on `Authors@R:`,
+  `Imports:`, `Suggests:`, Description block). `CONTRIBUTING.md` says don't restyle unrelated
+  lines. The author/affiliation additions themselves are fine.
+- **Dangling `.Rbuildignore` entry** — adds `^\.dir-locals\.el$` though no such file is tracked
+  (the Emacs config correctly did not land; the ignore line is now dead config).
+- Reimplements a pseudo-inverse via `eigen` rather than reusing `pseudo_inverse` in
+  `R/utils.R` — minor consistency point.
+
+---
+
+## 7. Tests
+
+- Helper coverage is good (incidence, concurrence, eigenvalue properties, correlation builders,
+  input validation).
+- **🟠 The two `speed()` integration tests don't validate spatial correctness** — they would
+  pass even with the §4.1 misalignment, since they only check rank/finiteness. A meaningful
+  test should assert the optimised design beats a fixed/random arrangement *under the intended*
+  $\Sigma$, and should exercise ordering explicitly.
+- **🟡 `block_col =` vs `block_column`** — tests at
+  [test:48](tests/testthat/test-objective_function_info.R#L48), 130, 143, 188 pass the wrong
+  argument name and survive only via R partial-argument matching (`compute_L_projection` has no
+  `...`). Fragile: breaks under `options(warnPartialMatchArgs = TRUE)` or if another `block_*`
+  argument is added. Use the real name.
+- No test for `criterion = "D"` *through* `speed()` (only direct calls).
+
+---
+
+## 8. Minor / style
+
+- `A_val = …` / `D_val = …` use `=` for assignment inside braces
+  ([objectives.R:288-289](R/objectives.R#L288)) — slipped past the "fix lintr" commit.
+- `match.arg(criterion)` is good; consider validating length-1 when called directly.
+
+---
+
+## 9. Pre-merge checklist
+
+- [ ] **(blocking)** Resolve `L`/`Σ` ordering vs. `speed()` row reorder (§4.1); align via
+      `spatial_cols`.
+- [ ] Add a test that actually proves spatial optimality (would catch §4.1).
+- [ ] Honour the incremental-update contract, or consciously document waiving it given the
+      performance cost (§5).
+- [ ] Fix the documented `L` formula (§4.3).
+- [ ] Add a `NEWS.md` bullet (§6).
+- [ ] Revert unrelated `DESCRIPTION` whitespace changes (§6).
+- [ ] Fix `block_col` → `block_column` in tests (§7).
+- [ ] Confirm scope of the 8 new exports with the author (§6).
diff --git a/PR-incidence-review.md b/PR-incidence-review.md
new file mode 100644
index 0000000..dd5054e
--- /dev/null
+++ b/PR-incidence-review.md
@@ -0,0 +1,579 @@
+# Review: `feature/incidence`
+
+**Reviewer perspective:** user experience — is this something a biometrician can pick up and use
+correctly without reading the source?
+
+**Stated objective:** functions to calculate and print (a) the incidence of treatments appearing
+next to, or in the same row/column as, another replicate of the same treatment, and (b) the
+frequency of pairs of treatments appearing together.
+
+**Scope reviewed:** `main...feature/incidence` — 14 files, +615/−11. New `R/incidence.R`
+(186 lines), new internal `build_design_matrix()` in `R/design_utils.R`, refactor of
+`calculate_adjacency_score()` and `objective_function_piepho()` onto it, 2 new test files
+(217 lines), 3 new `man/` pages, `NEWS.md`, `DESCRIPTION` 0.0.8 → 0.0.9.
+
+**Method:** findings were derived by reading the code, then **verified at the console** against
+R 4.6.1 with `pkgload::load_all()`. Numbers quoted below are measured, not inferred.
+
+**Recommendation: request changes.** The maths in `calculate_pair_incidence()` is correct and the
+edge enumeration is clean and vectorised. But as a *user-facing feature* it is not there yet: the
+happy path fails for any design whose columns aren't literally named `treatment`/`row`/`col`, it
+returns silently wrong answers for multi-site designs, the "print" half of the objective is
+unimplemented, and the "pairs of treatments appearing together" half is answered by a *different
+open PR* using a different naming convention.
+
+---
+
+## 1. Summary of findings
+
+| # | Severity | Finding |
+|---|----------|---------|
+| 2.1 | 🔴 | Silently wrong on multi-site (MET) designs; the two functions disagree with each other |
+| 2.2 | 🔴 | Passing a `design` object only works if columns are named `treatment`/`row`/`col` |
+| 2.3 | 🔴 | "Print" is not implemented — output is unusable at realistic trial sizes |
+| 2.4 | 🟠 | Objective only half met: no pairwise co-occurrence, no block support |
+| 2.5 | 🟠 | Naming/semantic collision with PR #91 (`calc_incidence_matrix`, `calc_concurrence_matrix`) |
+| 3.1 | 🟠 | "Incidence" conflicts with its established meaning in design theory |
+| 3.2 | 🟠 | Cannot match `calculate_adjacency_score()`'s neighbourhood; docs overclaim the link |
+| 3.3 | 🟠 | Undocumented bug fix: `objective_function_piepho()` was scoring a scrambled grid inside `speed()`; adjacency becomes robust for direct calls |
+| 3.4 | 🟡 | `add_buffers()` output pollutes the incidence matrix, inconsistently |
+| 3.5 | 🟡 | Hot-loop performance regression in `build_design_matrix()` |
+| 3.6 | 🟡 | Duplicates `calculate_nb()`'s neighbour enumeration, and disagrees with it on NA |
+| 3.7 | 🟡 | `calculate_position_incidence()` docs overclaim the balance-score link |
+| 4.x | 🟡/⚪ | `as_list`, error messages, column ordering, test gaps, `.gitignore`, NEWS |
+
+---
+
+## 2. Blocking issues
+
+### 2.1 🔴 Silently wrong on multi-site designs — and the two functions disagree
+
+`build_design_matrix()` ([R/design_utils.R:648](R/design_utils.R#L648)) places values by
+coordinate:
+
+```r
+m[cbind(as_numeric_factor(df[[row_column]]),
+        as_numeric_factor(df[[col_column]]))] <- as.character(df[[swap]])
+```
+
+If two data-frame rows share a `(row, col)` coordinate, **the later one silently overwrites the
+earlier one**. No warning, no error.
+
+This is not hypothetical — it is the MET example in `speed()`'s own documentation
+([R/speed.R:118-122](R/speed.R#L118-L122)):
+
+```r
+df_site   <- initialise_design_df(1, 28, 5, 14, 5)      # row 1:28, col 1:5
+df_initial <- rbind(df_site, df_site, df_site, df_site, df_site)
+df_initial$site <- rep(c("a", "b", "c", "d", "e"), each = 140)
+```
+
+Five sites, each reusing `row` 1:28 and `col` 1:5. The design disambiguates spatially via
+`site_row`/`site_col`. A user who runs that documented example and then asks the obvious
+follow-up question —
+
+```r
+calculate_pair_incidence(result, swap = "lines")
+```
+
+— gets a 28 × 5 grid containing **site "e" only** (the last 140 rows to be assigned), reported as
+if it were the whole 700-plot trial. Every count is ~1/5 of the truth. For the exact question the
+user wants answered — "did any line end up next to itself?" — this returns a confidently wrong,
+reassuringly low number.
+
+Worse, the sibling function does something *different*: `calculate_position_incidence()` uses
+`table(trt_fac, df[[row_column]])` over the full data frame
+([R/incidence.R:171](R/incidence.R#L171)), so it **pools across all five sites** — row "1" mixes
+five different physical rows. So on the same input, one function silently drops 80% of the data and
+the other silently aggregates over a factor it doesn't know exists. Neither warns.
+
+**Suggested fix.** Two parts:
+
+1. `build_design_matrix()` should refuse to guess. `if (anyDuplicated(cbind(r, c))) stop(...)` with
+   a message naming the offending coordinates — a hard error here is much kinder than a wrong
+   matrix.
+2. Give both functions a `by` / `within` argument (e.g. `by = "site"`) that splits the design and
+   returns per-group results, or at minimum have them detect that `nrow(df) > max(row) * max(col)`
+   and error with a pointer to that argument.
+
+### 2.2 🔴 Passing a `design` object only works by coincidence
+
+Both functions accept a `design` object, which reads as "this integrates with `speed()`":
+
+```r
+df <- if (inherits(design, "design")) design$design_df else design
+```
+
+But they then fall back to `swap = "treatment"`, `row_column = "row"`, `col_column = "col"`. The
+`design` object records **neither** the swap column name nor the grid factors — its documented
+fields are `design_df`, `score`, `scores`, `temperatures`, `iterations_run`, `stopped_early`,
+`treatments` (the treatment *values*, not the column name) and `seed`
+([R/speed.R:46-58](R/speed.R#L46-L58)). So:
+
+```r
+result <- speed(df, swap = "variety", seed = 42)
+calculate_pair_incidence(result)
+#> Error: Column(s) not found in data: treatment
+```
+
+The user did everything right. The object knows the answer. They still get an error, and the error
+doesn't tell them that `swap =` is the argument to set.
+
+It's worse for the grid columns, because `speed()` is *actively more capable* here. `infer_row_col()`
+([R/design_utils.R](R/design_utils.R)) matches `^(col(umn|)|range)(s|)$` case-insensitively — so a
+field trial with a `range` column optimises fine, and even prints `row and range are used as row and
+column, respectively`. Then `calculate_pair_incidence(result)` errors on the missing `col`. The
+diagnostic tool is less capable than the optimiser it diagnoses.
+
+**Suggested fix.**
+
+- Add `swap` (the column name) and the resolved `grid_factors` to the `design` object in `speed()`,
+  and have both functions default to them when `inherits(design, "design")`. Then
+  `calculate_pair_incidence(result)` just works — which is surely the intended headline call.
+- Fall back to `infer_row_col(df, quiet = TRUE)` for plain data frames so `range` is supported
+  consistently.
+- Note this is invisible to the current tests because both "design object accepted as input" tests
+  use a hand-built `structure(list(design_df = df), class = "design")` *and* pass `swap = "trt"`
+  explicitly ([test-pair-incidence.R:86-97](tests/testthat/test-pair-incidence.R#L86-L97)). No test
+  ever calls either function against a real `speed()` result, and no test ever exercises the
+  defaults.
+
+### 2.3 🔴 The "print" half of the objective is not implemented
+
+The objective says *calculate and print*. There is no `print` method, no `summary` method, no
+`autoplot` method, and neither return value carries a class — `NAMESPACE` still has exactly two S3
+methods (`print.design`, `autoplot.design`).
+
+That matters because the raw containers don't scale. The MET example has 100 treatments:
+`calculate_pair_incidence()` returns a 100 × 100 integer matrix — 10,000 numbers, overwhelmingly
+zero — dumped to the console, wrapping across dozens of screens. `calculate_position_incidence()`
+returns a 100 × 28 and a 100 × 5 matrix, and to answer "which treatments share a row with
+themselves?" the user must visually scan 2,800 cells for values `> 1`.
+
+The functions produce the raw ingredients but leave the user to do the actual analysis. For the
+stated use case the interesting output is tiny:
+
+```
+Adjacency incidence: 3 of 100 treatments are self-adjacent
+  T14 (2 occurrences), T57 (1), T88 (1)
+Row incidence: 4 treatments appear more than once in a row
+  T3 (row 7, ×2), ...
+Most frequent neighbour pairs: T2–T44 (3), T9–T17 (3), ...
+```
+
+**Suggested fix.** Give the return values classes (`speed_pair_incidence`,
+`speed_position_incidence`) wrapping the matrices, with `print` methods that surface the
+violations and summarise the rest, plus a `summary`/`as.data.frame` for programmatic use. Keep the
+matrix accessible underneath so nothing is lost. An `autoplot` (heatmap for the pair matrix) would
+be a natural third step and fits the package's existing ggplot2 tooling.
+
+### 2.4 🟠 The objective is only half met
+
+Mapping the objective onto what shipped:
+
+| Stated goal | Delivered? |
+|---|---|
+| Treatments next to another rep of themselves | ✅ `diag(calculate_pair_incidence(...))` |
+| Treatments in the same **row/column** as another rep of themselves | ⚠️ Only indirectly — marginal counts, user must spot `> 1` |
+| Frequency of **pairs of treatments appearing together** | ❌ Only *adjacent* pairs |
+
+Two real gaps:
+
+- **No pairwise co-occurrence.** "Appearing together" in design language means co-occurring in the
+  same row / column / block — the concurrence matrix. `calculate_pair_incidence()` only counts
+  pairs that are *physically adjacent*; `calculate_position_incidence()` gives per-treatment
+  marginals, not pairs. So *nothing in this branch answers "how often do treatments i and j share a
+  row?"* — and that appears to be one of the three things asked for.
+- **No block support.** `calculate_position_incidence()` hardcodes exactly two dimensions, `row`
+  and `col`, while `speed()` optimises over arbitrary `spatial_factors` (the MET example uses
+  `~ site_col + site_block`). A user optimising `~ row + col + block` can get a report on two of
+  their three spatial factors. Blocks are the most common structure in the package's target
+  audience.
+
+**Suggested fix.** Replace the `row_column`/`col_column` pair with a `spatial_factors` argument
+mirroring `speed()` (accepting the same `~ row + col` formula), returning one matrix per factor.
+That fixes the block gap and the naming asymmetry in one move, and makes the function's output
+actually align with what was optimised.
+
+### 2.5 🟠 Head-on collision with PR #91
+
+`info-objective` (reviewed in `PR-91-info-objective-review.md`) exports
+`calc_incidence_matrix()` and `calc_concurrence_matrix()`
+([`R/objectives.R:316`, `:344` on that branch](R/objectives.R)):
+
+```r
+calc_incidence_matrix(layout_df, treatment_column = "treatment", block_column = "block")  # N
+calc_concurrence_matrix(layout_df, treatment_column = "treatment", block_column = "block") # N Nᵀ
+```
+
+If both branches merge, users get **four** overlapping functions across **two** naming conventions
+and **two** meanings of "incidence":
+
+| Function | Branch | Prefix | Treatment arg | Meaning of "incidence" |
+|---|---|---|---|---|
+| `calc_incidence_matrix` | #91 | `calc_` | `treatment_column` | treatment × block (classical `N`) |
+| `calc_concurrence_matrix` | #91 | `calc_` | `treatment_column` | pairs sharing a block |
+| `calculate_pair_incidence` | this | `calculate_` | `swap` | pairs that are adjacent |
+| `calculate_position_incidence` | this | `calculate_` | `swap` | treatment × row/col |
+
+Note that **`calc_concurrence_matrix()` already delivers §2.4's missing piece** — "pairs of
+treatments appearing together" — for blocks. And the two branches even sort treatments
+differently: #91 uses `factor(...)` (lexical: `T1, T10, T100, T2`) while this branch uses
+`stringi::stri_sort(..., numeric = TRUE)` (`T1, T2, ..., T10`). The same treatment set will print
+in different row orders depending on which function you call. That is a genuinely confusing user
+experience, and it's the kind of thing that is very hard to fix after release.
+
+**Suggested fix.** Reconcile the two PRs before either merges. Concretely: one prefix
+(`calculate_`, matching the existing `calculate_*` family), one treatment-column argument name
+(`swap`, matching the rest of the package), one sort helper, and an explicit split of
+responsibilities — e.g. `calculate_incidence()` for treatment × factor, `calculate_concurrence()`
+for pairs sharing a factor level, and `calculate_adjacency_concurrence()` for pairs that are
+spatially adjacent. This branch's `calculate_pair_incidence()` is the odd one out on naming and
+should probably be the one to move.
+
+---
+
+## 3. Significant issues
+
+### 3.1 🟠 "Incidence" means something else to this audience
+
+In experimental design, the **incidence matrix** is conventionally the treatment × block matrix
+`N`, and pairwise co-occurrence is the **concurrence matrix** `N Nᵀ`. For a GRDC-funded package out
+of a biometry group, a user seeing `calculate_pair_incidence()` will reasonably expect the
+classical concurrence matrix. They get adjacency counts instead. `calculate_position_incidence()`
+is closer to correct usage (treatment × row and treatment × column *are* incidence matrices) —
+which makes the inconsistency between the two names worse, not better.
+
+`calculate_pair_incidence()` is, precisely, an *adjacency concurrence* matrix. Naming it so — or at
+minimum stating the relationship to the standard terms in `@description` — would prevent a
+predictable misreading. See also §2.5.
+
+### 3.2 🟠 Can't reproduce the neighbourhood that was actually optimised
+
+`calculate_adjacency_score()` is parameterised by `ring_dists`, `ring_weights`, `ring_type`
+(`"manhattan"`/`"chebyshev"`) and `relationship`
+([R/calculate_adjacency_score.R:242-254](R/calculate_adjacency_score.R#L242-L254)), and
+`objective_function()` forwards all four through `...`
+([R/metrics.R:58-69](R/metrics.R#L58-L69)). `calculate_pair_incidence()` supports **none** of them —
+it hardcodes rook adjacency at distance 1.
+
+So a user who optimises with `ring_dists = c(1, 2)`, or `ring_type = "chebyshev"`, or a
+`relationship` matrix, gets an incidence report that describes a **different neighbourhood than the
+one the design was optimised for** — with nothing to indicate the mismatch. Same for
+`adj_weight = 0` (as in the MET example): the function will happily report adjacency counts for a
+design where adjacency was never part of the objective.
+
+This makes the doc claim overreach ([R/incidence.R:9-11](R/incidence.R#L9-L11)):
+
+> The diagonal holds self-adjacency counts — the quantity that `calculate_adjacency_score()`
+> penalises.
+
+True only at the defaults. Measured on one 3×3 design:
+
+| Quantity | Value |
+|---|---|
+| `sum(diag(M))` | 3 |
+| `calculate_adjacency_score()` (defaults) | 3 ✅ |
+| `calculate_adjacency_score(ring_type = "chebyshev")` | 5 |
+| `calculate_adjacency_score(ring_dists = c(1, 2), ring_weights = c(1, 1))` | 6 |
+
+So the relationship does hold exactly at the defaults — and the incidence report understates the
+optimised objective by up to 2× as soon as the user changes the neighbourhood, with nothing to
+signal it.
+
+The relationship is also **never tested**. There is no assertion anywhere that
+`sum(diag(M)) == calculate_adjacency_score(df, swap)`. That is the single highest-value test this
+branch could add: I confirmed it passes today, it pins the documented claim, and it would catch
+future drift between the two neighbour enumerations (§3.6).
+
+**Suggested fix.** Accept and forward `ring_dists`/`ring_weights`/`ring_type`, or scope the docs
+honestly to rook-distance-1 and say so in the description rather than implying equivalence.
+
+Related, pre-existing on `main` and worth fixing if these args get forwarded:
+`calculate_adjacency_score(d, "trt", ring_dists = c(1, 2))` **errors** —
+`adjacency_score_vec()` asserts `length(dists) == length(weights)` but `ring_weights` defaults to
+scalar `1`, so the documented default is unusable with multi-ring `ring_dists`. It should recycle.
+
+### 3.3 🟠 This branch silently fixes a real bug — in `objective_function_piepho()`, not adjacency
+
+> **Correction.** An earlier revision of this review claimed adjacency scoring was broken for all
+> non-square designs produced by `speed()`. That was wrong. `speed()` sorts its input at
+> [R/speed.R:188-192](R/speed.R#L188-L192) —
+> `data[do.call(order, data[c(row_column, col_column)]), ]`, with the comment *"Sort the data frame
+> to start with to ensure consistency in calculating the adjacency later"* — which produces
+> **row-major** order, exactly what `byrow = TRUE` expects. Verified: `matrix(byrow = TRUE)` on
+> sorted data is `identical()` to coordinate placement for 3×3, 4×3, 3×4, 6×4, 4×4 and 12×5. Designs
+> produced by `speed()` with numeric row/col columns were correct.
+
+The real story is more interesting. Two functions assumed **opposite** data orderings, and neither
+read coordinates:
+
+| Function | Fill | Assumes | Correct inside `speed()` (row-major)? | Correct on raw `initialise_design_df()` (column-major)? |
+|---|---|---|---|---|
+| `calculate_adjacency_score()` | `matrix(..., byrow = TRUE)` | row-major | ✅ | ❌ |
+| `objective_function_piepho()` | `matrix(...)` | column-major | ❌ | ✅ |
+
+Each is wrong precisely where the other is right. So the genuine bug this branch fixes is in
+**piepho**, which receives `speed()`'s row-major data and fills it column-major. Verified wrong for
+**every** shape tested — including square grids, so the transpose-invariance let-off doesn't apply.
+On a 4×3 design:
+
+```
+piepho old grid          correct grid
+  T1  T2  T3               T1  T1  T1
+  T1  T2  T4               T2  T2  T2
+  T1  T3  T4               T3  T3  T3
+  T2  T3  T4               T4  T4  T4
+```
+
+The multiset of treatments is preserved (so replication counts look right), but the spatial
+arrangement is scrambled — which is the entire input to `calculate_ed()`. Anyone using
+`obj_function = objective_function_piepho` has designs optimised against a layout that isn't theirs.
+
+Two narrower issues the fix also resolves:
+
+- **`calculate_adjacency_score()` is exported and doesn't compose with `initialise_design_df()`.**
+  The latter produces column-major data ([R/design_utils.R:304](R/design_utils.R#L304)), so a direct
+  call `calculate_adjacency_score(initialise_design_df(...), "treatment")` is wrong on `main` — a 4×3
+  design scores 0 when the answer is 8. The function's own examples happen to use hand-written
+  row-major data, so they pass. Two exported functions in the same package that silently disagree.
+- **The sort is applied to factors, so lexical level order can defeat it.** `to_factor()` runs at
+  [R/speed.R:185](R/speed.R#L185) *before* the sort, so a **character** row column with ≥10 rows gets
+  levels `1, 10, 11, 2, …` and `order()` follows them. Demonstrated: an 11×1 design with a genuine
+  like-pair at rows 1-2 scores **0** on `main` and **1** correctly. Integer row columns are
+  unaffected (`as.factor(1:12)` sorts numerically).
+
+Remaining process points, unchanged:
+
+- `NEWS.md` says nothing about any of this — two "Major Changes" bullets for the new diagnostics, no
+  "Bug Fixes" heading.
+- `build_design_matrix()` has **zero tests**, despite now underpinning `calculate_adjacency_score()`,
+  `objective_function_piepho()` and `calculate_pair_incidence()`.
+
+Also confirmed harmless: `matrix()` on a factor already returns a **character** matrix in R 4.6.1,
+so `calculate_ed()`'s type-sensitive
+`design_matrix[!(design_matrix %in% swapped_items)] <- NA`
+([R/metrics.R:388](R/metrics.R#L388)) is unaffected by the change.
+
+**Suggested fix.** See `PLAN-adjacency-orientation-fix.md` for the full plan. In short: a `NEWS.md`
+"Bug Fixes" entry naming `objective_function_piepho()` as the affected code path (piepho users need
+to regenerate designs; default-objective users do not), plus `test-build_design_matrix.R` and
+direct-call regression tests for `calculate_adjacency_score()`.
+
+### 3.4 🟡 Buffers pollute the incidence matrix, and inconsistently
+
+`create_buffers()` fills treatment columns with the literal value `"buffer"`
+([R/buffers.R:112](R/buffers.R#L112)). So on a buffered design,
+`calculate_pair_incidence()` reports `"buffer"` as a treatment, and `M["buffer", "buffer"]` is a
+large number — buffer plots are mostly adjacent to each other. This immediately falsifies the
+documented reassurance that "a well-optimised design will have zeros (or near-zeros) on the
+diagonal", and there is no argument to exclude anything.
+
+It's inconsistent, too. `add_buffers()` only fills a column named exactly `"treatment"` when
+`treatment_cols` is `NULL` ([R/buffers.R:99-107](R/buffers.R#L99-L107)); other columns are left
+`NA`. So a design with a `treatment` column gets a polluted incidence matrix, while the same design
+with a `variety` column gets `NA` buffers that the function correctly drops. Identical workflows,
+different answers, depending on a column name.
+
+**Suggested fix.** An `exclude` argument (defaulting to `"buffer"`, since that value is
+package-generated) is the cheap fix. Better: have `add_buffers()` record which plots are buffers in
+a dedicated column so downstream functions can filter reliably rather than string-matching.
+
+### 3.5 🟡 Hot-loop performance regression
+
+`build_design_matrix()` now runs once per SA iteration (default 10,000, per level, times
+`random_initialisation` restarts). Compare the coercion work:
+
+```r
+# before: 2 passes
+max(as_numeric_factor(df[[row_column]]))   # 1
+max(as_numeric_factor(df[[col_column]]))   # 2
+
+# after: 4 passes + an extra as.character + cbind alloc + scattered subassignment
+nr <- max(as_numeric_factor(df[[row_column]]))   # 1
+nc <- max(as_numeric_factor(df[[col_column]]))   # 2
+m[cbind(as_numeric_factor(df[[row_column]]),     # 3
+        as_numeric_factor(df[[col_column]]))] <- # 4
+  as.character(df[[swap]])
+```
+
+`as_numeric_factor` is `as.numeric(as.character(x))` ([R/utils.R:301](R/utils.R#L301)) — two full
+allocating passes each. That's roughly double the per-iteration coercion cost, plus a scatter-assign
+instead of a contiguous fill, in the innermost loop of a package called `speed`. And it is entirely
+avoidable: **the row and column vectors never change across iterations** — only the `swap` column
+does.
+
+**Suggested fix.** Hoist the invariant part. Compute `nr`, `nc` and the `cbind(r, c)` index matrix
+once per level (outside the SA loop) and pass them in, or memoise on the row/col columns. A
+low-level `build_design_matrix(df, swap, index = precomputed)` variant keeps the public behaviour
+while removing the per-iteration cost. Worth benchmarking a realistic case (700 plots × 10,000
+iterations) before and after — the fix is cheap enough that it's worth doing regardless.
+
+### 3.6 🟡 Duplicates `calculate_nb()`, and disagrees with it
+
+`calculate_nb()` already enumerates rook-adjacent pairs from a design matrix
+([R/metrics.R:274-285](R/metrics.R#L274-L285)):
+
+```r
+lefts <- design_matrix[, -ncol(design_matrix)]; rights <- design_matrix[, -1]
+tops  <- design_matrix[-nrow(design_matrix), ]; bottoms <- design_matrix[-1, ]
+lr_pairs <- paste(lefts, rights, sep = ","); tb_pairs <- paste(tops, bottoms, sep = ",")
+nb <- table(sorted_pairs)
+```
+
+That is the same computation as `calculate_pair_incidence()`, in a different container. The package
+now has two independent implementations of "enumerate adjacent treatment pairs" that can drift, and
+they **already disagree on NA**: `calculate_nb()` uses `paste()`, which stringifies `NA` into the
+literal pair `"NA,A"`, while `calculate_pair_incidence()` explicitly drops NA edges
+([R/incidence.R:96-97](R/incidence.R#L96-L97)). So `calculate_nb()$nb` and
+`calculate_pair_incidence()` give different answers for the same design with missing plots — and
+missing plots are exactly what buffered and irregular trials have.
+
+**Suggested fix.** Express one in terms of the other, or factor the edge enumeration into a single
+internal helper (`.adjacent_pairs(m)`) that both call, and settle the NA question once. Add a test
+asserting the two agree.
+
+### 3.7 🟡 `calculate_position_incidence()` docs overclaim
+
+> This is the human-readable decomposition of what `calculate_balance_score()` collapses to a
+> scalar.
+
+Three ways that's not quite right. `calculate_balance_score()`
+([R/metrics.R:164-175](R/metrics.R#L164-L175)) is `sum(rowVars(table(...)))` — a sum of **variances**
+of counts, not the counts; a user can't recover the score from these matrices without knowing that.
+It operates over arbitrary `spatial_cols`, not just row and col (§2.4). And it tabulates
+`table(layout_df[[el]], layout_df[[swap]])` — spatial as rows, treatment as columns — the
+**transpose** of what this function returns.
+
+**Suggested fix.** Either soften to "related to" and state the actual relationship, or return the
+per-factor variances alongside the counts so the claim becomes literally true.
+
+---
+
+## 4. Polish
+
+### 4.1 🟡 `as_list` is redundant API surface
+
+The parameter's own documentation concedes it
+([R/incidence.R:22-24](R/incidence.R#L22-L24)): "This is the matrix split by row — identical data,
+different container." A `t × t` matrix already supports `M["A", ]`, so `as_list = TRUE` adds a
+parameter to document, test and maintain forever for no capability. It's also asymmetric — only
+`calculate_pair_incidence()` has it, and its return type now depends on an argument value, which
+complicates every downstream `print`/`summary` method added for §2.3. Recommend dropping it before
+release; it's much easier to add later than remove.
+
+### 4.2 🟡 Error messages are thin and duplicated
+
+```r
+stop("Column(s) not found in data: ", paste(missing_cols, collapse = ", "), call. = FALSE)
+```
+
+`call. = FALSE` removes the only clue about *where* the error came from, so the user sees a bare
+`Error: Column(s) not found in data: treatment` with no function name, no indication that `swap =`
+is the argument to set, and no list of what *is* available. The package already imports `rlang`,
+which makes `cli`-style messages with a bullet list of available columns cheap. The 8-line block is
+also copy-pasted verbatim into both functions — worth a `.verify_incidence_cols()` helper alongside
+the existing `.verify_*` family in `R/verify_utils.R`.
+
+### 4.3 🟡 Position columns aren't sorted with the same care as treatments
+
+`calculate_pair_incidence()` carefully uses `stri_sort(..., numeric = TRUE)` for treatment levels,
+but `calculate_position_incidence()` passes the position column straight to `table()`
+([R/incidence.R:171-172](R/incidence.R#L171-L172)). If `row` is character (or a factor with default
+lexical levels), the columns come out `1, 10, 11, 12, 2, 3, ...` — so `res$row[, 2]` is row 10, not
+row 2. Given the branch commit "Working on using factors for spatial columns" this seems in scope.
+Apply `stri_sort(numeric = TRUE)` to the position levels too.
+
+### 4.4 🟡 Test gaps
+
+The 17 tests that exist are focused and correct — I verified the 3×3 Latin square case by hand and
+the expected all-pairs-`= 4` result is right. But the coverage has holes that map onto every
+blocking issue above:
+
+- No test asserts `sum(diag(M)) == calculate_adjacency_score(df, swap)` — the documented
+  relationship (§3.2). I verified it holds (both `3` on a 3×3 test case), so this is a free test.
+- No test uses a real `speed()` result; both "design object" tests use a hand-built fake **and**
+  pass `swap` explicitly, so the defaults and the `design`-object integration are untested (§2.2).
+- No test for duplicate coordinates / multi-site (§2.1), buffered designs (§3.4), or
+  `build_design_matrix()` at all (§3.3).
+- The "missing column gives informative error" tests use `expect_error(..., "col")` — which matches
+  the word "Column" in the message regardless of whether the missing column is named. The test
+  cannot fail for the reason it claims to check. Match on the actual column name.
+- No test with a factor `swap` column carrying unused levels: `unique(as.character(...))` drops
+  them, so a treatment with zero plots silently vanishes from the matrix, while `speed()`'s
+  `table()`-based balance score would have counted it.
+
+### 4.5 ⚪ `.gitignore` scope creep and risk
+
+```gitignore
+*\.RData
+*\.rds
+*\.txt
+*\.xlsx
+*\.csv
+```
+
+Two problems. The backslashes are regex escaping, not gitignore glob syntax — they happen to work
+(`\.` collapses to a literal `.`) but they're misleading and will get copied. More substantively,
+blanket-ignoring `*.csv`, `*.txt` and `*.rds` across an R package is risky: it will silently
+swallow `inst/extdata/` fixtures, test data, and `data-raw/` inputs, and the failure mode is a
+contributor whose file "didn't get committed" with no error. Scope these to the directories where
+scratch files actually accumulate, and drop them from this branch — they're unrelated to the
+feature.
+
+### 4.6 ⚪ NEWS, versioning, and docs integration
+
+- The two new functions are filed under **Major Changes**. Two additive diagnostic helpers that
+  change no existing behaviour read as **Minor Changes**; the genuinely major item in this branch is
+  the undocumented adjacency fix (§3.3), which isn't listed at all.
+- This branch is the one that added to `CLAUDE.md`: *"Link to an issue number or PR such as (#999)
+  where relevant."* Neither new bullet has a link.
+- Neither function is mentioned in `README.Rmd` or any vignette, and `?speed` doesn't `@seealso`
+  them. They're discoverable only via the pkgdown reference index (where `starts_with("calculate")`
+  does pick them up correctly — no `_pkgdown.yml` change needed). For a feature whose whole purpose
+  is *inspecting* a design, a few lines in the main vignette right after the `autoplot()` example is
+  where users will actually find it.
+- `CLAUDE.md`'s "Source layout (R/)" section isn't updated for the new `R/incidence.R`.
+- `setNames()` is used bare at [R/incidence.R:111](R/incidence.R#L111); it resolves only because
+  `R/buffers.R` and `R/metrics.R` happen to carry `@importFrom stats setNames`. `R/buffers.R` uses
+  the explicit `stats::setNames()`. Match that, or add the tag to this file.
+
+---
+
+## 5. What's good
+
+Worth saying plainly, because the critique above is long:
+
+- **The pair-counting logic is correct and well constructed.** `raw + t(raw)` with
+  `diag(M) <- diag(raw)` is the right way to fold ordered edge counts into a symmetric matrix while
+  keeping self-adjacency counted once, and it's fully vectorised rather than looped. The `nc >= 2` /
+  `nr >= 2` guards correctly handle single-row and single-column designs, and both are tested.
+- **`build_design_matrix()` is the right abstraction** and fixes a real orientation bug (§3.3). It
+  just needs tests, a duplicate-coordinate guard, and the invariant work hoisted out of the loop.
+- **NA handling in `calculate_pair_incidence()` is deliberate and tested** — dropping edges rather
+  than counting `"NA"` as a treatment is the correct choice (and the one `calculate_nb()` gets
+  wrong).
+- **Roxygen is genuinely good** — `@seealso` cross-links, runnable examples, `@inheritParams` to
+  avoid duplication, and the return values are documented in enough detail to use. The overclaims
+  in §3.2/§3.7 are the flip side of docs that bother to explain *why* the function exists, which is
+  more than most packages manage.
+
+---
+
+## 6. Suggested path to merge
+
+Minimum to unblock:
+
+1. **Document and test the adjacency orientation fix (§3.3).** This is the most consequential change
+   in the branch and currently the least documented. Arguably it should be split into its own PR and
+   merged first, so the fix isn't gated on the review of a new feature.
+2. Guard duplicate coordinates in `build_design_matrix()` and decide the multi-site story (§2.1).
+3. Store `swap` + resolved grid factors on the `design` object; default to them; fall back to
+   `infer_row_col()` (§2.2).
+4. Add `print` methods so the functions satisfy "print" and are usable at 100 treatments (§2.3).
+5. Reconcile naming and scope with PR #91 before either merges (§2.5, §3.1) — see
+   `PR-91-alignment-and-integration.md`.
+6. Add the `sum(diag(M)) == calculate_adjacency_score(...)` cross-check test (§3.2).
+
+Then, ideally in the same PR: `spatial_factors` instead of hardcoded row/col (§2.4), hoist the
+hot-loop coercions (§3.5), drop `as_list` (§4.1), and revert the `.gitignore` churn (§4.5).