diff --git a/CLAUDE.md b/CLAUDE.md
index 70177764..0a82a3c5 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -58,6 +58,8 @@ 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`. NEWS entries should be kept concise, with just a short 1-2 sentence summary of the changes, not paragraphs of explanation. If a change is related to a GitHub issue, it can be referenced just by number in parentheses e.g. (#1), but do not reference GitHub issues unless it's certain that they are related.
 
+Within a release, sections must appear in the order **Major Changes**, **Minor Changes**, **Bug Fixes** - regardless of how many entries each holds. Omit a section entirely if it is empty.
+
 ## Architecture
 
 ### Entry point and control flow
diff --git a/NEWS.md b/NEWS.md
index dbe25ce0..ca4b5c0c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,12 +3,34 @@
 ## Major Changes
 
 - Added a `summary()` method for `"design"` objects, reporting structure and replication, a
-  decomposed optimisation score, and design-quality diagnostics (connectedness, concurrence,
-  replicate spans and spread across blocks, neighbour balance, and opt-in efficiency).
+  decomposed optimisation score, and design-quality diagnostics.
   ([#73](https://github.com/biometryhub/speed/issues/73))
+- `grid_factors` gains an optional `by` element naming the column that separates a design into
+  several grids, e.g. `list(dim1 = "row", dim2 = "col", by = "site")` for a multi-environment trial.
+  Each grid is scored on its own.
+
+## Minor Changes
+
+- Designs whose `row`/`col` columns are not numeric, or where two plots share a coordinate, now fail
+  with a message naming the problem.
 
 ## Bug Fixes
 
+- Design metrics are now built from each plot's `row`/`col` coordinates rather than the order of the
+  rows in the data frame. Designs generated with `objective_function_piepho()` should be regenerated.
+- Multi-site designs are no longer scored as one pooled grid, which discarded plots whose coordinates
+  collided and counted adjacencies between sites. Use `grid_factors$by` to name the grouping column.
+- `objective_function_piepho()` now scores evenness of distribution per grid and reports each grid
+  separately. A grid with no treatment replicated within it contributes `0` rather than `Inf`.
+- `calculate_efficiency_factor()` now errors for a design whose treatment contrasts are not
+  estimable, instead of returning an impossible value above 1. The row-column model gained an
+  intercept, which does not change results that were already valid.
+- `summary()` no longer errors on designs that cannot be placed on a single grid; the affected
+  diagnostics report why they are unavailable instead.
+- `calculate_nb()` no longer errors on designs with missing plots when `pair_mapping` is not
+  supplied.
+- `calculate_adjacency_score()` now recycles a single `ring_weights` value across every entry of
+  `ring_dists`, so the default is usable with more than one ring.
 - `swap_all = TRUE` no longer changes the replication of a design when an earlier level has
   unbalanced a swap group mid-search. Only treatments with matching replication are exchanged.
 
diff --git a/R/calculate_adjacency_score.R b/R/calculate_adjacency_score.R
index 137f4210..1bdb7c44 100755
--- a/R/calculate_adjacency_score.R
+++ b/R/calculate_adjacency_score.R
@@ -137,6 +137,11 @@ adjacency_score_vec <- function(
   relationship = NULL
 ) {
   ring_type <- match.arg(ring_type)
+  # A single weight applies to every ring, so the scalar default stays usable
+  # with a multi-ring `dists`. Any other length mismatch is still an error.
+  if (length(weights) == 1L) {
+    weights <- rep(weights, length(dists))
+  }
   stopifnot(length(dists) == length(weights))
   nr <- nrow(design_matrix)
   nc <- ncol(design_matrix)
@@ -207,6 +212,16 @@ adjacency_score_vec <- function(
 #'   to `NULL`, which keeps the strict identity match. Pass the raw matrix
 #'   through `prep_relationship()` first; the score functions consume only
 #'   the prepped form.
+#' @param by Optional column name grouping plots into separate grids (e.g.
+#'   `"site"` for a multi-environment trial). Each grid is scored on its own and
+#'   the counts summed, so no adjacency is counted between plots at different
+#'   sites. `NULL` (default) treats the design as a single grid, which errors if
+#'   two plots share a coordinate.
+#' @param grid_index Optional pre-built list of indices from [grid_indices()],
+#'   passed to [build_design_matrix()] to skip coordinate validation. `speed()`
+#'   supplies one so the annealing loop does not revalidate every iteration;
+#'   leave it `NULL` for a one-off call. Supplying it ignores `by`, which the
+#'   indices already encode.
 #'
 #' @return A non-negative numeric value: the number of like-treatment edges
 #'   in the row/column adjacency graph.
@@ -248,23 +263,38 @@ calculate_adjacency_score <- function(
   ring_dists = 1,
   ring_weights = 1,
   ring_type = c("manhattan", "chebyshev"),
-  relationship = NULL
+  relationship = NULL,
+  by = NULL,
+  grid_index = 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
-  )
+  if (is.null(grid_index)) {
+    grid_index <- grid_indices(layout_df, row_column, col_column, by = by)
+  }
 
-  per_cell <- adjacency_score_vec(
-    design_matrix,
-    dists = ring_dists,
-    weights = ring_weights,
-    ring_type = ring_type,
-    relationship = relationship
+  # Adjacency counts edges, and no edge crosses a grid boundary, so summing per
+  # grid is exact rather than an approximation.
+  totals <- vapply(
+    grid_index,
+    function(g) {
+      design_matrix <- build_design_matrix(
+        layout_df[g$rows, , drop = FALSE],
+        swap,
+        row_column = row_column,
+        col_column = col_column,
+        index = g$index
+      )
+      per_cell <- adjacency_score_vec(
+        design_matrix,
+        dists = ring_dists,
+        weights = ring_weights,
+        ring_type = ring_type,
+        relationship = relationship
+      )
+      return(sum(per_cell) / 2)
+    },
+    numeric(1)
   )
-  return(sum(per_cell) / 2)
+  return(sum(totals))
 }
diff --git a/R/design_utils.R b/R/design_utils.R
index be2547e2..4130bb20 100644
--- a/R/design_utils.R
+++ b/R/design_utils.R
@@ -889,3 +889,236 @@ random_initialise <- function(design, optimise, seed = NULL, ...) {
 #' @rdname initialise_design_df
 #' @export
 initialize_design_df <- initialise_design_df
+
+#' Signal a Coordinate Problem with a Classed Condition
+#'
+#' @description
+#' Carries two phrasings of the same problem: `message`, for someone calling a
+#' metric directly, and `reason`, a fragment `summary()` reports in place of a
+#' metric it cannot compute. Both are defined at the throw site so they cannot
+#' drift, and the class lets callers dispatch without matching on message text.
+#'
+#' @param class Condition subclass naming the specific problem.
+#' @param reason Short phrase for a `summary()` field.
+#' @param ... Pasted together to form the message.
+#'
+#' @keywords internal
+.grid_stop <- function(class, reason, ...) {
+  stop(structure(
+    class = c(class, "speed_grid_error", "error", "condition"),
+    list(message = paste0(...), reason = reason, call = NULL)
+  ))
+}
+
+#' Validate a Design's Coordinates and Build its Grid Index
+#'
+#' @description
+#' Coerces and validates the `row_column`/`col_column` coordinates, returning
+#' everything [build_design_matrix()] needs to place plots on a grid: the
+#' two-column matrix index and the grid's dimensions.
+#'
+#' Split out from [build_design_matrix()] because it is the expensive half and
+#' the *invariant* half: during annealing only the treatment column changes, so
+#' the index can be built once per `speed()` run and reused every iteration.
+#'
+#' @param df A data frame with columns named by `row_column` and `col_column`.
+#' @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 list with `idx` (an `nrow(df)` x 2 integer matrix of grid
+#'   positions), `nrow` and `ncol` (the grid's dimensions), and `n` (the number
+#'   of plots the index was built for, used to detect a stale index).
+#'
+#' @keywords internal
+grid_index <- function(df, row_column = "row", col_column = "col") {
+  # Checked before coercion: absent columns would otherwise reach max() as empty
+  # vectors and yield -Inf dimensions with a warning, rather than saying what is
+  # wrong.
+  missing_cols <- setdiff(c(row_column, col_column), names(df))
+  if (length(missing_cols)) {
+    .grid_stop(
+      "speed_grid_missing",
+      "no row/column factors",
+      "Cannot place the design on a grid: no ",
+      paste0("`", missing_cols, "`", collapse = " or "),
+      " column."
+    )
+  }
+  # Coercion of non-numeric labels warns; the check below reports it properly.
+  rows <- suppressWarnings(as_numeric_factor(df[[row_column]]))
+  cols <- suppressWarnings(as_numeric_factor(df[[col_column]]))
+
+  if (anyNA(rows) || anyNA(cols)) {
+    .grid_stop(
+      "speed_grid_nonnumeric",
+      sprintf("`%s`/`%s` labels are not numeric", row_column, col_column),
+      "Cannot place the design on a grid: `",
+      row_column,
+      "` and `",
+      col_column,
+      "` must be numeric, or coercible to numeric."
+    )
+  }
+  # Used directly as matrix indices, so they must be positive whole numbers.
+  if (
+    any(rows < 1 | cols < 1) ||
+      any(rows != trunc(rows) | cols != trunc(cols))
+  ) {
+    .grid_stop(
+      "speed_grid_notinteger",
+      sprintf(
+        "`%s`/`%s` are not positive whole numbers",
+        row_column,
+        col_column
+      ),
+      "`",
+      row_column,
+      "` and `",
+      col_column,
+      "` must be positive whole numbers to index a grid."
+    )
+  }
+  idx <- cbind(rows, cols)
+  # Duplicated coordinates would silently overwrite each other. Multi-site
+  # designs reuse row/col per site, so they must be split before scoring.
+  if (anyDuplicated(idx)) {
+    .grid_stop(
+      "speed_grid_duplicate",
+      sprintf(
+        "duplicate `%s`/`%s` coordinates (e.g. a multi-site design)",
+        row_column,
+        col_column
+      ),
+      "Duplicate (",
+      row_column,
+      ", ",
+      col_column,
+      ") coordinates: the design cannot be placed on a single grid. ",
+      "Split multi-site designs by site first."
+    )
+  }
+
+  return(list(
+    idx = idx,
+    nrow = max(rows),
+    ncol = max(cols),
+    n = nrow(df)
+  ))
+}
+
+#' Split a Design into One Grid Index per Grid
+#'
+#' @description
+#' The multi-grid counterpart of [grid_index()]. A multi-environment trial is
+#' several grids that share a treatment set and never share an edge, so it
+#' cannot be one matrix: sites reuse `row`/`col`, and pooling them either
+#' silently overwrites plots or invents adjacencies between sites.
+#'
+#' With `by = NULL` this returns a one-element list, so callers have a single
+#' code path whether or not the design spans grids.
+#'
+#' @inheritParams grid_index
+#' @param by Optional column name grouping plots into grids (e.g. `"site"`).
+#'   `NULL` treats the design as one grid.
+#'
+#' @return A named list with one element per grid, each a list of `rows` (the
+#'   positions in `df` belonging to that grid) and `index` (that grid's
+#'   [grid_index()]). Named `"1"` when `by` is `NULL`.
+#'
+#' @keywords internal
+grid_indices <- function(
+  df,
+  row_column = "row",
+  col_column = "col",
+  by = NULL
+) {
+  if (is.null(by)) {
+    return(list("1" = list(
+      rows = seq_len(nrow(df)),
+      index = grid_index(df, row_column, col_column)
+    )))
+  }
+  if (!by %in% names(df)) {
+    .grid_stop(
+      "speed_grid_missing_by",
+      sprintf("no `%s` column to group grids by", by),
+      "Cannot split the design into grids: no `",
+      by,
+      "` column."
+    )
+  }
+  # drop = TRUE so an unused factor level does not produce an empty grid, which
+  # grid_index() would then reject for having no coordinates.
+  groups <- split(seq_len(nrow(df)), df[[by]], drop = TRUE)
+  out <- lapply(groups, function(rows) {
+    return(list(
+      rows = rows,
+      index = grid_index(df[rows, , drop = FALSE], row_column, col_column)
+    ))
+  })
+  # Carried so consumers can label per-grid output without being told separately
+  # which column the grids came from.
+  attr(out, "by") <- by
+  return(out)
+}
+
+#' 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`.
+#'
+#' Each plot's position comes from its own coordinates, so the row ordering of
+#' `df` is irrelevant, as is the level order of factor coordinate columns.
+#'
+#' Coordinates are used as-is, never renumbered: a gap in the coordinates is a
+#' real gap in the field (a missing plot, or a buffer that was removed), so
+#' collapsing it would make non-adjacent plots into neighbours. Callers must
+#' therefore cope with `NA` cells.
+#'
+#' @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"`).
+#' @param index Optional pre-built index from [grid_index()]. Supplying one skips
+#'   coordinate coercion and validation, which is the bulk of the work and is
+#'   invariant during annealing. `speed()` builds one per run; anything calling
+#'   this once should leave it `NULL`.
+#'
+#' @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",
+  index = NULL
+) {
+  if (is.null(index)) {
+    index <- grid_index(df, row_column, col_column)
+  } else if (!identical(index$n, nrow(df))) {
+    # Only catches a mismatched index when the plot count differs; callers still
+    # own keeping index and design in step.
+    stop(
+      "`index` was built for ",
+      index$n,
+      " plots but `df` has ",
+      nrow(df),
+      ". Rebuild it with `grid_index()`.",
+      call. = FALSE
+    )
+  }
+
+  design_matrix <- matrix(
+    NA_character_,
+    nrow = index$nrow,
+    ncol = index$ncol
+  )
+  design_matrix[index$idx] <- as.character(df[[swap]])
+  return(design_matrix)
+}
diff --git a/R/metrics.R b/R/metrics.R
index 02c20faf..03da17fa 100644
--- a/R/metrics.R
+++ b/R/metrics.R
@@ -58,7 +58,14 @@ objective_function <- function(layout_df,
   ring_args <- list(...)
   ring_args <- ring_args[intersect(
     names(ring_args),
-    c("ring_dists", "ring_weights", "ring_type", "relationship")
+    c(
+      "ring_dists",
+      "ring_weights",
+      "ring_type",
+      "relationship",
+      "by",
+      "grid_index"
+    )
   )]
   adj_score <- ifelse(adj_weight != 0,
     do.call(
@@ -195,6 +202,11 @@ calculate_balance_score <- function(layout_df, swap, spatial_cols) {
 #' @inheritParams calculate_ed
 #' @param design A data frame representing the spatial information of the design
 #' @param current_score_obj A named list containing the current score
+#' @param by,grid_index Optional grouping of plots into separate grids; see
+#'   [calculate_adjacency_score()]. Neighbour balance sums across grids, while
+#'   evenness of distribution is scored **per grid** and the scores summed,
+#'   reported both in total and per grid. A grid with no treatment replicated
+#'   inside it contributes `0`.
 #'
 #' @examples
 #' design_df <- initialise_design_df(
@@ -230,35 +242,98 @@ objective_function_piepho <- function(design,
                                       pair_mapping = NULL,
                                       row_column = "row",
                                       col_column = "col",
+                                      by = NULL,
+                                      grid_index = NULL,
                                       ...) {
-  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)
-  )
+  # `by`/`grid_index` are documented on calculate_adjacency_score()
+  if (is.null(grid_index)) {
+    grid_index <- grid_indices(design, row_column, col_column, by = by)
+  }
+
+  # Each grid is scored on its own: there is no distance between plots at
+  # different sites, so a pooled spanning tree would be meaningless.
+  ed <- list()
+  ed_scores <- setNames(numeric(length(grid_index)), names(grid_index))
+  nb_counts <- list()
+
+  for (nm in names(grid_index)) {
+    g <- grid_index[[nm]]
+    design_matrix <- build_design_matrix(
+      design[g$rows, , drop = FALSE],
+      swap,
+      row_column = row_column,
+      col_column = col_column,
+      index = g$index
+    )
+    ed[[nm]] <- calculate_ed(
+      design_matrix,
+      current_score_obj$ed[[nm]],
+      swapped_items
+    )
+    # A grid with nothing replicated has no spanning tree, so it contributes 0
+    # rather than 1/0, which would make the whole score `Inf`.
+    ed_scores[[nm]] <- if (length(ed[[nm]]) == 0L) {
+      0
+    } else {
+      1 /
+        sum(vapply(ed[[nm]], function(ed_rep) return(ed_rep$min_mst), numeric(1)))
+    }
+    nb_counts[[nm]] <- unlist(calculate_nb(design_matrix, pair_mapping)$nb)
+  }
+
+  # Summed rather than pooled into one reciprocal, so ED scales with adjacency -
+  # which also sums per grid - instead of shrinking as sites are added.
+  ed_score <- sum(ed_scores)
 
-  ed <- calculate_ed(design_matrix, current_score_obj$ed, swapped_items)
-  # sum(1/) or 1/sum
-  ed_score <- 1 / sum(vapply(ed, function(ed_rep) ed_rep$min_mst, numeric(1)))
-  nb <- calculate_nb(design_matrix, pair_mapping)
+  # Neighbour balance counts edges and no edge crosses a grid boundary, so the
+  # counts sum across grids before the variance is taken.
+  all_pairs <- unique(unlist(lapply(nb_counts, names)))
+  totals <- setNames(numeric(length(all_pairs)), all_pairs)
+  for (counts in nb_counts) {
+    totals[names(counts)] <- totals[names(counts)] + counts
+  }
+  nb <- list(
+    nb = totals,
+    max_nb = max(totals),
+    max_pairs = names(totals)[totals == max(totals)],
+    var = stats::var(totals)
+  )
   nb_score <- nb$var
 
-  design[[swap]] <- as.factor(design_matrix)
+  # Balance and adjacency read `design` directly: the treatment column and the
+  # coordinates must stay aligned, so neither takes a flattened `design_matrix`.
   bal_score <- calculate_balance_score(design, swap, spatial_cols)
-  adj_score <- calculate_adjacency_score(design, swap, row_column, col_column)
+  adj_score <- calculate_adjacency_score(
+    design,
+    swap,
+    row_column,
+    col_column,
+    grid_index = grid_index
+  )
+
+  # Reported alongside the total, not instead of it: per-grid values are not
+  # comparable between grids, so they are never ranked or averaged.
+  components <- c(
+    neighbour_balance = nb_score,
+    even_distribution = ed_score,
+    balance           = bal_score,
+    adjacency         = adj_score
+  )
+  if (length(grid_index) > 1) {
+    components <- c(
+      components,
+      setNames(ed_scores, paste0("even_distribution_", names(ed_scores)))
+    )
+  }
 
   return(list(
     score = round(nb_score + ed_score + bal_score + adj_score, 10),
     ed = ed,
+    ed_per_grid = ed_scores,
     bal = bal_score,
     adj = adj_score,
     nb = nb,
-    components = c(
-      neighbour_balance = nb_score,
-      even_distribution = ed_score,
-      balance           = bal_score,
-      adjacency         = adj_score
-    )
+    components = components
   ))
 }
 
@@ -337,26 +412,35 @@ calculate_nb <- function(design_matrix, pair_mapping = NULL) {
   for (row_ in 1:n_rows) {
     for (col_ in 1:n_cols) {
       node <- design_matrix[row_, col_]
+      # Empty cells (a missing plot, or a removed buffer) have no pairs to
+      # contribute, matching the pair_mapping path.
+      if (is.na(node)) {
+        next
+      }
       if (row_ < n_rows) {
         bottom <- design_matrix[row_ + 1, col_]
-        if (node < bottom) {
-          pair_str <- paste0(node, ",", bottom)
-        } else {
-          pair_str <- paste0(bottom, ",", node)
+        if (!is.na(bottom)) {
+          if (node < bottom) {
+            pair_str <- paste0(node, ",", bottom)
+          } else {
+            pair_str <- paste0(bottom, ",", node)
+          }
+
+          env_add_one(nb, pair_str)
         }
-
-        env_add_one(nb, pair_str)
       }
 
       if (col_ < n_cols) {
         right <- design_matrix[row_, col_ + 1]
-        if (node < right) {
-          pair_str <- paste0(node, ",", right)
-        } else {
-          pair_str <- paste0(right, ",", node)
+        if (!is.na(right)) {
+          if (node < right) {
+            pair_str <- paste0(node, ",", right)
+          } else {
+            pair_str <- paste0(right, ",", node)
+          }
+
+          env_add_one(nb, pair_str)
         }
-
-        env_add_one(nb, pair_str)
       }
     }
   }
@@ -654,15 +738,38 @@ create_pair_mapping <- function(items) {
 #' @param col_column Name of the column giving the column of the design (default: "col")
 #'
 #' @examples
+#' # `initialise_design_df()` fills `items` down columns, so the literal below is
+#' # column-major; the grid it produces is
+#' #   a b d c
+#' #   e a f b
+#' #   c f e d
 #' df_design <- initialise_design_df(c(
-#'   "a", "b", "d", "c",
-#'   "e", "a", "f", "b",
-#'   "c", "f", "e", "d"
+#'   "a", "e", "c",
+#'   "b", "a", "f",
+#'   "d", "f", "e",
+#'   "c", "b", "d"
 #' ), 3, 4)
 #'
 #' calculate_efficiency_factor(df_design, "treatment")
 #'
-#' @return A numeric value representing the efficiency factor of the design. Higher values indicate more efficient designs.
+#' # Not every design can support the estimate. Here each treatment fills one
+#' # grid row, so the treatment differences cannot be separated from the row
+#' # effects and there is no efficiency factor to report:
+#' #   a a a a
+#' #   b b b b
+#' #   c c c c
+#' confounded <- initialise_design_df(rep(c("a", "b", "c"), 4), 3, 4)
+#' try(calculate_efficiency_factor(confounded, "treatment"))
+#'
+#' @return A numeric value representing the efficiency factor of the design,
+#'   between 0 and 1. Higher values indicate more efficient designs.
+#'
+#'   Errors with a `speed_efficiency_rank` condition if the design cannot support
+#'   the estimate - that is, if some treatment contrast is not estimable once row
+#'   and column effects are eliminated, whether because too few residual degrees
+#'   of freedom remain or because a treatment is confounded with a row or column.
+#'   Such a design has no efficiency factor; before this check the formula
+#'   returned a plausible-looking value, usually above 1.
 #'
 #' @references Piepho, H. P., Williams, E., & Michel, V. (2015). Nonresolvable Row-Column Designs with an Even
 #'   Distribution of Treatment Replications. Journal of Agricultural, Biological, and Environmental Statistics,
@@ -677,38 +784,40 @@ calculate_efficiency_factor <- function(
 ) {
   item <- as.character(substitute(item))
 
+  # An efficiency factor is a property of one experiment, and several cannot be
+  # combined. Validated explicitly because pooled sites otherwise return a value
+  # above 1 rather than erroring; `summary()` reports one value per site.
+  grid_index(design_df, row_column, col_column)
+
   # Design parameters
   encoded_items <- as.integer(as.factor(design_df[[item]]))
   n_treatments <- length(unique(encoded_items))
-  n_rows <- max(as_numeric_factor(design_df[[row_column]]), na.rm = TRUE)
-  n_cols <- max(as_numeric_factor(design_df[[col_column]]), na.rm = TRUE)
+  rows <- as_numeric_factor(design_df[[row_column]])
+  cols <- as_numeric_factor(design_df[[col_column]])
+  n_rows <- max(rows, na.rm = TRUE)
+  n_cols <- max(cols, na.rm = TRUE)
   n_plots <- nrow(design_df)
 
   # Create design matrix X for treatments
   X <- matrix(0, nrow = n_plots, ncol = n_treatments)
-  for (i in 1:n_plots) {
-    X[i, encoded_items[i]] <- 1
-  }
+  X[cbind(seq_len(n_plots), encoded_items)] <- 1
 
-  # Create design matrix Z for rows and columns
-  # Row and col effects (excluding last row and col to avoid singularity)
+  # Create design matrix Z for rows and columns, indexed by each plot's own
+  # coordinates so the row ordering of `design_df` does not matter. Row and col
+  # effects exclude the last row and col to avoid singularity.
   Z_row <- matrix(0, nrow = n_plots, ncol = n_rows - 1)
   Z_col <- matrix(0, nrow = n_plots, ncol = n_cols - 1)
-  plot_index <- 1
-  for (i in 1:n_rows) {
-    for (j in 1:n_cols) {
-      if (i < n_rows) {
-        Z_row[plot_index, i] <- 1
-      }
-      if (j < n_cols) {
-        Z_col[plot_index, j] <- 1
-      }
-      plot_index <- plot_index + 1
-    }
-  }
-
-  # Combine row and column design matrices
-  Z <- cbind(Z_row, Z_col)
+  in_row <- which(rows < n_rows)
+  in_col <- which(cols < n_cols)
+  Z_row[cbind(in_row, rows[in_row])] <- 1
+  Z_col[cbind(in_col, cols[in_col])] <- 1
+
+  # Intercept, then row and column design matrices. The row-column model has a
+  # mean, and including it is what makes the estimability test below exact:
+  # without it the mean stays inside the treatment term (X's rows sum to 1) and
+  # no rank test on `A_RC` can separate estimable contrasts from inestimable
+  # ones. Reported values for estimable designs are unchanged.
+  Z <- cbind(1, Z_row, Z_col)
 
   # Check if Z^TZ is invertible
   ZtZ <- t(Z) %*% Z
@@ -727,6 +836,33 @@ calculate_efficiency_factor <- function(
   I_n <- diag(n_plots)
   A_RC <- t(X) %*% (I_n - P_Z) %*% X
 
+  # Rank n_treatments - 1 means every treatment contrast is estimable. Without
+  # this gate pseudo_inverse() drops the null directions and returns a
+  # plausible-looking value above 1 instead of failing. Rank catches both
+  # aliasing and too few residual degrees of freedom; counting degrees of
+  # freedom catches only the latter. The tolerance matches pseudo_inverse()'s so
+  # the two cannot disagree, and `qr()` is unusable here - its relative default
+  # reports full rank for eigenvalues 2, 9e-16, 6e-16.
+  if (sum(svd(A_RC)$d > 1e-10) != n_treatments - 1) {
+    stop(structure(
+      class = c(
+        "speed_efficiency_rank",
+        "speed_efficiency_error",
+        "error",
+        "condition"
+      ),
+      list(
+        message = paste0(
+          "Not all treatment contrasts are estimable after eliminating ",
+          "`", row_column, "` and `", col_column, "` effects, so this design ",
+          "cannot support an efficiency factor."
+        ),
+        reason = "treatment contrasts not estimable given row + col",
+        call = NULL
+      )
+    ))
+  }
+
   # Calculate Moore-Penrose inverse of A_RC, variance matrix
   V <- pseudo_inverse(A_RC)
 
diff --git a/R/speed.R b/R/speed.R
index 1b35bd5c..16ac2139 100644
--- a/R/speed.R
+++ b/R/speed.R
@@ -23,6 +23,14 @@
 #' @param grid_factors A named list specifying grid factors to construct a
 #'   matrix for calculating adjacency score, `dim1` for row and `dim2` for
 #'   column. (default: `list(dim1 = "row", dim2 = "col")`).
+#'
+#'   An optional third element, `by`, names a column that groups plots into
+#'   *separate* grids - a multi-environment trial, where each site reuses the
+#'   same `row`/`col` numbering. Each grid is then scored on its own and the
+#'   adjacency counts summed, so no adjacency is counted between plots at
+#'   different sites, e.g.
+#'   `list(dim1 = "row", dim2 = "col", by = "site")`. Without it, a design whose
+#'   sites share coordinates is refused rather than silently pooled.
 #' @param iterations Maximum number of iterations for the simulated annealing
 #'   algorithm (default: 10000). For hierarchical designs, can be a named list
 #'   with names matching `swap`.
@@ -183,6 +191,10 @@ speed <- function(data,
     }
   }
 
+  # `by` groups plots into separate grids (a multi-environment trial)
+  .verify_grid_by(data, grid_factors)
+  grid_by <- grid_factors$by
+
   # Infer row and column columns
   inferred <- infer_row_col(data, grid_factors, quiet)
   row_column <- inferred$row
@@ -193,7 +205,8 @@ speed <- function(data,
   data <- factored$df
 
   if (inferred$inferred) {
-    # Sort the data frame to start with to ensure consistency in calculating the adjacency later
+    # Metrics are built from each plot's coordinates now, but neighbour
+    # generation and plotting may still rely on row order, so the sort stays.
     data <- data[do.call(order, data[c(row_column, col_column)]), ]
     # Only reset row labels for base data frames; tibbles are positional and
     # warn on `rownames<-`, and nothing downstream reads the design's row names.
@@ -229,7 +242,8 @@ speed <- function(data,
 
   design <- do.call(speed_hierarchical, c(
     list(data = data, optimise = optimise, quiet = quiet, seed = seed,
-         row_column = row_column, col_column = col_column),
+         row_column = row_column, col_column = col_column,
+         grid_by = grid_by),
     dots
   ))
   # Set here, not passed through do.call(): do.call would evaluate a language
@@ -259,6 +273,20 @@ speed_hierarchical <- function(data, optimise, quiet, seed, ...) {
   current_design <- layout_df
   best_design <- current_design
 
+  # Only the treatment column moves during annealing, so build the index once.
+  # `NULL` on failure defers to build_design_matrix(), so a design that cannot
+  # form a grid still runs if its objective never needs one.
+  dots <- list(...)
+  grid_idx <- tryCatch(
+    grid_indices(
+      current_design,
+      dots$row_column %||% "row",
+      dots$col_column %||% "col",
+      by = dots$grid_by
+    ),
+    speed_grid_error = function(e) return(NULL)
+  )
+
   # Sequential optimisation for each hierarchy level
   all_scores <- list()
   all_temperatures <- list()
@@ -279,7 +307,7 @@ speed_hierarchical <- function(data, optimise, quiet, seed, ...) {
 
     # Calculate initial score for this level
     current_score_obj <- opt$obj_function(current_design, opt$swap, spatial_cols, adj_weight = adj_weight,
-                                          bal_weight = bal_weight, ...)
+                                          bal_weight = bal_weight, grid_index = grid_idx, ...)
     current_score <- current_score_obj$score
 
     if (!is.numeric(current_score)) {
@@ -313,7 +341,7 @@ speed_hierarchical <- function(data, optimise, quiet, seed, ...) {
       # Calculate new score
       new_score_obj <- opt$obj_function(new_design$design,opt$swap, spatial_cols, adj_weight = adj_weight,
                                         bal_weight = bal_weight, current_score_obj = current_score_obj,
-                                        swapped_items = new_design$swapped_items, ...)
+                                        swapped_items = new_design$swapped_items, grid_index = grid_idx, ...)
       new_score <- new_score_obj$score
 
       # Decide whether to accept the new design
@@ -379,7 +407,7 @@ speed_hierarchical <- function(data, optimise, quiet, seed, ...) {
     # report faithful score components.
     treatments[[level]] <- stringi::stri_sort(unique(as.vector(best_design[[opt$swap]])), numeric = TRUE)
     score_obj <- opt$obj_function(best_design, opt$swap, spatial_cols, adj_weight = adj_weight,
-                                  bal_weight = bal_weight, ...)
+                                  bal_weight = bal_weight, grid_index = grid_idx, ...)
     level_scores[level] <- score_obj$score
     per_level_meta[[level]] <- list(
       swap             = opt$swap,
@@ -402,6 +430,8 @@ speed_hierarchical <- function(data, optimise, quiet, seed, ...) {
     levels     = hierarchy_levels,
     row_column = .dots$row_column %||% "row",
     col_column = .dots$col_column %||% "col",
+    # NULL for a single-grid design, so summary() need not guess the grouping
+    grid_by    = .dots$grid_by,
     per_level  = per_level_meta
   )
 
diff --git a/R/summary.R b/R/summary.R
index 9f2dd415..212b33de 100644
--- a/R/summary.R
+++ b/R/summary.R
@@ -73,7 +73,12 @@
 #'
 #' @returns A list of class `"summary.design"`:
 #' - **hierarchical** - `TRUE` for a multi-level (e.g. split-plot) design.
-#' - **layout** - `n_plots`, `nrow`, `ncol`, `row_column`, `col_column`, `has_grid`.
+#' - **layout** - `n_plots`, `nrow`, `ncol`, `row_column`, `col_column`,
+#'   `has_grid` (`TRUE` when the design is reportable as a single grid), and
+#'   `grid_reason` (why not, or `NA`). `nrow`/`ncol` count the rows and columns
+#'   the design *occupies*, so a design with a gap in its coordinates (a missing
+#'   plot, or a removed buffer) reports fewer than the coordinates span. Both are
+#'   `NA` unless `has_grid`.
 #' - **levels** - character vector of level names (e.g. `"wp"`/`"sp"`; a single
 #'   name for a simple design).
 #' - **per_level** - one element per level (named by `levels`), each a list with:
@@ -135,14 +140,35 @@ summary.design <- function(
 
   want_neighbour <- is.null(neighbour) || isTRUE(neighbour)
 
-  has_grid <- all(c(rc, cc) %in% names(df))
+  # The one coordinate validation for the whole summary: `grid` is either a
+  # `grid_index()` list or the reason there isn't one.
+  #
+  # `has_grid` means "reportable as one grid", not merely "row/col columns
+  # exist": duplicate coordinates span several grids of differing shapes.
+  grid <- tryCatch(
+    grid_indices(df, row_column = rc, col_column = cc, by = meta$grid_by),
+    speed_grid_error = function(e) return(e$reason)
+  )
+  # A multi-environment trial occupies several grids that share a treatment set
+  # and never share an edge, so no single `nrow` x `ncol` describes it.
+  n_grids <- if (is.character(grid)) NA_integer_ else length(grid)
+  has_grid <- !is.character(grid) && n_grids == 1L
   layout <- list(
     n_plots = nrow(df),
     nrow = if (has_grid) length(unique(df[[rc]])) else NA_integer_,
     ncol = if (has_grid) length(unique(df[[cc]])) else NA_integer_,
     row_column = rc,
     col_column = cc,
-    has_grid = has_grid
+    has_grid = has_grid,
+    n_grids = n_grids,
+    grid_by = meta$grid_by,
+    grid_reason = if (has_grid) {
+      NA_character_
+    } else if (is.character(grid)) {
+      grid
+    } else {
+      sprintf("%d grids, grouped by `%s`", n_grids, meta$grid_by)
+    }
   )
 
   per_level <- lapply(levels, function(lv) {
@@ -190,7 +216,7 @@ summary.design <- function(
 
     # --- Evaluation metrics ---
     evaluation <- list(
-      replicate_span = .replicate_spans(df, swap, rc, cc),
+      replicate_span = .replicate_spans(df, swap, rc, cc, grid),
       connectedness = if (isFALSE(connectedness)) {
         list(
           available = FALSE,
@@ -222,19 +248,19 @@ summary.design <- function(
         .block_spread(df, swap, block)
       },
       efficiency = if (isTRUE(efficiency)) {
-        .efficiency_factor(df, swap, rc, cc)
+        .efficiency_factor(df, swap, rc, cc, grid)
       } else {
         list(
           available = FALSE,
           reason = "not requested (set efficiency = TRUE)"
         )
       },
+      # A design that cannot be gridded reports `grid`'s reason rather than
+      # erroring.
       neighbour = if (!want_neighbour) {
         list(available = FALSE, reason = "not requested (neighbour = FALSE)")
-      } else if (!has_grid) {
-        list(available = FALSE, reason = "no row/column factors")
       } else {
-        .neighbour_balance(df, swap, layout$nrow, layout$ncol)
+        .neighbour_balance(df, swap, rc, cc, grid)
       }
     )
 
@@ -537,9 +563,30 @@ print.summary.design <- function(x, ...) {
     cat(lab("Repl. span:"), rs$reason, "\n", sep = "")
   }
 
-  # Efficiency
+  # Efficiency. A multi-grid design reports one value per grid and no total:
+  # there is no meaningful way to combine them (see `.efficiency_factor()`).
   ef <- e$efficiency
-  if (isTRUE(ef$available)) {
+  if (!is.null(ef$per_grid)) {
+    cat(
+      lab("Efficiency:"),
+      sprintf(
+        "per %s (A-efficiency, row-column model)\n",
+        if (is.null(ef$grid_by)) "grid" else paste0("`", ef$grid_by, "`")
+      ),
+      sep = ""
+    )
+    for (nm in names(ef$per_grid)) {
+      one <- ef$per_grid[[nm]]
+      cat(
+        "    ",
+        format(nm, width = max(nchar(names(ef$per_grid)))),
+        "  ",
+        if (isTRUE(one$available)) fmt_num(one$value) else one$reason,
+        "\n",
+        sep = ""
+      )
+    }
+  } else if (isTRUE(ef$available)) {
     cat(
       lab("Efficiency:"),
       fmt_num(ef$value),
@@ -677,10 +724,23 @@ print.summary.design <- function(x, ...) {
 #' @param df Design data frame.
 #' @param swap Treatment column name.
 #' @param rc,cc Row and column column names.
+#' @param grid A [grid_index()] list, or a character reason there is no grid.
 #' @keywords internal
-.replicate_spans <- function(df, swap, rc, cc) {
-  if (!all(c(rc, cc) %in% names(df))) {
-    return(list(available = FALSE, reason = "no row/column factors"))
+.replicate_spans <- function(df, swap, rc, cc, grid) {
+  # Spans are distances within one grid; across grids they are meaningless
+  # (two sites' row 3 are not one plot apart), so refuse rather than pool.
+  if (is.character(grid)) {
+    return(list(available = FALSE, reason = grid))
+  }
+  if (length(grid) > 1L) {
+    return(list(
+      available = FALSE,
+      reason = sprintf(
+        "design spans %d grids (grouped by `%s`)",
+        length(grid),
+        attr(grid, "by")
+      )
+    ))
   }
   span1 <- function(x) {
     if (length(x) < 2) {
@@ -898,30 +958,57 @@ print.summary.design <- function(x, ...) {
 #' reason rather than erroring when its assumptions are not met.
 #'
 #' @param rc,cc Row and column column names.
+#' @param grid A [grid_index()] list, or a character reason there is no grid.
 #' @keywords internal
-.efficiency_factor <- function(df, swap, rc, cc) {
-  if (!all(c(rc, cc) %in% names(df))) {
-    return(list(available = FALSE, reason = "requires a row/column grid"))
+.efficiency_factor <- function(df, swap, rc, cc, grid) {
+  # Not just a guard against erroring: on duplicate coordinates
+  # calculate_efficiency_factor() pools the grids and returns an impossible >1.
+  if (is.character(grid)) {
+    return(list(available = FALSE, reason = grid))
   }
   if (length(unique(df[[swap]])) < 3) {
     return(list(available = FALSE, reason = "requires >= 3 treatments"))
   }
-  ef <- tryCatch(
-    eval(bquote(calculate_efficiency_factor(
-      df,
-      .(as.name(swap)),
-      row_column = rc,
-      col_column = cc
-    ))),
-    error = function(e) return(NULL)
-  )
-  if (is.null(ef) || !is.finite(ef)) {
-    return(list(
-      available = FALSE,
-      reason = "could not be computed for this design"
-    ))
+
+  one <- function(sub) {
+    ef <- tryCatch(
+      eval(bquote(calculate_efficiency_factor(
+        sub,
+        .(as.name(swap)),
+        row_column = rc,
+        col_column = cc
+      ))),
+      # A rank failure is a property of the design, so it carries its own
+      # reason; anything else is unexpected and gets the generic one.
+      speed_efficiency_error = function(e) return(e$reason),
+      error = function(e) return(NULL)
+    )
+    if (is.character(ef)) {
+      return(list(available = FALSE, reason = ef))
+    }
+    if (is.null(ef) || !is.finite(ef)) {
+      return(list(
+        available = FALSE,
+        reason = "could not be computed for this design"
+      ))
+    }
+    return(list(available = TRUE, value = ef))
+  }
+
+  if (length(grid) == 1L) {
+    return(one(df))
   }
-  return(list(available = TRUE, value = ef))
+
+  # One value per grid, never summed or averaged: averaging gives a different
+  # quantity from the combined analysis, which is not identified at design time
+  # anyway. Each grid is gated on its own rank, so one unreplicated site reports
+  # its reason without withholding the others.
+  per_grid <- lapply(grid, function(g) return(one(df[g$rows, , drop = FALSE])))
+  return(list(
+    available = any(vapply(per_grid, function(x) x$available, logical(1))),
+    per_grid = per_grid,
+    grid_by = attr(grid, "by")
+  ))
 }
 
 #' Neighbour-balance diagnostics
@@ -938,22 +1025,40 @@ print.summary.design <- function(x, ...) {
 #' whereas a distinct pair that never neighbours is an imbalance. Lumping them
 #' together hides self-adjacency behind the same `min 0` as the harmless case.
 #'
-#' Takes `nrow`/`ncol` from the caller's `layout` (counted via
-#' `length(unique(...))`) rather than deriving them from `max(row)`/`max(col)`:
-#' buffer plots (`add_buffers()`) can shift row/col numbering so it no longer
-#' starts at 1, which would otherwise reshape the grid with the wrong
-#' dimensions. Assumes `rc`/`cc` are present in `df`; callers should check
-#' `has_grid` first (see `summary.design()`).
+#' The grid comes from [build_design_matrix()], which places each plot at its own
+#' `rc`/`cc` coordinates, so the counts describe the layout whatever order `df`
+#' is in. Coordinates are read as-is, so plots separated by a buffer row or
+#' column (`add_buffers()` offsets and scales them) keep that separation and are
+#' not counted as neighbours.
 #'
+#' A design that cannot be placed on one grid is reported as unavailable rather
+#' than propagating [grid_index()]'s error out of `summary()`.
+#'
+#' @param rc,cc Row and column column names.
+#' @param grid A [grid_index()] list, reused as the [build_design_matrix()]
+#'   index, or a character reason there is no grid.
 #' @keywords internal
-.neighbour_balance <- function(df, swap, nrow, ncol) {
-  dm <- matrix(df[[swap]], nrow = nrow, ncol = ncol)
+.neighbour_balance <- function(df, swap, rc, cc, grid) {
+  if (is.character(grid)) {
+    return(list(available = FALSE, reason = grid))
+  }
+  # One pair mapping for the whole design, so a pair absent from one site counts
+  # as zero rather than being dropped. Counts sum across grids: no edge crosses
+  # a grid boundary.
   pair_mapping <- create_pair_mapping(df[[swap]])
-  nb <- calculate_nb(dm, pair_mapping)
-
   all_pairs <- unique(pair_mapping)
   counts <- setNames(rep(0L, length(all_pairs)), all_pairs)
-  counts[names(nb$nb)] <- nb$nb
+  for (g in grid) {
+    dm <- build_design_matrix(
+      df[g$rows, , drop = FALSE],
+      swap,
+      row_column = rc,
+      col_column = cc,
+      index = g$index
+    )
+    nb <- calculate_nb(dm, pair_mapping)
+    counts[names(nb$nb)] <- counts[names(nb$nb)] + nb$nb
+  }
 
   # create_pair_mapping() keys are "trt1,trt2"; a self-pair repeats the level.
   parts <- strsplit(names(counts), ",", fixed = TRUE)
diff --git a/R/verify_utils.R b/R/verify_utils.R
index 7dc4518e..4befa148 100644
--- a/R/verify_utils.R
+++ b/R/verify_utils.R
@@ -196,6 +196,29 @@
   }
 }
 
+#' Verify the `by` element of `grid_factors`
+#'
+#' @description
+#' `grid_factors` is a plain list, so a mistyped `by` would be ignored and every
+#' grid silently pooled. Checked before any optimisation happens.
+#'
+#' @inheritParams speed
+#'
+#' @rdname verify
+#'
+#' @keywords internal
+.verify_grid_by <- function(data, grid_factors) {
+  grid_by <- grid_factors$by
+  if (!is.null(grid_by)) {
+    if (!is.character(grid_by) || length(grid_by) != 1) {
+      data_type_error("grid_factors$by", "a single column name")
+    }
+    verify_column_exists(grid_by, data, "grid grouping column")
+  }
+
+  return(invisible(NULL))
+}
+
 
 # Other functions for verifying
 
diff --git a/man/build_design_matrix.Rd b/man/build_design_matrix.Rd
new file mode 100644
index 00000000..4a4d8ea0
--- /dev/null
+++ b/man/build_design_matrix.Rd
@@ -0,0 +1,47 @@
+% 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_column = "row",
+  col_column = "col",
+  index = NULL
+)
+}
+\arguments{
+\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_column}{Column name of the row position variable (default \code{"row"}).}
+
+\item{col_column}{Column name of the column position variable
+(default \code{"col"}).}
+
+\item{index}{Optional pre-built index from \code{\link[=grid_index]{grid_index()}}. Supplying one skips
+coordinate coercion and validation, which is the bulk of the work and is
+invariant during annealing. \code{speed()} builds one per run; anything calling
+this once should leave it \code{NULL}.}
+}
+\value{
+A character matrix of dimensions \code{max(row)} by \code{max(col)}.
+}
+\description{
+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
+\code{max(row)} by \code{max(col)}. Cells with no corresponding row in \code{df} are \code{NA}.
+
+Each plot's position comes from its own coordinates, so the row ordering of
+\code{df} is irrelevant, as is the level order of factor coordinate columns.
+
+Coordinates are used as-is, never renumbered: a gap in the coordinates is a
+real gap in the field (a missing plot, or a buffer that was removed), so
+collapsing it would make non-adjacent plots into neighbours. Callers must
+therefore cope with \code{NA} cells.
+}
+\keyword{internal}
diff --git a/man/calculate_adjacency_score.Rd b/man/calculate_adjacency_score.Rd
index 0eac7010..4dd73db3 100644
--- a/man/calculate_adjacency_score.Rd
+++ b/man/calculate_adjacency_score.Rd
@@ -12,7 +12,9 @@ calculate_adjacency_score(
   ring_dists = 1,
   ring_weights = 1,
   ring_type = c("manhattan", "chebyshev"),
-  relationship = NULL
+  relationship = NULL,
+  by = NULL,
+  grid_index = NULL
 )
 }
 \arguments{
@@ -41,6 +43,18 @@ otherwise. NA-padded cells off the design edge contribute \code{0}. Defaults
 to \code{NULL}, which keeps the strict identity match. Pass the raw matrix
 through \code{prep_relationship()} first; the score functions consume only
 the prepped form.}
+
+\item{by}{Optional column name grouping plots into separate grids (e.g.
+\code{"site"} for a multi-environment trial). Each grid is scored on its own and
+the counts summed, so no adjacency is counted between plots at different
+sites. \code{NULL} (default) treats the design as a single grid, which errors if
+two plots share a coordinate.}
+
+\item{grid_index}{Optional pre-built list of indices from \code{\link[=grid_indices]{grid_indices()}},
+passed to \code{\link[=build_design_matrix]{build_design_matrix()}} to skip coordinate validation. \code{speed()}
+supplies one so the annealing loop does not revalidate every iteration;
+leave it \code{NULL} for a one-off call. Supplying it ignores \code{by}, which the
+indices already encode.}
 }
 \value{
 A non-negative numeric value: the number of like-treatment edges
diff --git a/man/calculate_efficiency_factor.Rd b/man/calculate_efficiency_factor.Rd
index 06beb08d..b2591d43 100644
--- a/man/calculate_efficiency_factor.Rd
+++ b/man/calculate_efficiency_factor.Rd
@@ -21,20 +21,43 @@ calculate_efficiency_factor(
 \item{col_column}{Name of the column giving the column of the design (default: "col")}
 }
 \value{
-A numeric value representing the efficiency factor of the design. Higher values indicate more efficient designs.
+A numeric value representing the efficiency factor of the design,
+between 0 and 1. Higher values indicate more efficient designs.
+
+Errors with a \code{speed_efficiency_rank} condition if the design cannot support
+the estimate - that is, if some treatment contrast is not estimable once row
+and column effects are eliminated, whether because too few residual degrees
+of freedom remain or because a treatment is confounded with a row or column.
+Such a design has no efficiency factor; before this check the formula
+returned a plausible-looking value, usually above 1.
 }
 \description{
 Calculates an efficiency factor of a design according to Piepho 2015.
 }
 \examples{
+# `initialise_design_df()` fills `items` down columns, so the literal below is
+# column-major; the grid it produces is
+#   a b d c
+#   e a f b
+#   c f e d
 df_design <- initialise_design_df(c(
-  "a", "b", "d", "c",
-  "e", "a", "f", "b",
-  "c", "f", "e", "d"
+  "a", "e", "c",
+  "b", "a", "f",
+  "d", "f", "e",
+  "c", "b", "d"
 ), 3, 4)
 
 calculate_efficiency_factor(df_design, "treatment")
 
+# Not every design can support the estimate. Here each treatment fills one
+# grid row, so the treatment differences cannot be separated from the row
+# effects and there is no efficiency factor to report:
+#   a a a a
+#   b b b b
+#   c c c c
+confounded <- initialise_design_df(rep(c("a", "b", "c"), 4), 3, 4)
+try(calculate_efficiency_factor(confounded, "treatment"))
+
 }
 \references{
 Piepho, H. P., Williams, E., & Michel, V. (2015). Nonresolvable Row-Column Designs with an Even
diff --git a/man/create_speed_input.Rd b/man/create_speed_input.Rd
index b67f0b36..ef1fcc4c 100644
--- a/man/create_speed_input.Rd
+++ b/man/create_speed_input.Rd
@@ -37,7 +37,15 @@ consider for balance (default: \code{~row + col}).}
 
 \item{grid_factors}{A named list specifying grid factors to construct a
 matrix for calculating adjacency score, \code{dim1} for row and \code{dim2} for
-column. (default: \code{list(dim1 = "row", dim2 = "col")}).}
+column. (default: \code{list(dim1 = "row", dim2 = "col")}).
+
+An optional third element, \code{by}, names a column that groups plots into
+\emph{separate} grids - a multi-environment trial, where each site reuses the
+same \code{row}/\code{col} numbering. Each grid is then scored on its own and the
+adjacency counts summed, so no adjacency is counted between plots at
+different sites, e.g.
+\code{list(dim1 = "row", dim2 = "col", by = "site")}. Without it, a design whose
+sites share coordinates is refused rather than silently pooled.}
 
 \item{iterations}{Maximum number of iterations for the simulated annealing
 algorithm (default: 10000). For hierarchical designs, can be a named list
diff --git a/man/dot-efficiency_factor.Rd b/man/dot-efficiency_factor.Rd
index 5f5d8bbd..aebb671b 100644
--- a/man/dot-efficiency_factor.Rd
+++ b/man/dot-efficiency_factor.Rd
@@ -4,10 +4,12 @@
 \alias{.efficiency_factor}
 \title{A-efficiency factor (opt-in wrapper)}
 \usage{
-.efficiency_factor(df, swap, rc, cc)
+.efficiency_factor(df, swap, rc, cc, grid)
 }
 \arguments{
 \item{rc, cc}{Row and column column names.}
+
+\item{grid}{A \code{\link[=grid_index]{grid_index()}} list, or a character reason there is no grid.}
 }
 \description{
 Thin guarded wrapper over \code{\link[=calculate_efficiency_factor]{calculate_efficiency_factor()}} (a row--column model
diff --git a/man/dot-grid_stop.Rd b/man/dot-grid_stop.Rd
new file mode 100644
index 00000000..63fa6921
--- /dev/null
+++ b/man/dot-grid_stop.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/design_utils.R
+\name{.grid_stop}
+\alias{.grid_stop}
+\title{Signal a Coordinate Problem with a Classed Condition}
+\usage{
+.grid_stop(class, reason, ...)
+}
+\arguments{
+\item{class}{Condition subclass naming the specific problem.}
+
+\item{reason}{Short phrase for a \code{summary()} field.}
+
+\item{...}{Pasted together to form the message.}
+}
+\description{
+Carries two phrasings of the same problem: \code{message}, for someone calling a
+metric directly, and \code{reason}, a fragment \code{summary()} reports in place of a
+metric it cannot compute. Both are defined at the throw site so they cannot
+drift, and the class lets callers dispatch without matching on message text.
+}
+\keyword{internal}
diff --git a/man/dot-neighbour_balance.Rd b/man/dot-neighbour_balance.Rd
index f155f651..4ef83142 100644
--- a/man/dot-neighbour_balance.Rd
+++ b/man/dot-neighbour_balance.Rd
@@ -4,7 +4,13 @@
 \alias{.neighbour_balance}
 \title{Neighbour-balance diagnostics}
 \usage{
-.neighbour_balance(df, swap, nrow, ncol)
+.neighbour_balance(df, swap, rc, cc, grid)
+}
+\arguments{
+\item{rc, cc}{Row and column column names.}
+
+\item{grid}{A \code{\link[=grid_index]{grid_index()}} list, reused as the \code{\link[=build_design_matrix]{build_design_matrix()}}
+index, or a character reason there is no grid.}
 }
 \description{
 Builds the treatment grid and counts how often each treatment pair ends up
@@ -20,11 +26,13 @@ zero self-adjacency is the desirable outcome the optimiser works towards,
 whereas a distinct pair that never neighbours is an imbalance. Lumping them
 together hides self-adjacency behind the same \verb{min 0} as the harmless case.
 
-Takes \code{nrow}/\code{ncol} from the caller's \code{layout} (counted via
-\code{length(unique(...))}) rather than deriving them from \code{max(row)}/\code{max(col)}:
-buffer plots (\code{add_buffers()}) can shift row/col numbering so it no longer
-starts at 1, which would otherwise reshape the grid with the wrong
-dimensions. Assumes \code{rc}/\code{cc} are present in \code{df}; callers should check
-\code{has_grid} first (see \code{summary.design()}).
+The grid comes from \code{\link[=build_design_matrix]{build_design_matrix()}}, which places each plot at its own
+\code{rc}/\code{cc} coordinates, so the counts describe the layout whatever order \code{df}
+is in. Coordinates are read as-is, so plots separated by a buffer row or
+column (\code{add_buffers()} offsets and scales them) keep that separation and are
+not counted as neighbours.
+
+A design that cannot be placed on one grid is reported as unavailable rather
+than propagating \code{\link[=grid_index]{grid_index()}}'s error out of \code{summary()}.
 }
 \keyword{internal}
diff --git a/man/dot-replicate_spans.Rd b/man/dot-replicate_spans.Rd
index 6374a0ce..9b8f6916 100644
--- a/man/dot-replicate_spans.Rd
+++ b/man/dot-replicate_spans.Rd
@@ -4,7 +4,7 @@
 \alias{.replicate_spans}
 \title{Replicate spatial spans}
 \usage{
-.replicate_spans(df, swap, rc, cc)
+.replicate_spans(df, swap, rc, cc, grid)
 }
 \arguments{
 \item{df}{Design data frame.}
@@ -12,6 +12,8 @@
 \item{swap}{Treatment column name.}
 
 \item{rc, cc}{Row and column column names.}
+
+\item{grid}{A \code{\link[=grid_index]{grid_index()}} list, or a character reason there is no grid.}
 }
 \description{
 For each treatment, the minimum Manhattan separation between its replicate
diff --git a/man/grid_index.Rd b/man/grid_index.Rd
new file mode 100644
index 00000000..ab205444
--- /dev/null
+++ b/man/grid_index.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/design_utils.R
+\name{grid_index}
+\alias{grid_index}
+\title{Validate a Design's Coordinates and Build its Grid Index}
+\usage{
+grid_index(df, row_column = "row", col_column = "col")
+}
+\arguments{
+\item{df}{A data frame with columns named by \code{row_column} and \code{col_column}.}
+
+\item{row_column}{Column name of the row position variable (default \code{"row"}).}
+
+\item{col_column}{Column name of the column position variable
+(default \code{"col"}).}
+}
+\value{
+A list with \code{idx} (an \code{nrow(df)} x 2 integer matrix of grid
+positions), \code{nrow} and \code{ncol} (the grid's dimensions), and \code{n} (the number
+of plots the index was built for, used to detect a stale index).
+}
+\description{
+Coerces and validates the \code{row_column}/\code{col_column} coordinates, returning
+everything \code{\link[=build_design_matrix]{build_design_matrix()}} needs to place plots on a grid: the
+two-column matrix index and the grid's dimensions.
+
+Split out from \code{\link[=build_design_matrix]{build_design_matrix()}} because it is the expensive half and
+the \emph{invariant} half: during annealing only the treatment column changes, so
+the index can be built once per \code{speed()} run and reused every iteration.
+}
+\keyword{internal}
diff --git a/man/grid_indices.Rd b/man/grid_indices.Rd
new file mode 100644
index 00000000..8c4ec23c
--- /dev/null
+++ b/man/grid_indices.Rd
@@ -0,0 +1,34 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/design_utils.R
+\name{grid_indices}
+\alias{grid_indices}
+\title{Split a Design into One Grid Index per Grid}
+\usage{
+grid_indices(df, row_column = "row", col_column = "col", by = NULL)
+}
+\arguments{
+\item{df}{A data frame with columns named by \code{row_column} and \code{col_column}.}
+
+\item{row_column}{Column name of the row position variable (default \code{"row"}).}
+
+\item{col_column}{Column name of the column position variable
+(default \code{"col"}).}
+
+\item{by}{Optional column name grouping plots into grids (e.g. \code{"site"}).
+\code{NULL} treats the design as one grid.}
+}
+\value{
+A named list with one element per grid, each a list of \code{rows} (the
+positions in \code{df} belonging to that grid) and \code{index} (that grid's
+\code{\link[=grid_index]{grid_index()}}). Named \code{"1"} when \code{by} is \code{NULL}.
+}
+\description{
+The multi-grid counterpart of \code{\link[=grid_index]{grid_index()}}. A multi-environment trial is
+several grids that share a treatment set and never share an edge, so it
+cannot be one matrix: sites reuse \code{row}/\code{col}, and pooling them either
+silently overwrites plots or invents adjacencies between sites.
+
+With \code{by = NULL} this returns a one-element list, so callers have a single
+code path whether or not the design spans grids.
+}
+\keyword{internal}
diff --git a/man/objective_function_piepho.Rd b/man/objective_function_piepho.Rd
index 699f5090..288ac3f6 100644
--- a/man/objective_function_piepho.Rd
+++ b/man/objective_function_piepho.Rd
@@ -13,6 +13,8 @@ objective_function_piepho(
   pair_mapping = NULL,
   row_column = "row",
   col_column = "col",
+  by = NULL,
+  grid_index = NULL,
   ...
 )
 }
@@ -33,6 +35,12 @@ objective_function_piepho(
 
 \item{col_column}{Name of column representing the column of the design (default: "col")}
 
+\item{by, grid_index}{Optional grouping of plots into separate grids; see
+\code{\link[=calculate_adjacency_score]{calculate_adjacency_score()}}. Neighbour balance sums across grids, while
+evenness of distribution is scored \strong{per grid} and the scores summed,
+reported both in total and per grid. A grid with no treatment replicated
+inside it contributes \code{0}.}
+
 \item{...}{Extra parameters passed from \link{speed}}
 }
 \value{
diff --git a/man/speed.Rd b/man/speed.Rd
index daf356b9..d266441f 100644
--- a/man/speed.Rd
+++ b/man/speed.Rd
@@ -43,7 +43,15 @@ consider for balance (default: \code{~row + col}).}
 
 \item{grid_factors}{A named list specifying grid factors to construct a
 matrix for calculating adjacency score, \code{dim1} for row and \code{dim2} for
-column. (default: \code{list(dim1 = "row", dim2 = "col")}).}
+column. (default: \code{list(dim1 = "row", dim2 = "col")}).
+
+An optional third element, \code{by}, names a column that groups plots into
+\emph{separate} grids - a multi-environment trial, where each site reuses the
+same \code{row}/\code{col} numbering. Each grid is then scored on its own and the
+adjacency counts summed, so no adjacency is counted between plots at
+different sites, e.g.
+\code{list(dim1 = "row", dim2 = "col", by = "site")}. Without it, a design whose
+sites share coordinates is refused rather than silently pooled.}
 
 \item{iterations}{Maximum number of iterations for the simulated annealing
 algorithm (default: 10000). For hierarchical designs, can be a named list
diff --git a/man/summary.design.Rd b/man/summary.design.Rd
index 49f1f30e..42412e29 100644
--- a/man/summary.design.Rd
+++ b/man/summary.design.Rd
@@ -41,7 +41,12 @@ See Details for more information.}
 A list of class \code{"summary.design"}:
 \itemize{
 \item \strong{hierarchical} - \code{TRUE} for a multi-level (e.g. split-plot) design.
-\item \strong{layout} - \code{n_plots}, \code{nrow}, \code{ncol}, \code{row_column}, \code{col_column}, \code{has_grid}.
+\item \strong{layout} - \code{n_plots}, \code{nrow}, \code{ncol}, \code{row_column}, \code{col_column},
+\code{has_grid} (\code{TRUE} when the design is reportable as a single grid), and
+\code{grid_reason} (why not, or \code{NA}). \code{nrow}/\code{ncol} count the rows and columns
+the design \emph{occupies}, so a design with a gap in its coordinates (a missing
+plot, or a removed buffer) reports fewer than the coordinates span. Both are
+\code{NA} unless \code{has_grid}.
 \item \strong{levels} - character vector of level names (e.g. \code{"wp"}/\code{"sp"}; a single
 name for a simple design).
 \item \strong{per_level} - one element per level (named by \code{levels}), each a list with:
diff --git a/man/verify.Rd b/man/verify.Rd
index ac0175cf..5a0e9d61 100644
--- a/man/verify.Rd
+++ b/man/verify.Rd
@@ -5,6 +5,7 @@
 \alias{.verify_hierarchical_inputs}
 \alias{.verify_optim_params}
 \alias{.verify_swap_all_replication}
+\alias{.verify_grid_by}
 \title{Verify Inputs for \code{speed}}
 \usage{
 .verify_speed_inputs(
@@ -42,6 +43,8 @@
 )
 
 .verify_swap_all_replication(data, optimise, dummy_group = NULL)
+
+.verify_grid_by(data, grid_factors)
 }
 \arguments{
 \item{data}{A data frame containing the experimental design with spatial
@@ -85,6 +88,18 @@ see more in example.}
 
 \item{dummy_group}{Name of the internal placeholder column used for a level
 with no \code{swap_within} boundary, so it can be described as the whole design.}
+
+\item{grid_factors}{A named list specifying grid factors to construct a
+matrix for calculating adjacency score, \code{dim1} for row and \code{dim2} for
+column. (default: \code{list(dim1 = "row", dim2 = "col")}).
+
+An optional third element, \code{by}, names a column that groups plots into
+\emph{separate} grids - a multi-environment trial, where each site reuses the
+same \code{row}/\code{col} numbering. Each grid is then scored on its own and the
+adjacency counts summed, so no adjacency is counted between plots at
+different sites, e.g.
+\code{list(dim1 = "row", dim2 = "col", by = "site")}. Without it, a design whose
+sites share coordinates is refused rather than silently pooled.}
 }
 \description{
 Verify inputs for the \code{speed} function.
@@ -99,5 +114,8 @@ any optimisation happens rather than silently altering replication.
 Called on the resolved \code{optimise} list, so it covers simple, legacy
 hierarchical and \verb{optimise = } calls alike, including levels that set
 \code{swap_all} individually.
+
+\code{grid_factors} is a plain list, so a mistyped \code{by} would be ignored and every
+grid silently pooled. Checked before any optimisation happens.
 }
 \keyword{internal}
diff --git a/tests/testthat/test-build_design_matrix.R b/tests/testthat/test-build_design_matrix.R
new file mode 100644
index 00000000..b73fbced
--- /dev/null
+++ b/tests/testthat/test-build_design_matrix.R
@@ -0,0 +1,434 @@
+# build_design_matrix() places treatments by their (row, col) coordinates rather
+# than assuming a data frame ordering. These tests pin that it is correct for
+# *both* orderings the package produces: speed() sorts row-major, while
+# initialise_design_df() emits column-major.
+
+test_that("row-major and column-major frames give the same grid", {
+  # Same physical 2x3 layout, described two different ways.
+  #   col 1  col 2  col 3
+  # row 1: a      c      b
+  # row 2: b      a      c
+  row_major <- data.frame(
+    row = c(1, 1, 1, 2, 2, 2),
+    col = c(1, 2, 3, 1, 2, 3),
+    trt = c("a", "c", "b", "b", "a", "c")
+  )
+  col_major <- data.frame(
+    row = c(1, 2, 1, 2, 1, 2),
+    col = c(1, 1, 2, 2, 3, 3),
+    trt = c("a", "b", "c", "a", "b", "c")
+  )
+
+  expected <- matrix(
+    c("a", "b", "c", "a", "b", "c"),
+    nrow = 2,
+    ncol = 3
+  )
+  expect_equal(build_design_matrix(row_major, "trt"), expected)
+  expect_equal(build_design_matrix(col_major, "trt"), expected)
+})
+
+test_that("shuffled row order does not change the grid", {
+  df <- data.frame(
+    row = rep(1:3, each = 4),
+    col = rep(1:4, times = 3),
+    trt = LETTERS[1:12]
+  )
+  expect_equal(
+    build_design_matrix(df[sample.int(12), ], "trt"),
+    build_design_matrix(df, "trt")
+  )
+})
+
+test_that("non-square grids keep their orientation", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2, 3, 3),
+    col = c(1, 2, 1, 2, 1, 2),
+    trt = c("a", "b", "c", "d", "e", "f")
+  )
+  m <- build_design_matrix(df, "trt")
+
+  expect_equal(dim(m), c(3L, 2L))
+  expect_equal(m[3, 1], "e")
+  expect_equal(m[1, 2], "b")
+})
+
+test_that("factor coordinates with lexical level order are handled", {
+  # as.factor() on characters orders levels 1, 10, 11, 2, ... A fill that
+  # relied on order() following the levels would build a permuted grid.
+  df <- data.frame(
+    row = factor(as.character(1:11)),
+    col = factor(rep("1", 11)),
+    trt = LETTERS[1:11]
+  )
+  m <- build_design_matrix(df, "trt")
+
+  expect_equal(dim(m), c(11L, 1L))
+  expect_equal(as.vector(m), LETTERS[1:11])
+})
+
+test_that("gaps in the coordinates become NA cells, not a collapsed grid", {
+  # Rows 1, 2, 4, 5 exist; row 3 does not (e.g. a road). Plots in rows 2 and 4
+  # must not become neighbours.
+  df <- data.frame(
+    row = rep(c(1, 2, 4, 5), each = 2),
+    col = rep(1:2, times = 4),
+    trt = c("A", "B", "C", "C", "C", "C", "A", "B")
+  )
+  m <- build_design_matrix(df, "trt")
+
+  expect_equal(dim(m), c(5L, 2L))
+  expect_true(all(is.na(m[3, ])))
+  expect_equal(m[2, ], c("C", "C"))
+  expect_equal(m[4, ], c("C", "C"))
+})
+
+test_that("coordinates that do not start at 1 are preserved, not shifted", {
+  # add_buffers() offsets the real design's coordinates; dropping the buffer
+  # rows leaves them starting above 1.
+  df <- data.frame(
+    row = c(2, 2, 3, 3),
+    col = c(2, 3, 2, 3),
+    trt = c("A", "B", "B", "A")
+  )
+  m <- build_design_matrix(df, "trt")
+
+  expect_equal(dim(m), c(3L, 3L))
+  expect_true(all(is.na(m[1, ])))
+  expect_true(all(is.na(m[, 1])))
+  expect_equal(m[2, 2], "A")
+})
+
+test_that("NA treatments are placed as NA", {
+  df <- data.frame(
+    row = c(1, 1, 2, 2),
+    col = c(1, 2, 1, 2),
+    trt = c("A", NA, "B", "A")
+  )
+  expect_true(is.na(build_design_matrix(df, "trt")[1, 2]))
+})
+
+test_that("non-numeric coordinates give an informative error", {
+  df <- data.frame(
+    row = c("R1", "R1"),
+    col = c("C1", "C2"),
+    trt = c("A", "B")
+  )
+  expect_error(
+    build_design_matrix(df, "trt"),
+    "must be numeric, or coercible to numeric"
+  )
+})
+
+test_that("non-positive or fractional coordinates give an informative error", {
+  zero_based <- data.frame(
+    row = c(0, 0, 1, 1),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "B", "A")
+  )
+  expect_error(
+    build_design_matrix(zero_based, "trt"),
+    "positive whole numbers"
+  )
+
+  fractional <- data.frame(
+    row = c(1, 1, 1.5, 1.5),
+    col = c(1, 2, 1, 2),
+    trt = c("A", "B", "B", "A")
+  )
+  expect_error(
+    build_design_matrix(fractional, "trt"),
+    "positive whole numbers"
+  )
+})
+
+test_that("duplicate coordinates error rather than silently overwriting", {
+  # Multi-site designs reuse row/col per site. Placing them on one grid would
+  # keep only the last site written.
+  met <- data.frame(
+    site = rep(c("s1", "s2"), each = 4),
+    row = rep(c(1, 1, 2, 2), times = 2),
+    col = rep(c(1, 2, 1, 2), times = 2),
+    trt = c("A", "B", "B", "A", "A", "A", "B", "B")
+  )
+  expect_error(
+    build_design_matrix(met, "trt"),
+    "Duplicate \\(row, col\\) coordinates"
+  )
+})
+
+test_that("column names are reported in errors", {
+  df <- data.frame(
+    range = c("R1", "R1"),
+    bed = c("C1", "C2"),
+    variety = c("A", "B")
+  )
+  expect_error(
+    build_design_matrix(
+      df,
+      "variety",
+      row_column = "range",
+      col_column = "bed"
+    ),
+    "`range` and `bed`"
+  )
+})
+
+test_that("grid_index() returns the index and dimensions build_design_matrix needs", {
+  d <- initialise_design_df(rep(LETTERS[1:6], 4), 3, 8)
+  gi <- grid_index(d)
+
+  expect_named(gi, c("idx", "nrow", "ncol", "n"))
+  expect_equal(gi$nrow, 3)
+  expect_equal(gi$ncol, 8)
+  expect_equal(gi$n, nrow(d))
+  expect_equal(dim(gi$idx), c(nrow(d), 2L))
+})
+
+test_that("a supplied index gives the same grid as validating in place", {
+  # The index is the only thing hoisted out of the annealing loop, so the two
+  # paths must be indistinguishable.
+  for (dims in list(c(3, 8), c(8, 3), c(4, 4), c(2, 6))) {
+    d <- initialise_design_df(
+      rep(LETTERS[1:4], prod(dims) / 4),
+      dims[[1]],
+      dims[[2]]
+    )
+    expect_equal(
+      build_design_matrix(d, "treatment", index = grid_index(d)),
+      build_design_matrix(d, "treatment"),
+      info = paste(dims, collapse = "x")
+    )
+  }
+})
+
+test_that("a supplied index survives the treatment column being reshuffled", {
+  # What the annealing loop actually does: same coordinates, different treatments.
+  d <- initialise_design_df(rep(LETTERS[1:6], 4), 3, 8)
+  gi <- grid_index(d)
+  set.seed(1)
+  d$treatment <- sample(d$treatment)
+
+  expect_equal(
+    build_design_matrix(d, "treatment", index = gi),
+    build_design_matrix(d, "treatment")
+  )
+})
+
+test_that("a supplied index is checked against the design it is used with", {
+  # A stale index would place treatments at the wrong coordinates silently.
+  d <- initialise_design_df(rep(LETTERS[1:6], 4), 3, 8)
+  gi <- grid_index(d)
+
+  expect_error(
+    build_design_matrix(d[1:10, ], "treatment", index = gi),
+    "was built for 24 plots but `df` has 10"
+  )
+})
+
+test_that("calculate_adjacency_score() accepts a pre-built index", {
+  d <- initialise_design_df(rep(LETTERS[1:6], 4), 3, 8)
+  expect_equal(
+    calculate_adjacency_score(d, "treatment", grid_index = grid_indices(d)),
+    calculate_adjacency_score(d, "treatment")
+  )
+})
+
+test_that("calculate_adjacency_score() sums per grid and never across them", {
+  # Two sites sharing row/col numbering. Scoring them as one grid is impossible
+  # (duplicate coordinates); scoring them side by side in one grid would invent
+  # adjacencies between plots at different sites.
+  set.seed(1)
+  d <- rbind(
+    data.frame(
+      site = "a",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(rep(LETTERS[1:6], 2))
+    ),
+    data.frame(
+      site = "b",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(rep(LETTERS[1:6], 2))
+    )
+  )
+
+  expect_error(
+    calculate_adjacency_score(d, "treatment"),
+    class = "speed_grid_duplicate"
+  )
+
+  per_site <- sum(vapply(
+    split(d, d$site),
+    function(s) return(calculate_adjacency_score(s, "treatment")),
+    numeric(1)
+  ))
+  expect_equal(
+    calculate_adjacency_score(d, "treatment", by = "site"),
+    per_site
+  )
+
+  # Laid side by side the coordinates are unique, so nothing errors - but the
+  # score picks up adjacencies across the join, which is exactly what `by`
+  # prevents. This is the case no amount of coordinate validation can catch.
+  side_by_side <- d
+  side_by_side$col[side_by_side$site == "b"] <-
+    side_by_side$col[side_by_side$site == "b"] + 3
+  expect_gt(
+    calculate_adjacency_score(side_by_side, "treatment"),
+    calculate_adjacency_score(side_by_side, "treatment", by = "site")
+  )
+  expect_equal(
+    calculate_adjacency_score(side_by_side, "treatment", by = "site"),
+    per_site
+  )
+})
+
+test_that("grid_indices() returns one index per grid", {
+  d <- rbind(
+    data.frame(site = "a", expand.grid(row = 1:4, col = 1:3), treatment = "x"),
+    data.frame(site = "b", expand.grid(row = 1:4, col = 1:3), treatment = "y")
+  )
+
+  single <- grid_indices(d[d$site == "a", ])
+  expect_length(single, 1)
+  expect_equal(single[[1]]$index$n, 12)
+
+  split_grids <- grid_indices(d, by = "site")
+  expect_named(split_grids, c("a", "b"))
+  expect_equal(attr(split_grids, "by"), "site")
+  expect_equal(vapply(split_grids, function(g) return(g$index$n), numeric(1)),
+               c(a = 12, b = 12))
+  # The rows recorded for each grid are the rows of that site, so a caller can
+  # subset the design with them.
+  expect_equal(split_grids$b$rows, which(d$site == "b"))
+
+  expect_error(grid_indices(d, by = "nope"), class = "speed_grid_missing_by")
+})
+
+test_that("grid_index() reports why a design cannot be gridded", {
+  ok <- expand.grid(row = 1:2, col = 1:2)
+
+  nonnumeric <- ok
+  nonnumeric$row <- c("a", "b", "c", "d")
+  expect_error(grid_index(nonnumeric), class = "speed_grid_nonnumeric")
+
+  # Coordinates index a matrix directly, so both zero and fractional fail
+  below_one <- ok
+  below_one$row <- c(0L, 1L, 2L, 3L)
+  expect_error(grid_index(below_one), class = "speed_grid_notinteger")
+
+  fractional <- ok
+  fractional$col <- c(1.5, 2.5, 3.5, 4.5)
+  expect_error(grid_index(fractional), class = "speed_grid_notinteger")
+})
+
+test_that("every grid failure carries a `speed_grid_error` and a reason", {
+  # speed_hierarchical() and summary() both dispatch on the parent class rather
+  # than the specific one, so each condition must carry it and a `reason`.
+  ok <- expand.grid(row = 1:2, col = 1:2)
+
+  nonnumeric <- ok
+  nonnumeric$row <- c("a", "b", "c", "d")
+  fractional <- ok
+  fractional$col <- c(1.5, 2.5, 3.5, 4.5)
+  duplicated_coords <- rbind(ok, ok)
+
+  cases <- list(
+    missing = data.frame(lane = 1:4),
+    nonnumeric = nonnumeric,
+    notinteger = fractional,
+    duplicate = duplicated_coords
+  )
+
+  for (nm in names(cases)) {
+    err <- tryCatch(grid_index(cases[[nm]]), speed_grid_error = function(e) {
+      return(e)
+    })
+    expect_s3_class(err, "speed_grid_error")
+    expect_type(err$reason, "character")
+    expect_gt(nchar(err$reason), 0)
+  }
+
+  expect_error(
+    grid_indices(ok, by = "nope"),
+    class = "speed_grid_error"
+  )
+})
+
+test_that("speed() scores identically whether or not the index is hoisted", {
+  # The hoist is a performance change only; guard against it becoming a
+  # behaviour change. objective_function_piepho() builds a grid every iteration,
+  # so it is the strongest case.
+  d <- initialise_design_df(rep(LETTERS[1:6], 4), 3, 8)
+  args <- list(
+    swap = "treatment",
+    swap_within = "1",
+    spatial_factors = ~ row + col,
+    iterations = 200,
+    seed = 42,
+    quiet = TRUE
+  )
+  hoisted <- do.call(speed, c(list(d), args))
+  # Calling the objective directly takes the un-hoisted path (index = NULL).
+  direct <- objective_function(
+    hoisted$design_df,
+    "treatment",
+    c("row", "col")
+  )
+  expect_equal(direct$score, hoisted$score)
+
+  hoisted_p <- do.call(
+    speed,
+    c(list(d), args, list(obj_function = objective_function_piepho))
+  )
+  direct_p <- objective_function_piepho(
+    hoisted_p$design_df,
+    "treatment",
+    c("row", "col")
+  )
+  expect_equal(direct_p$score, hoisted_p$score)
+})
+
+test_that("a design whose coordinates cannot form a grid still runs when no grid is needed", {
+  # The index is built lazily (tryCatch -> NULL) so that a design which never
+  # needs a grid is unaffected. Two sites reusing row/col have duplicate
+  # coordinates, which grid_index() rejects.
+  d <- data.frame(
+    site = rep(c("A", "B"), each = 12),
+    row = rep(rep(1:4, times = 3), 2),
+    col = rep(rep(1:3, each = 4), 2),
+    treatment = rep(rep(LETTERS[1:3], 4), 2)
+  )
+  expect_error(grid_index(d), "Duplicate")
+  expect_no_error(
+    r <- speed(
+      d,
+      swap = "treatment",
+      swap_within = "site",
+      spatial_factors = ~ row + col + site,
+      iterations = 100,
+      seed = 1,
+      quiet = TRUE,
+      optimise_params = optim_params(adj_weight = 0)
+    )
+  )
+  expect_equal(r$score, 4)
+})
+
+test_that("grid_index() names a missing coordinate column", {
+  # A design with no grid at all reaches here from speed(); without this check
+  # the absent columns reach max() as empty vectors and give -Inf dimensions.
+  d <- data.frame(a = 1:4, b = 1:4, treatment = LETTERS[1:4])
+
+  expect_error(
+    grid_index(d),
+    "no `row` or `col` column",
+    class = "speed_grid_missing"
+  )
+  # Only one of the pair missing is still a missing-column failure
+  expect_error(
+    grid_index(d, "row", "b"),
+    "no `row` column",
+    class = "speed_grid_missing"
+  )
+})
diff --git a/tests/testthat/test-calculate_adjacency_score.R b/tests/testthat/test-calculate_adjacency_score.R
index 0c390513..38906224 100644
--- a/tests/testthat/test-calculate_adjacency_score.R
+++ b/tests/testthat/test-calculate_adjacency_score.R
@@ -130,7 +130,14 @@ test_that("adjacency_score_vec applies per-ring weights", {
 })
 
 test_that("adjacency_score_vec rejects mismatched dists/weights", {
+  # A length-one `weights` is recycled across every ring, so the scalar default
+  # stays usable with a multi-ring `dists`. Any other mismatch is still an error.
   expect_error(adjacency_score_vec(
+    matrix(1, 2, 2),
+    dists = c(1, 2, 3),
+    weights = c(1, 2)
+  ))
+  expect_no_error(adjacency_score_vec(
     matrix(1, 2, 2),
     dists = c(1, 2),
     weights = 1
diff --git a/tests/testthat/test-calculate_efficiency_factor.R b/tests/testthat/test-calculate_efficiency_factor.R
index 3e246bf5..d3206874 100644
--- a/tests/testthat/test-calculate_efficiency_factor.R
+++ b/tests/testthat/test-calculate_efficiency_factor.R
@@ -1,6 +1,19 @@
+# The published designs below are written out visually, one grid row per source
+# line, so they can be read against the paper. `initialise_design_df()` assigns
+# `items` *down columns* - `expand.grid(row, col)` varies `row` fastest - so a
+# row-major literal has to be transposed before it is handed over, or the design
+# stored at those coordinates is not the one written here.
+by_row <- function(items, nrows, ncols) {
+  initialise_design_df(
+    as.vector(matrix(items, nrow = nrows, ncol = ncols, byrow = TRUE)),
+    nrows,
+    ncols
+  )
+}
+
 test_that("calculate_efficiency_factor provides the same results as the paper", {
   # fmt: skip
-  df_design1 <- initialise_design_df(c(
+  df_design1 <- by_row(c(
     7, 5, 6, 9, 4, 1, 3, 2, 8,
     5, 6, 3, 1, 7, 8, 2, 4, 9,
     8, 9, 5, 6, 3, 4, 1, 7, 2,
@@ -13,7 +26,7 @@ test_that("calculate_efficiency_factor provides the same results as the paper",
   )
 
   # fmt: skip
-  df_design2 <- initialise_design_df(c(
+  df_design2 <- by_row(c(
     8, 5, 7, 2, 4, 1, 6, 9, 3,
     1, 9, 8, 6, 3, 2, 4, 7, 5,
     7, 6, 4, 9, 5, 8, 3, 2, 1,
@@ -26,7 +39,7 @@ test_that("calculate_efficiency_factor provides the same results as the paper",
   )
 
   # fmt: skip
-  df_design3 <- initialise_design_df(c(
+  df_design3 <- by_row(c(
     9, 8, 1, 4, 3, 7, 5, 2, 6,
     7, 5, 6, 2, 9, 1, 3, 8, 4,
     2, 4, 3, 5, 6, 8, 9, 7, 1,
@@ -39,7 +52,7 @@ test_that("calculate_efficiency_factor provides the same results as the paper",
   )
 
   # fmt: skip
-  df_design4 <- initialise_design_df(c(
+  df_design4 <- by_row(c(
     47, 16, 43, 42, 37, 35,  1, 59, 24, 19,  4, 18, 40, 28, 51, 29, 54, 57,
     12, 25,  6, 57, 47, 32, 39, 17, 31, 50, 15,  5, 55, 51,  9, 54, 41,  3,
     23, 18, 45, 36, 49,  7,  8, 60, 41, 29,  3, 58, 26, 52,  2, 15, 28, 27,
@@ -54,15 +67,20 @@ test_that("calculate_efficiency_factor provides the same results as the paper",
 })
 
 test_that("calculate_efficiency_factor provides better result for an optimised design", {
+  # A poor but *estimable* comparator. The obvious "unoptimised" layout - each
+  # treatment filling one grid row - is not merely inefficient, it confounds
+  # treatment with row, so it has no efficiency factor at all and is refused
+  # (see the estimability tests below). Comparing against it would have been
+  # comparing against a number that does not exist.
   # fmt: skip
-  df_design_initial <- initialise_design_df(c(
-    1, 1, 2, 2,
-    3, 3, 4, 4,
-    5, 5, 6, 6
+  df_design_initial <- by_row(c(
+    1, 2, 6, 3,
+    4, 3, 5, 5,
+    1, 4, 2, 6
   ), 3, 4)
 
   # fmt: skip
-  df_design_optimised <- initialise_design_df(c(
+  df_design_optimised <- by_row(c(
     1, 2, 4, 3,
     5, 1, 6, 2,
     3, 6, 5, 4
@@ -72,25 +90,28 @@ test_that("calculate_efficiency_factor provides better result for an optimised d
     abs(1 - calculate_efficiency_factor(df_design_optimised, "treatment")),
     abs(1 - calculate_efficiency_factor(df_design_initial, "treatment"))
   )
+  # Both are genuine values in [0, 1], so the comparison above is meaningful.
+  expect_lte(calculate_efficiency_factor(df_design_optimised, "treatment"), 1)
+  expect_gt(calculate_efficiency_factor(df_design_initial, "treatment"), 0)
 })
 
 test_that("calculate_efficiency_factor provides same result for mathematically identical designs", {
   # fmt: skip
-  df_design1 <- initialise_design_df(c(
+  df_design1 <- by_row(c(
     1, 2, 4, 3,
     5, 1, 6, 2,
     3, 6, 5, 4
   ), 3, 4)
 
   # fmt: skip
-  df_design2 <- initialise_design_df(c(
+  df_design2 <- by_row(c(
     "a", "b", "d", "c",
     "e", "a", "f", "b",
     "c", "f", "e", "d"
   ), 3, 4)
 
   # fmt: skip
-  df_design3 <- initialise_design_df(c(
+  df_design3 <- by_row(c(
     "b", "a", "c", "d",
     "e", "b", "f", "a",
     "d", "f", "e", "c"
@@ -107,28 +128,42 @@ test_that("calculate_efficiency_factor provides same result for mathematically i
   )
 })
 
-test_that("calculate_efficiency_factor handles near-singular matrices using pseudoinverse", {
-  # Instead of trying to create a design that naturally triggers the pseudoinverse,
-  # let's test with a design that we know works and verify the function handles
-  # both the regular and pseudoinverse cases properly
-
-  # Use a simple but slightly unbalanced design
-  df_design_test <- data.frame(
+test_that("calculate_efficiency_factor refuses a design whose contrasts are not estimable", {
+  # C occupies row 3 entirely, so the C-vs-A contrast cannot be separated from
+  # the row 3 effect. Confirmed independently: fitting
+  # `y ~ factor(row) + factor(col) + treatment` with lm() aliases one treatment
+  # coefficient. There is no efficiency factor for such a design - before this
+  # check the formula returned a finite, plausible-looking value anyway.
+  df_confounded_row <- data.frame(
     row = c(1, 1, 2, 2, 3, 3),
     col = c(1, 2, 1, 2, 1, 2),
     treatment = c("A", "B", "A", "B", "C", "C")
   )
 
-  # This should complete without error regardless of which inversion method is used
-  expect_no_error({
-    result <- calculate_efficiency_factor(df_design_test, "treatment")
-  })
+  expect_error(
+    calculate_efficiency_factor(df_confounded_row, "treatment"),
+    class = "speed_efficiency_rank"
+  )
+  expect_error(
+    calculate_efficiency_factor(df_confounded_row, "treatment"),
+    "not all treatment contrasts are estimable",
+    ignore.case = TRUE
+  )
+})
 
-  result <- calculate_efficiency_factor(df_design_test, "treatment")
-  expect_type(result, "double")
-  expect_length(result, 1)
+test_that("calculate_efficiency_factor keeps designs with no residual df to spare", {
+  # The boundary the gate must not overshoot: zero residual degrees of freedom
+  # is still estimable, so these must return a value rather than be refused.
+  # fmt: skip
+  df_df0 <- by_row(c(
+    1, 2, 3,
+    3, 1, 2
+  ), 2, 3)
+
+  result <- calculate_efficiency_factor(df_df0, "treatment")
   expect_true(is.finite(result))
   expect_gt(result, 0)
+  expect_lte(result, 1)
 })
 
 test_that("calculate_efficiency_factor works with minimal design dimensions", {
@@ -149,25 +184,41 @@ test_that("calculate_efficiency_factor works with minimal design dimensions", {
   expect_gt(result, 0)
 })
 
-test_that("calculate_efficiency_factor handles designs with high treatment-space confounding", {
-  # Create a design where treatments are highly confounded with a spatial dimension
-  # This creates a scenario more likely to need pseudoinverse without being perfectly singular
+test_that("calculate_efficiency_factor refuses a single-column design with blocked treatments", {
+  # One plot per row means every row effect is a plot effect, so nothing is left
+  # to estimate a treatment difference with. Refusing is the only honest answer;
+  # the formula previously returned a finite value here.
   df_design_confounded <- data.frame(
     row = rep(1:8, each = 1),
     col = rep(1, times = 8),
     treatment = c("A", "A", "A", "A", "B", "B", "B", "B") # Single column, treatments in blocks
   )
 
-  # This single-column design with blocked treatments should be numerically challenging
-  # but still solvable
-  expect_no_error({
-    result <- calculate_efficiency_factor(df_design_confounded, "treatment")
-  })
+  expect_error(
+    calculate_efficiency_factor(df_design_confounded, "treatment"),
+    class = "speed_efficiency_rank"
+  )
+})
 
-  result <- calculate_efficiency_factor(df_design_confounded, "treatment")
-  expect_type(result, "double")
-  expect_true(is.finite(result))
-  expect_gt(result, 0)
+test_that("calculate_efficiency_factor refuses unreplicated and degenerate designs", {
+  # Every route to an impossible value, pinned together. Each was measured
+  # returning > 1 before the rank check existed.
+  unreplicated <- initialise_design_df(as.character(1:12), 4, 3)
+  single_row <- data.frame(
+    row = rep(1, 6),
+    col = 1:6,
+    treatment = rep(c("A", "B", "C"), 2)
+  )
+  # Column-major storage, so this is one treatment per grid row - `each = 4`
+  # would give a diagonal pattern, which is estimable.
+  aliased_with_row <- initialise_design_df(rep(c("a", "b", "c"), 4), 3, 4)
+
+  for (d in list(unreplicated, single_row, aliased_with_row)) {
+    expect_error(
+      calculate_efficiency_factor(d, "treatment"),
+      class = "speed_efficiency_rank"
+    )
+  }
 })
 
 test_that("calculate_efficiency_factor uses pseudoinverse for matrices with high condition numbers", {
@@ -205,24 +256,12 @@ test_that("calculate_efficiency_factor uses pseudoinverse for matrices with high
 })
 
 test_that("calculate_efficiency_factor honours custom row/col column names", {
-  df <- initialise_design_df(
-    c(
-      "a",
-      "b",
-      "d",
-      "c",
-      "e",
-      "a",
-      "f",
-      "b",
-      "c",
-      "f",
-      "e",
-      "d"
-    ),
-    3,
-    4
-  )
+  # fmt: skip
+  df <- by_row(c(
+    "a", "b", "d", "c",
+    "e", "a", "f", "b",
+    "c", "f", "e", "d"
+  ), 3, 4)
   base <- calculate_efficiency_factor(df, "treatment")
 
   # Same design, grid columns renamed - must give the same efficiency.
diff --git a/tests/testthat/test-grid-orientation.R b/tests/testthat/test-grid-orientation.R
new file mode 100644
index 00000000..4f565ab5
--- /dev/null
+++ b/tests/testthat/test-grid-orientation.R
@@ -0,0 +1,347 @@
+# Grid metrics must describe the layout the coordinates define, whatever order
+# the data frame is in. Every expected value here is derived from the (row, col)
+# coordinates by hand, never by reshaping the treatment column.
+
+# --- calculate_adjacency_score() composability -------------------------------
+
+test_that("adjacency score is the same for either input ordering", {
+  # Same physical 2x3 layout described row-major and column-major.
+  row_major <- data.frame(
+    row = c(1, 1, 1, 2, 2, 2),
+    col = c(1, 2, 3, 1, 2, 3),
+    trt = c("a", "c", "b", "b", "a", "c")
+  )
+  col_major <- data.frame(
+    row = c(1, 2, 1, 2, 1, 2),
+    col = c(1, 1, 2, 2, 3, 3),
+    trt = c("a", "b", "c", "a", "b", "c")
+  )
+
+  expect_equal(
+    calculate_adjacency_score(row_major, "trt"),
+    calculate_adjacency_score(col_major, "trt")
+  )
+})
+
+test_that("calculate_adjacency_score() composes with initialise_design_df()", {
+  # initialise_design_df() emits column-major data, and no treatment neighbours
+  # itself in this layout, so the score is 0.
+  #   col 1  col 2  col 3
+  # row 1: a      c      b
+  # row 2: b      a      c
+  df <- initialise_design_df(
+    items = rep(c("a", "b", "c"), 2),
+    nrows = 2,
+    ncols = 3
+  )
+  expect_equal(calculate_adjacency_score(df, "treatment"), 0)
+})
+
+test_that("adjacency score counts like-treatment edges on a known layout", {
+  # row 1: A A A
+  # row 2: B B B
+  # row 3: A A A
+  # Horizontal like-pairs: 2 per row x 3 rows = 6. No vertical like-pairs.
+  df <- data.frame(
+    row = rep(1:3, each = 3),
+    col = rep(1:3, times = 3),
+    treatment = c("A", "A", "A", "B", "B", "B", "A", "A", "A")
+  )
+  expect_equal(calculate_adjacency_score(df, "treatment"), 6)
+})
+
+# --- objective_function_piepho() ---------------------------------------------
+
+test_that("piepho score does not depend on the input row ordering", {
+  # This is the property the whole change buys: the score describes the layout,
+  # not the order the rows happen to be in.
+  col_major <- initialise_design_df(
+    items = rep(c("a", "b", "c"), 4),
+    nrows = 2,
+    ncols = 6
+  )
+  row_major <- col_major[order(col_major$row, col_major$col), ]
+  rownames(row_major) <- NULL
+  pm <- create_pair_mapping(col_major$treatment)
+
+  from_col_major <- objective_function_piepho(
+    col_major,
+    "treatment",
+    c("row", "col"),
+    pair_mapping = pm
+  )
+  from_row_major <- objective_function_piepho(
+    row_major,
+    "treatment",
+    c("row", "col"),
+    pair_mapping = pm
+  )
+
+  expect_equal(from_row_major$score, from_col_major$score)
+  expect_equal(from_row_major$components, from_col_major$components)
+})
+
+test_that("piepho components match hand-derived values on a non-square grid", {
+  # 2x6 layout from initialise_design_df(rep(c("a","b","c"), 4), 2, 6):
+  #   col 1  2  3  4  5  6
+  # row 1:  a  c  b  a  c  b
+  # row 2:  b  a  c  b  a  c
+  # Every treatment appears 4 times, twice per row and never twice in a column,
+  # so the balance score is the row term only: 3 treatments x var(c(2,2)) = 0
+  # across columns, and rowVars over the 2 rows gives 2 in total.
+  df <- initialise_design_df(
+    items = rep(c("a", "b", "c"), 4),
+    nrows = 2,
+    ncols = 6
+  )
+  df <- df[order(df$row, df$col), ]
+  pm <- create_pair_mapping(df$treatment)
+
+  res <- objective_function_piepho(
+    df,
+    "treatment",
+    c("row", "col"),
+    pair_mapping = pm
+  )
+
+  # No like-treatment neighbours anywhere in this layout.
+  expect_equal(res$components[["adjacency"]], 0)
+  expect_equal(
+    res$components[["adjacency"]],
+    calculate_adjacency_score(df, "treatment")
+  )
+  # Balance must be computed on the real treatment column.
+  expect_equal(
+    res$components[["balance"]],
+    calculate_balance_score(df, "treatment", c("row", "col"))
+  )
+})
+
+test_that("piepho does not overwrite the treatment column it was given", {
+  # Balance and adjacency must come out identical to computing them on `design`
+  # directly, i.e. the treatments stay aligned with their coordinates.
+  df <- initialise_design_df(
+    items = rep(LETTERS[1:4], 6),
+    nrows = 4,
+    ncols = 6
+  )
+  df <- df[order(df$row, df$col), ]
+  pm <- create_pair_mapping(df$treatment)
+
+  res <- objective_function_piepho(
+    df,
+    "treatment",
+    c("row", "col"),
+    pair_mapping = pm
+  )
+
+  expect_equal(
+    res$components[["balance"]],
+    calculate_balance_score(df, "treatment", c("row", "col"))
+  )
+  expect_equal(
+    res$components[["adjacency"]],
+    calculate_adjacency_score(df, "treatment")
+  )
+})
+
+# --- sparse grids ------------------------------------------------------------
+
+test_that("neighbour balance skips empty cells, with and without a mapping", {
+  # Rows 1, 2, 4, 5 exist; row 3 is a gap. Pairs, by coordinate:
+  # horizontal A-B, C-C, C-C, A-B; vertical A-C, B-C (rows 1-2), A-C, B-C
+  # (rows 4-5). Nothing crosses the gap.
+  df <- data.frame(
+    row = rep(c(1, 2, 4, 5), each = 2),
+    col = rep(1:2, times = 4),
+    trt = c("A", "B", "C", "C", "C", "C", "A", "B")
+  )
+  m <- build_design_matrix(df, "trt")
+
+  with_mapping <- calculate_nb(m, create_pair_mapping(df$trt))
+  without_mapping <- calculate_nb(m)
+
+  expect_equal(as.integer(with_mapping$nb[["C,C"]]), 2L)
+  expect_equal(as.integer(without_mapping$nb[["C,C"]]), 2L)
+  # No pair may involve an empty cell.
+  expect_false(any(grepl("NA", names(with_mapping$nb), fixed = TRUE)))
+  expect_false(any(grepl("NA", names(without_mapping$nb), fixed = TRUE)))
+})
+
+test_that("piepho runs on a design with missing plots", {
+  df <- data.frame(
+    row = rep(c(1, 2, 4, 5), each = 2),
+    col = rep(1:2, times = 4),
+    treatment = c("A", "B", "C", "D", "D", "C", "A", "B")
+  )
+  # Both the default (NULL) and supplied pair_mapping paths must work.
+  expect_no_error(
+    objective_function_piepho(df, "treatment", c("row", "col"))
+  )
+  expect_no_error(
+    objective_function_piepho(
+      df,
+      "treatment",
+      c("row", "col"),
+      pair_mapping = create_pair_mapping(df$treatment)
+    )
+  )
+})
+
+test_that("adjacency score ignores gaps rather than closing them", {
+  # C-C pairs exist within rows 2 and 4 but not between them: the gap at row 3
+  # must not be collapsed.
+  df <- data.frame(
+    row = rep(c(1, 2, 4, 5), each = 2),
+    col = rep(1:2, times = 4),
+    trt = c("A", "B", "C", "C", "C", "C", "A", "B")
+  )
+  expect_equal(calculate_adjacency_score(df, "trt"), 2)
+})
+
+# --- ring_weights recycling --------------------------------------------------
+
+test_that("a single ring_weight applies to every ring_dist", {
+  # ring_weights defaults to the scalar 1, so the documented default was
+  # unusable with a multi-ring ring_dists: adjacency_score_vec() asserted the
+  # two were the same length.
+  df <- data.frame(
+    row = rep(1:3, each = 3),
+    col = rep(1:3, times = 3),
+    trt = c("A", "B", "C", "B", "C", "A", "C", "A", "B")
+  )
+
+  expect_no_error(
+    scalar_weight <- calculate_adjacency_score(df, "trt", ring_dists = c(1, 2))
+  )
+  expect_equal(
+    scalar_weight,
+    calculate_adjacency_score(
+      df,
+      "trt",
+      ring_dists = c(1, 2),
+      ring_weights = c(1, 1)
+    )
+  )
+})
+
+test_that("a genuine ring_weights length mismatch is still an error", {
+  df <- data.frame(
+    row = rep(1:3, each = 3),
+    col = rep(1:3, times = 3),
+    trt = c("A", "B", "C", "B", "C", "A", "C", "A", "B")
+  )
+  expect_error(
+    calculate_adjacency_score(
+      df,
+      "trt",
+      ring_dists = c(1, 2, 3),
+      ring_weights = c(1, 2)
+    )
+  )
+})
+
+# --- calculate_efficiency_factor() -------------------------------------------
+
+test_that("efficiency factor does not depend on the input row ordering", {
+  # initialise_design_df() emits column-major data, so this is the ordering a
+  # caller is most likely to arrive with.
+  col_major <- initialise_design_df(
+    items = rep(LETTERS[1:3], 4),
+    nrows = 4,
+    ncols = 3
+  )
+  row_major <- col_major[order(col_major$row, col_major$col), ]
+  rownames(row_major) <- NULL
+
+  expect_equal(
+    calculate_efficiency_factor(col_major, treatment),
+    calculate_efficiency_factor(row_major, treatment)
+  )
+  # Pinned so the shared value cannot drift: an efficiency factor is bounded
+  # above by 1, and a value over it means a rank-deficient information matrix.
+  expect_equal(calculate_efficiency_factor(row_major, treatment), 0.9375)
+})
+
+test_that("efficiency factor is unchanged by shuffling a real design", {
+  d <- speed(
+    initialise_design_df(items = rep(LETTERS[1:6], 4), nrows = 3, ncols = 8),
+    swap = "treatment",
+    iterations = 500,
+    seed = 11,
+    quiet = TRUE
+  )
+  ordered <- calculate_efficiency_factor(d$design_df, treatment)
+  shuffled <- calculate_efficiency_factor(
+    d$design_df[rev(seq_len(nrow(d$design_df))), ],
+    treatment
+  )
+
+  expect_equal(ordered, shuffled)
+  expect_lte(ordered, 1)
+})
+
+test_that("efficiency factor is unaffected by an offset coordinate origin", {
+  # add_buffers() displaces the real plots to make room and never undoes it, so
+  # a de-buffered design arrives with coordinates that neither start at 1 nor
+  # run consecutively. The gaps leave empty indicator columns in Z, making ZtZ
+  # singular; the kappa() check routes to pseudo_inverse() and they contribute
+  # nothing, so the value is unchanged.
+  base <- initialise_design_df(
+    items = rep(LETTERS[1:3], 4),
+    nrows = 4,
+    ncols = 3
+  )
+  expected <- calculate_efficiency_factor(base, treatment)
+
+  # The displacement each add_buffers() type applies (R/buffers.R), reproduced
+  # here rather than by calling add_buffers(): the assertion is about the metric
+  # tolerating offset coordinates, which outlives buffers moving out of speed.
+  displaced <- list(
+    edge = transform(base, row = row + 1, col = col + 1),
+    row = transform(base, row = row * 2),
+    col = transform(base, col = col * 2),
+    `double row` = transform(base, row = (3 * row) - 1),
+    `double col` = transform(base, col = (3 * col) - 1),
+    # Stacked edge then row: rows 4, 6, 8, 10 and cols 2-4.
+    `edge + row` = transform(base, row = 2 * (row + 1), col = col + 1)
+  )
+
+  # Buffers are not part of the statistical design, so every one of these must
+  # return the unbuffered design's value.
+  for (type in names(displaced)) {
+    expect_equal(
+      calculate_efficiency_factor(displaced[[type]], treatment),
+      expected
+    )
+  }
+})
+
+test_that("efficiency factor computes on a grid with a genuine hole", {
+  # A plot missing mid-grid is a real input class - an irregular trial edge, a
+  # road - and is not a buffer, so it must be scored on the coordinates it has
+  # rather than closed up. Here Z has no empty columns, so this is a separate
+  # path from the offset case above: solve(), not pseudo_inverse().
+  # The treatments are shuffled rather than laid out systematically:
+  # `rep(LETTERS[1:6], 4)` down a 4 x 6 grid puts each treatment in only two of
+  # the four rows, which confounds treatment with row and leaves the design with
+  # no efficiency factor to report at all.
+  set.seed(2)
+  full <- initialise_design_df(
+    items = sample(rep(LETTERS[1:6], 4)),
+    nrows = 4,
+    ncols = 6
+  )
+  holed <- full[!(full$row == 2 & full$col == 3), ]
+
+  expect_no_error(holed_ef <- calculate_efficiency_factor(holed, treatment))
+  expect_lte(holed_ef, 1)
+  # Both values measured at the coordinate-indexed fix and pinned here, so a
+  # silent change in either is caught.
+  expect_equal(
+    calculate_efficiency_factor(full, treatment),
+    0.5949797,
+    tolerance = 1e-6
+  )
+  expect_equal(holed_ef, 0.5764686, tolerance = 1e-6)
+})
diff --git a/tests/testthat/test-objective_functions.R b/tests/testthat/test-objective_functions.R
index c8077599..c0c86753 100644
--- a/tests/testthat/test-objective_functions.R
+++ b/tests/testthat/test-objective_functions.R
@@ -276,7 +276,7 @@ test_that("objective_function_piepho works with basic design", {
   )
 
   expect_type(result, "list")
-  expect_named(result, c("score", "ed", "bal", "adj", "nb", "components"))
+  expect_named(result, c("score", "ed", "ed_per_grid", "bal", "adj", "nb", "components"))
   expect_type(result$score, "double")
   expect_length(result$score, 1)
   expect_type(result$ed, "list")
@@ -335,10 +335,12 @@ test_that("objective_function_piepho handles different treatment replication pat
   expect_type(result_3_reps, "list")
   expect_type(result_4_reps, "list")
 
-  # Check that ed structure reflects different replication patterns
-  expect_true("2" %in% names(result_2_reps$ed))
-  expect_true("3" %in% names(result_3_reps$ed))
-  expect_true("4" %in% names(result_4_reps$ed))
+  # `ed` is keyed by grid, then by replication class within that grid. A design
+  # with no `by` is a single grid named "1".
+  expect_named(result_2_reps$ed, "1")
+  expect_true("2" %in% names(result_2_reps$ed[["1"]]))
+  expect_true("3" %in% names(result_3_reps$ed[["1"]]))
+  expect_true("4" %in% names(result_4_reps$ed[["1"]]))
 
   # Scores should be different for different replication patterns
   expect_false(identical(result_2_reps$score, result_3_reps$score))
@@ -376,7 +378,7 @@ test_that("objective_function_piepho handles incremental calculation with curren
   expect_type(incremental_result, "list")
   expect_named(
     incremental_result,
-    c("score", "ed", "bal", "adj", "nb", "components")
+    c("score", "ed", "ed_per_grid", "bal", "adj", "nb", "components")
   )
 
   # Test that incremental calculation works differently from full calculation
@@ -405,7 +407,7 @@ test_that("objective_function_piepho works without pair_mapping", {
 
   result <- objective_function_piepho(design_df, "treatment", c("row", "col"))
   expect_type(result, "list")
-  expect_named(result, c("score", "ed", "bal", "adj", "nb", "components"))
+  expect_named(result, c("score", "ed", "ed_per_grid", "bal", "adj", "nb", "components"))
 })
 
 test_that("objective_function_piepho handles different spatial column configurations", {
@@ -484,7 +486,7 @@ test_that("objective_function_piepho uses custom row and column names", {
     col_column = "Column"
   )
   expect_type(result, "list")
-  expect_named(result, c("score", "ed", "bal", "adj", "nb", "components"))
+  expect_named(result, c("score", "ed", "ed_per_grid", "bal", "adj", "nb", "components"))
 })
 
 test_that("objective_function_piepho score is properly rounded to 10 decimal places", {
@@ -534,7 +536,7 @@ test_that("objective_function_piepho handles designs with missing values", {
     pair_mapping = pair_mapping
   )
   expect_type(result, "list")
-  expect_named(result, c("score", "ed", "bal", "adj", "nb", "components"))
+  expect_named(result, c("score", "ed", "ed_per_grid", "bal", "adj", "nb", "components"))
 })
 
 test_that("objective_function_piepho errors on single treatment design", {
@@ -575,13 +577,132 @@ test_that("objective_function_piepho individual components are reasonable", {
   expect_true(is.finite(result$nb$var))
   expect_gte(result$nb$var, 0) # Variance should be non-negative
 
-  # Check ed component structure
+  # Check ed component structure: one entry per grid, each keyed by replication
+  # class.
   expect_type(result$ed, "list")
-  for (ed_rep in result$ed) {
-    expect_named(ed_rep, c("msts", "min_mst", "min_items"))
-    expect_true(is.finite(ed_rep$min_mst))
-    expect_gte(ed_rep$min_mst, 0) # MST should be non-negative
+  for (grid_ed in result$ed) {
+    for (ed_rep in grid_ed) {
+      expect_named(ed_rep, c("msts", "min_mst", "min_items"))
+      expect_true(is.finite(ed_rep$min_mst))
+      expect_gte(ed_rep$min_mst, 0) # MST should be non-negative
+    }
   }
+
+  # One evenness score per grid, summing to the reported component.
+  expect_length(result$ed_per_grid, length(result$ed))
+  expect_equal(
+    sum(result$ed_per_grid),
+    result$components[["even_distribution"]]
+  )
+})
+
+test_that("objective_function_piepho scores evenness per grid and sums them", {
+  # A p-rep multi-environment trial: partial replication *within* each site,
+  # overlapping but unequal subsets across them.
+  set.seed(7)
+  d <- rbind(
+    data.frame(
+      site = "a",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(c(
+        rep(c("A", "B", "C", "D"), 2),
+        "E", "F", "G", "H"
+      ))
+    ),
+    data.frame(
+      site = "b",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(c(rep(c("E", "F", "G"), 3), "A", "B", "I"))
+    )
+  )
+
+  res <- objective_function_piepho(
+    d,
+    "treatment",
+    c("row", "col", "site"),
+    by = "site"
+  )
+
+  expect_named(res$ed_per_grid, c("a", "b"))
+  expect_equal(
+    sum(res$ed_per_grid),
+    res$components[["even_distribution"]]
+  )
+  # Reported per grid as well as in total, so a site can be acted on.
+  expect_equal(
+    res$components[["even_distribution_a"]],
+    res$ed_per_grid[["a"]]
+  )
+
+  # Each grid's evenness is exactly what that grid scores on its own: a
+  # spanning tree measures distance between plots, and there is none between
+  # plots at different sites.
+  alone <- objective_function_piepho(
+    d[d$site == "a", ],
+    "treatment",
+    c("row", "col")
+  )
+  expect_equal(
+    res$ed_per_grid[["a"]],
+    alone$components[["even_distribution"]]
+  )
+
+  # A swap inside one grid cannot change another grid's evenness, which is why
+  # summing the per-grid scores optimises each grid independently.
+  moved <- d
+  i <- which(moved$site == "a")[c(1, 5)]
+  moved$treatment[i] <- moved$treatment[rev(i)]
+  after <- objective_function_piepho(
+    moved,
+    "treatment",
+    c("row", "col", "site"),
+    by = "site"
+  )
+  expect_equal(after$ed_per_grid[["b"]], res$ed_per_grid[["b"]])
+})
+
+test_that("a grid with nothing replicated inside it contributes no evenness", {
+  # Evenness measures how far apart a treatment's replicates are, so a grid
+  # with no replicated treatment has nothing to measure. It must contribute 0
+  # rather than 1/0, which would make the whole score Inf and leave the
+  # optimiser with nothing to compare.
+  set.seed(7)
+  d <- rbind(
+    data.frame(
+      site = "a",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(c(
+        rep(c("A", "B", "C", "D"), 2),
+        "E", "F", "G", "H"
+      ))
+    ),
+    data.frame(
+      site = "b",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = LETTERS[1:12] # every entry once at this site
+    )
+  )
+
+  res <- objective_function_piepho(
+    d,
+    "treatment",
+    c("row", "col", "site"),
+    by = "site"
+  )
+  expect_equal(res$ed_per_grid[["b"]], 0)
+  expect_gt(res$ed_per_grid[["a"]], 0)
+  expect_true(is.finite(res$score))
+
+  # The same rule covers a single grid with no replication at all, which
+  # previously scored Inf.
+  unreplicated <- initialise_design_df(as.character(1:12), 4, 3)
+  single <- objective_function_piepho(
+    unreplicated,
+    "treatment",
+    c("row", "col")
+  )
+  expect_equal(single$components[["even_distribution"]], 0)
+  expect_true(is.finite(single$score))
 })
 
 test_that("objective_function_piepho handles extra parameters via ...", {
@@ -606,7 +727,7 @@ test_that("objective_function_piepho handles extra parameters via ...", {
   })
 
   expect_type(result, "list")
-  expect_named(result, c("score", "ed", "bal", "adj", "nb", "components"))
+  expect_named(result, c("score", "ed", "ed_per_grid", "bal", "adj", "nb", "components"))
 })
 
 test_that("objective_function_factorial works", {
diff --git a/tests/testthat/test-speed.R b/tests/testthat/test-speed.R
index 270ce854..a0cb564c 100644
--- a/tests/testthat/test-speed.R
+++ b/tests/testthat/test-speed.R
@@ -913,39 +913,28 @@ test_that("autoplot handles factor columns with custom column names", {
   )
 })
 
-test_that("autoplot fails with factor columns with character levels", {
-  # Sample data with factor columns that have character levels
-  # Too hard to predict plot layout with character levels
+test_that("speed reports non-numeric row/col labels clearly", {
+  # Row/col labels like "R1"/"C1" cannot index a grid. This used to surface as
+  # five coercion warnings followed by a cryptic "invalid 'nrow' value (too
+  # large or NA)" from matrix(); build_design_matrix() now names the actual
+  # problem once, with no warnings.
   test_data <- data.frame(
     row = factor(rep(paste0("R", 1:5), times = 4)),
     col = factor(rep(paste0("C", 1:4), each = 5)),
     treatment = rep(LETTERS[1:4], 5)
   )
 
-  expect_warning(
-    expect_warning(
-      expect_warning(
-        expect_warning(
-          expect_warning(
-            expect_error(
-              speed(
-                data = test_data,
-                swap = "treatment",
-                iterations = 100,
-                seed = 42,
-                quiet = TRUE
-              ),
-              "invalid 'nrow' value \\(too large or NA\\)"
-            ),
-            "NAs introduced by coercion"
-          ),
-          "no non-missing arguments to max; returning -Inf"
-        ),
-        "no non-missing arguments to max; returning -Inf"
+  expect_no_warning(
+    expect_error(
+      speed(
+        data = test_data,
+        swap = "treatment",
+        iterations = 100,
+        seed = 42,
+        quiet = TRUE
       ),
-      "NAs introduced by coercion to integer range"
-    ),
-    "NAs introduced by coercion"
+      "must be numeric, or coercible to numeric"
+    )
   )
 })
 
@@ -1428,6 +1417,35 @@ test_that("speed runs with grid_factors", {
   expect_equal(result$score, 1)
 })
 
+test_that("`grid_factors$by` is checked before optimising", {
+  test_data <- data.frame(
+    lane = rep(1:5, times = 4),
+    position = rep(1:4, each = 5),
+    treatment = rep(LETTERS[1:4], 5)
+  )
+
+  speed_by <- function(by) {
+    return(speed(
+      data = test_data,
+      swap = "treatment",
+      spatial_factors = ~ lane + position,
+      grid_factors = list(dim1 = "lane", dim2 = "position", by = by),
+      iterations = 10,
+      seed = 42,
+      quiet = TRUE
+    ))
+  }
+
+  # A mistyped column would otherwise be ignored and every grid pooled
+  expect_error(speed_by("site"), "not found in", fixed = TRUE)
+  expect_error(
+    speed_by(c("site", "block")),
+    "must be a single column name",
+    fixed = TRUE
+  )
+  expect_error(speed_by(1), "must be a single column name", fixed = TRUE)
+})
+
 test_that("speed handles MET", {
   # 5 sites, 100 treatments, 7 total reps
   # 5x28x5
diff --git a/tests/testthat/test-summary.R b/tests/testthat/test-summary.R
index ad8bb1a9..f3b9066c 100644
--- a/tests/testthat/test-summary.R
+++ b/tests/testthat/test-summary.R
@@ -1,6 +1,16 @@
 # Tests for summary.design / print.summary.design (Phase 3: Structure +
 # Optimisation + flags). Evaluation metrics are covered separately.
 
+# The `grid` argument the evaluation helpers take: a grid_indices() list - one
+# entry per grid - or the reason there is no grid at all. summary.design() builds
+# it inline, once per call; tests calling a helper directly need the same value.
+grid_or_reason <- function(df, rc = "row", cc = "col", by = NULL) {
+  return(tryCatch(
+    grid_indices(df, row_column = rc, col_column = cc, by = by),
+    speed_grid_error = function(e) return(e$reason)
+  ))
+}
+
 simple_design <- function(iterations = 200, seed = 42) {
   d <- data.frame(
     row = rep(1:4, times = 3),
@@ -119,8 +129,7 @@ test_that("score components are faithful for non-default objectives (piepho)", {
     quiet = TRUE
   )
   sc <- summary(r)$per_level[[1]]$score
-  # Piepho exposes four additive components that sum to its score - the bug this
-  # fixes was the old adjacency+balance recompute not matching the piepho score.
+  # Piepho exposes four additive components, which must sum to its score.
   expect_named(
     sc$components,
     c("neighbour_balance", "even_distribution", "balance", "adjacency")
@@ -296,9 +305,12 @@ test_that("buffers are excluded from summary() and print() entirely", {
 })
 
 test_that("neighbour balance is unaffected by buffers (KNOWN_ISSUES.md #1a)", {
-  # add_buffers("edge") shifts row/col by 1; .neighbour_balance() must use the
-  # (offset-invariant) layout$nrow/ncol rather than max(row)/max(col), or it
-  # reshapes the treatment grid with the wrong dimensions.
+  # add_buffers("edge") shifts row/col by 1. The grid is built from the
+  # coordinates, so the offset just leaves an empty leading row and column,
+  # which contribute no pairs. "edge" only offsets - it inserts no gap between
+  # the design's own plots - so the counts must match the unbuffered design
+  # exactly. (A "row" or "block" buffer does insert a gap, and there the counts
+  # are expected to differ: plots either side of a buffer are not neighbours.)
   d <- data.frame(
     row = rep(1:4, times = 3),
     col = rep(1:3, each = 4),
@@ -640,16 +652,12 @@ test_that("efficiency is opt-in and guarded", {
   expect_true(is.finite(on$value))
 
   # Guard: < 3 treatments returns NA with a reason rather than erroring.
-  two <- .efficiency_factor(
-    data.frame(
-      row = rep(1:2, 2),
-      col = rep(1:2, each = 2),
-      treatment = rep(c("A", "B"), 2)
-    ),
-    "treatment",
-    "row",
-    "col"
+  d2 <- data.frame(
+    row = rep(1:2, 2),
+    col = rep(1:2, each = 2),
+    treatment = rep(c("A", "B"), 2)
   )
+  two <- .efficiency_factor(d2, "treatment", "row", "col", grid_or_reason(d2))
   expect_false(two$available)
 })
 
@@ -727,21 +735,44 @@ test_that("neighbour balance separates self-adjacency from distinct-pair counts"
   expect_true(nb$min_pair_count >= 0)
   expect_true(nb$self_adjacent >= 0)
 
-  # Cross-check directly against create_pair_mapping()/calculate_nb() over all
-  # 6 possible pairs for this 3-treatment design.
-  dm <- matrix(r$design_df$treatment, nrow = 4, ncol = 3)
-  pm <- create_pair_mapping(r$design_df$treatment)
-  raw <- calculate_nb(dm, pm)
-  all_counts <- setNames(rep(0L, length(unique(pm))), unique(pm))
-  all_counts[names(raw$nb)] <- raw$nb
-  self <- c("A,A", "B,B", "C,C")
-  distinct <- setdiff(names(all_counts), self)
-
-  expect_equal(nb$self_adjacent, sum(all_counts[self]))
-  expect_equal(nb$min_pair_count, min(all_counts[distinct]))
-  expect_equal(nb$max_pair_count, max(all_counts[distinct]))
-  expect_equal(nb$pair_var, var(all_counts[distinct]))
-  expect_equal(nb$n_zero_pairs, sum(all_counts[distinct] == 0))
+  # Cross-check by walking the (row, col) coordinates directly. Deliberately NOT
+  # matrix(treatment, nrow, ncol): reshaping is what the implementation must
+  # avoid, so an expectation built that way would only confirm itself.
+  coords <- r$design_df
+  rr <- as.numeric(as.character(coords$row))
+  cc <- as.numeric(as.character(coords$col))
+  trt <- as.character(coords$treatment)
+  at <- function(i, j) {
+    k <- which(rr == i & cc == j)
+    if (length(k) == 1) trt[k] else NA_character_
+  }
+  pm_levels <- sort(unique(trt))
+  counts <- list()
+  for (k in seq_along(trt)) {
+    for (nb_ij in list(c(rr[k], cc[k] + 1), c(rr[k] + 1, cc[k]))) {
+      other <- at(nb_ij[1], nb_ij[2])
+      if (!is.na(other)) {
+        key <- paste(sort(c(trt[k], other)), collapse = ",")
+        prev <- counts[[key]]
+        counts[[key]] <- if (is.null(prev)) 1L else prev + 1L
+      }
+    }
+  }
+  self_keys <- paste(pm_levels, pm_levels, sep = ",")
+  all_pairs <- combn(pm_levels, 2, function(p) paste(p, collapse = ","))
+  get <- function(k) if (is.null(counts[[k]])) 0L else counts[[k]]
+  self_total <- sum(vapply(self_keys, get, integer(1)))
+  pair_counts <- vapply(all_pairs, get, integer(1))
+
+  expect_equal(nb$self_adjacent, self_total)
+  expect_equal(nb$min_pair_count, min(pair_counts))
+  expect_equal(nb$max_pair_count, max(pair_counts))
+  expect_equal(nb$pair_var, var(pair_counts))
+  expect_equal(nb$n_zero_pairs, sum(pair_counts == 0))
+
+  # The design is 4x3, so a column-major reshape would disagree - which is the
+  # whole point of the coordinate-based grid.
+  expect_equal(nb$self_adjacent, calculate_adjacency_score(coords, "treatment"))
 })
 
 test_that("a non-zero self-adjacency count is highlighted in the printed output", {
@@ -892,6 +923,285 @@ test_that("designs without a row/column grid summarise and print without a grid"
   expect_match(out, "Neighbour:\\s+no row/column factors")
 })
 
+test_that("grid_index() conditions carry the reason summary() reports", {
+  # summary() reports `reason` in place of a metric it cannot compute, so every
+  # condition class must supply one.
+  ok <- data.frame(row = c(1, 1, 2, 2), col = c(1, 2, 1, 2), treatment = "A")
+  expect_false(is.character(grid_or_reason(ok)))
+
+  # A genuine hole (missing plot) is still a single grid - coordinates stay
+  # unique, so the metrics must not refuse it.
+  expect_false(is.character(grid_or_reason(ok[-2, ])))
+
+  # Factor coordinates whose level order is lexical are still fine.
+  lex <- data.frame(
+    row = factor(c(1, 2, 10)),
+    col = factor(c(1, 1, 1)),
+    treatment = "A"
+  )
+  expect_false(is.character(grid_or_reason(lex)))
+
+  expect_equal(grid_or_reason(ok, "nope"), "no row/column factors")
+  expect_match(
+    grid_or_reason(data.frame(row = c("A", "B"), col = c(1, 1))),
+    "labels are not numeric"
+  )
+  expect_match(
+    grid_or_reason(data.frame(row = c(0, 1), col = c(1, 1))),
+    "not positive whole numbers"
+  )
+  expect_match(
+    grid_or_reason(ok[c(1, 1, 2), ]),
+    "duplicate .row./.col. coordinates"
+  )
+})
+
+test_that("multi-site (MET) designs summarise instead of erroring", {
+  # initialise_design_df(designs = ) reuses row/col per site, so a MET design has
+  # duplicate coordinates and cannot be placed on one grid. Every grid metric is
+  # wrong for it (they pool sites that share no edge), so all must report a
+  # reason. Before this gate, summary() errored outright.
+  met <- initialise_design_df(
+    items = rep(LETTERS[1:3], 6),
+    designs = list(
+      a = list(nrows = 3, ncols = 2),
+      b = list(nrows = 3, ncols = 4)
+    )
+  )
+  expect_true(anyDuplicated(met[c("row", "col")]) > 0)
+
+  r <- speed(
+    met,
+    swap = "treatment",
+    swap_within = "site",
+    spatial_factors = ~ row + col + site,
+    optimise_params = optim_params(adj_weight = 0),
+    iterations = 50,
+    seed = 1,
+    quiet = TRUE
+  )
+  s <- expect_no_error(summary(r, efficiency = TRUE))
+
+  # No single nrow x ncol describes two grids of different shapes.
+  expect_false(s$layout$has_grid)
+  expect_match(s$layout$grid_reason, "multi-site")
+  expect_true(is.na(s$layout$nrow))
+  expect_equal(s$layout$n_plots, 18)
+
+  e <- s$per_level[[1]]$evaluation
+  for (metric in c("neighbour", "replicate_span", "efficiency")) {
+    expect_false(e[[metric]]$available, info = metric)
+    expect_match(e[[metric]]$reason, "duplicate", info = metric)
+  }
+
+  # Non-grid metrics still work - only the grid ones are withheld.
+  expect_equal(s$per_level[[1]]$n_treatments, 3)
+
+  out <- capture_output(print(s))
+  expect_match(out, "Layout:\\s+18 plots")
+  expect_no_match(out, "rows x")
+  expect_match(out, "Neighbour:\\s+duplicate")
+})
+
+test_that("efficiency is withheld rather than reported above 1 for MET designs", {
+  # An efficiency factor is a property of one experiment's information matrix,
+  # so a multi-site frame has no single answer. Nothing about the arithmetic
+  # stops it being computed - duplicate coordinates just pool the sites into one
+  # row/column model - so the refusal has to be explicit. Before it, this
+  # returned a value above 1, impossible for an efficiency factor.
+  met <- initialise_design_df(
+    items = rep(LETTERS[1:3], 6),
+    designs = list(
+      a = list(nrows = 3, ncols = 2),
+      b = list(nrows = 3, ncols = 4)
+    )
+  )
+  expect_error(
+    calculate_efficiency_factor(met, treatment),
+    class = "speed_grid_duplicate"
+  )
+  gate <- .efficiency_factor(
+    met,
+    "treatment",
+    "row",
+    "col",
+    grid_or_reason(met)
+  )
+  expect_false(gate$available)
+  expect_match(gate$reason, "duplicate")
+})
+
+test_that("a MET design reports one efficiency per site and never a pooled one", {
+  set.seed(1)
+  d <- rbind(
+    data.frame(
+      site = "a",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(rep(LETTERS[1:6], 2))
+    ),
+    data.frame(
+      site = "b",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(rep(LETTERS[1:6], 2))
+    )
+  )
+  r <- speed(
+    d,
+    swap = "treatment",
+    swap_within = "site",
+    spatial_factors = ~ row + col + site,
+    grid_factors = list(dim1 = "row", dim2 = "col", by = "site"),
+    iterations = 100,
+    seed = 42,
+    quiet = TRUE
+  )
+  s <- summary(r, efficiency = TRUE)
+
+  # The grouping is recorded, not guessed from the column name.
+  expect_equal(r$metadata$grid_by, "site")
+  expect_false(s$layout$has_grid)
+  expect_equal(s$layout$n_grids, 2L)
+  expect_match(s$layout$grid_reason, "2 grids")
+  # nrow/ncol stay NA: no single shape describes a design spanning two grids.
+  expect_true(is.na(s$layout$nrow))
+
+  ef <- s$per_level[[1]]$evaluation$efficiency
+  expect_named(ef$per_grid, c("a", "b"))
+  expect_equal(ef$grid_by, "site")
+  expect_null(ef$value) # no pooled value, ever
+  for (one in ef$per_grid) {
+    expect_true(one$available)
+    expect_lte(one$value, 1)
+  }
+  # Each site's value is that site scored on its own.
+  expect_equal(
+    ef$per_grid$a$value,
+    calculate_efficiency_factor(
+      r$design_df[r$design_df$site == "a", ],
+      treatment
+    )
+  )
+
+  out <- capture_output(print(s))
+  expect_match(out, "per `site`")
+
+  # Spans pool sites if computed naively - two sites' row 3 are not one plot
+  # apart - so they are withheld with a reason instead.
+  spans <- s$per_level[[1]]$evaluation$replicate_span
+  expect_false(spans$available)
+  expect_match(spans$reason, "spans 2 grids")
+})
+
+test_that("a site that cannot support an efficiency factor reports only for itself", {
+  # Site b puts each treatment in one column, so its treatment contrasts are
+  # confounded with the column effects. That must not withhold site a's value.
+  d <- rbind(
+    data.frame(
+      site = "a",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = c("A", "B", "C", "D", "C", "D", "A", "B", "B", "A", "D", "C")
+    ),
+    data.frame(
+      site = "b",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = rep(c("A", "B", "C"), each = 4)
+    )
+  )
+  grid <- grid_or_reason(d, by = "site")
+  ef <- .efficiency_factor(d, "treatment", "row", "col", grid)
+
+  expect_true(ef$per_grid$a$available)
+  expect_false(ef$per_grid$b$available)
+  expect_match(ef$per_grid$b$reason, "not estimable")
+  # `available` is TRUE overall because at least one site reported.
+  expect_true(ef$available)
+})
+
+test_that("neighbour balance sums across grids rather than pooling them", {
+  set.seed(3)
+  d <- rbind(
+    data.frame(
+      site = "a",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(rep(LETTERS[1:6], 2))
+    ),
+    data.frame(
+      site = "b",
+      expand.grid(row = 1:4, col = 1:3),
+      treatment = sample(rep(LETTERS[1:6], 2))
+    )
+  )
+  both <- .neighbour_balance(
+    d,
+    "treatment",
+    "row",
+    "col",
+    grid_or_reason(d, by = "site")
+  )
+  expect_true(both$available)
+
+  # The self-adjacency count is an edge count, so it is exactly the sum of the
+  # per-site counts - no edge crosses a site boundary.
+  per_site <- vapply(
+    split(d, d$site),
+    function(s) {
+      return(.neighbour_balance(
+        s,
+        "treatment",
+        "row",
+        "col",
+        grid_or_reason(s)
+      )$self_adjacent)
+    },
+    numeric(1)
+  )
+  expect_equal(both$self_adjacent, sum(per_site))
+})
+
+test_that("non-numeric row/col labels are reported, not coerced silently", {
+  d <- data.frame(
+    row = rep(c("A", "B", "C"), each = 4),
+    col = rep(c("w", "x", "y", "z"), 3),
+    treatment = rep(LETTERS[1:4], 3)
+  )
+  r <- speed(
+    d,
+    swap = "treatment",
+    spatial_factors = ~ row + col,
+    optimise_params = optim_params(adj_weight = 0),
+    iterations = 50,
+    seed = 1,
+    quiet = TRUE
+  )
+
+  # Previously this errored from build_design_matrix(), and .replicate_spans()
+  # leaked two "NAs introduced by coercion" warnings on the way.
+  s <- expect_no_warning(expect_no_error(summary(r)))
+
+  expect_false(s$layout$has_grid)
+  e <- s$per_level[[1]]$evaluation
+  expect_match(e$neighbour$reason, "labels are not numeric")
+  expect_match(e$replicate_span$reason, "labels are not numeric")
+})
+
+test_that("split-plot designs keep their grid metrics at every level", {
+  # The gate must not withhold metrics from a legitimate hierarchical design:
+  # a split-plot shares one grid, and both levels' plots have unique
+  # coordinates, so nothing here spans multiple grids.
+  s <- summary(split_plot_design(), efficiency = TRUE)
+
+  expect_true(s$layout$has_grid)
+  expect_true(is.na(s$layout$grid_reason))
+  expect_equal(c(s$layout$nrow, s$layout$ncol), c(6, 4))
+
+  expect_named(s$per_level, c("wp", "sp"))
+  for (lv in c("wp", "sp")) {
+    e <- s$per_level[[lv]]$evaluation
+    expect_true(e$neighbour$available, info = lv)
+    expect_true(e$replicate_span$available, info = lv)
+  }
+})
+
 test_that("the disconnected flag names the affected levels for hierarchical designs", {
   s <- summary(split_plot_design())
   s$flags$disconnected <- c("wp", "sp")
@@ -938,13 +1248,14 @@ test_that("buffer removal drops the unused factor level from a factor swap colum
 test_that("evaluation helpers return a reason instead of erroring on unmet assumptions", {
   no_grid <- data.frame(a = 1:3, treatment = c("A", "B", "C"))
 
-  rs <- .replicate_spans(no_grid, "treatment", "row", "col")
+  ng <- grid_or_reason(no_grid)
+  rs <- .replicate_spans(no_grid, "treatment", "row", "col", ng)
   expect_false(rs$available)
   expect_equal(rs$reason, "no row/column factors")
 
-  ef <- .efficiency_factor(no_grid, "treatment", "row", "col")
+  ef <- .efficiency_factor(no_grid, "treatment", "row", "col", ng)
   expect_false(ef$available)
-  expect_equal(ef$reason, "requires a row/column grid")
+  expect_equal(ef$reason, "no row/column factors")
 
   # A contrast needs two treatments to be a contrast.
   one <- .design_connectedness(
@@ -982,15 +1293,54 @@ test_that(".design_connectedness is trivially connected with no nuisance factors
 })
 
 test_that(".efficiency_factor reports a reason when the computation fails", {
-  # A 1x1 "grid" leaves no row/column effects to fit, so
-  # calculate_efficiency_factor() errors; the wrapper must absorb that.
-  degenerate <- data.frame(
-    row = rep(1, 3),
-    col = rep(1, 3),
-    treatment = c("A", "B", "C")
+  # The wrapper must absorb an error from the underlying metric rather than
+  # propagate it. Mocked rather than provoked with a degenerate design, because a
+  # design degenerate enough to fail the metric (e.g. a 1x1 grid, three plots
+  # sharing one coordinate) is rejected by grid_index() first, leaving this
+  # backstop untested.
+  d <- data.frame(
+    row = rep(1:2, 3),
+    col = rep(1:3, each = 2),
+    treatment = rep(c("A", "B", "C"), 2)
+  )
+  grid <- grid_or_reason(d)
+  expect_false(is.character(grid))
+
+  local_mocked_bindings(
+    calculate_efficiency_factor = function(...) stop("cannot compute")
   )
-  ef <- .efficiency_factor(degenerate, "treatment", "row", "col")
+  ef <- .efficiency_factor(d, "treatment", "row", "col", grid)
 
   expect_false(ef$available)
   expect_equal(ef$reason, "could not be computed for this design")
 })
+
+test_that("neighbour balance reads coordinates, not the data frame order", {
+  # A 3x8 design optimised to a genuine zero self-adjacency. A column-major
+  # reshape of a row-major frame scrambles it and invents adjacencies, so this
+  # is the regression test for that class of bug: the reported figure must
+  # agree with the objective's own adjacency score for the same design.
+  d <- speed(
+    initialise_design_df(items = rep(LETTERS[1:6], 4), nrows = 3, ncols = 8),
+    swap = "treatment",
+    iterations = 3000,
+    seed = 11,
+    quiet = TRUE
+  )
+  nb <- summary(d)$per_level[[1]]$evaluation$neighbour
+
+  expect_equal(nb$self_adjacent, 0)
+  expect_equal(
+    nb$self_adjacent,
+    calculate_adjacency_score(d$design_df, "treatment")
+  )
+
+  # Row order must not matter: the same design with its rows reversed is the
+  # same field layout, and must report the same neighbour balance.
+  reversed <- d
+  reversed$design_df <- d$design_df[rev(seq_len(nrow(d$design_df))), ]
+  expect_equal(
+    summary(reversed)$per_level[[1]]$evaluation$neighbour,
+    nb
+  )
+})
diff --git a/vignettes/met.qmd b/vignettes/met.qmd
index 1b60ab36..6a5df19e 100644
--- a/vignettes/met.qmd
+++ b/vignettes/met.qmd
@@ -99,6 +99,10 @@ plot_layout(met_design, "block")
 For MET designs, we use lists of named arguments to specify the hierarchical structure. The `optimise` parameter
 defines what to optimise and constraints at each level.
 
+Each site numbers its rows and columns from 1, so several plots share the same `row`/`col` pair. A MET is not
+one grid but several, which share a treatment set and never share an edge. Tell `speed()` which column
+separates them with the `by` element of `grid_factors`:
+
 ```{r met-example}
 optimise <- list(
   connectivity = list(spatial_factors = ~site),
@@ -110,6 +114,7 @@ met_result <- speed(
   swap = "treatment",
   early_stop_iterations = 5000,
   optimise = optimise,
+  grid_factors = list(dim1 = "row", dim2 = "col", by = "site"),
   optimise_params = optim_params(random_initialisation = TRUE, adj_weight = 0),
   seed = 112
 )
@@ -117,6 +122,10 @@ met_result <- speed(
 met_result
 ```
 
+Each site is then scored on its own and the results combined, so no two plots at different sites are ever
+treated as neighbours. Without `by`, a design whose sites reuse coordinates is refused rather than silently
+pooled - the grid-based diagnostics would otherwise describe a layout that does not exist.
+
 ## Output of the Optimisation
 
 The output shows optimisation results for the design. The score and iterations are combined for the entire
@@ -203,6 +212,7 @@ met_result <- speed(
   early_stop_iterations = 10000,
   iterations = 50000,
   optimise = optimise,
+  grid_factors = list(dim1 = "row", dim2 = "col", by = "site"),
   optimise_params = optim_params(random_initialisation = 10, adj_weight = 0),
   seed = 112
 )
@@ -330,6 +340,7 @@ met_result <- speed(
   early_stop_iterations = 8000,
   iterations = 50000,
   optimise = optimise,
+  grid_factors = list(dim1 = "row", dim2 = "col", by = "site"),
   optimise_params = optim_params(random_initialisation = 30, adj_weight = 0),
   seed = 112
 )