diff --git a/.Rbuildignore b/.Rbuildignore
index d361182..c83037c 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -21,6 +21,7 @@ codemeta.json
 ^Meta$
 ^CODE_OF_CONDUCT\.md$
 ^doc$
+^\.dir-locals\.el$
 \.[Rr][Dd]ata$
 \.[Rr]ds$
 \.[Xx]lsx$
diff --git a/DESCRIPTION b/DESCRIPTION
index d0b2fda..c4ad802 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Type: Package
 Package: speed
 Title: Generate Spatially Efficient Experimental Designs
 Version: 0.0.11
-Authors@R: 
+Authors@R:
     c(person(given = "Sam",
              family = "Rogers",
              role = c("aut", "cre"),
@@ -15,18 +15,26 @@ Authors@R:
              family = "Taylor",
              role = "aut",
              email = "julian.taylor@adelaide.edu.au"),
+      person(given = "Kai",
+             family = "Bagley",
+             email = "kai.bagley@curtin.edu.au",
+             role = "aut",
+             comment = c(ORCID = "0009-0004-6579-6959")),
       person(given = "University of Adelaide",
              role = c("cph", "fnd"),
              comment = "https://adelaide.edu.au/"),
       person(given = "Grains Research and Development Corporation",
              role = c("cph", "fnd"),
-             comment = "https://grdc.com.au/"))
-Description: The Speed package optimises spatial experimental designs by rearranging 
-    treatments to improve statistical efficiency while maintaining statistical validity. 
-    It employs customisable optimisation metrics, with a default approach that balances 
-    minimising treatment adjacency and maintaining spatial balance. Users can supply their 
-    own optimisation metrics to tailor designs to specific experimental requirements. 
-    The package also includes visualisation tools that allow users to inspect 
+             comment = "https://grdc.com.au/"),
+      person(given = "Curtin University",
+             role = c("cph"),
+             comment = c("http://www.curtin.edu.au/", ROR = "02n415q13")))
+Description: The Speed package optimises spatial experimental designs by rearranging
+    treatments to improve statistical efficiency while maintaining statistical validity.
+    It employs customisable optimisation metrics, with a default approach that balances
+    minimising treatment adjacency and maintaining spatial balance. Users can supply their
+    own optimisation metrics to tailor designs to specific experimental requirements.
+    The package also includes visualisation tools that allow users to inspect
     optimised designs.
 License: MIT + file LICENSE
 URL: https://biometryhub.github.io/speed/
diff --git a/NAMESPACE b/NAMESPACE
index 64c9265..96d2559 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -11,7 +11,11 @@ export(calculate_balance_score)
 export(calculate_ed)
 export(calculate_efficiency_factor)
 export(calculate_nb)
+export(compute_L_projection)
+export(cor_ar1)
+export(cor_ar1_ar1)
 export(create_pair_mapping)
+export(factorial_contrast_matrix)
 export(get_edges)
 export(get_vertices)
 export(initialise_design_df)
@@ -19,6 +23,7 @@ export(initialise_split_design_df)
 export(initialize_design_df)
 export(objective_function)
 export(objective_function_factorial)
+export(objective_function_info)
 export(objective_function_piepho)
 export(objective_function_signature)
 export(optim_params)
