From d91aeef0a8a09645cbfdee7fa8955e9a121a395c Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Fri, 8 May 2026 00:17:48 +0930 Subject: [PATCH 1/5] feat(relationship): add relationship matrix lookup for adjacency scoring Introduces an optional `relationship` argument to `calculate_adjacency_score()` and `speed()`, letting neighbour pairs contribute a graded similarity from a user-supplied matrix instead of a strict identity match. `speed()` validates and pre-flattens the matrix once before the SA loop via the new exported `prep_relationship()` helper. Co-Authored-By: Claude Opus 4.7 (1M context) --- NAMESPACE | 1 + R/calculate_adjacency_score.R | 106 +++++++++++++- R/metrics.R | 2 +- R/speed.R | 16 ++- man/adjacency_score_vec.Rd | 11 +- man/calculate_adjacency_score.Rd | 23 +++- man/prep_relationship.Rd | 27 ++++ .../testthat/test-calculate_adjacency_score.R | 130 ++++++++++++++++++ 8 files changed, 304 insertions(+), 12 deletions(-) create mode 100644 man/prep_relationship.Rd diff --git a/NAMESPACE b/NAMESPACE index 889d8fa..923ae6e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -20,6 +20,7 @@ export(objective_function_piepho) export(objective_function_signature) export(optim_params) export(plot_progress) +export(prep_relationship) export(speed) importFrom(farver,decode_colour) importFrom(ggplot2,aes) diff --git a/R/calculate_adjacency_score.R b/R/calculate_adjacency_score.R index f3a8c3c..b53e64e 100755 --- a/R/calculate_adjacency_score.R +++ b/R/calculate_adjacency_score.R @@ -42,6 +42,10 @@ shift_pad <- function(m, dx, dy, fill = NA) { nr <- nrow(m) nc <- ncol(m) out <- matrix(fill, nr, nc) + # entire source falls off the grid when the shift meets or exceeds the + # corresponding dimension; `seq.int(1, nr - dy)` would otherwise produce + # a descending sequence and out-of-bounds indices + if (abs(dx) >= nc || abs(dy) >= nr) return(out) r_src <- if (dy >= 0) seq.int(1, nr - dy) else seq.int(1 - dy, nr) c_src <- if (dx >= 0) seq.int(1, nc - dx) else seq.int(1 - dx, nc) @@ -49,6 +53,64 @@ shift_pad <- function(m, dx, dy, fill = NA) { return(out) } +#' Validate and Flatten a Relationship Matrix for Lookup +#' +#' @description +#' Validates `relationship` and packs it into a flat numeric vector plus +#' row/column level vectors that [adjacency_score_vec()] can index into in +#' constant time per cell. [speed()] calls this once before the SA loop so +#' the per-iteration cost stays in the lookup itself. Direct callers of +#' [calculate_adjacency_score()] must wrap their matrix with this function +#' before passing — the score functions consume the prepped form only. +#' +#' @param relationship A numeric matrix with rownames and colnames covering +#' every treatment value to be scored. +#' @param treatments Optional character vector of treatment values that +#' `relationship` must cover. Skipped when `NULL`. +#' +#' @return A list with `flat` (column-major flattened matrix), `row_levels`, +#' `col_levels`, and `n_row`. +#' +#' @export +prep_relationship <- function(relationship, treatments = NULL) { + if (!is.matrix(relationship) || !is.numeric(relationship)) { + stop("`relationship` must be a numeric matrix.", call. = FALSE) + } + + if (is.null(rownames(relationship)) || is.null(colnames(relationship))) { + stop( + "`relationship` must have rownames and colnames covering all ", + "treatments to be scored.", + call. = FALSE + ) + } + + if (!is.null(treatments)) { + treatments <- as.character(treatments) + treatments <- unique(treatments[!is.na(treatments)]) + missing <- unique(c( + setdiff(treatments, rownames(relationship)), + setdiff(treatments, colnames(relationship)) + )) + + if (length(missing)) { + stop( + "`relationship` is missing entries for treatments: ", + paste(missing, collapse = ", "), + ".", + call. = FALSE + ) + } + } + + return(list( + flat = as.vector(relationship), + row_levels = rownames(relationship), + col_levels = colnames(relationship), + n_row = nrow(relationship) + )) +} + #' Per-Cell Weighted Adjacency Score #' #' @description @@ -70,7 +132,8 @@ shift_pad <- function(m, dx, dy, fill = NA) { adjacency_score_vec <- function(design_matrix, dists = c(1, 2), weights = c(1, 2), - ring_type = c("manhattan", "chebyshev")) { + ring_type = c("manhattan", "chebyshev"), + relationship = NULL) { ring_type <- match.arg(ring_type) stopifnot(length(dists) == length(weights)) nr <- nrow(design_matrix) @@ -91,10 +154,19 @@ adjacency_score_vec <- function(design_matrix, # force 3D even when n_offsets == 1, where simplify2array returns a matrix dim(shifted_stack) <- c(nr, nc, n_offsets) - matches <- sweep(shifted_stack, c(1, 2), design_matrix, FUN = "==") - matches[is.na(matches)] <- FALSE + if (is.null(relationship)) { + pair_values <- sweep(shifted_stack, c(1, 2), design_matrix, FUN = "==") + pair_values[is.na(pair_values)] <- FALSE + } else { + orig_idx <- match(as.character(design_matrix), relationship$row_levels) + shifted_idx <- match(as.character(shifted_stack), relationship$col_levels) + flat_idx <- rep(orig_idx, n_offsets) + relationship$n_row * (shifted_idx - 1L) + pair_values <- relationship$flat[flat_idx] + pair_values[is.na(pair_values)] <- 0 + dim(pair_values) <- c(nr, nc, n_offsets) + } - weighted <- sweep(matches, 3, offset_weights, FUN = "*") + weighted <- sweep(pair_values, 3, offset_weights, FUN = "*") score <- rowSums(weighted, dims = 2) return(score) } @@ -113,6 +185,9 @@ adjacency_score_vec <- function(design_matrix, #' matches at larger ring radii. The per-cell score is summed and halved #' because each adjacency is counted once from each endpoint. #' +#' Pass a `relationship` matrix to score neighbour pairs by a graded +#' similarity (e.g. genetic relatedness) instead of a strict identity match. +#' #' @param layout_df A data frame containing the design. #' @param swap Column name of the treatments to be scored. #' @param row_column Name of the column representing rows (default `"row"`). @@ -124,6 +199,13 @@ adjacency_score_vec <- function(design_matrix, #' `1`). #' @param ring_type Ring shape: `"manhattan"` (default; diamond ring) #' or `"chebyshev"` (square ring). See [ring_offsets()]. +#' @param relationship Optional pairwise-relationship lookup produced by +#' [prep_relationship()]. When supplied, each neighbour pair contributes +#' `relationship[cell, neighbour]` rather than `1` for matches and `0` +#' otherwise. NA-padded cells off the design edge contribute `0`. Defaults +#' to `NULL`, which keeps the strict identity match. Pass the raw matrix +#' through `prep_relationship()` first; the score functions consume only +#' the prepped form. #' #' @return A non-negative numeric value: the number of like-treatment edges #' in the row/column adjacency graph. @@ -145,6 +227,15 @@ adjacency_score_vec <- function(design_matrix, #' ) #' calculate_adjacency_score(design_with_adj, "treatment") # 6 #' +#' # Example 3: graded relationship between A and B +#' rel <- prep_relationship(matrix( +#' c(1, 0.3, 0.3, 1), +#' nrow = 2, +#' dimnames = list(c("A", "B"), c("A", "B")) +#' )) +#' calculate_adjacency_score(design_no_adj, "treatment", relationship = rel) +#' # 3.6: each of the 12 A-B edges contributes 0.3 +#' #' @seealso [adjacency_score_vec()] #' #' @export @@ -154,8 +245,10 @@ calculate_adjacency_score <- function(layout_df, col_column = "col", ring_dists = 1, ring_weights = 1, - ring_type = c("manhattan", "chebyshev")) { + ring_type = c("manhattan", "chebyshev"), + relationship = NULL) { ring_type <- match.arg(ring_type) + design_matrix <- matrix( layout_df[[swap]], nrow = max(as_numeric_factor(layout_df[[row_column]]), na.rm = TRUE), @@ -167,7 +260,8 @@ calculate_adjacency_score <- function(layout_df, design_matrix, dists = ring_dists, weights = ring_weights, - ring_type = ring_type + ring_type = ring_type, + relationship = relationship ) return(sum(per_cell) / 2) } diff --git a/R/metrics.R b/R/metrics.R index f272a68..0479893 100644 --- a/R/metrics.R +++ b/R/metrics.R @@ -58,7 +58,7 @@ objective_function <- function(layout_df, ring_args <- list(...) ring_args <- ring_args[intersect( names(ring_args), - c("ring_dists", "ring_weights", "ring_type") + c("ring_dists", "ring_weights", "ring_type", "relationship") )] adj_score <- ifelse(adj_weight != 0, do.call( diff --git a/R/speed.R b/R/speed.R index b4f81b7..0bd9ea9 100644 --- a/R/speed.R +++ b/R/speed.R @@ -208,8 +208,18 @@ speed <- function(data, } } - design <- speed_hierarchical(data, optimise, quiet, seed, row_column = row_column, - col_column = col_column, ...) + dots <- list(...) + if (!is.null(dots$relationship)) { + swap_cols <- unique(vapply(optimise, function(o) o$swap, character(1))) + treatments <- unlist(lapply(swap_cols, function(s) as.character(data[[s]]))) + dots$relationship <- prep_relationship(dots$relationship, treatments) + } + + design <- do.call(speed_hierarchical, c( + list(data = data, optimise = optimise, quiet = quiet, seed = seed, + row_column = row_column, col_column = col_column), + dots + )) design$design_df[[dummy_group]] <- NULL design$design_df <- to_types(design$design_df, factored$input_types) @@ -288,7 +298,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, ...) new_score <- new_score_obj$score # Decide whether to accept the new design diff --git a/man/adjacency_score_vec.Rd b/man/adjacency_score_vec.Rd index 8631e31..7cb6627 100644 --- a/man/adjacency_score_vec.Rd +++ b/man/adjacency_score_vec.Rd @@ -8,7 +8,8 @@ adjacency_score_vec( design_matrix, dists = c(1, 2), weights = c(1, 2), - ring_type = c("manhattan", "chebyshev") + ring_type = c("manhattan", "chebyshev"), + relationship = NULL ) } \arguments{ @@ -20,6 +21,14 @@ adjacency_score_vec( \item{ring_type}{Ring shape: \code{"manhattan"} (default; diamond ring) or \code{"chebyshev"} (square ring). See \code{\link[=ring_offsets]{ring_offsets()}}.} + +\item{relationship}{Optional pairwise-relationship lookup produced by +\code{\link[=prep_relationship]{prep_relationship()}}. When supplied, each neighbour pair contributes +\code{relationship[cell, neighbour]} rather than \code{1} for matches and \code{0} +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.} } \value{ Matrix of the same dimensions as \code{design_matrix}, where diff --git a/man/calculate_adjacency_score.Rd b/man/calculate_adjacency_score.Rd index cebbb1e..391fc6a 100644 --- a/man/calculate_adjacency_score.Rd +++ b/man/calculate_adjacency_score.Rd @@ -11,7 +11,8 @@ calculate_adjacency_score( col_column = "col", ring_dists = 1, ring_weights = 1, - ring_type = c("manhattan", "chebyshev") + ring_type = c("manhattan", "chebyshev"), + relationship = NULL ) } \arguments{ @@ -32,6 +33,14 @@ calculate_adjacency_score( \item{ring_type}{Ring shape: \code{"manhattan"} (default; diamond ring) or \code{"chebyshev"} (square ring). See \code{\link[=ring_offsets]{ring_offsets()}}.} + +\item{relationship}{Optional pairwise-relationship lookup produced by +\code{\link[=prep_relationship]{prep_relationship()}}. When supplied, each neighbour pair contributes +\code{relationship[cell, neighbour]} rather than \code{1} for matches and \code{0} +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.} } \value{ A non-negative numeric value: the number of like-treatment edges @@ -48,6 +57,9 @@ neighbourhood the simulated-annealing loop in \code{\link[=speed]{speed()}} mini Pass longer \code{ring_dists} / \code{ring_weights} to also penalise like-treatment matches at larger ring radii. The per-cell score is summed and halved because each adjacency is counted once from each endpoint. + +Pass a \code{relationship} matrix to score neighbour pairs by a graded +similarity (e.g. genetic relatedness) instead of a strict identity match. } \examples{ # Example 1: design with no like-treatment adjacencies @@ -66,6 +78,15 @@ design_with_adj <- data.frame( ) calculate_adjacency_score(design_with_adj, "treatment") # 6 +# Example 3: graded relationship between A and B +rel <- prep_relationship(matrix( + c(1, 0.3, 0.3, 1), + nrow = 2, + dimnames = list(c("A", "B"), c("A", "B")) +)) +calculate_adjacency_score(design_no_adj, "treatment", relationship = rel) +# 3.6: each of the 12 A-B edges contributes 0.3 + } \seealso{ \code{\link[=adjacency_score_vec]{adjacency_score_vec()}} diff --git a/man/prep_relationship.Rd b/man/prep_relationship.Rd new file mode 100644 index 0000000..e90df54 --- /dev/null +++ b/man/prep_relationship.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calculate_adjacency_score.R +\name{prep_relationship} +\alias{prep_relationship} +\title{Validate and Flatten a Relationship Matrix for Lookup} +\usage{ +prep_relationship(relationship, treatments = NULL) +} +\arguments{ +\item{relationship}{A numeric matrix with rownames and colnames covering +every treatment value to be scored.} + +\item{treatments}{Optional character vector of treatment values that +\code{relationship} must cover. Skipped when \code{NULL}.} +} +\value{ +A list with \code{flat} (column-major flattened matrix), \code{row_levels}, +\code{col_levels}, and \code{n_row}. +} +\description{ +Validates \code{relationship} and packs it into a flat numeric vector plus +row/column level vectors that \code{\link[=adjacency_score_vec]{adjacency_score_vec()}} can index into in +constant time per cell. \code{\link[=speed]{speed()}} calls this once before the SA loop so +the per-iteration cost stays in the lookup itself. Direct callers of +\code{\link[=calculate_adjacency_score]{calculate_adjacency_score()}} must wrap their matrix with this function +before passing — the score functions consume the prepped form only. +} diff --git a/tests/testthat/test-calculate_adjacency_score.R b/tests/testthat/test-calculate_adjacency_score.R index 23616ef..bfe143f 100644 --- a/tests/testthat/test-calculate_adjacency_score.R +++ b/tests/testthat/test-calculate_adjacency_score.R @@ -64,6 +64,14 @@ test_that("shift_pad shifts and NA-pads correctly", { expect_equal(filled[1, ], rep(0L, 3)) }) +test_that("shift_pad returns all-fill when shift exceeds matrix bounds", { + m <- matrix(1:2, nrow = 1, ncol = 2) + expect_true(all(is.na(shift_pad(m, dx = 0, dy = 1)))) + expect_true(all(is.na(shift_pad(m, dx = 0, dy = -1)))) + expect_true(all(is.na(shift_pad(m, dx = 2, dy = 0)))) + expect_true(all(is.na(shift_pad(m, dx = -2, dy = 0)))) +}) + test_that("adjacency_score_vec scores per-cell matches", { # 1 1 # 2 2 @@ -122,3 +130,125 @@ test_that("adjacency_score_vec rejects mismatched dists/weights", { ) ) }) + +test_that("calculate_adjacency_score uses relationship for graded similarity", { + # A B A + # B A B + # A B A -> 12 A-B edges, no like-treatment adjacencies + design <- data.frame( + row = rep(1:3, each = 3), + col = rep(1:3, times = 3), + swap = c("A", "B", "A", "B", "A", "B", "A", "B", "A") + ) + expect_equal(calculate_adjacency_score(design, "swap"), 0) + + rel <- prep_relationship(matrix( + c(1.0, 0.3, 0.3, 1.0), + nrow = 2, + dimnames = list(c("A", "B"), c("A", "B")) + )) + expect_equal( + calculate_adjacency_score(design, "swap", relationship = rel), + 12 * 0.3 + ) +}) + +test_that("calculate_adjacency_score with identity relationship matches default", { + # A A A + # B B B + # A A A -> 6 like-treatment edges, 6 mixed edges + design <- data.frame( + row = rep(1:3, each = 3), + col = rep(1:3, times = 3), + swap = c("A", "A", "A", "B", "B", "B", "A", "A", "A") + ) + rel_mat <- diag(2) + dimnames(rel_mat) <- list(c("A", "B"), c("A", "B")) + + expect_equal( + calculate_adjacency_score(design, "swap", relationship = prep_relationship(rel_mat)), + calculate_adjacency_score(design, "swap") + ) +}) + +test_that("prep_relationship validates input", { + expect_error(prep_relationship(list()), "numeric matrix") + + expect_error( + prep_relationship(matrix("x", 2, 2, dimnames = list(c("A", "B"), c("A", "B")))), + "numeric matrix" + ) + + expect_error( + prep_relationship(matrix(0, nrow = 2, ncol = 2)), + "rownames and colnames" + ) + + expect_error( + prep_relationship( + matrix(0, 1, 1, dimnames = list("A", "A")), + treatments = c("A", "B") + ), + "missing entries for treatments" + ) +}) + +test_that("speed preps a raw relationship matrix once before the SA loop", { + df <- data.frame( + row = rep(1:3, times = 3), + col = rep(1:3, each = 3), + treatment = rep(c("A", "B", "C"), times = 3) + ) + rel <- matrix( + c(1, 0.5, 0.2, 0.5, 1, 0.5, 0.2, 0.5, 1), + nrow = 3, + dimnames = list(c("A", "B", "C"), c("A", "B", "C")) + ) + + # capture what reaches the objective function — must be the prepped form, + # not the raw matrix the user passed + seen <- new.env(parent = emptyenv()) + spy_obj <- function(layout_df, swap, spatial_cols, ...) { + dots <- list(...) + if (!is.null(dots$relationship)) { + seen$relationship <- dots$relationship + } + objective_function(layout_df, swap, spatial_cols, ...) + } + + result <- speed( + df, + swap = "treatment", + iterations = 5, + early_stop_iterations = 5, + obj_function = spy_obj, + relationship = rel, + quiet = TRUE, + seed = 1 + ) + + expect_named(seen$relationship, c("flat", "row_levels", "col_levels", "n_row")) + expect_s3_class(result, "design") +}) + +test_that("adjacency_score_vec respects relationship asymmetry per cell", { + # A B + # B A + m <- matrix(c("A", "B", "B", "A"), nrow = 2, byrow = TRUE) + rel <- prep_relationship(matrix( + c(0, 0.1, 0.7, 0), + nrow = 2, + dimnames = list(c("A", "B"), c("A", "B")) + )) + # rel["A","B"] = 0.7, rel["B","A"] = 0.1 (column-major fill) + # A cells see 2 B neighbours -> 2 * 0.7 = 1.4 + # B cells see 2 A neighbours -> 2 * 0.1 = 0.2 + score <- adjacency_score_vec( + m, + dists = 1, + weights = 1, + ring_type = "manhattan", + relationship = rel + ) + expect_equal(score, matrix(c(1.4, 0.2, 0.2, 1.4), nrow = 2, byrow = TRUE)) +}) From e4416e90dc747ced3cea9ffe817c18bca2e73106 Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Fri, 8 May 2026 00:18:02 +0930 Subject: [PATCH 2/5] docs(vignettes): add relationship-aware adjacency vignette Walks through a 16-treatment, 4-group blocked design and contrasts the default identity-only optimisation against `speed()` driven by a relationship matrix. Co-Authored-By: Claude Opus 4.7 (1M context) --- vignettes/genetic_relationship.qmd | 275 +++++++++++++++++++++++++++++ 1 file changed, 275 insertions(+) create mode 100644 vignettes/genetic_relationship.qmd diff --git a/vignettes/genetic_relationship.qmd b/vignettes/genetic_relationship.qmd new file mode 100644 index 0000000..e20a4ba --- /dev/null +++ b/vignettes/genetic_relationship.qmd @@ -0,0 +1,275 @@ +--- +title: "Relationship-Aware Adjacency in speed" +author: "Wasin Pipattungsakul" +date: "`r Sys.Date()`" +output: + html_document: + toc: true + toc_depth: 3 + toc_float: true + theme: flatly + highlight: tango +vignette: > + %\VignetteIndexEntry{Relationship-Aware Adjacency in speed} + %\VignetteEngine{quarto::html} + %\VignetteEncoding{UTF-8} +--- + +# Design with Relationship Matrix + +## Overview + +The default `speed` adjacency score asks *"are these two cells the same +treatment?"*. Pairs that match contribute `1`; pairs that differ +contribute `0`. Sometimes that's too coarse: two distinct treatments may +still be "close" to each other, and placing them next to each other +carries some of the cost of a true repeat. + +The `relationship` argument lets the simulated-annealing loop look up a +graded similarity from a user-supplied matrix instead of hard-coding `1` +for matches and `0` for non-matches. This vignette walks through a +blocked design with sixteen treatments split into four groups whose +within-group relatedness varies from "tightly related" to "unrelated". + +```{r load-package} +library(speed) +library(patchwork) +``` + +## When to Use + +- Trials where treatments fall into clusters of close relatives — for + example plant-breeding lines that share parents (encoded via a pedigree + A-matrix or marker-based G-matrix), variety entries from the same + maturity group, or scaled mixtures of a base formulation +- Anywhere "minimise neighbours of the same kind" is too coarse and you + would rather "minimise neighbours of related kind" + +# Example + +## Setting Up the Design + +Sixteen treatments `T01`–`T16` are organised into four groups of four: + +| Group | Treatments | Within-group similarity | +|------:|:-----------|------------------------:| +| `A` | `T01–T04` | `0.8` (tightly related) | +| `B` | `T05–T08` | `0.4` (moderate) | +| `C` | `T09–T12` | `0.1` (loose) | +| `D` | `T13–T16` | `0.0` (unrelated) | + +The trial uses an `8 × 8` grid divided into four `4 × 4` blocks. + +```{r design-df} +items <- paste0("T", 1:16) +groups <- setNames(rep(c("A", "B", "C", "D"), each = 4), items) + +design_df <- initialise_design_df( + items = rep(items, 4), + nrows = 8, + ncols = 8, + block_nrows = 4, + block_ncols = 4 +) +design_df$group <- groups[design_df$treatment] + +head(design_df) +``` + +Plotting these factors shows the initial layout of the treatments and their group. + +```{r} +#| code-fold: true +#| code-summary: "Code" +#| label: fig-init +#| fig-cap: "Initial systematic layout: treatments and groups, with block boundaries." +class(design_df) <- c("design", class(design_df)) +p_treat <- autoplot(design_df, treatments = "treatment") + + ggplot2::labs(title = "Treatments") +p_group <- autoplot(design_df, treatments = "group") + + ggplot2::labs(title = "Groups") +p_treat + p_group + plot_layout(ncol = 2) +``` + +## Building the Relationship Matrix + +The relationship matrix is a numeric matrix with rownames and colnames +covering every value of the swap column. Each entry is the contribution +of having that pair adjacent. The diagonal is typically `1` +(self-pair = identity). Off-diagonal entries are usually in `[0, 1]`; +here we encode four within-group levels and zero between-group +similarity. + +```{r relationship} +n <- length(items) +rel <- matrix(0, n, n, dimnames = list(items, items)) +within <- c(A = 0.8, B = 0.4, C = 0.1, D = 0.0) +for (g in names(within)) { + ids <- items[groups == g] + rel[ids, ids] <- within[[g]] +} +diag(rel) <- 1 +``` + +```{r} +#| code-fold: true +#| code-summary: "Code" +#| label: fig-relationship +#| fig-cap: "Relationship matrix as a heatmap. Diagonal cells are self-pairs (`1`); each `4 × 4` block on the diagonal corresponds to a group, with brightness encoding within-group similarity." +#| fig-width: 6 +#| fig-height: 5 +rel_long <- as.data.frame.table(rel, responseName = "value") +names(rel_long)[1:2] <- c("row", "col") + +ggplot2::ggplot(rel_long, ggplot2::aes(x = col, y = row, fill = value)) + + ggplot2::geom_tile(colour = "white") + + ggplot2::geom_text( + ggplot2::aes(label = sprintf("%.1f", value)), + size = 2.5 + ) + + ggplot2::scale_x_discrete(position = "top") + + ggplot2::scale_y_discrete(limits = rev(rownames(rel))) + + ggplot2::scale_fill_gradient( + low = "white", high = "steelblue4", limits = c(0, 1), + name = "similarity" + ) + + ggplot2::coord_equal() + + ggplot2::labs(x = NULL, y = NULL) + + ggplot2::theme_minimal(base_size = 10) + + ggplot2::theme(panel.grid = ggplot2::element_blank()) +``` + +## Running the Optimisation + +We run `speed()` twice with the same seed: first using the default +identity-only adjacency, then with the relationship matrix (`relationship = rel`). + +```{r both-runs} +default_result <- speed( + design_df, + swap = "treatment", + swap_within = "block", + seed = 42, + quiet = TRUE +) + +related_result <- speed( + design_df, + swap = "treatment", + swap_within = "block", + relationship = rel, + seed = 42, + quiet = TRUE +) + +default_result +related_result +``` + +## Visualising the Output + +```{r} +#| code-fold: true +#| code-summary: "Code" +#| label: fig-final-group +#| fig-cap: "Group view: default identity-only adjacency vs. relationship-aware adjacency." +#| fig-width: 10 +default_result$design_df$group <- groups[default_result$design_df$treatment] +related_result$design_df$group <- groups[related_result$design_df$treatment] + +p_def_g <- autoplot(default_result, treatments = "group") + + ggplot2::labs(title = "Default speed: groups") +p_rel_g <- autoplot(related_result, treatments = "group") + + ggplot2::labs(title = "speed + relationship: groups") +p_def_g + p_rel_g + plot_layout(ncol = 2) +``` + +```{r} +#| code-fold: true +#| code-summary: "Code" +#| label: fig-final-treat +#| fig-cap: "Treatment view: default identity-only adjacency vs. relationship-aware adjacency." +#| fig-width: 10 +p_def_t <- autoplot(default_result, treatments = "treatment") + + ggplot2::labs(title = "Default speed: treatments") +p_rel_t <- autoplot(related_result, treatments = "treatment") + + ggplot2::labs(title = "speed + relationship: treatments") +p_def_t + p_rel_t + plot_layout(ncol = 2) +``` + +The default optimisation only avoids same-treatment adjacency; related +treatments from the same group are still free to sit next to each other. +The relationship-aware optimisation pushes whole groups apart, with the +strongest pressure on group `A` (within-similarity `0.8`) and no +pressure at all on group `D` (`0.0`). + + + +# Related Vignettes + +- [Common Agricultural Experimental Designs with speed](speed.html) — + foundational designs and the default optimisation flow. +- [Custom Objective Functions and Advanced Options for Optimisation in speed](custom_objective_functions.html) — + for cases where weighted-relationship adjacency is one of several + competing objectives. +- [Multi-Environment Trials with speed](met.html) — how the same workflow + extends to hierarchical layouts. + +--- + +*This vignette demonstrates the `relationship` argument added to +`calculate_adjacency_score()` and `speed()`. For other features of the +package, consult the documentation and additional vignettes.* From 750d1fb33cfcfdcfd75ad9f2a7ec526eff534847 Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Fri, 8 May 2026 00:18:13 +0930 Subject: [PATCH 3/5] chore: ignore scratch and local-config files .Rbuildignore: exclude .claude/ and Rplots.pdf from package builds. .gitignore: exclude *.local.json. Co-Authored-By: Claude Opus 4.7 (1M context) --- .Rbuildignore | 2 ++ .gitignore | 1 + 2 files changed, 3 insertions(+) diff --git a/.Rbuildignore b/.Rbuildignore index 3565385..489ce00 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -22,3 +22,5 @@ codemeta.json ^CODE_OF_CONDUCT\.md$ ^doc$ CLAUDE.md +.claude +Rplots.pdf diff --git a/.gitignore b/.gitignore index ca11726..6a2325c 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ inst/doc /Meta/ Rplots.pdf +*.local.json From 889f052d2b19597f52ad45f25a0a902a0fcca4d9 Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Fri, 8 May 2026 00:40:18 +0930 Subject: [PATCH 4/5] chore: add prep relationship --- _pkgdown.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index f2a3ae6..bd777ab 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -30,6 +30,7 @@ reference: - create_pair_mapping - initialise_design_df - add_buffers + - prep_relationship - title: Advanced options - contents: - optim_params From b37316250947e04ada413dff96c3e15a8d7604d1 Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Sat, 9 May 2026 20:15:01 +0930 Subject: [PATCH 5/5] fix(speed): reject adj_weight/bal_weight passed via `...` speed() forwards these explicitly from optim_params() while also splatting `...` to the objective function, so a direct `speed(..., bal_weight = 2)` triggered "matched by multiple actual arguments". `speed()` now stops with a message pointing the caller at `optimise_params = optim_params(...)`. Factor the dots guard and the pre-loop relationship prep into `.reject_optim_params_in_dots()` and `.prep_dots()` helpers. Vignette `genetic_relationship.qmd` gains a "Per-Pair Similarities" example with closely-related-pair highlighting. Co-Authored-By: Claude Opus 4.7 (1M context) --- R/calculate_adjacency_score.R | 14 +- R/speed.R | 43 +++++- man/adjacency_score_vec.Rd | 2 +- man/calculate_adjacency_score.Rd | 2 +- man/dot-prep_dots.Rd | 24 +++ man/dot-reject_optim_params_in_dots.Rd | 19 +++ man/prep_relationship.Rd | 2 +- man/ring_offsets.Rd | 4 +- tests/testthat/test-speed.R | 23 +++ vignettes/genetic_relationship.qmd | 196 ++++++++++++++++++++++++- 10 files changed, 305 insertions(+), 24 deletions(-) create mode 100644 man/dot-prep_dots.Rd create mode 100644 man/dot-reject_optim_params_in_dots.Rd diff --git a/R/calculate_adjacency_score.R b/R/calculate_adjacency_score.R index b53e64e..5062254 100755 --- a/R/calculate_adjacency_score.R +++ b/R/calculate_adjacency_score.R @@ -7,8 +7,8 @@ #' @inheritParams calculate_adjacency_score #' @param d Positive integer ring radius. #' -#' @return Integer matrix with two columns — column 1 holds `dx`, column 2 -#' holds `dy` — one row per cell on the ring. +#' @return Integer matrix with two columns - column 1 holds `dx`, column 2 +#' holds `dy` - one row per cell on the ring. #' #' @keywords internal ring_offsets <- function(d, ring_type = c("manhattan", "chebyshev")) { @@ -45,7 +45,9 @@ shift_pad <- function(m, dx, dy, fill = NA) { # entire source falls off the grid when the shift meets or exceeds the # corresponding dimension; `seq.int(1, nr - dy)` would otherwise produce # a descending sequence and out-of-bounds indices - if (abs(dx) >= nc || abs(dy) >= nr) return(out) + if (abs(dx) >= nc || abs(dy) >= nr) { + return(out) + } r_src <- if (dy >= 0) seq.int(1, nr - dy) else seq.int(1 - dy, nr) c_src <- if (dx >= 0) seq.int(1, nc - dx) else seq.int(1 - dx, nc) @@ -61,7 +63,7 @@ shift_pad <- function(m, dx, dy, fill = NA) { #' constant time per cell. [speed()] calls this once before the SA loop so #' the per-iteration cost stays in the lookup itself. Direct callers of #' [calculate_adjacency_score()] must wrap their matrix with this function -#' before passing — the score functions consume the prepped form only. +#' before passing - the score functions consume the prepped form only. #' #' @param relationship A numeric matrix with rownames and colnames covering #' every treatment value to be scored. @@ -117,7 +119,7 @@ prep_relationship <- function(relationship, treatments = NULL) { #' For each cell of `design_matrix`, counts neighbours on rings of radius #' `dists` whose value equals the cell's own, weighted by `weights`. #' Implemented by stacking shifted copies of the matrix once per offset and -#' reducing across the stack — avoids a per-offset loop when scoring many +#' reducing across the stack - avoids a per-offset loop when scoring many #' rings. #' #' @inheritParams calculate_adjacency_score @@ -174,7 +176,7 @@ adjacency_score_vec <- function(design_matrix, #' Calculate Adjacency Score for Design #' #' @description -#' Counts adjacent plots — immediate horizontal and vertical neighbours — +#' Counts adjacent plots - immediate horizontal and vertical neighbours - #' that share the same treatment. Lower scores indicate better separation. #' #' Internally this is a thin wrapper around [adjacency_score_vec()]. With diff --git a/R/speed.R b/R/speed.R index 0bd9ea9..4f6da78 100644 --- a/R/speed.R +++ b/R/speed.R @@ -209,11 +209,8 @@ speed <- function(data, } dots <- list(...) - if (!is.null(dots$relationship)) { - swap_cols <- unique(vapply(optimise, function(o) o$swap, character(1))) - treatments <- unlist(lapply(swap_cols, function(s) as.character(data[[s]]))) - dots$relationship <- prep_relationship(dots$relationship, treatments) - } + .reject_optim_params_in_dots(dots) + dots <- .prep_dots(dots, optimise, data) design <- do.call(speed_hierarchical, c( list(data = data, optimise = optimise, quiet = quiet, seed = seed, @@ -434,3 +431,39 @@ print.design <- function(x, ...) { return(invisible(x)) } + +#' Reject `...` arguments that must travel through [optim_params()]. +#' +#' Stops with a message pointing the user at `optimise_params = optim_params(...)` +#' when any of the listed names is found in `dots`. +#' +#' @param dots A named list captured from `...`. +#' @return `NULL`, invisibly. Called for its side effect. +#' @keywords internal +.reject_optim_params_in_dots <- function(dots) { + forbidden <- intersect(names(dots), c("adj_weight", "bal_weight")) + if (length(forbidden) == 0) return(invisible(NULL)) + stop( + "Argument(s) ", paste(sprintf("`%s`", forbidden), collapse = ", "), + " must be passed via `optim_params()`, not directly to `speed()`. ", + "For example: `optimise_params = optim_params(", + paste0(forbidden[1], " = ..."), ")`.", + call. = FALSE + ) +} + +#' Prep `dots$relationship` once with the union of treatments seen at +#' every swap level. +#' +#' @param dots A named list captured from `...`. +#' @param optimise Per-level `optimise` list as built by [create_speed_input()]. +#' @param data The (factor-converted) design data frame. +#' @return `dots`, with `relationship` replaced by the prepped form when present. +#' @keywords internal +.prep_dots <- function(dots, optimise, data) { + if (is.null(dots$relationship)) return(dots) + swap_cols <- unique(vapply(optimise, function(o) o$swap, character(1))) + treatments <- unlist(lapply(swap_cols, function(s) as.character(data[[s]]))) + dots$relationship <- prep_relationship(dots$relationship, treatments) + dots +} diff --git a/man/adjacency_score_vec.Rd b/man/adjacency_score_vec.Rd index 7cb6627..979d4aa 100644 --- a/man/adjacency_score_vec.Rd +++ b/man/adjacency_score_vec.Rd @@ -38,7 +38,7 @@ each cell holds its weighted matching-neighbour count. For each cell of \code{design_matrix}, counts neighbours on rings of radius \code{dists} whose value equals the cell's own, weighted by \code{weights}. Implemented by stacking shifted copies of the matrix once per offset and -reducing across the stack — avoids a per-offset loop when scoring many +reducing across the stack - avoids a per-offset loop when scoring many rings. } \keyword{internal} diff --git a/man/calculate_adjacency_score.Rd b/man/calculate_adjacency_score.Rd index 391fc6a..4cadd43 100644 --- a/man/calculate_adjacency_score.Rd +++ b/man/calculate_adjacency_score.Rd @@ -47,7 +47,7 @@ A non-negative numeric value: the number of like-treatment edges in the row/column adjacency graph. } \description{ -Counts adjacent plots — immediate horizontal and vertical neighbours — +Counts adjacent plots - immediate horizontal and vertical neighbours - that share the same treatment. Lower scores indicate better separation. Internally this is a thin wrapper around \code{\link[=adjacency_score_vec]{adjacency_score_vec()}}. With diff --git a/man/dot-prep_dots.Rd b/man/dot-prep_dots.Rd new file mode 100644 index 0000000..dbc9eed --- /dev/null +++ b/man/dot-prep_dots.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/speed.R +\name{.prep_dots} +\alias{.prep_dots} +\title{Prep \code{dots$relationship} once with the union of treatments seen at +every swap level.} +\usage{ +.prep_dots(dots, optimise, data) +} +\arguments{ +\item{dots}{A named list captured from \code{...}.} + +\item{optimise}{Per-level \code{optimise} list as built by \code{\link[=create_speed_input]{create_speed_input()}}.} + +\item{data}{The (factor-converted) design data frame.} +} +\value{ +\code{dots}, with \code{relationship} replaced by the prepped form when present. +} +\description{ +Prep \code{dots$relationship} once with the union of treatments seen at +every swap level. +} +\keyword{internal} diff --git a/man/dot-reject_optim_params_in_dots.Rd b/man/dot-reject_optim_params_in_dots.Rd new file mode 100644 index 0000000..29538b9 --- /dev/null +++ b/man/dot-reject_optim_params_in_dots.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/speed.R +\name{.reject_optim_params_in_dots} +\alias{.reject_optim_params_in_dots} +\title{Reject \code{...} arguments that must travel through \code{\link[=optim_params]{optim_params()}}.} +\usage{ +.reject_optim_params_in_dots(dots) +} +\arguments{ +\item{dots}{A named list captured from \code{...}.} +} +\value{ +\code{NULL}, invisibly. Called for its side effect. +} +\description{ +Stops with a message pointing the user at \code{optimise_params = optim_params(...)} +when any of the listed names is found in \code{dots}. +} +\keyword{internal} diff --git a/man/prep_relationship.Rd b/man/prep_relationship.Rd index e90df54..4ab8d6b 100644 --- a/man/prep_relationship.Rd +++ b/man/prep_relationship.Rd @@ -23,5 +23,5 @@ row/column level vectors that \code{\link[=adjacency_score_vec]{adjacency_score_ constant time per cell. \code{\link[=speed]{speed()}} calls this once before the SA loop so the per-iteration cost stays in the lookup itself. Direct callers of \code{\link[=calculate_adjacency_score]{calculate_adjacency_score()}} must wrap their matrix with this function -before passing — the score functions consume the prepped form only. +before passing - the score functions consume the prepped form only. } diff --git a/man/ring_offsets.Rd b/man/ring_offsets.Rd index 0beb847..86c84a8 100644 --- a/man/ring_offsets.Rd +++ b/man/ring_offsets.Rd @@ -13,8 +13,8 @@ ring_offsets(d, ring_type = c("manhattan", "chebyshev")) or \code{"chebyshev"} (square ring). See \code{\link[=ring_offsets]{ring_offsets()}}.} } \value{ -Integer matrix with two columns — column 1 holds \code{dx}, column 2 -holds \code{dy} — one row per cell on the ring. +Integer matrix with two columns - column 1 holds \code{dx}, column 2 +holds \code{dy} - one row per cell on the ring. } \description{ Returns the \verb{(dx, dy)} offsets of cells at exact distance \code{d} from the diff --git a/tests/testthat/test-speed.R b/tests/testthat/test-speed.R index dae3c20..c7dd43c 100644 --- a/tests/testthat/test-speed.R +++ b/tests/testthat/test-speed.R @@ -2184,6 +2184,29 @@ test_that("speed honours ring_dists and ring_weights via ...", { expect_false(isTRUE(all.equal(base$score, with_rings$score))) }) +test_that("speed errors when adj_weight/bal_weight go through ...", { + test_data <- data.frame( + row = rep(1:3, times = 3), + col = rep(1:3, each = 3), + treatment = rep(LETTERS[1:3], 3) + ) + + expect_error( + speed( + data = test_data, swap = "treatment", iterations = 5, + seed = 1, quiet = TRUE, bal_weight = 2 + ), + "must be passed via `optim_params\\(\\)`" + ) + expect_error( + speed( + data = test_data, swap = "treatment", iterations = 5, + seed = 1, quiet = TRUE, adj_weight = 0.5 + ), + "must be passed via `optim_params\\(\\)`" + ) +}) + # TODO: Test cases to add/update # - Add more detailed checking of current designs # - NSE diff --git a/vignettes/genetic_relationship.qmd b/vignettes/genetic_relationship.qmd index e20a4ba..0e8be32 100644 --- a/vignettes/genetic_relationship.qmd +++ b/vignettes/genetic_relationship.qmd @@ -38,14 +38,14 @@ library(patchwork) ## When to Use -- Trials where treatments fall into clusters of close relatives — for +- Trials where treatments fall into clusters of close relatives - for example plant-breeding lines that share parents (encoded via a pedigree A-matrix or marker-based G-matrix), variety entries from the same maturity group, or scaled mixtures of a base formulation - Anywhere "minimise neighbours of the same kind" is too coarse and you would rather "minimise neighbours of related kind" -# Example +# Grouped Treatments ## Setting Up the Design @@ -142,7 +142,8 @@ ggplot2::ggplot(rel_long, ggplot2::aes(x = col, y = row, fill = value)) + ## Running the Optimisation We run `speed()` twice with the same seed: first using the default -identity-only adjacency, then with the relationship matrix (`relationship = rel`). +identity-only adjacency, then with the relationship matrix +(`relationship = rel`). ```{r both-runs} default_result <- speed( @@ -203,6 +204,185 @@ The relationship-aware optimisation pushes whole groups apart, with the strongest pressure on group `A` (within-similarity `0.8`) and no pressure at all on group `D` (`0.0`). +# Per-Pair Similarities + +## Setting Up the Design + +Eight varieties `V01`–`V08` are evaluated on a `4 × 8` grid divided +into two `4 × 4` blocks, with two replicates of every variety per +block. Unlike the previous example, there is no group structure - +every pair of varieties carries its own similarity coefficient. + +```{r design-df-cor} +items2 <- paste0("V", 1:8) +design_df2 <- initialise_design_df( + items = rep(items2, 4), + nrows = 4, + ncols = 8, + block_nrows = 4, + block_ncols = 4 +) + +head(design_df2) +``` + +## Building the Relationship Matrix + +We populate the upper triangle with random similarities in `[0, 0.9]`, +mirror it for symmetry, and pin the diagonal to `1`. The result is a +relationship matrix in which every off-diagonal entry is distinct - +there is no plateau for the optimiser to exploit, so each pair +contributes its own adjacency cost. + +```{r relationship-cor} +n_v <- length(items2) + +rel2 <- matrix(0, n_v, n_v, dimnames = list(items2, items2)) +set.seed(112) +rel2[upper.tri(rel2)] <- runif(n_v * (n_v - 1) / 2, min = 0, max = 0.9) +rel2 <- rel2 + t(rel2) +diag(rel2) <- 1 +``` + +```{r} +#| code-fold: true +#| code-summary: "Code" +#| label: fig-relationship-cor +#| fig-cap: "Per-pair relationship matrix as a heatmap. Every off-diagonal cell is unique; the diagonal is fixed at `1`." +#| fig-width: 6 +#| fig-height: 5 +rel2_long <- as.data.frame.table(rel2, responseName = "value") +names(rel2_long)[1:2] <- c("row", "col") + +ggplot2::ggplot(rel2_long, ggplot2::aes(x = col, y = row, fill = value)) + + ggplot2::geom_tile(colour = "white") + + ggplot2::geom_text( + ggplot2::aes(label = sprintf("%.2f", value)), + size = 2.5 + ) + + ggplot2::scale_x_discrete(position = "top") + + ggplot2::scale_y_discrete(limits = rev(rownames(rel2))) + + ggplot2::scale_fill_gradient( + low = "white", high = "steelblue4", limits = c(0, 1), + name = "similarity" + ) + + ggplot2::coord_equal() + + ggplot2::labs(x = NULL, y = NULL) + + ggplot2::theme_minimal(base_size = 10) + + ggplot2::theme(panel.grid = ggplot2::element_blank()) +``` + +## Running the Optimisation + +We run `speed()` twice with the same seed: first using the default +identity-only adjacency, then with the per-pair relationship matrix +(`relationship = rel2`) and a stronger balance weight set through +`optim_params()`. + +```{r both-runs-cor} +default_result2 <- speed( + design_df2, + swap = "treatment", + swap_within = "block", + seed = 112, + quiet = TRUE +) + +related_result2 <- speed( + design_df2, + swap = "treatment", + swap_within = "block", + relationship = rel2, + optimise_params = optim_params(bal_weight = 2), + seed = 112, + quiet = TRUE +) + +default_result2 +related_result2 +``` + +## Visualising the Output + +```{r} +#| code-fold: true +#| code-summary: "Code" +#| label: fig-final-treat-cor +#| fig-cap: "Treatment view: default identity-only adjacency vs. relationship-aware adjacency under per-pair similarities." +#| fig-width: 10 +p_def_t2 <- autoplot(default_result2, treatments = "treatment") + + ggplot2::labs(title = "Default speed: treatments") +p_rel_t2 <- autoplot(related_result2, treatments = "treatment") + + ggplot2::labs(title = "speed + relationship: treatments") +p_def_t2 + p_rel_t2 + plot_layout(ncol = 2) +``` + +The default optimisation treats every non-self adjacency as equally costly. With a relationship matrix, the +algorithm can separate the similar varieties - neighbours of any given variety tend to be the varieties with the +smallest entries in its row of `rel2`. + +The placements of the closely related pairs of treatments: + +```{r} +#| code-fold: true +#| code-summary: "Code" +#| label: fig-close-treatments-cor +#| fig-cap: "Optimised layout. Only varieties involved in the high-similarity pairs are shown." +# Plot a layout where cells with NA treatments are greyed and the rest +# are labelled in red. +plot_highlighted_layout <- function(layout_df) { + x_breaks <- 1:max(layout_df$col) + y_breaks <- 1:max(layout_df$row) + ggplot2::ggplot(layout_df, ggplot2::aes(col, row, fill = treatment)) + + ggplot2::geom_tile(color = "black") + + ggplot2::scale_fill_viridis_d(na.value = "grey") + + ggplot2::scale_x_continuous(expand = c(0, 0), breaks = x_breaks) + + ggplot2::scale_y_continuous(expand = c(0, 0), breaks = y_breaks, trans = scales::reverse_trans()) + + ggplot2::geom_text( + ggplot2::aes(label = treatment), + size = 3, + color = "red", + na.rm = TRUE, + fontface = "bold" + ) + + ggplot2::theme_bw() + + ggplot2::theme( + panel.grid.major = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank(), + legend.position = "none" + ) +} + +ut <- which(upper.tri(rel2), arr.ind = TRUE) +top_orders <- order(rel2[upper.tri(rel2)], decreasing = TRUE) +plots <- NULL +for (i in 1:3) { + top <- ut[top_orders[i], , drop = FALSE] + close_treatments <- unique(c( + rownames(rel2)[top[, "row"]], + colnames(rel2)[top[, "col"]] + )) + + df_def <- default_result2$design_df + df_def$treatment <- ifelse(df_def$treatment %in% close_treatments, as.character(df_def$treatment), NA) + + df <- related_result2$design_df + df$treatment <- ifelse(df$treatment %in% close_treatments, as.character(df$treatment), NA) + + hl_def <- plot_highlighted_layout(df_def) + hl_rel <- plot_highlighted_layout(df) + if (is.null(plots)) { + plots <- hl_def + hl_rel + } else { + plots <- plots + hl_def + hl_rel + } +} +plots + plot_layout(ncol = 2) +``` + +Without a relationship matrix, these related treatments can be placed next to each other. This can be optimised +by providing such matrix. +