Build design grids from plot coordinates, and support designs spanning multiple grids

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

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.
 

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))
 }

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)
+}