diff --git a/R/objective_fisher.R b/R/objective_fisher.R
new file mode 100644
index 0000000..76dc528
--- /dev/null
+++ b/R/objective_fisher.R
@@ -0,0 +1,580 @@
+#' Objective function using the treatment information matrix
+#'
+#' Creates an objective function that optimises experimental designs using the
+#' Fisher information for treatment contrasts after adjusting for nuisance
+#' effects, and takes an optional spatial covariance structure for spatial
+#' optimisation.
+#'
+#' @param layout_df Data frame representing the spatial layout of the
+#'   experiment.
+#' @param swap Column name to swap, usually the treatment.
+#' @param spatial_cols Column name of the spatial factors.
+#' @param criterion Either \code{"A"} or \code{"D"}, representing A or D
+#'   optimality.
+#'     - A-optimality: Minimises \eqn{\mathrm{tr} \left( \mathcal{I}^-
+#'       \right)}{tr(I⁻)}.
+#'     - D-optimality: Minimises \eqn{-\log \left| \mathcal{I} \right|}{-
+#'       log(|I|)}
+#' @param L_matrix Optional precomputed projection matrix. This is an advanced
+#'   alternative to `Sigma`; the two cannot be supplied together.
+#' @param Sigma Optional covariance matrix whose rows and columns follow the row
+#'   order of `layout_df`. The projection is computed on the first objective
+#'   call and reused during optimisation. If both `Sigma` and `L_matrix` are
+#'   `NULL`, identity covariance is assumed.
+#' @param contrast_matrix Numeric contrast matrix with one column per treatment
+#'   level. Column names, when supplied, must match the treatment levels. Only
+#'   these contrasts are required to be estimable and contribute to the
+#'   optimality criterion. The default targets the complete treatment-contrast
+#'   space.
+#' @param block_column Column name of the design's block factor. Used when
+#'   \code{L_matrix} isn't supplied.
+#' @param nuisance_formula Optional one-sided formula describing nuisance
+#'   effects. Used when \code{L_matrix} is not supplied. If omitted,
+#'   \code{block_column} is treated as a factor.
+#' @param tolerance Relative numerical tolerance for rank and estimability
+#'   checks.
+#' @param current_score_obj The current objective result supplied by [speed()]
+#'   after the initial call. Used to reuse invariant matrix calculations.
+#' @param swapped_items The treatments swapped by [speed()]. Changed plot rows
+#'   are identified from the cached and proposed allocations so repeated
+#'   treatments and multi-swaps are handled exactly.
+#' @param ... Extra parameters passed from \code{speed}.
+#'
+#' @details
+#' This function computes the treatment information matrix:
+#' \deqn{I = X_1^\intercal L X_1}{I = X₁ᵀ L X₁}
+#' Where \eqn{X_1}{X₁} is the treatment design matrix, and \eqn{L} is the
+#' projection that removes nuisance fixed effects from the GLS's inverse
+#' covariance:
+#' \deqn{L = \Sigma^{-1} - \Sigma^{-1} X_2 \left( X_2^\intercal \Sigma^{-1} X_2
+#' \right)^{-1} X_2^\intercal \Sigma^{-1}}{L = Σ⁻¹ - Σ⁻¹ X₂
+#' (X₂ᵀ Σ⁻¹ X₂)⁻¹ X₂ᵀ Σ⁻¹}
+#' The user specifies the spatial covariance \eqn{\Sigma}{Σ} they intend to use
+#' at the time of analysis, e.g. lag-1 autoregressive \eqn{\text{AR}_1 \otimes
+#' \text{AR}_1}{AR₁ ⊗ AR₁}, and this objective will optimise under that
+#' structure.
+#'
+#' @examples
+#' # Small RCBD layout: 6 treatments in 4 blocks of 6 plots
+#' df <- initialise_design_df(
+#'   items = 6, nrows = 4, ncols = 6,
+#'   block_nrows = 1, block_ncols = 6
+#' )
+#'
+#' # Non-spatial: identity covariance
+#' result <- speed(
+#'   df,
+#'   swap = "treatment",
+#'   swap_within = "block",
+#'   obj_function = objective_function_info,
+#'   criterion = "A",
+#'   seed = 42,
+#'   quiet = TRUE
+#' )
+#'
+#' # Spatial: AR(1) x AR(1)
+#' Sigma <- cor_ar1_ar1(
+#'   n_rows = 4, n_cols = 6,
+#'   rho_row = 0.6, rho_col = 0.3
+#' )
+#' result_spatial <- speed(
+#'   df,
+#'   swap = "treatment",
+#'   swap_within = "block",
+#'   obj_function = objective_function_info,
+#'   Sigma = Sigma,
+#'   criterion = "A",
+#'   optimise_params = optim_params(random_initialisation = TRUE),
+#'   seed = 42,
+#'   quiet = TRUE
+#' )
+#'
+#' @return A named list with \code{score}, \code{info_matrix},
+#'   \code{eigenvalues}, contrast diagnostics, and \code{criterion}.
+#'
+#' @export
+objective_function_info <- function(
+  layout_df,
+  swap,
+  spatial_cols,
+  criterion = c("A", "D"),
+  L_matrix = NULL,
+  Sigma = NULL,
+  contrast_matrix = NULL,
+  block_column = "block",
+  nuisance_formula = NULL,
+  tolerance = 1e-9,
+  current_score_obj = NULL,
+  swapped_items = NULL,
+  ...
+) {
+  criterion <- match.arg(criterion)
+
+  fisher_cache <- NULL
+  if (!is.null(swapped_items)) {
+    fisher_cache <- attr(current_score_obj, "fisher_cache", exact = TRUE)
+  }
+  ci <- .compute_info(
+    layout_df,
+    swap,
+    Sigma,
+    L_matrix,
+    block_column,
+    nuisance_formula,
+    tolerance,
+    fisher_cache
+  )
+  diagnostics <- .information_criteria(
+    ci$info, ci$trt_levels, contrast_matrix, tolerance
+  )
+  score <- diagnostics[[paste0(criterion, "_value")]]
+  if (!is.finite(score)) score <- 1e10
+
+  result <- list(
+    score = score,
+    info_matrix = ci$info,
+    eigenvalues = diagnostics$eigenvalues,
+    contrast_matrix = diagnostics$contrast_matrix,
+    contrast_covariance = diagnostics$contrast_covariance,
+    estimable = diagnostics$estimable,
+    rank = diagnostics$rank,
+    criterion = criterion,
+    A_value = diagnostics$A_value,
+    D_value = diagnostics$D_value
+  )
+  attr(result, "fisher_cache") <- ci$cache
+  return(result)
+}
+
+#' Build the treatment indicator matrix X1
+#' @noRd
+.build_treatment_matrix <- function(treatments, trt_levels, n, v) {
+  X1 <- matrix(0, n, v)
+  trt_idx <- match(treatments, trt_levels)
+  if (anyNA(trt_idx)) stop("Treatment allocations cannot contain missing values.", call. = FALSE)
+  X1[cbind(seq_len(n), trt_idx)] <- 1
+  return(X1)
+}
+
+#' Build an orthonormal factorial contrast matrix
+#'
+#' Constructs a contrast basis for selected factorial effects while treatments
+#' remain encoded as atomic combinations during optimisation. The returned rows
+#' are orthonormal, so A-optimality is invariant to the particular contrast
+#' coding used by [stats::model.matrix()].
+#'
+#' @param treatment_df One row per treatment combination, containing the
+#'   component factor columns.
+#' @param formula One-sided formula selecting factorial effects, for example
+#'   \code{~ (stage + cultivar + inoculum)^2}.
+#' @param treatment_column Column containing the atomic treatment labels.
+#' @param tolerance Numerical rank tolerance.
+#'
+#' @return A numeric contrast matrix whose columns are named by treatment level.
+#'
+#' @export
+factorial_contrast_matrix <- function(
+  treatment_df,
+  formula,
+  treatment_column = "treatment",
+  tolerance = 1e-10
+) {
+  if (!is.data.frame(treatment_df)) {
+    stop("`treatment_df` must be a data frame.", call. = FALSE)
+  }
+  if (!inherits(formula, "formula") || length(formula) != 2L) {
+    stop("`formula` must be a one-sided formula.", call. = FALSE)
+  }
+  if (!treatment_column %in% names(treatment_df)) {
+    stop("Treatment column `", treatment_column, "` not found.", call. = FALSE)
+  }
+
+  treatment_levels <- as.character(treatment_df[[treatment_column]])
+  if (anyNA(treatment_levels) || anyDuplicated(treatment_levels)) {
+    stop("`treatment_df` must contain one non-missing row per treatment.", call. = FALSE)
+  }
+
+  formula_vars <- all.vars(formula)
+  missing_vars <- setdiff(formula_vars, names(treatment_df))
+  if (length(missing_vars) > 0L) {
+    stop(
+      "Formula variable(s) not found: ", paste(missing_vars, collapse = ", "),
+      call. = FALSE
+    )
+  }
+
+  model_data <- treatment_df
+  for (column in formula_vars) {
+    if (is.character(model_data[[column]])) {
+      model_data[[column]] <- factor(model_data[[column]])
+    }
+  }
+
+  effect_matrix <- stats::model.matrix(formula, model_data)
+  effect_matrix <- sweep(effect_matrix, 2L, colMeans(effect_matrix), "-")
+  decomposition <- qr(effect_matrix, tol = tolerance)
+  effect_rank <- decomposition$rank
+  if (effect_rank == 0L) {
+    stop("`formula` does not define any treatment contrasts.", call. = FALSE)
+  }
+
+  basis <- qr.Q(decomposition, complete = FALSE)[, seq_len(effect_rank), drop = FALSE]
+  contrasts <- t(basis)
+  rownames(contrasts) <- paste0("contrast_", seq_len(nrow(contrasts)))
+  colnames(contrasts) <- treatment_levels
+  return(contrasts)
+}
+
+#' Prepare a requested contrast space
+#' @noRd
+.prepare_contrast_matrix <- function(contrast_matrix, trt_levels, tolerance) {
+  v <- length(trt_levels)
+  if (is.null(contrast_matrix)) {
+    if (v < 2L) {
+      return(matrix(numeric(0), nrow = 0L, ncol = v))
+    }
+    basis <- qr.Q(qr(stats::contr.helmert(v)), complete = FALSE)
+    contrast_matrix <- t(basis)
+    colnames(contrast_matrix) <- trt_levels
+  } else {
+    if (!is.matrix(contrast_matrix) || !is.numeric(contrast_matrix)) {
+      stop("`contrast_matrix` must be a numeric matrix.", call. = FALSE)
+    }
+    if (any(!is.finite(contrast_matrix))) {
+      stop("`contrast_matrix` must contain only finite values.", call. = FALSE)
+    }
+    if (is.null(colnames(contrast_matrix))) {
+      if (ncol(contrast_matrix) != v) {
+        stop("`contrast_matrix` must have one column per treatment.", call. = FALSE)
+      }
+      colnames(contrast_matrix) <- trt_levels
+    } else {
+      missing <- setdiff(trt_levels, colnames(contrast_matrix))
+      extra <- setdiff(colnames(contrast_matrix), trt_levels)
+      if (length(missing) > 0L || length(extra) > 0L) {
+        stop("Contrast columns must match the treatment levels.", call. = FALSE)
+      }
+      contrast_matrix <- contrast_matrix[, trt_levels, drop = FALSE]
+    }
+  }
+
+  if (nrow(contrast_matrix) > 0L && qr(t(contrast_matrix), tol = tolerance)$rank < nrow(contrast_matrix)) {
+    stop("Rows of `contrast_matrix` must be linearly independent.", call. = FALSE)
+  }
+  return(contrast_matrix)
+}
+
+#' Calculate information criteria for a requested contrast space
+#' @noRd
+.information_criteria <- function(info, trt_levels, contrast_matrix, tolerance) {
+  info <- (info + t(info)) / 2
+  eig <- eigen(info, symmetric = TRUE)
+  scale <- max(abs(eig$values))
+  keep <- if (scale == 0) rep(FALSE, length(eig$values)) else eig$values > scale * tolerance
+  rank <- sum(keep)
+  pos_eig <- sort(eig$values[keep], decreasing = TRUE)
+
+  info_inverse <- matrix(0, nrow(info), ncol(info))
+  if (rank > 0L) {
+    vectors <- eig$vectors[, keep, drop = FALSE]
+    info_inverse <- vectors %*% diag(1 / eig$values[keep], nrow = rank) %*% t(vectors)
+  }
+
+  contrasts <- .prepare_contrast_matrix(contrast_matrix, trt_levels, tolerance)
+  projection <- info_inverse %*% info
+  residual <- contrasts - contrasts %*% projection
+  estimable <- nrow(contrasts) == 0L ||
+    max(abs(residual)) <= tolerance * (1 + max(abs(contrasts)))
+
+  contrast_covariance <- matrix(numeric(0), nrow = 0L, ncol = 0L)
+  A_value <- Inf
+  D_value <- Inf
+  if (estimable) {
+    contrast_covariance <- contrasts %*% info_inverse %*% t(contrasts)
+    contrast_covariance <- (contrast_covariance + t(contrast_covariance)) / 2
+    if (nrow(contrasts) == 0L) {
+      A_value <- 0
+      D_value <- 0
+    } else {
+      contrast_eig <- eigen(contrast_covariance, symmetric = TRUE, only.values = TRUE)$values
+      contrast_scale <- max(abs(contrast_eig))
+      if (all(contrast_eig > contrast_scale * tolerance)) {
+        A_value <- sum(diag(contrast_covariance))
+        D_value <- sum(log(contrast_eig))
+      }
+    }
+  }
+
+  return(list(
+    eigenvalues = pos_eig,
+    rank = rank,
+    contrast_matrix = contrasts,
+    contrast_covariance = contrast_covariance,
+    estimable = estimable,
+    A_value = A_value,
+    D_value = D_value
+  ))
+}
+
+#' Function to precompute \eqn{L}
+#'
+#' \eqn{L} is the projection matrix that removes nuisance effects under a
+#' spatial correlation structure, given an \eqn{n \times n}{n * n} covariance
+#' \eqn{\Sigma}{Sigma} and nuisance effects design \eqn{X_2}.
+#'
+#' This only depends on lambda and the block structure, so precompute this.
+#' @noRd
+.build_nuisance_matrix <- function(layout_df, block_column, nuisance_formula) {
+  if (is.null(nuisance_formula)) {
+    if (!block_column %in% names(layout_df)) {
+      return(matrix(1, nrow(layout_df), 1L))
+    }
+    return(stats::model.matrix(~ 0 + block, data.frame(block = factor(layout_df[[block_column]]))))
+  }
+  if (!inherits(nuisance_formula, "formula") || length(nuisance_formula) != 2L) {
+    stop("`nuisance_formula` must be a one-sided formula.", call. = FALSE)
+  }
+  return(stats::model.matrix(nuisance_formula, layout_df))
+}
+
+#' Build a projection under identity covariance
+#' @noRd
+.build_L_from_df <- function(layout_df, block_column, n, nuisance_formula = NULL, tolerance = 1e-10) {
+  X2 <- .build_nuisance_matrix(layout_df, block_column, nuisance_formula)
+  crossprod_inverse <- pseudo_inverse(crossprod(X2), tolerance)
+  return(diag(n) - X2 %*% crossprod_inverse %*% t(X2))
+}
+
+#' Check whether a cached projection matches the requested projection
+#' @noRd
+.fisher_projection_matches <- function(
+  cache,
+  source,
+  input,
+  block_column,
+  nuisance_formula,
+  tolerance,
+  n
+) {
+  matches <- is.list(cache) &&
+    identical(cache$projection_source, source) &&
+    identical(cache$projection_input, input) &&
+    identical(cache$block_column, block_column) &&
+    identical(cache$nuisance_formula, nuisance_formula) &&
+    identical(cache$tolerance, tolerance) &&
+    is.matrix(cache$L_matrix) &&
+    identical(dim(cache$L_matrix), c(n, n))
+  return(isTRUE(matches))
+}
+
+#' Check whether an information-matrix cache can be updated
+#' @noRd
+.fisher_cache_is_valid <- function(cache, swap, trt_levels, L_matrix, n, v) {
+  valid <- is.list(cache) &&
+    identical(cache$swap, swap) &&
+    identical(cache$trt_levels, trt_levels) &&
+    length(cache$treatments) == n &&
+    all(cache$treatments %in% trt_levels) &&
+    is.matrix(cache$LX_matrix) &&
+    identical(dim(cache$LX_matrix), c(n, v)) &&
+    identical(cache$L_matrix, L_matrix)
+  return(isTRUE(valid))
+}
+
+#' Function for computing and caching the information matrix
+#' @noRd
+.compute_info <- function(
+  layout_df,
+  swap,
+  Sigma,
+  L_matrix,
+  block_column,
+  nuisance_formula = NULL,
+  tolerance = 1e-10,
+  cache = NULL
+) {
+  treatments <- layout_df[[swap]]
+  trt_levels <- levels(factor(treatments))
+  v <- length(trt_levels)
+  n <- nrow(layout_df)
+  treatment_labels <- as.character(treatments)
+
+  if (!is.null(Sigma) && !is.null(L_matrix)) {
+    stop("Only one of `Sigma` and `L_matrix` can be supplied.", call. = FALSE)
+  }
+
+  if (!is.null(L_matrix)) {
+    projection_source <- "L_matrix"
+    projection_input <- L_matrix
+  } else if (!is.null(Sigma)) {
+    projection_source <- "Sigma"
+    projection_input <- Sigma
+    if (.fisher_projection_matches(
+      cache,
+      projection_source,
+      projection_input,
+      block_column,
+      nuisance_formula,
+      tolerance,
+      n
+    )) {
+      L_matrix <- cache$L_matrix
+    } else {
+      L_matrix <- compute_L_projection(
+        layout_df,
+        Sigma,
+        block_column,
+        nuisance_formula,
+        tolerance
+      )
+    }
+  } else {
+    projection_source <- "identity"
+    projection_input <- NULL
+    if (.fisher_projection_matches(
+      cache,
+      projection_source,
+      projection_input,
+      block_column,
+      nuisance_formula,
+      tolerance,
+      n
+    )) {
+      L_matrix <- cache$L_matrix
+    } else {
+      L_matrix <- .build_L_from_df(
+        layout_df,
+        block_column,
+        n,
+        nuisance_formula,
+        tolerance
+      )
+    }
+  }
+
+  if (!is.matrix(L_matrix) || !identical(dim(L_matrix), c(n, n))) {
+    stop("`L_matrix` must be an n by n matrix for `layout_df`.", call. = FALSE)
+  }
+
+  X1 <- .build_treatment_matrix(treatment_labels, trt_levels, n, v)
+  if (.fisher_cache_is_valid(cache, swap, trt_levels, L_matrix, n, v)) {
+    changed <- which(cache$treatments != treatment_labels)
+    LX <- cache$LX_matrix
+    if (length(changed) > 0L) {
+      old_index <- match(cache$treatments[changed], trt_levels)
+      new_index <- match(treatment_labels[changed], trt_levels)
+      delta <- matrix(0, nrow = length(changed), ncol = v)
+      delta[cbind(seq_along(changed), old_index)] <- -1
+      delta[cbind(seq_along(changed), new_index)] <- 1
+      LX <- LX + L_matrix[, changed, drop = FALSE] %*% delta
+    }
+  } else {
+    LX <- L_matrix %*% X1
+  }
+
+  info <- crossprod(X1, LX)
+  info <- (info + t(info)) / 2
+  dimnames(info) <- list(trt_levels, trt_levels)
+
+  cache <- list(
+    swap = swap,
+    treatments = treatment_labels,
+    trt_levels = trt_levels,
+    L_matrix = L_matrix,
+    LX_matrix = LX,
+    projection_source = projection_source,
+    projection_input = projection_input,
+    block_column = block_column,
+    nuisance_formula = nuisance_formula,
+    tolerance = tolerance
+  )
+
+  return(list(
+    info = info,
+    v = v,
+    trt_levels = trt_levels,
+    cache = cache
+  ))
+}
+
+# Covariance utilities
+
+#' Compute L projection
+#'
+#' This is the projection matrix that removes nuisance effects from the GLS's
+#' covariance
+#'
+#' @param layout_df Data frame representing the spatial layout of the
+#'   experiment.
+#' @param Sigma Covariance structure to use. Its rows and columns must follow
+#'   the row order of `layout_df`.
+#' @param block_column Column name of the design's block factor in
+#'   \code{layout_df}.
+#' @param nuisance_formula Optional one-sided formula describing nuisance
+#'   effects. If omitted, \code{block_column} is treated as a factor.
+#' @param tolerance Numerical tolerance used for generalized inverses.
+#'
+#' @return \eqn{(n \times n)} numeric matrix.
+#'
+#' @export
+compute_L_projection <- function(
+  layout_df,
+  Sigma,
+  block_column = "block",
+  nuisance_formula = NULL,
+  tolerance = 1e-10
+) {
+  n <- nrow(layout_df)
+
+  if (!is.matrix(Sigma)) {
+    stop("Sigma must be a matrix")
+  }
+  if (!identical(dim(Sigma), c(n, n))) {
+    stop(
+      "Sigma matrix must have dimension equal ",
+      "to dimension of layout_df dataframe"
+    )
+  }
+  if (!isSymmetric(Sigma, tol = 1e-8)) stop("Sigma must be symmetric")
+  X2 <- .build_nuisance_matrix(layout_df, block_column, nuisance_formula)
+
+  Sigma_inv <- solve(Sigma)
+  SiX2 <- Sigma_inv %*% X2
+  nuisance_inverse <- pseudo_inverse(crossprod(X2, SiX2), tolerance)
+  L <- Sigma_inv - SiX2 %*% nuisance_inverse %*% t(SiX2)
+  return((L + t(L)) / 2)
+}
+
+
+# Convenience correlation structure constructors
+
+#' Construct a 1 dimensional lag-1 autoregressive covariance
+#'
+#' @param n Size of the covariance vector to generate
+#' @param rho Correlation parameter
+#'
+#' @return A symmetric \eqn{(n \times n)} correlation matrix.
+#'
+#' @export
+cor_ar1 <- function(n, rho) {
+  return(rho^abs(outer(seq_len(n), seq_len(n), "-")))
+}
+
+
+#' Construct a 2 dimensional lag-1 autoregressive covariance
+#'
+#' @param n_rows Number of rows in the design.
+#' @param n_cols Number of columns in the design.
+#' @param rho_row Correlation parameter along the row direction.
+#' @param rho_col Correlation parameter along the column direction.
+#'
+#' @return A symmetric \eqn{(\text{n_rows} \cdot \text{n_cols}) \times
+#'   (\text{n_rows} \cdot \text{n_cols})} correlation matrix.
+#'
+#' @export
+cor_ar1_ar1 <- function(n_rows, n_cols, rho_row, rho_col) {
+  return(kronecker(
+    cor_ar1(n_rows, rho_row),
+    cor_ar1(n_cols, rho_col)
+  ))
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 87ebe5e..18955f1 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -20,6 +20,13 @@ reference:
       - get_edges
       - get_vertices
       - prep_relationship
+  - subtitle: Fisher information utilities
+    desc: Functions for constructing contrast and covariance inputs for Fisher information objectives
+  - contents:
+      - factorial_contrast_matrix
+      - compute_L_projection
+      - cor_ar1
+      - cor_ar1_ar1
   - title: Output
     desc: Functions for checking designs.
   - contents:
diff --git a/man/compute_L_projection.Rd b/man/compute_L_projection.Rd
new file mode 100644
index 0000000..9923f9c
--- /dev/null
+++ b/man/compute_L_projection.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/objective_fisher.R
+\name{compute_L_projection}
+\alias{compute_L_projection}
+\title{Compute L projection}
+\usage{
+compute_L_projection(
+  layout_df,
+  Sigma,
+  block_column = "block",
+  nuisance_formula = NULL,
+  tolerance = 1e-10
+)
+}
+\arguments{
+\item{layout_df}{Data frame representing the spatial layout of the
+experiment.}
+
+\item{Sigma}{Covariance structure to use. Its rows and columns must follow
+the row order of \code{layout_df}.}
+
+\item{block_column}{Column name of the design's block factor in
+\code{layout_df}.}
+
+\item{nuisance_formula}{Optional one-sided formula describing nuisance
+effects. If omitted, \code{block_column} is treated as a factor.}
+
+\item{tolerance}{Numerical tolerance used for generalized inverses.}
+}
+\value{
+\eqn{(n \times n)} numeric matrix.
+}
+\description{
+This is the projection matrix that removes nuisance effects from the GLS's
+covariance
+}
diff --git a/man/cor_ar1.Rd b/man/cor_ar1.Rd
new file mode 100644
index 0000000..4dd478b
--- /dev/null
+++ b/man/cor_ar1.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/objective_fisher.R
+\name{cor_ar1}
+\alias{cor_ar1}
+\title{Construct a 1 dimensional lag-1 autoregressive covariance}
+\usage{
+cor_ar1(n, rho)
+}
+\arguments{
+\item{n}{Size of the covariance vector to generate}
+
+\item{rho}{Correlation parameter}
+}
+\value{
+A symmetric \eqn{(n \times n)} correlation matrix.
+}
+\description{
+Construct a 1 dimensional lag-1 autoregressive covariance
+}
diff --git a/man/cor_ar1_ar1.Rd b/man/cor_ar1_ar1.Rd
new file mode 100644
index 0000000..bbb4b2b
--- /dev/null
+++ b/man/cor_ar1_ar1.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/objective_fisher.R
+\name{cor_ar1_ar1}
+\alias{cor_ar1_ar1}
+\title{Construct a 2 dimensional lag-1 autoregressive covariance}
+\usage{
+cor_ar1_ar1(n_rows, n_cols, rho_row, rho_col)
+}
+\arguments{
+\item{n_rows}{Number of rows in the design.}
+
+\item{n_cols}{Number of columns in the design.}
+
+\item{rho_row}{Correlation parameter along the row direction.}
+
+\item{rho_col}{Correlation parameter along the column direction.}
+}
+\value{
+A symmetric \eqn{(\text{n_rows} \cdot \text{n_cols}) \times
+  (\text{n_rows} \cdot \text{n_cols})} correlation matrix.
+}
+\description{
+Construct a 2 dimensional lag-1 autoregressive covariance
+}
diff --git a/man/factorial_contrast_matrix.Rd b/man/factorial_contrast_matrix.Rd
new file mode 100644
index 0000000..b64c792
--- /dev/null
+++ b/man/factorial_contrast_matrix.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/objective_fisher.R
+\name{factorial_contrast_matrix}
+\alias{factorial_contrast_matrix}
+\title{Build an orthonormal factorial contrast matrix}
+\usage{
+factorial_contrast_matrix(
+  treatment_df,
+  formula,
+  treatment_column = "treatment",
+  tolerance = 1e-10
+)
+}
+\arguments{
+\item{treatment_df}{One row per treatment combination, containing the
+component factor columns.}
+
+\item{formula}{One-sided formula selecting factorial effects, for example
+\code{~ (stage + cultivar + inoculum)^2}.}
+
+\item{treatment_column}{Column containing the atomic treatment labels.}
+
+\item{tolerance}{Numerical rank tolerance.}
+}
+\value{
+A numeric contrast matrix whose columns are named by treatment level.
+}
+\description{
+Constructs a contrast basis for selected factorial effects while treatments
+remain encoded as atomic combinations during optimisation. The returned rows
+are orthonormal, so A-optimality is invariant to the particular contrast
+coding used by \code{\link[stats:model.matrix]{stats::model.matrix()}}.
+}
diff --git a/man/objective_function_info.Rd b/man/objective_function_info.Rd
new file mode 100644
index 0000000..db4acc0
--- /dev/null
+++ b/man/objective_function_info.Rd
@@ -0,0 +1,130 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/objective_fisher.R
+\name{objective_function_info}
+\alias{objective_function_info}
+\title{Objective function using the treatment information matrix}
+\usage{
+objective_function_info(
+  layout_df,
+  swap,
+  spatial_cols,
+  criterion = c("A", "D"),
+  L_matrix = NULL,
+  Sigma = NULL,
+  contrast_matrix = NULL,
+  block_column = "block",
+  nuisance_formula = NULL,
+  tolerance = 1e-09,
+  current_score_obj = NULL,
+  swapped_items = NULL,
+  ...
+)
+}
+\arguments{
+\item{layout_df}{Data frame representing the spatial layout of the
+experiment.}
+
+\item{swap}{Column name to swap, usually the treatment.}
+
+\item{spatial_cols}{Column name of the spatial factors.}
+
+\item{criterion}{Either \code{"A"} or \code{"D"}, representing A or D
+optimality.
+- A-optimality: Minimises \eqn{\mathrm{tr} \left( \mathcal{I}^-
+      \right)}{tr(I⁻)}.
+    - D-optimality: Minimises \eqn{-\log \left| \mathcal{I} \right|}{-
+      log(|I|)}}
+
+\item{L_matrix}{Optional precomputed projection matrix. This is an advanced
+alternative to \code{Sigma}; the two cannot be supplied together.}
+
+\item{Sigma}{Optional covariance matrix whose rows and columns follow the row
+order of \code{layout_df}. The projection is computed on the first objective
+call and reused during optimisation. If both \code{Sigma} and \code{L_matrix} are
+\code{NULL}, identity covariance is assumed.}
+
+\item{contrast_matrix}{Numeric contrast matrix with one column per treatment
+level. Column names, when supplied, must match the treatment levels. Only
+these contrasts are required to be estimable and contribute to the
+optimality criterion. The default targets the complete treatment-contrast
+space.}
+
+\item{block_column}{Column name of the design's block factor. Used when
+\code{L_matrix} isn't supplied.}
+
+\item{nuisance_formula}{Optional one-sided formula describing nuisance
+effects. Used when \code{L_matrix} is not supplied. If omitted,
+\code{block_column} is treated as a factor.}
+
+\item{tolerance}{Relative numerical tolerance for rank and estimability
+checks.}
+
+\item{current_score_obj}{The current objective result supplied by \code{\link[=speed]{speed()}}
+after the initial call. Used to reuse invariant matrix calculations.}
+
+\item{swapped_items}{The treatments swapped by \code{\link[=speed]{speed()}}. Changed plot rows
+are identified from the cached and proposed allocations so repeated
+treatments and multi-swaps are handled exactly.}
+
+\item{...}{Extra parameters passed from \code{speed}.}
+}
+\value{
+A named list with \code{score}, \code{info_matrix},
+\code{eigenvalues}, contrast diagnostics, and \code{criterion}.
+}
+\description{
+Creates an objective function that optimises experimental designs using the
+Fisher information for treatment contrasts after adjusting for nuisance
+effects, and takes an optional spatial covariance structure for spatial
+optimisation.
+}
+\details{
+This function computes the treatment information matrix:
+\deqn{I = X_1^\intercal L X_1}{I = X₁ᵀ L X₁}
+Where \eqn{X_1}{X₁} is the treatment design matrix, and \eqn{L} is the
+projection that removes nuisance fixed effects from the GLS's inverse
+covariance:
+\deqn{L = \Sigma^{-1} - \Sigma^{-1} X_2 \left( X_2^\intercal \Sigma^{-1} X_2
+\right)^{-1} X_2^\intercal \Sigma^{-1}}{L = Σ⁻¹ - Σ⁻¹ X₂
+(X₂ᵀ Σ⁻¹ X₂)⁻¹ X₂ᵀ Σ⁻¹}
+The user specifies the spatial covariance \eqn{\Sigma}{Σ} they intend to use
+at the time of analysis, e.g. lag-1 autoregressive \eqn{\text{AR}_1 \otimes
+\text{AR}_1}{AR₁ ⊗ AR₁}, and this objective will optimise under that
+structure.
+}
+\examples{
+# Small RCBD layout: 6 treatments in 4 blocks of 6 plots
+df <- initialise_design_df(
+  items = 6, nrows = 4, ncols = 6,
+  block_nrows = 1, block_ncols = 6
+)
+
+# Non-spatial: identity covariance
+result <- speed(
+  df,
+  swap = "treatment",
+  swap_within = "block",
+  obj_function = objective_function_info,
+  criterion = "A",
+  seed = 42,
+  quiet = TRUE
+)
+
+# Spatial: AR(1) x AR(1)
+Sigma <- cor_ar1_ar1(
+  n_rows = 4, n_cols = 6,
+  rho_row = 0.6, rho_col = 0.3
+)
+result_spatial <- speed(
+  df,
+  swap = "treatment",
+  swap_within = "block",
+  obj_function = objective_function_info,
+  Sigma = Sigma,
+  criterion = "A",
+  optimise_params = optim_params(random_initialisation = TRUE),
+  seed = 42,
+  quiet = TRUE
+)
+
+}
diff --git a/man/speed-package.Rd b/man/speed-package.Rd
index 8e018f8..0264306 100644
--- a/man/speed-package.Rd
+++ b/man/speed-package.Rd
@@ -23,12 +23,14 @@ Authors:
   \item Sam Rogers \email{sam.rogers@adelaide.edu.au}
   \item Wasin Pipattungsakul \email{wasin.pipattungsakul@adelaide.edu.au}
   \item Julian Taylor \email{julian.taylor@adelaide.edu.au}
+  \item Kai Bagley \email{kai.bagley@curtin.edu.au} (\href{https://orcid.org/0009-0004-6579-6959}{ORCID})
 }
 
 Other contributors:
 \itemize{
   \item University of Adelaide (https://adelaide.edu.au/) [copyright holder, funder]
   \item Grains Research and Development Corporation (https://grdc.com.au/) [copyright holder, funder]
+  \item Curtin University (\href{https://ror.org/02n415q13}{ROR}) (http://www.curtin.edu.au/) [copyright holder]
 }
 
 }
diff --git a/tests/testthat/test-compute_L_projection.R b/tests/testthat/test-compute_L_projection.R
new file mode 100644
index 0000000..ff3e58a
--- /dev/null
+++ b/tests/testthat/test-compute_L_projection.R
@@ -0,0 +1,48 @@
+test_that("compute_L_projection validates covariance inputs", {
+  df <- data.frame(
+    row = 1:4,
+    col = 1:4,
+    block = c(1, 1, 2, 2),
+    treatment = c("A", "B", "A", "B")
+  )
+
+  expect_error(compute_L_projection(df, 1:4), "matrix")
+  expect_error(compute_L_projection(df, diag(3)), "dimension")
+  expect_error(compute_L_projection(df, matrix(1:16, 4, 4)), "symmetric")
+})
+
+test_that("identity covariance matches the direct projection", {
+  df <- data.frame(
+    row       = rep(1:3, each = 4),
+    col       = rep(1:4, times = 3),
+    block     = rep(1:3, each = 4),
+    treatment = rep(paste0("T", 1:4), 3)
+  )
+
+  spatial <- compute_L_projection(df, diag(12), block_column = "block")
+  direct <- .build_L_from_df(df, "block", 12)
+
+  expect_equal(unname(spatial), unname(direct), tolerance = 1e-10)
+})
+
+test_that("compute_L_projection supports nuisance formulas", {
+  df <- data.frame(
+    row       = rep(1:3, each = 4),
+    col       = rep(1:4, times = 3),
+    treatment = rep(paste0("T", 1:4), 3)
+  )
+
+  spatial <- compute_L_projection(
+    df,
+    diag(12),
+    nuisance_formula = ~ row + col
+  )
+  direct <- .build_L_from_df(
+    df,
+    "block",
+    12,
+    nuisance_formula = ~ row + col
+  )
+
+  expect_equal(unname(spatial), unname(direct), tolerance = 1e-10)
+})
diff --git a/tests/testthat/test-cor_ar1.R b/tests/testthat/test-cor_ar1.R
new file mode 100644
index 0000000..96a8e0c
--- /dev/null
+++ b/tests/testthat/test-cor_ar1.R
@@ -0,0 +1,15 @@
+test_that("cor_ar1 produces a valid correlation matrix", {
+  correlation <- cor_ar1(5, 0.7)
+
+  expect_equal(dim(correlation), c(5, 5))
+  expect_true(isSymmetric(correlation))
+  expect_true(all(diag(correlation) == 1))
+  expect_true(all(eigen(correlation, only.values = TRUE)$values > 0))
+  expect_equal(correlation[1, 2], 0.7)
+  expect_equal(correlation[1, 3], 0.7^2)
+  expect_equal(correlation[1, 5], 0.7^4)
+})
+
+test_that("cor_ar1 with zero correlation gives identity", {
+  expect_equal(cor_ar1(4, 0), diag(4))
+})
diff --git a/tests/testthat/test-cor_ar1_ar1.R b/tests/testthat/test-cor_ar1_ar1.R
new file mode 100644
index 0000000..6079926
--- /dev/null
+++ b/tests/testthat/test-cor_ar1_ar1.R
@@ -0,0 +1,14 @@
+test_that("cor_ar1_ar1 produces a valid separable correlation matrix", {
+  correlation <- cor_ar1_ar1(3, 5, 0.6, 0.4)
+
+  expect_equal(dim(correlation), c(15, 15))
+  expect_true(isSymmetric(correlation))
+  expect_true(all(diag(correlation) == 1))
+  expect_true(all(eigen(correlation, only.values = TRUE)$values > 0))
+})
+
+test_that("cor_ar1_ar1 is the Kronecker product of marginal correlations", {
+  expected <- kronecker(cor_ar1(3, 0.6), cor_ar1(5, 0.4))
+
+  expect_equal(cor_ar1_ar1(3, 5, 0.6, 0.4), expected)
+})
diff --git a/tests/testthat/test-factorial_contrast_matrix.R b/tests/testthat/test-factorial_contrast_matrix.R
new file mode 100644
index 0000000..8f37241
--- /dev/null
+++ b/tests/testthat/test-factorial_contrast_matrix.R
@@ -0,0 +1,50 @@
+test_that("factorial_contrast_matrix builds an orthonormal effect space", {
+  treatments <- expand.grid(
+    stage = c("early", "late"),
+    cultivar = c("A", "B", "C")
+  )
+  treatments$treatment <- with(treatments, paste(stage, cultivar, sep = ":"))
+
+  contrasts <- factorial_contrast_matrix(
+    treatments,
+    ~ stage + cultivar
+  )
+
+  expect_equal(dim(contrasts), c(3, 6))
+  expect_equal(colnames(contrasts), treatments$treatment)
+  expect_equal(unname(rowSums(contrasts)), rep(0, 3), tolerance = 1e-12)
+  expect_equal(unname(tcrossprod(contrasts)), diag(3), tolerance = 1e-12)
+})
+
+test_that("factorial_contrast_matrix supports interactions", {
+  treatments <- expand.grid(
+    stage = c("early", "late"),
+    cultivar = c("A", "B", "C")
+  )
+  treatments$treatment <- with(treatments, paste(stage, cultivar, sep = ":"))
+
+  contrasts <- factorial_contrast_matrix(
+    treatments,
+    ~ stage * cultivar
+  )
+
+  expect_equal(nrow(contrasts), nrow(treatments) - 1)
+  expect_equal(unname(tcrossprod(contrasts)), diag(5), tolerance = 1e-12)
+})
+
+test_that("factorial_contrast_matrix validates its inputs", {
+  treatments <- data.frame(
+    treatment = c("A", "B"),
+    factor_a = c("low", "high")
+  )
+
+  expect_error(factorial_contrast_matrix(1:2, ~factor_a), "data frame")
+  expect_error(factorial_contrast_matrix(treatments, factor_a ~ 1), "one-sided")
+  expect_error(factorial_contrast_matrix(treatments, ~missing), "not found")
+
+  duplicated <- rbind(treatments, treatments[1, ])
+  expect_error(
+    factorial_contrast_matrix(duplicated, ~factor_a),
+    "one non-missing row"
+  )
+})
diff --git a/tests/testthat/test-objective_fisher.R b/tests/testthat/test-objective_fisher.R
new file mode 100644
index 0000000..6cbb22e
--- /dev/null
+++ b/tests/testthat/test-objective_fisher.R
@@ -0,0 +1,327 @@
+# tests/testthat/test-objective_fisher.R
+
+# ---- objective_function_info
+
+# all RCBDs are equally optimal
+test_that("all valid RCBD randomisations produce the same score under identity covariance", {
+  set.seed(1)
+  n_treatments <- 4
+  n_blocks <- 3
+
+  df <- data.frame(
+    row       = rep(1:n_blocks, each = n_treatments),
+    col       = rep(1:n_treatments, times = n_blocks),
+    block     = rep(1:n_blocks, each = n_treatments),
+    treatment = NA_character_
+  )
+
+  scores <- numeric(20)
+  for (i in seq_len(20)) {
+    for (b in 1:n_blocks) {
+      idx <- df$block == b
+      df$treatment[idx] <- paste0("T", sample(n_treatments))
+    }
+    result <- objective_function_info(
+      df,
+      swap = "treatment", spatial_cols = c("row", "col"),
+      criterion = "A"
+    )
+    scores[i] <- result$score
+  }
+
+  expect_true(max(scores) - min(scores) < 1e-10)
+})
+
+# when there IS spatial correlation, two different arrangements of the same
+# treatments must give different scores. df1 should be worse, since treatments
+# are always in the same column position, but in df2 they're "randomised".
+test_that("spatial designs differentiate between arrangements", {
+  n_treatments <- 4
+  n_blocks <- 3
+
+  df <- data.frame(
+    row   = rep(1:n_blocks, each = n_treatments),
+    col   = rep(1:n_treatments, times = n_blocks),
+    block = rep(1:n_blocks, each = n_treatments)
+  )
+
+  Sigma <- cor_ar1_ar1(n_blocks, n_treatments, rho_row = 0.7, rho_col = 0.5)
+
+  # Same treatment ordering in every block
+  df1 <- df
+  df1$treatment <- rep(paste0("T", 1:n_treatments), n_blocks)
+
+  # Alternated ordering
+  df2 <- df
+  df2$treatment <- c(
+    paste0("T", c(1, 2, 3, 4)),
+    paste0("T", c(3, 4, 1, 2)),
+    paste0("T", c(4, 3, 2, 1))
+  )
+
+  s1 <- objective_function_info(
+    df1, "treatment", c("row", "col"),
+    criterion = "A", Sigma = Sigma
+  )$score
+
+  s2 <- objective_function_info(
+    df2, "treatment", c("row", "col"),
+    criterion = "A", Sigma = Sigma
+  )$score
+
+  expect_false(abs(s1 - s2) < 1e-10)
+})
+
+# just make sure this works lel
+test_that("D-optimality criterion works", {
+  df <- data.frame(
+    row       = rep(1:3, each = 4),
+    col       = rep(1:4, times = 3),
+    block     = rep(1:3, each = 4),
+    treatment = rep(paste0("T", 1:4), 3)
+  )
+
+  result <- objective_function_info(
+    df, "treatment", c("row", "col"),
+    criterion = "D"
+  )
+
+  expect_true(is.finite(result$score))
+  expect_equal(result$criterion, "D")
+  expect_equal(length(result$eigenvalues), 3) # v - 1
+})
+
+# if treatments never exist in the same block together, you cant compare them,
+# so the design is fully broken
+test_that("disconnected designs have massive penalty", {
+  df <- data.frame(
+    row       = 1:6,
+    col       = c(1, 2, 1, 2, 1, 2),
+    block     = c(1, 1, 2, 2, 3, 3),
+    treatment = c("A", "A", "B", "B", "C", "C")
+  )
+
+  result <- objective_function_info(
+    df, "treatment", c("row", "col"),
+    criterion = "A"
+  )
+
+  expect_gte(result$score, 1e9)
+  expect_true(is.finite(result$score))
+  expect_equal(result$A_value, Inf)
+  expect_equal(result$D_value, Inf)
+})
+
+test_that("information matrix for an RCBD has the expected rank", {
+  df <- data.frame(
+    row       = rep(1:3, each = 5),
+    col       = rep(1:5, times = 3),
+    block     = rep(1:3, each = 5),
+    treatment = rep(paste0("T", 1:5), 3)
+  )
+
+  result <- objective_function_info(df, "treatment", c("row", "col"))
+
+  expect_equal(result$rank, 4)
+  expect_equal(result$eigenvalues, rep(3, 4), tolerance = 1e-10)
+  expect_equal(dim(result$info_matrix), c(5, 5))
+  expect_equal(rownames(result$info_matrix), paste0("T", 1:5))
+  expect_equal(colnames(result$info_matrix), paste0("T", 1:5))
+})
+
+test_that("a requested estimable contrast does not require full connectedness", {
+  df <- data.frame(
+    row       = 1:6,
+    col       = rep(1:2, 3),
+    block     = rep(1:3, each = 2),
+    treatment = c("A", "B", "A", "B", "C", "C")
+  )
+  contrast <- matrix(c(1, -1, 0), nrow = 1)
+  colnames(contrast) <- c("A", "B", "C")
+
+  full <- objective_function_info(df, "treatment", c("row", "col"))
+  targeted <- objective_function_info(
+    df,
+    "treatment",
+    c("row", "col"),
+    contrast_matrix = contrast
+  )
+
+  expect_false(full$estimable)
+  expect_true(targeted$estimable)
+  expect_true(is.finite(targeted$score))
+})
+
+test_that("incremental calculation matches full recalculation", {
+  df <- data.frame(
+    row       = rep(1:3, each = 4),
+    col       = rep(1:4, times = 3),
+    block     = rep(1:3, each = 4),
+    treatment = rep(LETTERS[1:4], 3)
+  )
+  Sigma <- cor_ar1_ar1(3, 4, rho_row = 0.6, rho_col = 0.3)
+  current <- objective_function_info(
+    df,
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma
+  )
+
+  proposal <- df
+  proposal$treatment[c(1, 2, 7, 8)] <- proposal$treatment[c(2, 1, 8, 7)]
+  incremental <- objective_function_info(
+    proposal,
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma,
+    current_score_obj = current,
+    swapped_items = c("A", "B", "C", "D")
+  )
+  full <- objective_function_info(
+    proposal,
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma
+  )
+
+  expect_equal(incremental$info_matrix, full$info_matrix, tolerance = 1e-12)
+  expect_equal(incremental$score, full$score, tolerance = 1e-12)
+
+  # Reuse the original score object as if the first proposal was rejected.
+  second_proposal <- df
+  second_proposal$treatment[c(3, 4)] <- second_proposal$treatment[c(4, 3)]
+  second_incremental <- objective_function_info(
+    second_proposal,
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma,
+    current_score_obj = current,
+    swapped_items = c("C", "D")
+  )
+  second_full <- objective_function_info(
+    second_proposal,
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma
+  )
+
+  expect_equal(
+    second_incremental$info_matrix,
+    second_full$info_matrix,
+    tolerance = 1e-12
+  )
+})
+
+test_that("covariance matrices remain aligned when layout rows are reordered", {
+  df <- data.frame(
+    row       = rep(1:3, each = 4),
+    col       = rep(1:4, times = 3),
+    block     = rep(1:3, each = 4),
+    treatment = rep(LETTERS[1:4], 3)
+  )
+  Sigma <- cor_ar1_ar1(3, 4, rho_row = 0.6, rho_col = 0.3)
+  original <- objective_function_info(
+    df,
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma
+  )
+
+  row_order <- c(5:12, 1:4)
+  reordered <- objective_function_info(
+    df[row_order, ],
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma[row_order, row_order]
+  )
+
+  expect_equal(reordered$info_matrix, original$info_matrix, tolerance = 1e-12)
+  expect_equal(reordered$score, original$score, tolerance = 1e-12)
+})
+
+test_that("Sigma and L_matrix are mutually exclusive", {
+  df <- data.frame(
+    row = 1:4,
+    col = 1:4,
+    treatment = rep(c("A", "B"), 2)
+  )
+
+  expect_error(
+    objective_function_info(
+      df,
+      "treatment",
+      c("row", "col"),
+      Sigma = diag(4),
+      L_matrix = diag(4)
+    ),
+    "Only one"
+  )
+})
+
+
+# ---- Test with speed function
+
+# Actually test by passing to speed
+test_that("speed() optimises an RCBD with the info objective under AR1xAR1", {
+  skip_if_not_installed("speed")
+
+  df <- initialise_design_df(
+    items = 6, nrows = 4, ncols = 6,
+    block_nrows = 1, block_ncols = 6
+  )
+
+  Sigma <- cor_ar1_ar1(4, 6, rho_row = 0.6, rho_col = 0.3)
+
+  result <- speed(
+    df,
+    swap = "treatment",
+    swap_within = "block",
+    spatial_factors = ~ block + row + col,
+    obj_function = objective_function_info,
+    Sigma = Sigma,
+    criterion = "A",
+    optimise_params = optim_params(random_initialisation = TRUE),
+    seed = 42,
+    quiet = TRUE
+  )
+
+  expect_s3_class(result, "design")
+  expect_lt(result$score, Inf)
+
+  info <- objective_function_info(
+    result$design_df,
+    "treatment",
+    c("row", "col"),
+    Sigma = Sigma
+  )
+  expect_equal(info$rank, 5)
+  expect_true(is.finite(info$A_value))
+})
+
+test_that("speed() optimises a BIBD with the info objective", {
+  skip_if_not_installed("speed")
+
+  df <- initialise_design_df(
+    items = 5, nrows = 3, ncols = 10,
+    block_nrows = 3, block_ncols = 1
+  )
+
+  result <- speed(
+    df,
+    swap = "treatment",
+    spatial_factors = ~ block + row + col,
+    obj_function = objective_function_info,
+    criterion = "A",
+    optimise_params = optim_params(random_initialisation = TRUE),
+    seed = 42,
+    quiet = TRUE
+  )
+
+  info <- objective_function_info(
+    result$design_df,
+    "treatment",
+    c("row", "col")
+  )
+  expect_equal(info$rank, 4)
+  expect_true(is.finite(info$A_value))
+})