From c94bf63716f97ffdc7d4617e000cddef64bc1cb1 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 11:42:22 +0100
Subject: [PATCH 01/10] `plot.ArrayTimeBySpeciesBySize` now also works on
 single-species model

---
 R/ArrayTimeBySpeciesBySize-class.R             |  4 +++-
 tests/testthat/test-ArrayTimeBySpeciesBySize.R | 15 +++++++++++++++
 2 files changed, 18 insertions(+), 1 deletion(-)

diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index 27cb2498c..0e6a14d65 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -158,7 +158,9 @@ plot.ArrayTimeBySpeciesBySize <- function(x, species = NULL, time = NULL,
     }
 
     arr <- unclass(x)
-    slice <- arr[tidx, , ]
+    slice <- matrix(arr[tidx, , , drop = FALSE],
+                    nrow = dim(arr)[2],
+                    dimnames = dimnames(arr)[2:3])
     slice <- ArraySpeciesBySize(slice, value_name = value_name,
                                 units = units, params = params)
 
diff --git a/tests/testthat/test-ArrayTimeBySpeciesBySize.R b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
index eef55c0d6..641287ec9 100644
--- a/tests/testthat/test-ArrayTimeBySpeciesBySize.R
+++ b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
@@ -86,6 +86,21 @@ test_that("plot.ArrayTimeBySpeciesBySize time argument selects correct slice", {
     expect_true(is.data.frame(p))
 })
 
+test_that("plot.ArrayTimeBySpeciesBySize preserves single species dimension", {
+    arr <- array(seq_len(6), dim = c(2, 1, 3),
+                 dimnames = list(time = c("2000", "2001"),
+                                 sp = "Cod",
+                                 w = c("1", "10", "100")))
+    rate <- ArrayTimeBySpeciesBySize(arr, value_name = "Test rate",
+                                     units = "1/year")
+
+    p <- plot(rate, time = 2001, return_data = TRUE)
+
+    expect_true(is.data.frame(p))
+    expect_identical(p$Species, rep("Cod", 3))
+    expect_equal(p$value, unname(arr["2001", "Cod", ]))
+})
+
 test_that("as.data.frame.ArrayTimeBySpeciesBySize returns correct structure", {
     fmort <- getFMort(NS_sim)
     df <- as.data.frame(fmort)

From 2520dce193d1699463eda4bd0a137255a1d95165 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 12:09:36 +0100
Subject: [PATCH 02/10] Enhance ArrayTimeBySpeciesBySize subsetting to return
 appropriate object types for single time and size selections

---
 NEWS.md                                       |  4 ++++
 R/ArrayTimeBySpeciesBySize-class.R            | 16 +++++++++++++++
 .../testthat/test-ArrayTimeBySpeciesBySize.R  | 20 ++++++++++++++++++-
 3 files changed, 39 insertions(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index a1624af0b..1d74cc871 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -112,6 +112,10 @@ individual variability in growth to be modelled.
   `ArrayTimeBySpeciesBySize`. An `animate()` method allows interactive
   playback.
 
+- Subsetting an `ArrayTimeBySpeciesBySize` object returns an
+  `ArraySpeciesBySize` object when a single time is selected, and an
+  `ArrayTimeBySpecies` object when a single size is selected.
+
 - New `ggplotly()` methods for `ArraySpeciesBySize`, `ArrayTimeBySpecies`, and
   `ArrayTimeBySpeciesBySize` convert the `ggplot2` output of `plot()` into
   interactive plotly figures.
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index 0e6a14d65..72bfc8a92 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -296,6 +296,22 @@ as.data.frame.ArrayTimeBySpeciesBySize <- function(x, row.names = NULL,
         attr(result, "units") <- attr(x, "units")
         attr(result, "params") <- attr(x, "params")
         class(result) <- c("ArrayTimeBySpeciesBySize", "array")
+    } else if (is.matrix(result)) {
+        dim_names <- names(dimnames(result))
+        attrs <- list(value_name = attr(x, "value_name"),
+                      units = attr(x, "units"),
+                      params = attr(x, "params"))
+        if (identical(dim_names, c("sp", "w"))) {
+            result <- ArraySpeciesBySize(result,
+                                         value_name = attrs$value_name,
+                                         units = attrs$units,
+                                         params = attrs$params)
+        } else if (identical(dim_names, c("time", "sp"))) {
+            result <- ArrayTimeBySpecies(result,
+                                         value_name = attrs$value_name,
+                                         units = attrs$units,
+                                         params = attrs$params)
+        }
     }
     result
 }
diff --git a/tests/testthat/test-ArrayTimeBySpeciesBySize.R b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
index 641287ec9..98967a4a8 100644
--- a/tests/testthat/test-ArrayTimeBySpeciesBySize.R
+++ b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
@@ -117,10 +117,28 @@ test_that("[.ArrayTimeBySpeciesBySize preserves class for 3D result", {
     expect_identical(attr(sub, "value_name"), "Fishing mortality")
 })
 
-test_that("[.ArrayTimeBySpeciesBySize drops class for 2D result", {
+test_that("[.ArrayTimeBySpeciesBySize returns ArraySpeciesBySize when time is dropped", {
     fmort <- getFMort(NS_sim)
     slice <- fmort[1, , ]
     expect_false(is.ArrayTimeBySpeciesBySize(slice))
+    expect_true(is.ArraySpeciesBySize(slice))
+    expect_identical(attr(slice, "value_name"), "Fishing mortality")
+})
+
+test_that("[.ArrayTimeBySpeciesBySize returns ArrayTimeBySpecies when size is dropped", {
+    fmort <- getFMort(NS_sim)
+    slice <- fmort[, , 1]
+    expect_false(is.ArrayTimeBySpeciesBySize(slice))
+    expect_true(is.ArrayTimeBySpecies(slice))
+    expect_identical(attr(slice, "value_name"), "Fishing mortality")
+})
+
+test_that("[.ArrayTimeBySpeciesBySize leaves time by size matrices plain", {
+    fmort <- getFMort(NS_sim)
+    slice <- fmort[, 1, ]
+    expect_false(is.ArrayTimeBySpeciesBySize(slice))
+    expect_false(is.ArraySpeciesBySize(slice))
+    expect_false(is.ArrayTimeBySpecies(slice))
     expect_true(is.matrix(slice))
 })
 

From b3f59204c272c8b02b642818e62e3695acec1a22 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 12:18:08 +0100
Subject: [PATCH 03/10] Add log scaling options to plotSpectra

---
 NEWS.md                     |  3 +++
 R/plots.R                   | 29 ++++++++++++++++++++++++-----
 tests/testthat/test-plots.R | 21 +++++++++++++++++++++
 3 files changed, 48 insertions(+), 5 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 1d74cc871..6349c9b0a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -129,6 +129,9 @@ individual variability in growth to be modelled.
 - New `addPlot()` generic with methods for adding `ArraySpeciesBySize` and
   `ArrayTimeBySpecies` values as extra lines on an existing compatible ggplot.
 
+- `plotSpectra()` now accepts `log_x`, `log_y`, and `log` arguments for
+  controlling axis scaling, matching the mizer array `plot()` methods.
+
 - The `plot()` and `summary()` methods for `MizerParams`, `MizerSim`, and the
   mizer array classes are now registered as S3 methods rather than S4 methods,
   so `plot()` and `summary()` remain plain S3 generics when mizer is loaded,
diff --git a/R/plots.R b/R/plots.R
index 849187adb..1c30a459c 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -627,6 +627,11 @@ plotlyYieldGear <- function(sim, species = NULL,
 #'   are included. Ignored if the model does not contain background species.
 #'   Default is TRUE.
 #' @param highlight Name or vector of names of the species to be highlighted.
+#' @param log_x If `TRUE` (default), use a log10 x-axis.
+#' @param log_y If `TRUE` (default), use a log10 y-axis.
+#' @param log Character string specifying which axes should use log10 scales,
+#'   in the same form as the base [plot()] argument. For example, `"x"`,
+#'   `"y"`, `"xy"` or `""`. If supplied, this overrides `log_x` and `log_y`.
 #' @param return_data A boolean value that determines whether the formatted data
 #' used for the plot is returned instead of the plot itself. Default value is FALSE
 #' @param ... Other arguments (currently unused)
@@ -666,11 +671,16 @@ plotSpectra.MizerSim <- function(object, species = NULL,
                         power = 1, biomass = TRUE,
                         total = FALSE, resource = TRUE,
                         background = TRUE,
-                        highlight = NULL, return_data = FALSE, ...) {
+                        highlight = NULL, log_x = TRUE, log_y = TRUE,
+                        log = NULL, return_data = FALSE, ...) {
     # to deal with old-type biomass argument
     if (missing(power)) {
         power <- as.numeric(biomass)
     }
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
+    log_x <- log_axes$log_x
+    log_y <- log_axes$log_y
+
     assert_that(is.flag(total), is.flag(resource),
                 is.flag(background),
                 is.number(power),
@@ -696,6 +706,7 @@ plotSpectra.MizerSim <- function(object, species = NULL,
                  species = species, wlim = wlim, ylim = ylim,
                  power = power, total = total, resource = resource,
                  background = background, highlight = highlight,
+                 log_x = log_x, log_y = log_y,
                  return_data = return_data)
 }
 
@@ -706,11 +717,16 @@ plotSpectra.MizerParams <- function(object, species = NULL,
                         power = 1, biomass = TRUE,
                         total = FALSE, resource = TRUE,
                         background = TRUE,
-                        highlight = NULL, return_data = FALSE, ...) {
+                        highlight = NULL, log_x = TRUE, log_y = TRUE,
+                        log = NULL, return_data = FALSE, ...) {
     # to deal with old-type biomass argument
     if (missing(power)) {
         power <- as.numeric(biomass)
     }
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
+    log_x <- log_axes$log_x
+    log_y <- log_axes$log_y
+
     assert_that(is.flag(total), is.flag(resource),
                 is.flag(background),
                 is.number(power),
@@ -725,6 +741,7 @@ plotSpectra.MizerParams <- function(object, species = NULL,
                  species = species, wlim = wlim, ylim = ylim,
                  power = power, total = total, resource = resource,
                  background = background, highlight = highlight,
+                 log_x = log_x, log_y = log_y,
                  return_data = return_data)
 }
 
@@ -732,7 +749,7 @@ plotSpectra.MizerParams <- function(object, species = NULL,
 plot_spectra <- function(params, n, n_pp,
                          species, wlim, ylim, power,
                          total, resource, background,
-                         highlight, return_data) {
+                         highlight, log_x, log_y, return_data) {
     params <- validParams(params)
     if (is.na(wlim[1])) {
         wlim[1] <- if (resource) min(params@w) / 100 else min(params@w)
@@ -813,7 +830,8 @@ plot_spectra <- function(params, n, n_pp,
     if (return_data) return(plot_dat)
 
     plotDataFrame(plot_dat, params, xlab = "Size [g]", ylab = y_label,
-                  xtrans = "log10", ytrans = "log10",
+                  xtrans = if (log_x) "log10" else "identity",
+                  ytrans = if (log_y) "log10" else "identity",
                   xlim = wlim, ylim = ylim,
                   highlight = highlight, legend_var = "Legend")
 }
@@ -826,7 +844,8 @@ plotlySpectra <- function(object, species = NULL,
                         power = 1, biomass = TRUE,
                         total = FALSE, resource = TRUE,
                         background = TRUE,
-                        highlight = NULL, ...) {
+                        highlight = NULL, log_x = TRUE, log_y = TRUE,
+                        log = NULL, ...) {
     argg <- as.list(environment())
     ggplotly(do.call("plotSpectra", argg),
              tooltip = c("Species", "w", "value"))
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 7db2b6256..141874cb2 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -203,6 +203,27 @@ test_that("plotSpectra validates empty selection and can return total only", {
     expect_true(all(df$Legend == "Total"))
 })
 
+test_that("plotSpectra supports base plot log argument", {
+    p_y <- plotSpectra(params, species = species, log = "y")
+    expect_identical(p_y$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_y$scales$get_scales("y")$trans$name, "log-10")
+
+    p_xy <- plotSpectra(params, species = species, log = "xy")
+    expect_identical(p_xy$scales$get_scales("x")$trans$name, "log-10")
+    expect_identical(p_xy$scales$get_scales("y")$trans$name, "log-10")
+
+    p_none <- plotSpectra(params, species = species, log = "")
+    expect_identical(p_none$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_none$scales$get_scales("y")$trans$name, "identity")
+
+    p_sim <- plotSpectra(sim, species = species, log_x = FALSE, log_y = FALSE)
+    expect_identical(p_sim$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_sim$scales$get_scales("y")$trans$name, "identity")
+
+    expect_error(plotSpectra(params, species = species, log = "z"),
+                 "`log` must be a character string")
+})
+
 test_that("plotPredMort and plotFMort trim to species size range by default", {
     pred_trimmed <- plotPredMort(params, species = 2, return_data = TRUE)
     pred_full <- plotPredMort(params, species = 2, all.sizes = TRUE, return_data = TRUE)

From a6ba4e9c40a35b5c622596b6c89f9bc51ce2e299 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 13:02:44 +0100
Subject: [PATCH 04/10] Add `plot2()` generic and related methods for comparing
 mizer array objects. Moved `plotSpectra2` from mizeExperimental to mizer

---
 NAMESPACE                                     |   5 +
 NEWS.md                                       |   5 +
 R/ArraySpeciesBySize-class.R                  |  89 ++++++++++++
 R/ArrayTimeBySpecies-class.R                  |  33 +++++
 R/ArrayTimeBySpeciesBySize-class.R            |  47 +++++++
 R/plots.R                                     | 129 ++++++++++++++++--
 pkgdown/_pkgdown.yml                          |   4 +-
 tests/testthat/test-ArraySpeciesBySize.R      |  21 +++
 tests/testthat/test-ArrayTimeBySpecies.R      |  32 +++++
 .../testthat/test-ArrayTimeBySpeciesBySize.R  |  17 +++
 tests/testthat/test-plots.R                   |  34 +++++
 11 files changed, 405 insertions(+), 11 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 286267939..28db4e328 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -125,6 +125,9 @@ S3method(plot,ArrayTimeBySpecies)
 S3method(plot,ArrayTimeBySpeciesBySize)
 S3method(plot,MizerParams)
 S3method(plot,MizerSim)
+S3method(plot2,ArraySpeciesBySize)
+S3method(plot2,ArrayTimeBySpecies)
+S3method(plot2,ArrayTimeBySpeciesBySize)
 S3method(plotBiomass,MizerSim)
 S3method(plotBiomassObservedVsModel,MizerParams)
 S3method(plotBiomassObservedVsModel,MizerSim)
@@ -396,6 +399,7 @@ export(newSingleSpeciesParams)
 export(newTraitParams)
 export(noRDD)
 export(other_params)
+export(plot2)
 export(plotBiomass)
 export(plotBiomassObservedVsModel)
 export(plotDataFrame)
@@ -406,6 +410,7 @@ export(plotGrowthCurves)
 export(plotM2)
 export(plotPredMort)
 export(plotSpectra)
+export(plotSpectra2)
 export(plotYield)
 export(plotYieldGear)
 export(plotYieldObservedVsModel)
diff --git a/NEWS.md b/NEWS.md
index 6349c9b0a..2a6e15c1f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -129,6 +129,11 @@ individual variability in growth to be modelled.
 - New `addPlot()` generic with methods for adding `ArraySpeciesBySize` and
   `ArrayTimeBySpecies` values as extra lines on an existing compatible ggplot.
 
+- New `plot2()` generic with methods for comparing two compatible mizer array
+  objects in one plot, with species or group shown by colour and model by
+  linetype. The `plotSpectra2()` helper has moved from `mizerExperimental` into
+  mizer for comparing two abundance spectra.
+
 - `plotSpectra()` now accepts `log_x`, `log_y`, and `log` arguments for
   controlling axis scaling, matching the mizer array `plot()` methods.
 
diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R
index d71bca1aa..a78307fff 100644
--- a/R/ArraySpeciesBySize-class.R
+++ b/R/ArraySpeciesBySize-class.R
@@ -228,6 +228,95 @@ parsePlotLog <- function(log, log_x = FALSE, log_y = FALSE) {
     )
 }
 
+#' Compare two mizer array objects in one plot
+#'
+#' `plot2()` compares two compatible mizer array objects in a single ggplot.
+#' Colours identify species or groups, and linetype identifies which object the
+#' values came from.
+#'
+#' @param x,y Two compatible mizer array objects of the same class.
+#' @param name1,name2 Labels for the two objects, used in the linetype legend.
+#' @inheritParams plot
+#'
+#' @return A ggplot2 object.
+#' @export
+#' @family plotting functions
+#'
+#' @examples
+#' \donttest{
+#' enc <- getEncounter(NS_params)
+#' plot2(enc, enc, name1 = "Original", name2 = "Changed")
+#' plot2(getBiomass(NS_sim), getBiomass(NS_sim), species = "Cod")
+#' }
+plot2 <- function(x, y, ...) {
+    UseMethod("plot2", x)
+}
+
+#' @rdname plot2
+#' @export
+plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
+                                     species = NULL, all.sizes = FALSE,
+                                     log_x = TRUE, log_y = FALSE, log = NULL,
+                                     wlim = c(NA, NA), ylim = c(NA, NA),
+                                     total = FALSE, background = TRUE,
+                                     y_ticks = 6, ...) {
+    check_plot2_compatible(x, y, "ArraySpeciesBySize")
+    compare_array_metadata(x, y)
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
+    log_x <- log_axes$log_x
+    log_y <- log_axes$log_y
+
+    params <- attr(x, "params")
+    y_label <- array_y_label(x, default = "Rate")
+    plot_dat1 <- prepare_ArraySpeciesBySize_plot_data(
+        x, species = species, all.sizes = all.sizes, wlim = wlim,
+        total = total, background = background)
+    plot_dat2 <- prepare_ArraySpeciesBySize_plot_data(
+        y, species = species, all.sizes = all.sizes, wlim = wlim,
+        total = total, background = background)
+
+    plotComparisonDataFrame(plot_dat1, plot_dat2, params,
+                            name1 = name1, name2 = name2,
+                            xlab = "Size [g]", ylab = y_label,
+                            xtrans = if (log_x) "log10" else "identity",
+                            ytrans = if (log_y) "log10" else "identity",
+                            xlim = wlim, ylim = ylim,
+                            y_ticks = y_ticks, legend_var = "Legend")
+}
+
+check_plot2_compatible <- function(x, y, class) {
+    if (!inherits(y, class)) {
+        stop("Both objects must be of class `", class, "`.")
+    }
+}
+
+compare_array_metadata <- function(x, y) {
+    value_name1 <- attr(x, "value_name")
+    value_name2 <- attr(y, "value_name")
+    if (!is.null(value_name1) && !is.null(value_name2) &&
+            !identical(value_name1, value_name2)) {
+        warning("The first array has value name `", value_name1,
+                "`, but the second array has value name `", value_name2, "`.")
+    }
+    units1 <- attr(x, "units")
+    units2 <- attr(y, "units")
+    if (!is.null(units1) && !is.null(units2) &&
+            nzchar(units1) && nzchar(units2) &&
+            !identical(units1, units2)) {
+        warning("The first array has y units `", units1,
+                "`, but the second array has y units `", units2, "`.")
+    }
+}
+
+array_y_label <- function(x, default = "Value") {
+    value_name <- attr(x, "value_name") %||% default
+    units_str <- attr(x, "units")
+    if (!is.null(units_str) && nzchar(units_str)) {
+        value_name <- paste0(value_name, " [", units_str, "]")
+    }
+    value_name
+}
+
 #' Add values to an existing plot
 #'
 #' `r lifecycle::badge("experimental")`
diff --git a/R/ArrayTimeBySpecies-class.R b/R/ArrayTimeBySpecies-class.R
index de1b695fc..b206005b2 100644
--- a/R/ArrayTimeBySpecies-class.R
+++ b/R/ArrayTimeBySpecies-class.R
@@ -234,6 +234,39 @@ addPlot.ArrayTimeBySpecies <- function(plot, x, species = NULL,
     plot + do.call(geom_line, layer_args)
 }
 
+#' @rdname plot2
+#' @export
+plot2.ArrayTimeBySpecies <- function(x, y, name1 = "First", name2 = "Second",
+                                     species = NULL,
+                                     start_time = NULL, end_time = NULL,
+                                     y_ticks = 6, ylim = c(NA, NA),
+                                     total = FALSE, background = TRUE,
+                                     log_x = FALSE, log_y = TRUE,
+                                     log = NULL, ...) {
+    check_plot2_compatible(x, y, "ArrayTimeBySpecies")
+    compare_array_metadata(x, y)
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
+    log_x <- log_axes$log_x
+    log_y <- log_axes$log_y
+
+    params <- attr(x, "params")
+    y_label <- array_y_label(x, default = "Value")
+    plot_dat1 <- prepare_ArrayTimeBySpecies_plot_data(
+        x, species = species, start_time = start_time, end_time = end_time,
+        ylim = ylim, total = total, background = background)
+    plot_dat2 <- prepare_ArrayTimeBySpecies_plot_data(
+        y, species = species, start_time = start_time, end_time = end_time,
+        ylim = ylim, total = total, background = background)
+
+    plotComparisonDataFrame(plot_dat1, plot_dat2, params,
+                            name1 = name1, name2 = name2,
+                            xlab = "Year", ylab = y_label,
+                            xtrans = if (log_x) "log10" else "identity",
+                            ytrans = if (log_y) "log10" else "identity",
+                            ylim = ylim, y_ticks = y_ticks,
+                            legend_var = "Legend")
+}
+
 prepare_ArrayTimeBySpecies_plot_data <- function(x, species = NULL,
                                                  start_time = NULL,
                                                  end_time = NULL,
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index 72bfc8a92..0c0e04640 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -171,6 +171,53 @@ plot.ArrayTimeBySpeciesBySize <- function(x, species = NULL, time = NULL,
                             background = background, y_ticks = y_ticks, ...)
 }
 
+#' @rdname plot2
+#'
+#' @param time The time to display. Default (`NULL`) is the final time step.
+#'   Only applies to `ArrayTimeBySpeciesBySize`.
+#' @export
+plot2.ArrayTimeBySpeciesBySize <- function(x, y, name1 = "First",
+                                           name2 = "Second",
+                                           species = NULL, time = NULL,
+                                           all.sizes = FALSE,
+                                           log_x = TRUE, log_y = FALSE,
+                                           log = NULL,
+                                           wlim = c(NA, NA),
+                                           ylim = c(NA, NA),
+                                           total = FALSE,
+                                           background = TRUE,
+                                           y_ticks = 6, ...) {
+    check_plot2_compatible(x, y, "ArrayTimeBySpeciesBySize")
+    slice1 <- ArrayTimeBySpeciesBySize_slice(x, time = time)
+    slice2 <- ArrayTimeBySpeciesBySize_slice(y, time = time)
+
+    plot2.ArraySpeciesBySize(slice1, slice2, name1 = name1, name2 = name2,
+                             species = species, all.sizes = all.sizes,
+                             log_x = log_x, log_y = log_y, log = log,
+                             wlim = wlim, ylim = ylim, total = total,
+                             background = background, y_ticks = y_ticks, ...)
+}
+
+ArrayTimeBySpeciesBySize_slice <- function(x, time = NULL) {
+    params <- attr(x, "params")
+    value_name <- attr(x, "value_name")
+    units <- attr(x, "units")
+
+    times <- as.numeric(dimnames(x)[[1]])
+    if (is.null(time)) {
+        tidx <- dim(x)[1]
+    } else {
+        tidx <- which.min(abs(times - time))
+    }
+
+    arr <- unclass(x)
+    slice <- matrix(arr[tidx, , , drop = FALSE],
+                    nrow = dim(arr)[2],
+                    dimnames = dimnames(arr)[2:3])
+    ArraySpeciesBySize(slice, value_name = value_name,
+                       units = units, params = params)
+}
+
 #' @rdname plot
 #' @exportS3Method plotly::ggplotly
 #' @examples
diff --git a/R/plots.R b/R/plots.R
index 1c30a459c..53d0092c8 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -41,8 +41,8 @@
 #' The same array objects can be passed to [ggplotly()] to produce interactive
 #' versions, for example `ggplotly(getBiomass(sim))` or
 #' `ggplotly(getEncounter(params))`. To add another compatible array to an
-#' existing ggplot, use [addPlot()]. This is useful for comparing two simulations
-#' or two parameter sets on the same axes. To visualise how spectra or rates
+#' existing ggplot, use [addPlot()]. To compare two compatible mizer arrays
+#' directly, use [plot2()]. To visualise how spectra or rates
 #' change through time, use [animate()] on a `MizerSim` or an
 #' `ArrayTimeBySpeciesBySize` object.
 #'
@@ -54,6 +54,7 @@
 #'   [plotYield()] \tab Plots the total yield of each species across all fishing gears against time. \cr
 #'   [plotYieldGear()] \tab Plots the total yield of each species by gear against time. \cr
 #'   [plotSpectra()] \tab Plots the abundance (biomass or numbers) spectra of each species and the background community. It is possible to specify a minimum size which is useful for truncating the plot. \cr
+#'   [plotSpectra2()] \tab Compares the spectra from two simulations or parameter objects in one plot. \cr
 #'   [plotFeedingLevel()] \tab Plots the feeding level of each species against size. \cr
 #'   [plotPredMort()] \tab Plots the predation mortality of each species against size. \cr
 #'   [plotFMort()] \tab Plots the total fishing mortality of each species against size. \cr
@@ -128,7 +129,8 @@ NULL
 utils::globalVariables(c("time", "value", "Species", "w", "gear", "Age",
                          "x", "y", "Year", "Yield", "Biomass", "Size",
                          "Proportion", "Prey", "Legend", "Type", "Gear",
-                         "Predator", "weight", "a", "b", "age", "w_max"))
+                         "Predator", "weight", "a", "b", "age", "w_max",
+                         "Model"))
 
 #' Make a plot from a data frame
 #'
@@ -250,6 +252,58 @@ plotDataFrame <- function(frame, params, style = "line", xlab = waiver(),
     p
 }
 
+plotComparisonDataFrame <- function(frame1, frame2, params,
+                                    name1 = "First", name2 = "Second",
+                                    xlab = waiver(), ylab = waiver(),
+                                    xtrans = "identity", ytrans = "identity",
+                                    xlim = c(NA, NA), ylim = c(NA, NA),
+                                    y_ticks = 6, legend_var = "Legend") {
+    assert_that(is.data.frame(frame1),
+                is.data.frame(frame2),
+                is(params, "MizerParams"))
+
+    names(frame2)[seq_len(min(3, ncol(frame2)))] <-
+        names(frame1)[seq_len(min(3, ncol(frame1)))]
+    frame1$Model <- name1
+    frame2$Model <- name2
+    frame <- rbind(frame1, frame2)
+    frame$Model <- factor(frame$Model, levels = c(name1, name2))
+
+    var_names <- names(frame)
+    x_var <- var_names[[1]]
+    y_var <- var_names[[2]]
+    group_var <- var_names[[3]]
+    if (!(legend_var %in% var_names)) {
+        stop("The `legend_var` argument must be the name of a variable ",
+             "in the data frame.")
+    }
+
+    legend_levels <- intersect(names(params@linecolour), frame[[legend_var]])
+    frame[[legend_var]] <- factor(frame[[legend_var]], levels = legend_levels)
+    if (sum(is.na(frame[[legend_var]]))) {
+        warning("missing legend in params@linecolour, some groups won't be displayed")
+    }
+    linecolour <- params@linecolour[legend_levels]
+
+    xbreaks <- waiver()
+    if (xtrans == "log10") xbreaks <- log_breaks()
+    ybreaks <- waiver()
+    if (ytrans == "log10") ybreaks <- log_breaks(n = y_ticks)
+
+    ggplot(frame,
+           aes(group = interaction(.data[[group_var]], .data[["Model"]]))) +
+        scale_y_continuous(trans = ytrans, breaks = ybreaks,
+                           labels = prettyNum, name = ylab,
+                           limits = ylim) +
+        scale_x_continuous(trans = xtrans, breaks = xbreaks, name = xlab,
+                           limits = xlim) +
+        geom_line(aes(x = .data[[x_var]], y = .data[[y_var]],
+                      colour = .data[[legend_var]],
+                      linetype = .data[["Model"]])) +
+        scale_colour_manual(values = linecolour) +
+        scale_linetype_discrete(drop = FALSE)
+}
+
 #' Helper function to produce nice breaks on logarithmic axes
 #'
 #' This is needed when the logarithmic y-axis spans less than one order of
@@ -767,12 +821,7 @@ plot_spectra <- function(params, n, n_pp,
     }
     species <- valid_species_arg(params, species)
     # Deal with power argument
-    if (power %in% c(0, 1, 2)) {
-        y_label <- c("Number density [1/g]", "Biomass density",
-                    "Biomass density [g]")[power + 1]
-    } else {
-        y_label <- paste0("Number density * w^", power)
-    }
+    y_label <- spectra_y_label(power)
     n <- sweep(n, 2, params@w^power, "*")
     # Select only the desired species
     spec_n <- n[as.character(dimnames(n)[[1]]) %in% species, , drop = FALSE]
@@ -836,6 +885,68 @@ plot_spectra <- function(params, n, n_pp,
                   highlight = highlight, legend_var = "Legend")
 }
 
+#' Compare two size spectra in the same plot
+#'
+#' `plotSpectra2()` compares the abundance spectra from two `MizerParams` or
+#' `MizerSim` objects in a single plot. Colours identify species or groups and
+#' linetype identifies the object.
+#'
+#' @param object1 First `MizerParams` or `MizerSim` object.
+#' @param object2 Second `MizerParams` or `MizerSim` object.
+#' @param name1,name2 Labels for the two objects, used in the linetype legend.
+#' @inheritParams plotSpectra
+#' @param log_x If `TRUE` (default), use a log10 x-axis.
+#' @param log_y If `TRUE` (default), use a log10 y-axis.
+#' @param log Character string specifying which axes should use log10 scales,
+#'   in the same form as the base [plot()] argument. For example, `"x"`,
+#'   `"y"`, `"xy"` or `""`. If supplied, this overrides `log_x` and `log_y`.
+#' @param ... Arguments passed to [plotSpectra()] for preparing the spectra
+#'   data, for example `species`, `time_range`, `wlim`, `ylim`, `resource`,
+#'   `background` or `total`.
+#'
+#' @return A ggplot2 object.
+#' @export
+#' @family plotting functions
+#'
+#' @examples
+#' \donttest{
+#' sim1 <- project(NS_params, t_max = 10, progress_bar = FALSE)
+#' sim2 <- project(NS_params, effort = 0.5, t_max = 10, progress_bar = FALSE)
+#' plotSpectra2(sim1, sim2, "Original", "Effort = 0.5")
+#' }
+plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
+                         power = 1, log_x = TRUE, log_y = TRUE,
+                         log = NULL, ...) {
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
+    log_x <- log_axes$log_x
+    log_y <- log_axes$log_y
+
+    args <- list(...)
+    wlim <- args$wlim %||% c(NA, NA)
+    ylim <- args$ylim %||% c(NA, NA)
+
+    sf1 <- plotSpectra(object1, power = power, return_data = TRUE, ...)
+    sf2 <- plotSpectra(object2, power = power, return_data = TRUE, ...)
+    params <- if (is(object1, "MizerSim")) object1@params else object1
+
+    plotComparisonDataFrame(sf1, sf2, validParams(params),
+                            name1 = name1, name2 = name2,
+                            xlab = "Size [g]",
+                            ylab = spectra_y_label(power),
+                            xtrans = if (log_x) "log10" else "identity",
+                            ytrans = if (log_y) "log10" else "identity",
+                            xlim = wlim, ylim = ylim,
+                            legend_var = "Legend")
+}
+
+spectra_y_label <- function(power) {
+    if (power %in% c(0, 1, 2)) {
+        return(c("Number density [1/g]", "Biomass density",
+                 "Biomass density [g]")[power + 1])
+    }
+    paste0("Number density * w^", power)
+}
+
 #' @rdname plotSpectra
 #' @export
 plotlySpectra <- function(object, species = NULL,
diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml
index 8c6a67517..90856e3ad 100644
--- a/pkgdown/_pkgdown.yml
+++ b/pkgdown/_pkgdown.yml
@@ -171,10 +171,10 @@ reference:
 - title: Plotting results
   contents:
   - plotting_functions
-  - has_concept("plotting functions")
-  - addPlot
   - starts_with("animate")
   - starts_with("plot")
+  - addPlot
+  - has_concept("plotting functions")
   - -plotM2
   - setColours
   - setLinetypes
diff --git a/tests/testthat/test-ArraySpeciesBySize.R b/tests/testthat/test-ArraySpeciesBySize.R
index e887bed4e..b53d6d17a 100644
--- a/tests/testthat/test-ArraySpeciesBySize.R
+++ b/tests/testthat/test-ArraySpeciesBySize.R
@@ -108,6 +108,27 @@ test_that("plot.ArraySpeciesBySize supports base plot log argument", {
     expect_error(plot(enc, log = "z"), "`log` must be a character string")
 })
 
+test_that("plot2.ArraySpeciesBySize compares compatible arrays", {
+    enc <- getEncounter(NS_params)
+
+    p <- plot2(enc, enc, name1 = "Original", name2 = "Changed",
+               species = "Cod", total = TRUE, background = FALSE,
+               wlim = c(1, NA), log = "xy")
+    expect_s3_class(p, "ggplot")
+    expect_identical(levels(p$data$Model), c("Original", "Changed"))
+    expect_true(all(p$data$Species %in% c("Cod", "Total")))
+    expect_true(all(p$data$w >= 1))
+    expect_identical(p$scales$get_scales("x")$trans$name, "log-10")
+    expect_identical(p$scales$get_scales("y")$trans$name, "log-10")
+
+    p_none <- plot2(enc, enc, species = "Cod", log = "")
+    expect_identical(p_none$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_none$scales$get_scales("y")$trans$name, "identity")
+
+    expect_error(plot2(enc, getBiomass(NS_sim)), "Both objects must be")
+    expect_error(plot2(enc, enc, log = "z"), "`log` must be a character string")
+})
+
 test_that("addPlot.ArraySpeciesBySize adds lines to an existing ggplot", {
     enc <- getEncounter(NS_params)
     pred_mort <- getPredMort(NS_params)
diff --git a/tests/testthat/test-ArrayTimeBySpecies.R b/tests/testthat/test-ArrayTimeBySpecies.R
index 004dc9fca..08924e656 100644
--- a/tests/testthat/test-ArrayTimeBySpecies.R
+++ b/tests/testthat/test-ArrayTimeBySpecies.R
@@ -78,6 +78,38 @@ test_that("plot.ArrayTimeBySpecies supports base plot log argument", {
     expect_error(plot(bio, log = TRUE), "`log` must be a character string")
 })
 
+test_that("plot2.ArrayTimeBySpecies compares compatible arrays", {
+    bio <- getBiomass(NS_sim)
+    years <- as.numeric(rownames(bio))
+
+    p <- plot2(bio, bio, name1 = "Original", name2 = "Changed",
+               species = "Cod", total = TRUE, start_time = years[2],
+               end_time = years[5], log = "xy")
+    expect_s3_class(p, "ggplot")
+    expect_identical(levels(p$data$Model), c("Original", "Changed"))
+    expect_true(all(p$data$Species %in% c("Cod", "Total")))
+    expect_true(all(p$data$Year >= years[2]))
+    expect_true(all(p$data$Year <= years[5]))
+    expect_identical(p$scales$get_scales("x")$trans$name, "log-10")
+    expect_identical(p$scales$get_scales("y")$trans$name, "log-10")
+
+    p_none <- plot2(bio, bio, species = "Cod", log = "")
+    expect_identical(p_none$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_none$scales$get_scales("y")$trans$name, "identity")
+
+    warnings <- character()
+    withCallingHandlers(
+        plot2(bio, getYield(NS_sim), species = "Cod"),
+        warning = function(w) {
+            warnings <<- c(warnings, conditionMessage(w))
+            invokeRestart("muffleWarning")
+        }
+    )
+    expect_true(any(grepl("value name", warnings)))
+    expect_true(any(grepl("y units", warnings)))
+    expect_error(plot2(bio, getEncounter(NS_params)), "Both objects must be")
+})
+
 test_that("addPlot.ArrayTimeBySpecies adds lines to an existing ggplot", {
     bio <- getBiomass(NS_sim)
     yield <- getYield(NS_sim)
diff --git a/tests/testthat/test-ArrayTimeBySpeciesBySize.R b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
index 98967a4a8..680a28602 100644
--- a/tests/testthat/test-ArrayTimeBySpeciesBySize.R
+++ b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
@@ -86,6 +86,23 @@ test_that("plot.ArrayTimeBySpeciesBySize time argument selects correct slice", {
     expect_true(is.data.frame(p))
 })
 
+test_that("plot2.ArrayTimeBySpeciesBySize compares selected time slices", {
+    fmort <- getFMort(NS_sim)
+    times <- as.numeric(dimnames(fmort)[[1]])
+
+    p <- plot2(fmort, fmort, name1 = "Original", name2 = "Changed",
+               species = "Cod", time = times[5], total = TRUE,
+               wlim = c(1, NA), log = "xy")
+    expect_s3_class(p, "ggplot")
+    expect_identical(levels(p$data$Model), c("Original", "Changed"))
+    expect_true(all(p$data$Species %in% c("Cod", "Total")))
+    expect_true(all(p$data$w >= 1))
+    expect_identical(p$scales$get_scales("x")$trans$name, "log-10")
+    expect_identical(p$scales$get_scales("y")$trans$name, "log-10")
+
+    expect_error(plot2(fmort, getBiomass(NS_sim)), "Both objects must be")
+})
+
 test_that("plot.ArrayTimeBySpeciesBySize preserves single species dimension", {
     arr <- array(seq_len(6), dim = c(2, 1, 3),
                  dimnames = list(time = c("2000", "2001"),
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 141874cb2..0850efc88 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -113,6 +113,40 @@ test_that("plotly wrappers return plotly objects for spectra and rate plots", {
                                        include_critical = TRUE), "plotly")
 })
 
+test_that("plotSpectra2 compares spectra from params and sims", {
+    p_params <- plotSpectra2(params, params, name1 = "Original",
+                             name2 = "Changed", species = species,
+                             total = TRUE)
+    expect_s3_class(p_params, "ggplot")
+    expect_identical(levels(p_params$data$Model), c("Original", "Changed"))
+    expect_true("Total" %in% p_params$data$Legend)
+
+    expect_s3_class(plotSpectra2(sim, sim0, species = species), "ggplot")
+    expect_s3_class(plotSpectra2(params, sim, species = species), "ggplot")
+})
+
+test_that("plotSpectra2 supports base plot log argument", {
+    p_y <- plotSpectra2(params, sim0, species = species, log = "y")
+    expect_identical(p_y$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_y$scales$get_scales("y")$trans$name, "log-10")
+
+    p_xy <- plotSpectra2(params, sim0, species = species, log = "xy")
+    expect_identical(p_xy$scales$get_scales("x")$trans$name, "log-10")
+    expect_identical(p_xy$scales$get_scales("y")$trans$name, "log-10")
+
+    p_none <- plotSpectra2(params, sim0, species = species, log = "")
+    expect_identical(p_none$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_none$scales$get_scales("y")$trans$name, "identity")
+
+    p_flags <- plotSpectra2(params, sim0, species = species,
+                            log_x = FALSE, log_y = FALSE)
+    expect_identical(p_flags$scales$get_scales("x")$trans$name, "identity")
+    expect_identical(p_flags$scales$get_scales("y")$trans$name, "identity")
+
+    expect_error(plotSpectra2(params, sim0, species = species, log = "z"),
+                 "`log` must be a character string")
+})
+
 test_that("yield plotting helpers validate comparison and gear selection", {
     sim_shifted <- sim
     dimnames(sim_shifted@n)$time <- as.character(10:13)

From 6bd66aa5aff67a5b6743336f7ae136ef97340692 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 15:54:05 +0100
Subject: [PATCH 05/10] Add `plotRelative()` generic and methods for symmetric
 relative difference plotting between mizer array objects. Move
 `plotSpectraRelative()` and `plotlySpectraRelative()` from mizerExperimental
 to mizer.

---
 NAMESPACE                                     |   6 +
 NEWS.md                                       |   5 +
 R/ArraySpeciesBySize-class.R                  |  50 ++++++++
 R/ArrayTimeBySpecies-class.R                  |  25 ++++
 R/ArrayTimeBySpeciesBySize-class.R            |  24 ++++
 R/plots.R                                     | 118 +++++++++++++++++-
 tests/testthat/test-ArraySpeciesBySize.R      |  19 +++
 tests/testthat/test-ArrayTimeBySpecies.R      |  21 ++++
 .../testthat/test-ArrayTimeBySpeciesBySize.R  |  17 +++
 tests/testthat/test-plots.R                   |  22 ++++
 10 files changed, 306 insertions(+), 1 deletion(-)

diff --git a/NAMESPACE b/NAMESPACE
index 28db4e328..5a5257e39 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -141,6 +141,9 @@ S3method(plotGrowthCurves,MizerParams)
 S3method(plotGrowthCurves,MizerSim)
 S3method(plotPredMort,MizerParams)
 S3method(plotPredMort,MizerSim)
+S3method(plotRelative,ArraySpeciesBySize)
+S3method(plotRelative,ArrayTimeBySpecies)
+S3method(plotRelative,ArrayTimeBySpeciesBySize)
 S3method(plotSpectra,MizerParams)
 S3method(plotSpectra,MizerSim)
 S3method(plotYield,MizerSim)
@@ -409,8 +412,10 @@ export(plotFeedingLevel)
 export(plotGrowthCurves)
 export(plotM2)
 export(plotPredMort)
+export(plotRelative)
 export(plotSpectra)
 export(plotSpectra2)
+export(plotSpectraRelative)
 export(plotYield)
 export(plotYieldGear)
 export(plotYieldObservedVsModel)
@@ -421,6 +426,7 @@ export(plotlyFeedingLevel)
 export(plotlyGrowthCurves)
 export(plotlyPredMort)
 export(plotlySpectra)
+export(plotlySpectraRelative)
 export(plotlyYield)
 export(plotlyYieldGear)
 export(plotlyYieldObservedVsModel)
diff --git a/NEWS.md b/NEWS.md
index 2a6e15c1f..991277012 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -134,6 +134,11 @@ individual variability in growth to be modelled.
   linetype. The `plotSpectra2()` helper has moved from `mizerExperimental` into
   mizer for comparing two abundance spectra.
 
+- New `plotRelative()` generic with methods for plotting the symmetric relative
+  difference between two compatible mizer array objects. The
+  `plotSpectraRelative()` and `plotlySpectraRelative()` helpers have moved from
+  `mizerExperimental` into mizer.
+
 - `plotSpectra()` now accepts `log_x`, `log_y`, and `log` arguments for
   controlling axis scaling, matching the mizer array `plot()` methods.
 
diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R
index a78307fff..4e48f3e01 100644
--- a/R/ArraySpeciesBySize-class.R
+++ b/R/ArraySpeciesBySize-class.R
@@ -284,6 +284,56 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
                             y_ticks = y_ticks, legend_var = "Legend")
 }
 
+#' Plot the relative difference between two mizer array objects
+#'
+#' `plotRelative()` plots the difference between two compatible mizer array
+#' objects relative to their average. If the values in the first object are
+#' \eqn{N_1} and the values in the second are \eqn{N_2}, it plots
+#' \deqn{2 (N_2 - N_1) / (N_1 + N_2).}
+#'
+#' @inheritParams plot2
+#' @param log_x If `TRUE`, use a log10 x-axis. Default is `TRUE` for size
+#'   spectra and `FALSE` for time series.
+#'
+#' @return A ggplot2 object.
+#' @export
+#' @family plotting functions
+#'
+#' @examples
+#' \donttest{
+#' enc <- getEncounter(NS_params)
+#' plotRelative(enc, enc, species = "Cod")
+#' plotRelative(getBiomass(NS_sim), getBiomass(NS_sim), species = "Cod")
+#' }
+plotRelative <- function(x, y, ...) {
+    UseMethod("plotRelative", x)
+}
+
+#' @rdname plotRelative
+#' @export
+plotRelative.ArraySpeciesBySize <- function(x, y, species = NULL,
+                                            all.sizes = FALSE,
+                                            log_x = TRUE,
+                                            wlim = c(NA, NA),
+                                            ylim = c(NA, NA),
+                                            total = FALSE,
+                                            background = TRUE, ...) {
+    check_plot2_compatible(x, y, "ArraySpeciesBySize")
+    compare_array_metadata(x, y)
+    params <- attr(x, "params")
+    plot_dat1 <- prepare_ArraySpeciesBySize_plot_data(
+        x, species = species, all.sizes = all.sizes, wlim = wlim,
+        total = total, background = background)
+    plot_dat2 <- prepare_ArraySpeciesBySize_plot_data(
+        y, species = species, all.sizes = all.sizes, wlim = wlim,
+        total = total, background = background)
+
+    plotRelativeDataFrame(plot_dat1, plot_dat2, params,
+                          xlab = "Size [g]",
+                          xtrans = if (log_x) "log10" else "identity",
+                          xlim = wlim, ylim = ylim, legend_var = "Legend")
+}
+
 check_plot2_compatible <- function(x, y, class) {
     if (!inherits(y, class)) {
         stop("Both objects must be of class `", class, "`.")
diff --git a/R/ArrayTimeBySpecies-class.R b/R/ArrayTimeBySpecies-class.R
index b206005b2..5a5e8884a 100644
--- a/R/ArrayTimeBySpecies-class.R
+++ b/R/ArrayTimeBySpecies-class.R
@@ -267,6 +267,31 @@ plot2.ArrayTimeBySpecies <- function(x, y, name1 = "First", name2 = "Second",
                             legend_var = "Legend")
 }
 
+#' @rdname plotRelative
+#' @export
+plotRelative.ArrayTimeBySpecies <- function(x, y, species = NULL,
+                                            start_time = NULL,
+                                            end_time = NULL,
+                                            ylim = c(NA, NA),
+                                            total = FALSE,
+                                            background = TRUE,
+                                            log_x = FALSE, ...) {
+    check_plot2_compatible(x, y, "ArrayTimeBySpecies")
+    compare_array_metadata(x, y)
+    params <- attr(x, "params")
+    plot_dat1 <- prepare_ArrayTimeBySpecies_plot_data(
+        x, species = species, start_time = start_time, end_time = end_time,
+        total = total, background = background)
+    plot_dat2 <- prepare_ArrayTimeBySpecies_plot_data(
+        y, species = species, start_time = start_time, end_time = end_time,
+        total = total, background = background)
+
+    plotRelativeDataFrame(plot_dat1, plot_dat2, params,
+                          xlab = "Year",
+                          xtrans = if (log_x) "log10" else "identity",
+                          ylim = ylim, legend_var = "Legend")
+}
+
 prepare_ArrayTimeBySpecies_plot_data <- function(x, species = NULL,
                                                  start_time = NULL,
                                                  end_time = NULL,
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index 0c0e04640..a326abee6 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -198,6 +198,30 @@ plot2.ArrayTimeBySpeciesBySize <- function(x, y, name1 = "First",
                              background = background, y_ticks = y_ticks, ...)
 }
 
+#' @rdname plotRelative
+#'
+#' @param time The time to display. Default (`NULL`) is the final time step.
+#'   Only applies to `ArrayTimeBySpeciesBySize`.
+#' @export
+plotRelative.ArrayTimeBySpeciesBySize <- function(x, y, species = NULL,
+                                                  time = NULL,
+                                                  all.sizes = FALSE,
+                                                  log_x = TRUE,
+                                                  wlim = c(NA, NA),
+                                                  ylim = c(NA, NA),
+                                                  total = FALSE,
+                                                  background = TRUE, ...) {
+    check_plot2_compatible(x, y, "ArrayTimeBySpeciesBySize")
+    slice1 <- ArrayTimeBySpeciesBySize_slice(x, time = time)
+    slice2 <- ArrayTimeBySpeciesBySize_slice(y, time = time)
+
+    plotRelative.ArraySpeciesBySize(slice1, slice2, species = species,
+                                    all.sizes = all.sizes, log_x = log_x,
+                                    wlim = wlim, ylim = ylim,
+                                    total = total, background = background,
+                                    ...)
+}
+
 ArrayTimeBySpeciesBySize_slice <- function(x, time = NULL) {
     params <- attr(x, "params")
     value_name <- attr(x, "value_name")
diff --git a/R/plots.R b/R/plots.R
index 53d0092c8..2069656c1 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -130,7 +130,7 @@ utils::globalVariables(c("time", "value", "Species", "w", "gear", "Age",
                          "x", "y", "Year", "Yield", "Biomass", "Size",
                          "Proportion", "Prey", "Legend", "Type", "Gear",
                          "Predator", "weight", "a", "b", "age", "w_max",
-                         "Model"))
+                         "Model", "rel_diff"))
 
 #' Make a plot from a data frame
 #'
@@ -304,6 +304,59 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
         scale_linetype_discrete(drop = FALSE)
 }
 
+plotRelativeDataFrame <- function(frame1, frame2, params,
+                                  xlab = waiver(),
+                                  xtrans = "identity",
+                                  xlim = c(NA, NA),
+                                  ylim = c(NA, NA),
+                                  legend_var = "Legend") {
+    assert_that(is.data.frame(frame1),
+                is.data.frame(frame2),
+                is(params, "MizerParams"))
+
+    names(frame2)[seq_len(min(3, ncol(frame2)))] <-
+        names(frame1)[seq_len(min(3, ncol(frame1)))]
+    var_names <- names(frame1)
+    x_var <- var_names[[1]]
+    y_var <- var_names[[2]]
+    group_var <- var_names[[3]]
+    if (!(legend_var %in% var_names)) {
+        stop("The `legend_var` argument must be the name of a variable ",
+             "in the data frame.")
+    }
+
+    by_vars <- c(x_var, group_var, legend_var)
+    frame <- dplyr::inner_join(frame1, frame2, by = by_vars,
+                               suffix = c(".x", ".y"))
+    frame$rel_diff <- relative_difference(frame[[paste0(y_var, ".x")]],
+                                          frame[[paste0(y_var, ".y")]])
+    frame <- frame[is.finite(frame$rel_diff), ]
+
+    legend_levels <- intersect(names(params@linecolour), frame[[legend_var]])
+    frame[[legend_var]] <- factor(frame[[legend_var]], levels = legend_levels)
+    if (sum(is.na(frame[[legend_var]]))) {
+        warning("missing legend in params@linecolour, some groups won't be displayed")
+    }
+    linecolour <- params@linecolour[legend_levels]
+
+    xbreaks <- waiver()
+    if (xtrans == "log10") xbreaks <- log_breaks()
+
+    ggplot(frame, aes(group = .data[[group_var]])) +
+        scale_y_continuous(name = "Relative difference", limits = ylim) +
+        scale_x_continuous(trans = xtrans, breaks = xbreaks, name = xlab,
+                           limits = xlim) +
+        geom_hline(yintercept = 0, linetype = 1,
+                   colour = "dark grey", linewidth = 0.75) +
+        geom_line(aes(x = .data[[x_var]], y = .data[["rel_diff"]],
+                      colour = .data[[legend_var]])) +
+        scale_colour_manual(values = linecolour)
+}
+
+relative_difference <- function(first, second) {
+    2 * (second - first) / (first + second)
+}
+
 #' Helper function to produce nice breaks on logarithmic axes
 #'
 #' This is needed when the logarithmic y-axis spans less than one order of
@@ -947,6 +1000,69 @@ spectra_y_label <- function(power) {
     paste0("Number density * w^", power)
 }
 
+#' Plot the relative difference between two spectra
+#'
+#' `plotSpectraRelative()` plots the difference between the spectra relative to
+#' their average. If we denote the number density from the first object as
+#' \eqn{N_1(w)} and that from the second object as \eqn{N_2(w)}, then this plot
+#' shows
+#' \deqn{2 (N_2(w) - N_1(w)) / (N_2(w) + N_1(w)).}
+#'
+#' The individual spectra are calculated by [plotSpectra()], to which all
+#' additional arguments are passed. For example, you can determine a time range
+#' over which to average simulation results via `time_range`. See
+#' [plotSpectra()] for more options.
+#'
+#' Note that it does not matter whether the relative difference is calculated
+#' for number density, biomass density, or biomass density in log weight,
+#' because the factors of \eqn{w} by which the densities differ cancel out in
+#' the relative difference.
+#'
+#' @param object1 First `MizerParams` or `MizerSim` object.
+#' @param object2 Second `MizerParams` or `MizerSim` object.
+#' @param log_x If `TRUE` (default), use a log10 x-axis.
+#' @param ylim A numeric vector of length two providing lower and upper limits
+#'   for the relative difference (y) axis. Use `NA` to refer to the existing
+#'   minimum or maximum.
+#' @param ... Arguments passed to [plotSpectra()] for preparing the spectra
+#'   data, for example `species`, `time_range`, `wlim`, `resource`,
+#'   `background` or `total`.
+#'
+#' @return A ggplot2 object.
+#' @export
+#' @family plotting functions
+#'
+#' @examples
+#' \donttest{
+#' sim1 <- project(NS_params, t_max = 10, progress_bar = FALSE)
+#' sim2 <- project(NS_params, effort = 0.5, t_max = 10, progress_bar = FALSE)
+#' plotSpectraRelative(sim1, sim2)
+#' }
+plotSpectraRelative <- function(object1, object2, log_x = TRUE,
+                                ylim = c(NA, NA), ...) {
+    args <- list(...)
+    wlim <- args$wlim %||% c(NA, NA)
+
+    sf1 <- plotSpectra(object1, return_data = TRUE, ...)
+    sf2 <- plotSpectra(object2, return_data = TRUE, ...)
+    params <- if (is(object1, "MizerSim")) object1@params else object1
+
+    plotRelativeDataFrame(sf1, sf2, validParams(params),
+                          xlab = "Size [g]",
+                          xtrans = if (log_x) "log10" else "identity",
+                          xlim = wlim, ylim = ylim,
+                          legend_var = "Legend")
+}
+
+#' @rdname plotSpectraRelative
+#' @export
+plotlySpectraRelative <- function(object1, object2, log_x = TRUE,
+                                  ylim = c(NA, NA), ...) {
+    ggplotly(plotSpectraRelative(object1, object2, log_x = log_x,
+                                  ylim = ylim, ...),
+             tooltip = c("Legend", "w", "rel_diff"))
+}
+
 #' @rdname plotSpectra
 #' @export
 plotlySpectra <- function(object, species = NULL,
diff --git a/tests/testthat/test-ArraySpeciesBySize.R b/tests/testthat/test-ArraySpeciesBySize.R
index b53d6d17a..4e779144a 100644
--- a/tests/testthat/test-ArraySpeciesBySize.R
+++ b/tests/testthat/test-ArraySpeciesBySize.R
@@ -129,6 +129,25 @@ test_that("plot2.ArraySpeciesBySize compares compatible arrays", {
     expect_error(plot2(enc, enc, log = "z"), "`log` must be a character string")
 })
 
+test_that("plotRelative.ArraySpeciesBySize plots symmetric relative difference", {
+    enc <- getEncounter(NS_params)
+    enc2 <- enc
+    enc2[] <- unclass(enc) * 2
+
+    p <- plotRelative(enc, enc2, species = "Cod", total = TRUE,
+                      background = FALSE, wlim = c(1, NA))
+    expect_s3_class(p, "ggplot")
+    expect_true(all(p$data$Species %in% c("Cod", "Total")))
+    expect_true(all(p$data$w >= 1))
+    expect_true(all(abs(p$data$rel_diff - 2 / 3) < 1e-12))
+    expect_identical(p$scales$get_scales("x")$trans$name, "log-10")
+
+    p_linear <- plotRelative(enc, enc2, species = "Cod", log_x = FALSE)
+    expect_identical(p_linear$scales$get_scales("x")$trans$name, "identity")
+
+    expect_error(plotRelative(enc, getBiomass(NS_sim)), "Both objects must be")
+})
+
 test_that("addPlot.ArraySpeciesBySize adds lines to an existing ggplot", {
     enc <- getEncounter(NS_params)
     pred_mort <- getPredMort(NS_params)
diff --git a/tests/testthat/test-ArrayTimeBySpecies.R b/tests/testthat/test-ArrayTimeBySpecies.R
index 08924e656..8361258dd 100644
--- a/tests/testthat/test-ArrayTimeBySpecies.R
+++ b/tests/testthat/test-ArrayTimeBySpecies.R
@@ -110,6 +110,27 @@ test_that("plot2.ArrayTimeBySpecies compares compatible arrays", {
     expect_error(plot2(bio, getEncounter(NS_params)), "Both objects must be")
 })
 
+test_that("plotRelative.ArrayTimeBySpecies plots symmetric relative difference", {
+    bio <- getBiomass(NS_sim)
+    bio2 <- bio
+    bio2[] <- unclass(bio) * 2
+    years <- as.numeric(rownames(bio))
+
+    p <- plotRelative(bio, bio2, species = "Cod", total = TRUE,
+                      start_time = years[2], end_time = years[5])
+    expect_s3_class(p, "ggplot")
+    expect_true(all(p$data$Species %in% c("Cod", "Total")))
+    expect_true(all(p$data$Year >= years[2]))
+    expect_true(all(p$data$Year <= years[5]))
+    expect_true(all(abs(p$data$rel_diff - 2 / 3) < 1e-12))
+    expect_identical(p$scales$get_scales("x")$trans$name, "identity")
+
+    p_log <- plotRelative(bio, bio2, species = "Cod", log_x = TRUE)
+    expect_identical(p_log$scales$get_scales("x")$trans$name, "log-10")
+
+    expect_error(plotRelative(bio, getEncounter(NS_params)), "Both objects must be")
+})
+
 test_that("addPlot.ArrayTimeBySpecies adds lines to an existing ggplot", {
     bio <- getBiomass(NS_sim)
     yield <- getYield(NS_sim)
diff --git a/tests/testthat/test-ArrayTimeBySpeciesBySize.R b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
index 680a28602..44a242e4c 100644
--- a/tests/testthat/test-ArrayTimeBySpeciesBySize.R
+++ b/tests/testthat/test-ArrayTimeBySpeciesBySize.R
@@ -103,6 +103,23 @@ test_that("plot2.ArrayTimeBySpeciesBySize compares selected time slices", {
     expect_error(plot2(fmort, getBiomass(NS_sim)), "Both objects must be")
 })
 
+test_that("plotRelative.ArrayTimeBySpeciesBySize compares selected time slices", {
+    fmort <- getFMort(NS_sim)
+    fmort2 <- fmort
+    fmort2[] <- unclass(fmort) * 2
+    times <- as.numeric(dimnames(fmort)[[1]])
+
+    p <- plotRelative(fmort, fmort2, species = "Cod", time = times[5],
+                      total = TRUE, wlim = c(1, NA))
+    expect_s3_class(p, "ggplot")
+    expect_true(all(p$data$Species %in% c("Cod", "Total")))
+    expect_true(all(p$data$w >= 1))
+    expect_true(all(abs(p$data$rel_diff - 2 / 3) < 1e-12))
+    expect_identical(p$scales$get_scales("x")$trans$name, "log-10")
+
+    expect_error(plotRelative(fmort, getBiomass(NS_sim)), "Both objects must be")
+})
+
 test_that("plot.ArrayTimeBySpeciesBySize preserves single species dimension", {
     arr <- array(seq_len(6), dim = c(2, 1, 3),
                  dimnames = list(time = c("2000", "2001"),
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 0850efc88..c31527d33 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -106,6 +106,8 @@ test_that("plotly functions do not throw error", {
 
 test_that("plotly wrappers return plotly objects for spectra and rate plots", {
     expect_s3_class(plotlySpectra(params, species = species), "plotly")
+    expect_s3_class(plotlySpectraRelative(params, params, species = species,
+                                          resource = FALSE), "plotly")
     expect_s3_class(plotlyPredMort(sim, species = species), "plotly")
     expect_s3_class(plotlyFMort(sim, species = species), "plotly")
     expect_s3_class(plotlyGrowthCurves(sim, species = species), "plotly")
@@ -147,6 +149,26 @@ test_that("plotSpectra2 supports base plot log argument", {
                  "`log` must be a character string")
 })
 
+test_that("plotSpectraRelative plots symmetric relative difference", {
+    params2 <- params
+    params2@initial_n[] <- params@initial_n * 2
+
+    p <- plotSpectraRelative(params, params2, species = species,
+                             resource = FALSE)
+    expect_s3_class(p, "ggplot")
+    expect_true(all(abs(p$data$rel_diff - 2 / 3) < 1e-12))
+    expect_identical(p$scales$get_scales("x")$trans$name, "log-10")
+
+    p_linear <- plotSpectraRelative(params, params2, species = species,
+                                    resource = FALSE, log_x = FALSE)
+    expect_identical(p_linear$scales$get_scales("x")$trans$name, "identity")
+
+    expect_s3_class(plotSpectraRelative(sim, sim0, species = species,
+                                        resource = FALSE), "ggplot")
+    expect_s3_class(plotSpectraRelative(params, sim, species = species,
+                                        resource = FALSE), "ggplot")
+})
+
 test_that("yield plotting helpers validate comparison and gear selection", {
     sim_shifted <- sim
     dimnames(sim_shifted@n)$time <- as.character(10:13)

From 989c8907fa2baf4ef1e9fd7c3a27d8a4ce2dd145 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 15:54:36 +0100
Subject: [PATCH 06/10] Updated vignette about cohort dynamics

---
 vignettes/cohort_dynamics_and_diffusion.Rmd | 283 ++++++--------------
 1 file changed, 80 insertions(+), 203 deletions(-)

diff --git a/vignettes/cohort_dynamics_and_diffusion.Rmd b/vignettes/cohort_dynamics_and_diffusion.Rmd
index 83459747c..bbe77dd02 100644
--- a/vignettes/cohort_dynamics_and_diffusion.Rmd
+++ b/vignettes/cohort_dynamics_and_diffusion.Rmd
@@ -43,11 +43,11 @@ library(plotly)
 
 We start by creating a single-species model using
 `newSingleSpeciesParams()`. This sets up a species embedded in a
-power-law background community.
+power-law background community, so that the encounter rate scales as $w^{3/4}$
+and the mortality rate scales as $w^{-1/4}.
 
 ```{r}
 params <- newSingleSpeciesParams(h = 10, no_w = 400)
-params <- steady(params)
 ```
 
 # Pulsed reproduction
@@ -58,15 +58,6 @@ custom density-dependent reproduction rate function (RDD function) that
 only allows reproduction during a brief window at the start of each
 year.
 
-First, we calculate the steady-state reproduction rate. This tells us
-the total egg production rate needed to maintain the population. We will
-use this value as the magnitude of our annual pulse.
-
-```{r}
-rdd_steady <- getRDD(params)
-cat("Steady-state RDD:", rdd_steady, "\n")
-```
-
 The RDD function receives the current time `t` as an argument. We use
 this to turn reproduction on only during a short window at the start of
 each year and off at all other times. To maintain the same total annual
@@ -80,22 +71,14 @@ pulse_width <- 0.1  # Reproduce during first 10% of each year
 annual_pulse_RDD <- function(rdi, species_params, t, ...) {
     frac <- t %% 1
     if (frac < pulse_width) {
-        rdd <- 1 - cos(frac / pulse_width * 2 * pi) * species_params$rdd_steady
-        # Scale up to maintain total annual reproduction
-        return(rdd / pulse_width)
+        rdd <- rdi # to get vector of right length
+        rdd[] <- 1 - cos(frac / pulse_width * 2 * pi)
+        return(rdd)
     } else {
         return(0 * rdi)
     }
 }
-```
-
-We store the steady-state RDD value in the species parameters so our
-function can access it, and then register the function:
-
-```{r}
-params_pulse <- params
-species_params(params_pulse)$rdd_steady <- rdd_steady
-params_pulse <- setRateFunction(params_pulse, "RDD", "annual_pulse_RDD")
+params <- setRateFunction(params, "RDD", "annual_pulse_RDD")
 ```
 
 Of course this is not a realistic way to model seasonal reproduction. A more
@@ -110,221 +93,112 @@ We start from an empty spectrum (no fish) and let the pulsed
 reproduction create cohorts from scratch.
 
 ```{r}
-params_empty <- params_pulse
-initialN(params_empty)[] <- 0
-
-sim_no_diff <- project(params_empty, t_max = 5, dt = 0.01,
-                       t_save = 0.1, progress_bar = FALSE,
-                       method = "predictor-corrector")
-
-animate(sim_no_diff, log_x = TRUE, log_y = FALSE, power = 2,
-        time_range = c(0, 5), resource = FALSE, interpolate = FALSE)
-```
-
-Let's visualise the size spectrum at different time points to see the
-cohorts:
-
-```{r fig.height=6}
-plot_cohort_evolution <- function(sim, times_to_plot, species = 1,
-                                  title = "Cohort evolution",
-                                  y_label = "Biomass density [g]",
-                                  interactive = TRUE) {
-    w <- sim@params@w
-    sim_times <- as.numeric(dimnames(sim@n)$time)
-    plot_list <- list()
-    labels <- character()
-
-    for (tt in times_to_plot) {
-        idx <- which.min(abs(sim_times - tt))
-        actual_time <- sim_times[idx]
-        n_at_t <- as.numeric(sim@n[idx, species, ])
-        pos <- n_at_t > 0
-        if (any(pos)) {
-            label <- paste0("t = ", actual_time)
-            labels <- c(labels, label)
-            plot_list[[length(plot_list) + 1]] <- data.frame(
-                w = w[pos],
-                n = n_at_t[pos] * w[pos]^2,
-                time = label
-            )
-        }
-    }
-
-    if (length(plot_list) > 0) {
-        plot_data <- do.call(rbind, plot_list)
-        plot_data$time <- factor(plot_data$time, levels = unique(labels))
-    } else {
-        plot_data <- data.frame(
-            w = numeric(),
-            n = numeric(),
-            time = factor(character())
-        )
-    }
-
-    p <- ggplot(plot_data, aes(x = w, y = n, colour = time)) +
-        geom_line(linewidth = 0.8) +
-        scale_x_log10(limits = c(1e-3, 100)) +
-        # scale_y_log10() +
-        labs(x = "Weight [g]", y = y_label,
-             title = title,
-             colour = "Time") +
-        theme_minimal(base_size = 14)
-
-    if (interactive) {
-        return(ggplotly(p))
-    } else {
-        return(p)
-    }
-}
-
-plot_cohort_evolution(
-    sim_no_diff,
-    times_to_plot = c(0.5, 1, 1.5, 2, 3, 4),
-    title = "Cohort evolution without diffusion"
-)
-```
-
-```{r}
+initialN(params)[] <- 0
 
-sim_corr <- project(params_empty, t_max = 5, dt = 0.02,
-                       t_save = 0.1, progress_bar = FALSE,
+sim_no_diff <- project(params, t_max = 5, dt = 0.01,
+                       t_save = 0.05, progress_bar = FALSE,
                        method = "predictor-corrector")
 
-plot_cohort_evolution(
-    sim_corr,
-    times_to_plot = c(0.5, 1, 1.5, 2, 3, 4),
-    title = "Cohort evolution without diffusion"
-)
-```
-
-```{r}
-animate(sim_corr, log_x = TRUE, log_y = FALSE, power = 2,
-        time_range = c(0.5, 5), resource = FALSE)
+animate(sim_no_diff, log_x = TRUE, log_y = FALSE, power = 2, resource = FALSE,
+        transition_duration = 0, frame_duration = 200)
 ```
 
 
-Without diffusion, each cohort appears as a relatively sharp peak that
-moves to the right (towards larger sizes) as the fish grow.
+Without diffusion, each cohort peak 
+moves to the right (towards larger sizes) as the fish grow. The width of the cohort stays constant on the logarithmic axis. The growth on
+the logarithmic axis slows down towards larger sizes and so eventually
+the older cohorts start to merge.
 
 # Adding predation diffusion
 
+As discussed elsewhere, some of the randomness in the growth rate comes from the randomness in the size of prey
+encountered by the predator. We refer to this as the predation diffusion, even though there are other sources
+of randomness associated with predation, arising for example from the patchiness in the spatial distribution of
+prey. So the predation diffusion can be seen as a lower bound on the amount of diffusion arising from randomness
+in growth.
+
 ```{r}
-params_pred_diff <- params_empty
+params_pred_diff <- params
 use_predation_diffusion(params_pred_diff) <- TRUE
-sim_pred_diff <- project(params_pred_diff, t_max = 5, dt = 0.02,
-                         t_save = 0.1, progress_bar = FALSE,
+sim_pred_diff <- project(params_pred_diff, t_max = 5, dt = 0.01,
+                         t_save = 0.05, progress_bar = FALSE,
                          method = "predictor-corrector")
 
+animate(sim_pred_diff, log_x = TRUE, log_y = FALSE, power = 2, resource = FALSE,
+        transition_duration = 0, frame_duration = 200)
+
 ```
 
-```{r}
-animate(sim_pred_diff, log_x = TRUE, log_y = FALSE, power = 2,
-        time_range = c(0.5, 5), resource = FALSE)
+We see that there is a slight broadening of the cohorts as they grow up, which is most noticeable at larger
+sizes where they merge together sooner. To see this more clearly we plot the biomass density at
+time $t=5$ for both cases:
 
+```{r}
+plotSpectra2(sim_no_diff, sim_pred_diff,
+             name1 = "No diffusion", name2 = "Predation diffusion",
+             resource = FALSE, power = 2, log = "x")
 ```
 
+The diffusion rate is a power law in $w$ with exponent $7/4$. The coefficient is
 ```{r}
-plot_cohort_evolution(
-    sim_pred_diff,
-    times_to_plot = c(0.5, 1, 1.5, 2, 3, 4),
-    title = "Cohort evolution with predation diffusion"
-)
+(getDiffusion(params_pred_diff) / w(params)^(7/4))[1]
 ```
 
-
 # Adding external diffusion
 
-Now let's add diffusion to the model. Diffusion is set as an array with
-dimensions species × size via `setExtDiffusion()`. We'll set a constant
-diffusion rate across all sizes.
-
-We define a helper function that runs the simulation for a given
-diffusion coefficient:
+The predation diffusion is the only part of the diffusion that is explicitly modelled in mizer. We refer to all non-modelled diffusion as "external" to the model. We assume that it follows the same power law but with an a priory unknown coefficient that will need to be determined by looking at the rate at which cohort size distributions widen
+in the real world. The coefficient of external diffusion is set with the species parameter `D_ext`.
 
-```{r}
-run_with_diffusion <- function(params_base, diff_coeff, diff_exp, t_max = 5) {
-    p <- params_base
-    w <- p@w
-    d <- p@ext_diffusion
-    d[] <- diff_coeff * w ^ diff_exp
-    p <- setExtDiffusion(p, ext_diffusion = d)
-    initialN(p)[] <- 0
-    sim <- project(p, t_max = t_max, dt = 0.02,
-                   t_save = 0.1, progress_bar = FALSE,
-                   method = "predictor-corrector")
-    return(sim)
-}
-```
 
 # Comparing different diffusion rates
 
-Let's compare simulations with no diffusion, low diffusion, and high
-diffusion:
+We now run the model with two different levels of external diffusion in order to
+compare the effects.
 
 ```{r}
-diff_exp <- params_pulse@species_params$n + 1
-sim_d0 <- run_with_diffusion(params_empty, diff_coeff = 0, diff_exp = diff_exp)
-sim_d_low <- run_with_diffusion(params_empty, diff_coeff = 0.1, diff_exp = diff_exp)
-sim_d_high <- run_with_diffusion(params_empty, diff_coeff = 0.5, diff_exp = diff_exp)
+species_params(params)$D_ext <- 0.1
+sim_medium_diff <- project(params, t_max = 5, dt = 0.01,
+                      t_save = 0.05, progress_bar = FALSE,
+                      method = "predictor-corrector")
+
+species_params(params)$D_ext <- 0.5
+sim_high_diff <- project(params, t_max = 5, dt = 0.01,
+                     t_save = 0.05, progress_bar = FALSE,
+                     method = "predictor-corrector")
 ```
 
-Now let's visualise the cohorts at a specific time point to see how
-diffusion affects their shape:
+Let's compare simulations with no diffusion, only predation diffusion, a medium level of diffusion and a level of diffusion at time $t=5$:
 
 ```{r fig.height=5, fig.width=9}
-snapshot_time <- 3
-
-build_snapshot <- function(sim, label) {
-    idx <- which.min(abs(as.numeric(dimnames(sim@n)$time) - snapshot_time))
-    n_at_t <- as.numeric(sim@n[idx, 1, ])
-    w <- sim@params@w
-    pos <- n_at_t > 0
-    if (any(pos)) {
-        data.frame(w = w[pos], n = n_at_t[pos], diffusion = label)
-    } else {
-        data.frame(w = numeric(0), n = numeric(0), diffusion = character(0))
-    }
-}
-
+w <- w(params)
 snapshot_data <- rbind(
-    build_snapshot(sim_d0, "D = 0 (no diffusion)"),
-    build_snapshot(sim_d_low, "D = 0.1 (low)"),
-    build_snapshot(sim_d_high, "D = 0.5 (high)"),
-    build_snapshot(sim_pred_diff, "Predation diffusion")
+    data.frame(x = w, y = finalN(sim_no_diff)[1, ] * w^2,
+               label = "D = 0 (no diffusion)"),
+    data.frame(x = w, y = finalN(sim_pred_diff)[1, ] * w^2,
+               label = "D = 0.02 (predation)"),
+    data.frame(x = w, y = finalN(sim_medium_diff)[1, ] * w^2,
+               label = "D = 0.1 (medium)"),
+    data.frame(x = w, y = finalN(sim_high_diff)[1, ] * w^2,
+               label = "D = 0.5 (high)")
 )
 
-p <- ggplot(snapshot_data, aes(x = w, y = n * w^2, colour = diffusion)) +
-    geom_line(linewidth = 0.8) +
-    scale_x_log10(limits = c(1e-3, 100)) +
-    # scale_y_log10() +
-    labs(x = "Weight [g]", y = "Biomass density [g]",
-         title = paste0("Effect of diffusion on cohorts at t = ",
-                        snapshot_time),
-         colour = "Diffusion rate") +
-    theme_minimal(base_size = 14)
+p <- ggplot(snapshot_data, aes(x = x, y = y, colour = label)) +
+     geom_line(linewidth = 0.8) +
+     scale_x_log10(limits = c(1e-3, 100)) +
+     labs(x = "Weight [g]", y = "Biomass density [g]",
+          title = "Effect of diffusion on cohorts at t = 5",
+          colour = "Diffusion rate") +
+     theme_minimal(base_size = 14)
 
 ggplotly(p)
 ```
 
-We can see that diffusion has a large effect on the speed at which the cohort
-peaks are moving but not so much effect on the broadening of the peaks.
-
-# Time evolution with diffusion
-
-Let's look at the full time evolution of the size spectrum with
-moderate diffusion to see how cohorts spread over time:
+We can see that diffusion has a large effect on how quickly the cohorts widen and
+merge into each other. Let us look at an animation of the high-diffusion case:
 
-```{r fig.height=6}
-diff_exp <- params_pulse@species_params$n + 1
-sim_d_med <- run_with_diffusion(params_pulse, diff_coeff = 0.05, diff_exp = diff_exp)
-
-plot_cohort_evolution(
-    sim_d_med,
-    times_to_plot = c(0.5, 1, 2, 3, 4, 5),
-    title = "Cohort evolution with moderate diffusion (D = 0.05)",
-    y_label = "Number density [1/g]",
-    interactive = FALSE
-)
+```{r}
+animate(sim_high_diff, log_x = TRUE, log_y = FALSE, power = 2, resource = FALSE,
+        transition_duration = 0, frame_duration = 200)
 ```
 
 As time progresses, we see that:
@@ -339,15 +213,18 @@ As time progresses, we see that:
 # Heatmap visualisation
 
 A heatmap provides a compact view of the entire dynamics, showing how
-the size spectrum evolves continuously over time:
+the size spectrum evolves continuously over time. Here we look at the case 
+of high diffusion:
 
 ```{r fig.height=5, fig.width=9}
-all_times <- as.numeric(dimnames(sim_d_med@n)$time)
-w <- sim_d_med@params@w
+sim <- sim_high_diff
+all_times <- getTimes(sim)
+w <- params@w
+n <- N(sim)
 
 heatmap_list <- list()
 for (i in seq_along(all_times)) {
-    n_at_t <- as.numeric(sim_d_med@n[i, 1, ])
+    n_at_t <- as.numeric(n[i, 1, ])
     pos <- n_at_t > 0
     if (any(pos)) {
         heatmap_list[[length(heatmap_list) + 1]] <- data.frame(
@@ -364,12 +241,12 @@ ggplot(heatmap_data, aes(x = time, y = w, fill = log_n)) +
     scale_y_log10() +
     scale_fill_viridis_c(name = expression(log[10](N))) +
     labs(x = "Time [years]", y = "Weight [g]",
-         title = "Size spectrum over time (D = 0.05)") +
+         title = "Biomass density over time (D = 0.5)") +
     theme_minimal(base_size = 14)
 ```
 
 In the heatmap, the diagonal bands represent individual cohorts growing
-through the size spectrum. The broadening of these bands with time is
+through the size spectrum. The curving of these bands is due to the slowing down and the broadening of these bands with time is
 the effect of diffusion.
 
 # Summary

From f25589f004d185e1a4e147c7c0875bf4b43e0c60 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 17:45:49 +0100
Subject: [PATCH 07/10] Add `plotCDF()` generic and methods for plotting
 cumulative distributions in `MizerParams` and `MizerSim`.

---
 NAMESPACE                   |   3 +
 NEWS.md                     |   3 +
 R/plots.R                   | 173 +++++++++++++++++++++++++++++++++++-
 tests/testthat/test-plots.R |  62 +++++++++++++
 4 files changed, 239 insertions(+), 2 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 5a5257e39..255d132ad 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -131,6 +131,8 @@ S3method(plot2,ArrayTimeBySpeciesBySize)
 S3method(plotBiomass,MizerSim)
 S3method(plotBiomassObservedVsModel,MizerParams)
 S3method(plotBiomassObservedVsModel,MizerSim)
+S3method(plotCDF,MizerParams)
+S3method(plotCDF,MizerSim)
 S3method(plotDiet,MizerParams)
 S3method(plotDiet,MizerSim)
 S3method(plotFMort,MizerParams)
@@ -405,6 +407,7 @@ export(other_params)
 export(plot2)
 export(plotBiomass)
 export(plotBiomassObservedVsModel)
+export(plotCDF)
 export(plotDataFrame)
 export(plotDiet)
 export(plotFMort)
diff --git a/NEWS.md b/NEWS.md
index 991277012..5fd73b7ee 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -139,6 +139,9 @@ individual variability in growth to be modelled.
   `plotSpectraRelative()` and `plotlySpectraRelative()` helpers have moved from
   `mizerExperimental` into mizer.
 
+- New `plotCDF()` generic for plotting cumulative abundance or biomass
+  distributions from `MizerParams` and `MizerSim` objects.
+
 - `plotSpectra()` now accepts `log_x`, `log_y`, and `log` arguments for
   controlling axis scaling, matching the mizer array `plot()` methods.
 
diff --git a/R/plots.R b/R/plots.R
index 2069656c1..4ba5c11f0 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -42,8 +42,9 @@
 #' versions, for example `ggplotly(getBiomass(sim))` or
 #' `ggplotly(getEncounter(params))`. To add another compatible array to an
 #' existing ggplot, use [addPlot()]. To compare two compatible mizer arrays
-#' directly, use [plot2()]. To visualise how spectra or rates
-#' change through time, use [animate()] on a `MizerSim` or an
+#' directly, use [plot2()]. To plot cumulative distributions over body size,
+#' use [plotCDF()]. To visualise how spectra or rates change through time, use
+#' [animate()] on a `MizerSim` or an
 #' `ArrayTimeBySpeciesBySize` object.
 #'
 #' The named plotting functions give more specialised control. This table shows
@@ -54,6 +55,7 @@
 #'   [plotYield()] \tab Plots the total yield of each species across all fishing gears against time. \cr
 #'   [plotYieldGear()] \tab Plots the total yield of each species by gear against time. \cr
 #'   [plotSpectra()] \tab Plots the abundance (biomass or numbers) spectra of each species and the background community. It is possible to specify a minimum size which is useful for truncating the plot. \cr
+#'   [plotCDF()] \tab Plots cumulative distributions of abundance or biomass over size. \cr
 #'   [plotSpectra2()] \tab Compares the spectra from two simulations or parameter objects in one plot. \cr
 #'   [plotFeedingLevel()] \tab Plots the feeding level of each species against size. \cr
 #'   [plotPredMort()] \tab Plots the predation mortality of each species against size. \cr
@@ -938,6 +940,173 @@ plot_spectra <- function(params, n, n_pp,
                   highlight = highlight, legend_var = "Legend")
 }
 
+#' Plot cumulative abundance or biomass distributions
+#'
+#' `plotCDF()` plots the cumulative distribution over body size from small to
+#' large sizes. It uses the same spectra data preparation as [plotSpectra()].
+#' The density is first multiplied by `w^power`, then integrated over size.
+#' With `normalise = TRUE`, each curve is divided by its final value so that it
+#' ends at 1.
+#'
+#' @inheritParams plotSpectra
+#' @param normalise If `TRUE` (default), plot the cumulative proportion. If
+#'   `FALSE`, plot the cumulative abundance, biomass, or other unnormalised
+#'   integral.
+#' @param log_x If `TRUE` (default), use a log10 x-axis.
+#' @param log Character string specifying whether the x-axis should use a log10
+#'   scale, in the same form as the base [plot()] argument. For `plotCDF()`,
+#'   only `"x"` and `""` are supported. If supplied, this overrides `log_x`.
+#'
+#' @return A ggplot2 object, unless `return_data = TRUE`, in which case a data
+#'   frame with the four variables 'w', 'value', 'Species', 'Legend' is
+#'   returned.
+#' @export
+#' @family plotting functions
+#' @seealso [plotSpectra()]
+#' @examples
+#' \donttest{
+#' plotCDF(NS_params, species = c("Cod", "Herring"))
+#' plotCDF(NS_sim, power = 0, normalise = FALSE)
+#' }
+plotCDF <- function(object, ...) {
+    UseMethod("plotCDF")
+}
+
+#' @rdname plotCDF
+#' @export
+plotCDF.MizerSim <- function(object, species = NULL,
+                             time_range,
+                             geometric_mean = FALSE,
+                             wlim = c(NA, NA), ylim = c(NA, NA),
+                             power = 1, biomass = TRUE,
+                             total = FALSE, resource = TRUE,
+                             background = TRUE,
+                             highlight = NULL, normalise = TRUE,
+                             log_x = TRUE, log = NULL,
+                             return_data = FALSE, ...) {
+    if (missing(power)) {
+        power <- as.numeric(biomass)
+    }
+    log_x <- parsePlotCDFLog(log, log_x)
+    assert_that(is.flag(total), is.flag(resource),
+                is.flag(background), is.flag(normalise),
+                is.number(power),
+                length(wlim) == 2,
+                length(ylim) == 2)
+
+    args <- list(object = object, species = species,
+                 geometric_mean = geometric_mean,
+                 wlim = wlim, ylim = c(NA, NA),
+                 power = power, total = total,
+                 resource = resource, background = background,
+                 return_data = TRUE)
+    if (!missing(time_range)) {
+        args$time_range <- time_range
+    }
+    plot_dat <- do.call(plotSpectra, args)
+    plot_cdf(plot_dat, object@params, power = power, normalise = normalise,
+             log_x = log_x, wlim = wlim, ylim = ylim,
+             highlight = highlight, return_data = return_data)
+}
+
+#' @rdname plotCDF
+#' @export
+plotCDF.MizerParams <- function(object, species = NULL,
+                                wlim = c(NA, NA), ylim = c(NA, NA),
+                                power = 1, biomass = TRUE,
+                                total = FALSE, resource = TRUE,
+                                background = TRUE,
+                                highlight = NULL, normalise = TRUE,
+                                log_x = TRUE, log = NULL,
+                                return_data = FALSE, ...) {
+    if (missing(power)) {
+        power <- as.numeric(biomass)
+    }
+    log_x <- parsePlotCDFLog(log, log_x)
+    assert_that(is.flag(total), is.flag(resource),
+                is.flag(background), is.flag(normalise),
+                is.number(power),
+                length(wlim) == 2,
+                length(ylim) == 2)
+
+    plot_dat <- plotSpectra(object, species = species,
+                            wlim = wlim, ylim = c(NA, NA),
+                            power = power, total = total,
+                            resource = resource, background = background,
+                            return_data = TRUE)
+    plot_cdf(plot_dat, object, power = power, normalise = normalise,
+             log_x = log_x, wlim = wlim, ylim = ylim,
+             highlight = highlight, return_data = return_data)
+}
+
+plot_cdf <- function(plot_dat, params, power, normalise, log_x, wlim, ylim,
+                     highlight, return_data) {
+    cdf_dat <- prepare_spectra_cdf_data(plot_dat, params,
+                                        normalise = normalise)
+    if (return_data) return(cdf_dat)
+
+    plotDataFrame(cdf_dat, validParams(params),
+                  xlab = "Size [g]", ylab = cdf_y_label(power, normalise),
+                  xtrans = if (log_x) "log10" else "identity",
+                  ytrans = "identity",
+                  xlim = wlim, ylim = ylim,
+                  highlight = highlight, legend_var = "Legend")
+}
+
+prepare_spectra_cdf_data <- function(plot_dat, params, normalise = TRUE) {
+    params <- validParams(params)
+    plot_dat <- plot_dat[order(plot_dat$Species, plot_dat$w), ]
+    plot_dat$value <- plot_dat$value * spectra_bin_width(plot_dat$w, params)
+    plot_dat$value <- ave(plot_dat$value, plot_dat$Species, FUN = cumsum)
+    if (normalise) {
+        totals <- ave(plot_dat$value, plot_dat$Species, FUN = max)
+        plot_dat$value <- plot_dat$value / totals
+    }
+    plot_dat
+}
+
+spectra_bin_width <- function(w, params) {
+    idx <- match(w, params@w_full)
+    if (anyNA(idx)) {
+        missing <- which(is.na(idx))
+        for (i in missing) {
+            idx[i] <- which.min(abs(params@w_full - w[i]))
+        }
+        if (!isTRUE(all.equal(w, params@w_full[idx], scale = 1))) {
+            stop("Could not determine size-bin widths for the spectra data.")
+        }
+    }
+    params@dw_full[idx]
+}
+
+cdf_y_label <- function(power, normalise) {
+    if (normalise) {
+        if (power == 0) {
+            return("Cumulative proportion of abundance")
+        }
+        if (power == 1) {
+            return("Cumulative proportion of biomass")
+        }
+        return("Cumulative proportion")
+    }
+    if (power == 0) {
+        return("Cumulative abundance")
+    }
+    if (power == 1) {
+        return("Cumulative biomass [g]")
+    }
+    paste0("Cumulative number density * w^", power)
+}
+
+parsePlotCDFLog <- function(log, log_x) {
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = FALSE)
+    if (log_axes$log_y) {
+        stop("`plotCDF()` only supports log scaling on the x axis. ",
+             "Use `log = \"x\"` or `log = \"\"`.")
+    }
+    log_axes$log_x
+}
+
 #' Compare two size spectra in the same plot
 #'
 #' `plotSpectra2()` compares the abundance spectra from two `MizerParams` or
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index c31527d33..2321bc8e8 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -169,6 +169,68 @@ test_that("plotSpectraRelative plots symmetric relative difference", {
                                         resource = FALSE), "ggplot")
 })
 
+test_that("plotCDF plots cumulative spectra from small to large sizes", {
+    p <- plotCDF(params, species = species, resource = FALSE, power = 0,
+                 wlim = c(1, NA), return_data = TRUE)
+    expect_true(all(p$w >= 1))
+    expect_true(all(p$value >= 0))
+    for (sp in unique(p$Species)) {
+        sp_dat <- p[p$Species == sp, ]
+        expect_equal(max(sp_dat$value), 1)
+        expect_true(all(diff(sp_dat$value) >= -1e-12))
+    }
+
+    spectra <- plotSpectra(params, species = species, resource = FALSE,
+                           power = 1, wlim = c(1, NA), return_data = TRUE)
+    cdf <- plotCDF(params, species = species, resource = FALSE,
+                   power = 1, wlim = c(1, NA), normalise = FALSE,
+                   return_data = TRUE)
+    widths <- params@dw_full[match(spectra$w, params@w_full)]
+    expected <- sum(spectra$value[spectra$Species == species[[1]]] *
+                        widths[spectra$Species == species[[1]]])
+    observed <- max(cdf$value[cdf$Species == species[[1]]])
+    expect_equal(observed, expected)
+
+    p_plot <- plotCDF(params, species = species, resource = FALSE)
+    expect_s3_class(p_plot, "ggplot")
+    expect_identical(p_plot$scales$get_scales("x")$trans$name, "log-10")
+    expect_match(p_plot$scales$get_scales("y")$name, "biomass",
+                 ignore.case = TRUE)
+
+    p_linear <- plotCDF(params, species = species, resource = FALSE,
+                        log_x = FALSE)
+    expect_identical(p_linear$scales$get_scales("x")$trans$name, "identity")
+
+    p_log_none <- plotCDF(params, species = species, resource = FALSE,
+                          log = "")
+    expect_identical(p_log_none$scales$get_scales("x")$trans$name,
+                     "identity")
+
+    p_log_x <- plotCDF(params, species = species, resource = FALSE,
+                       log = "x")
+    expect_identical(p_log_x$scales$get_scales("x")$trans$name, "log-10")
+    expect_error(plotCDF(params, species = species, resource = FALSE,
+                         log = "y"),
+                 "only supports log scaling on the x axis")
+
+    p_abundance <- plotCDF(params, species = species, resource = FALSE,
+                           power = 0)
+    expect_match(p_abundance$scales$get_scales("y")$name, "abundance",
+                 ignore.case = TRUE)
+})
+
+test_that("plotCDF supports simulations, resource, total and unnormalised output", {
+    p <- plotCDF(sim, species = species, time_range = 1:3,
+                 total = TRUE, resource = TRUE, normalise = FALSE,
+                 return_data = TRUE)
+    expect_true(all(c("Resource", "Total") %in% p$Legend))
+    expect_true(max(p$value) > 1)
+
+    p_plot <- plotCDF(sim, species = species, time_range = 1:3,
+                      total = TRUE, resource = TRUE, normalise = FALSE)
+    expect_s3_class(p_plot, "ggplot")
+})
+
 test_that("yield plotting helpers validate comparison and gear selection", {
     sim_shifted <- sim
     dimnames(sim_shifted@n)$time <- as.character(10:13)

From 5e60cc9e56d9f67201564e725567d9a3f2ddcc6e Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 17:55:17 +0100
Subject: [PATCH 08/10] Add `plotCDF2()` and `plotlyCDF2()` for comparing
 cumulative distributions

---
 NAMESPACE                   |  3 ++
 NEWS.md                     |  5 +-
 R/plots.R                   | 95 +++++++++++++++++++++++++++++++++++++
 tests/testthat/test-plots.R | 26 ++++++++++
 4 files changed, 127 insertions(+), 2 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 255d132ad..2f7b36479 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -408,6 +408,7 @@ export(plot2)
 export(plotBiomass)
 export(plotBiomassObservedVsModel)
 export(plotCDF)
+export(plotCDF2)
 export(plotDataFrame)
 export(plotDiet)
 export(plotFMort)
@@ -424,6 +425,8 @@ export(plotYieldGear)
 export(plotYieldObservedVsModel)
 export(plotlyBiomass)
 export(plotlyBiomassObservedVsModel)
+export(plotlyCDF)
+export(plotlyCDF2)
 export(plotlyFMort)
 export(plotlyFeedingLevel)
 export(plotlyGrowthCurves)
diff --git a/NEWS.md b/NEWS.md
index 5fd73b7ee..7d6822219 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -139,8 +139,9 @@ individual variability in growth to be modelled.
   `plotSpectraRelative()` and `plotlySpectraRelative()` helpers have moved from
   `mizerExperimental` into mizer.
 
-- New `plotCDF()` generic for plotting cumulative abundance or biomass
-  distributions from `MizerParams` and `MizerSim` objects.
+- New `plotCDF()` and `plotCDF2()` generics for plotting cumulative abundance
+  or biomass distributions from `MizerParams` and `MizerSim` objects, together
+  with `plotlyCDF()` and `plotlyCDF2()` wrappers.
 
 - `plotSpectra()` now accepts `log_x`, `log_y`, and `log` arguments for
   controlling axis scaling, matching the mizer array `plot()` methods.
diff --git a/R/plots.R b/R/plots.R
index 4ba5c11f0..229174341 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -56,6 +56,7 @@
 #'   [plotYieldGear()] \tab Plots the total yield of each species by gear against time. \cr
 #'   [plotSpectra()] \tab Plots the abundance (biomass or numbers) spectra of each species and the background community. It is possible to specify a minimum size which is useful for truncating the plot. \cr
 #'   [plotCDF()] \tab Plots cumulative distributions of abundance or biomass over size. \cr
+#'   [plotCDF2()] \tab Compares cumulative distributions from two simulations or parameter objects in one plot. \cr
 #'   [plotSpectra2()] \tab Compares the spectra from two simulations or parameter objects in one plot. \cr
 #'   [plotFeedingLevel()] \tab Plots the feeding level of each species against size. \cr
 #'   [plotPredMort()] \tab Plots the predation mortality of each species against size. \cr
@@ -1107,6 +1108,56 @@ parsePlotCDFLog <- function(log, log_x) {
     log_axes$log_x
 }
 
+#' Compare two cumulative abundance or biomass distributions
+#'
+#' `plotCDF2()` compares cumulative distributions from two `MizerParams` or
+#' `MizerSim` objects in a single plot. Colours identify species or groups and
+#' linetype identifies the object.
+#'
+#' @param object1 First `MizerParams` or `MizerSim` object.
+#' @param object2 Second `MizerParams` or `MizerSim` object.
+#' @param name1,name2 Labels for the two objects, used in the linetype legend.
+#' @inheritParams plotCDF
+#' @param ... Arguments passed to [plotCDF()] for preparing the cumulative
+#'   distribution data, for example `species`, `time_range`, `wlim`,
+#'   `resource`, `background` or `total`.
+#'
+#' @return A ggplot2 object.
+#' @export
+#' @family plotting functions
+#'
+#' @examples
+#' \donttest{
+#' sim1 <- project(NS_params, t_max = 10, progress_bar = FALSE)
+#' sim2 <- project(NS_params, effort = 0.5, t_max = 10, progress_bar = FALSE)
+#' plotCDF2(sim1, sim2, "Original", "Effort = 0.5")
+#' }
+plotCDF2 <- function(object1, object2, name1 = "First", name2 = "Second",
+                     power = 1, normalise = TRUE, log_x = TRUE, log = NULL,
+                     ...) {
+    log_x <- parsePlotCDFLog(log, log_x)
+    assert_that(is.number(power), is.flag(normalise))
+
+    args <- list(...)
+    wlim <- args$wlim %||% c(NA, NA)
+    ylim <- args$ylim %||% c(NA, NA)
+
+    cf1 <- plotCDF(object1, power = power, normalise = normalise,
+                   return_data = TRUE, ...)
+    cf2 <- plotCDF(object2, power = power, normalise = normalise,
+                   return_data = TRUE, ...)
+    params <- if (is(object1, "MizerSim")) object1@params else object1
+
+    plotComparisonDataFrame(cf1, cf2, validParams(params),
+                            name1 = name1, name2 = name2,
+                            xlab = "Size [g]",
+                            ylab = cdf_y_label(power, normalise),
+                            xtrans = if (log_x) "log10" else "identity",
+                            ytrans = "identity",
+                            xlim = wlim, ylim = ylim,
+                            legend_var = "Legend")
+}
+
 #' Compare two size spectra in the same plot
 #'
 #' `plotSpectra2()` compares the abundance spectra from two `MizerParams` or
@@ -1232,6 +1283,50 @@ plotlySpectraRelative <- function(object1, object2, log_x = TRUE,
              tooltip = c("Legend", "w", "rel_diff"))
 }
 
+#' @rdname plotCDF
+#' @return `plotlyCDF()` returns a plotly object.
+#' @export
+plotlyCDF <- function(object, species = NULL,
+                      time_range, geometric_mean = FALSE,
+                      wlim = c(NA, NA), ylim = c(NA, NA),
+                      power = 1, biomass = TRUE,
+                      total = FALSE, resource = TRUE,
+                      background = TRUE,
+                      highlight = NULL, normalise = TRUE,
+                      log_x = TRUE, log = NULL, ...) {
+    args <- list(object = object, species = species,
+                 geometric_mean = geometric_mean,
+                 wlim = wlim, ylim = ylim,
+                 biomass = biomass, total = total,
+                 resource = resource, background = background,
+                 highlight = highlight, normalise = normalise,
+                 log_x = log_x, log = log, ...)
+    if (!missing(time_range)) {
+        args$time_range <- time_range
+    }
+    if (!missing(power)) {
+        args$power <- power
+    }
+    ggplotly(do.call("plotCDF", args),
+             tooltip = c("Species", "w", "value"))
+}
+
+#' @rdname plotCDF2
+#' @return `plotlyCDF2()` returns a plotly object.
+#' @export
+plotlyCDF2 <- function(object1, object2, name1 = "First", name2 = "Second",
+                       power = 1, normalise = TRUE,
+                       log_x = TRUE, log = NULL, ...) {
+    args <- list(object1 = object1, object2 = object2,
+                 name1 = name1, name2 = name2,
+                 normalise = normalise, log_x = log_x, log = log, ...)
+    if (!missing(power)) {
+        args$power <- power
+    }
+    ggplotly(do.call("plotCDF2", args),
+             tooltip = c("Species", "w", "value", "Model"))
+}
+
 #' @rdname plotSpectra
 #' @export
 plotlySpectra <- function(object, species = NULL,
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 2321bc8e8..bef68abe7 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -106,6 +106,10 @@ test_that("plotly functions do not throw error", {
 
 test_that("plotly wrappers return plotly objects for spectra and rate plots", {
     expect_s3_class(plotlySpectra(params, species = species), "plotly")
+    expect_s3_class(plotlyCDF(params, species = species,
+                              resource = FALSE), "plotly")
+    expect_s3_class(plotlyCDF2(params, params, species = species,
+                               resource = FALSE), "plotly")
     expect_s3_class(plotlySpectraRelative(params, params, species = species,
                                           resource = FALSE), "plotly")
     expect_s3_class(plotlyPredMort(sim, species = species), "plotly")
@@ -231,6 +235,28 @@ test_that("plotCDF supports simulations, resource, total and unnormalised output
     expect_s3_class(p_plot, "ggplot")
 })
 
+test_that("plotCDF2 compares cumulative distributions", {
+    p <- plotCDF2(params, params, name1 = "Original", name2 = "Changed",
+                  species = species, total = TRUE, resource = FALSE,
+                  wlim = c(1, NA), normalise = FALSE, log = "")
+    expect_s3_class(p, "ggplot")
+    expect_identical(levels(p$data$Model), c("Original", "Changed"))
+    expect_true("Total" %in% p$data$Legend)
+    expect_true(all(p$data$w >= 1))
+    expect_identical(p$scales$get_scales("x")$trans$name, "identity")
+
+    p_log <- plotCDF2(params, sim, species = species, resource = FALSE,
+                      log = "x")
+    expect_s3_class(p_log, "ggplot")
+    expect_identical(p_log$scales$get_scales("x")$trans$name, "log-10")
+
+    expect_s3_class(plotCDF2(sim, sim0, species = species,
+                             time_range = 1:3, resource = FALSE),
+                    "ggplot")
+    expect_error(plotCDF2(params, sim, species = species, log = "y"),
+                 "only supports log scaling on the x axis")
+})
+
 test_that("yield plotting helpers validate comparison and gear selection", {
     sim_shifted <- sim
     dimnames(sim_shifted@n)$time <- as.character(10:13)

From 6ac290a5b97a068fdd7221c34a0de23369c49aa8 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 18:06:04 +0100
Subject: [PATCH 09/10] Add plotly wrappers for new plotting functions:
 plotlyDiet() and plotlySpectra2()

---
 NAMESPACE                   |  2 ++
 R/plots.R                   | 28 +++++++++++++++++++++++++---
 tests/testthat/test-plots.R |  4 ++++
 3 files changed, 31 insertions(+), 3 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 2f7b36479..04a614a7e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -427,11 +427,13 @@ export(plotlyBiomass)
 export(plotlyBiomassObservedVsModel)
 export(plotlyCDF)
 export(plotlyCDF2)
+export(plotlyDiet)
 export(plotlyFMort)
 export(plotlyFeedingLevel)
 export(plotlyGrowthCurves)
 export(plotlyPredMort)
 export(plotlySpectra)
+export(plotlySpectra2)
 export(plotlySpectraRelative)
 export(plotlyYield)
 export(plotlyYieldGear)
diff --git a/R/plots.R b/R/plots.R
index 229174341..9d61537e6 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -70,9 +70,11 @@
 #'
 #' The static plotting functions use ggplot2 and return a ggplot object. This
 #' means that you can manipulate the plot further after its creation using the
-#' ggplot grammar of graphics. Many named plot functions also have a plotly
-#' counterpart, for example [plotlyBiomass()] or [plotlySpectra()], for
-#' interactive exploration.
+#' ggplot grammar of graphics. The named high-level plot functions have plotly
+#' counterparts, for example [plotlyBiomass()] or [plotlySpectra()], for
+#' interactive exploration. Generic and compositional plotting APIs, such as
+#' [plot()], [plot2()], [plotRelative()] and [addPlot()], do not have separate
+#' plotly wrappers. Use [ggplotly()] on the ggplot object they return.
 #'
 #' While most plot functions take their data from a MizerSim object, some of
 #' those that make plots representing data at a single time can also take their
@@ -1220,6 +1222,18 @@ spectra_y_label <- function(power) {
     paste0("Number density * w^", power)
 }
 
+#' @rdname plotSpectra2
+#' @return `plotlySpectra2()` returns a plotly object.
+#' @export
+plotlySpectra2 <- function(object1, object2, name1 = "First",
+                           name2 = "Second", power = 1,
+                           log_x = TRUE, log_y = TRUE, log = NULL, ...) {
+    ggplotly(plotSpectra2(object1, object2, name1 = name1, name2 = name2,
+                          power = power, log_x = log_x, log_y = log_y,
+                          log = log, ...),
+             tooltip = c("Species", "w", "value", "Model"))
+}
+
 #' Plot the relative difference between two spectra
 #'
 #' `plotSpectraRelative()` plots the difference between the spectra relative to
@@ -2088,6 +2102,14 @@ plot_diet <- function(params, n, diet, species, return_data) {
     p
 }
 
+#' @rdname plotDiet
+#' @return `plotlyDiet()` returns a plotly object.
+#' @export
+plotlyDiet <- function(object, species = NULL, ...) {
+    ggplotly(plotDiet(object, species = species, ...),
+             tooltip = c("Predator", "w", "Proportion", "Prey"))
+}
+
 
 #### plot ####
 #' Summary plot for `MizerSim` objects
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index bef68abe7..cefda6250 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -98,14 +98,17 @@ test_that("plotly functions do not throw error", {
     expect_error(plotlyYield(sim, sim), NA)
     expect_error(plotlyYieldGear(sim, species = species), NA)
     expect_error(plotlySpectra(params, species = species), NA)
+    expect_error(plotlySpectra2(params, sim, species = species), NA)
     expect_error(plotlyPredMort(sim, species = species), NA)
     expect_error(plotlyFMort(sim, species = species), NA)
     expect_error(plotlyGrowthCurves(sim, species = species), NA)
     expect_error(plotlyGrowthCurves(params, species = species), NA)
+    expect_error(plotlyDiet(params, species = species[[1]]), NA)
 })
 
 test_that("plotly wrappers return plotly objects for spectra and rate plots", {
     expect_s3_class(plotlySpectra(params, species = species), "plotly")
+    expect_s3_class(plotlySpectra2(params, sim, species = species), "plotly")
     expect_s3_class(plotlyCDF(params, species = species,
                               resource = FALSE), "plotly")
     expect_s3_class(plotlyCDF2(params, params, species = species,
@@ -117,6 +120,7 @@ test_that("plotly wrappers return plotly objects for spectra and rate plots", {
     expect_s3_class(plotlyGrowthCurves(sim, species = species), "plotly")
     expect_s3_class(plotlyFeedingLevel(sim, species = species,
                                        include_critical = TRUE), "plotly")
+    expect_s3_class(plotlyDiet(params, species = species[[1]]), "plotly")
 })
 
 test_that("plotSpectra2 compares spectra from params and sims", {

From dfe08b8569d57eeadb1fdb773f340d10bcad76c0 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Tue, 12 May 2026 18:39:58 +0100
Subject: [PATCH 10/10] Add mizer_plot class and tooltip handling for ggplotly
 integration

---
 NAMESPACE                   |  1 +
 R/plots.R                   | 39 +++++++++++++++++++++++++++++++++----
 tests/testthat/test-plots.R | 16 +++++++++++++++
 3 files changed, 52 insertions(+), 4 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 04a614a7e..1070ed4d7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -155,6 +155,7 @@ S3method(plotYieldObservedVsModel,MizerSim)
 S3method(plotly::ggplotly,ArraySpeciesBySize)
 S3method(plotly::ggplotly,ArrayTimeBySpecies)
 S3method(plotly::ggplotly,ArrayTimeBySpeciesBySize)
+S3method(plotly::ggplotly,mizer_plot)
 S3method(print,ArraySpeciesBySize)
 S3method(print,ArrayTimeBySpecies)
 S3method(print,ArrayTimeBySpeciesBySize)
diff --git a/R/plots.R b/R/plots.R
index 9d61537e6..ac2379b7c 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -254,7 +254,34 @@ plotDataFrame <- function(frame, params, style = "line", xlab = waiver(),
         p <- p + facet_wrap(wrap_var, scales = wrap_scale)
     }
 
-    p
+    make_mizer_plot(p, mizer_tooltip_vars(frame, group_var, x_var, y_var,
+                                          legend_var))
+}
+
+make_mizer_plot <- function(plot, tooltip) {
+    attr(plot, "mizer_tooltip") <- tooltip
+    class(plot) <- unique(c("mizer_plot", class(plot)))
+    plot
+}
+
+mizer_tooltip_vars <- function(frame, group_var, x_var, y_var,
+                               legend_var = NULL, extra = NULL) {
+    tooltip <- c(group_var, x_var, y_var)
+    if (!is.null(legend_var) && legend_var %in% names(frame) &&
+            !identical(legend_var, group_var) &&
+            any(as.character(frame[[legend_var]]) !=
+                    as.character(frame[[group_var]]), na.rm = TRUE)) {
+        tooltip <- c(tooltip, legend_var)
+    }
+    unique(c(tooltip, extra))
+}
+
+#' @exportS3Method plotly::ggplotly
+ggplotly.mizer_plot <- function(p = ggplot2::last_plot(), ...,
+                                tooltip = attr(p, "mizer_tooltip") %||%
+                                    "all") {
+    class(p) <- setdiff(class(p), "mizer_plot")
+    ggplotly(p, ..., tooltip = tooltip)
 }
 
 plotComparisonDataFrame <- function(frame1, frame2, params,
@@ -295,8 +322,8 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
     ybreaks <- waiver()
     if (ytrans == "log10") ybreaks <- log_breaks(n = y_ticks)
 
-    ggplot(frame,
-           aes(group = interaction(.data[[group_var]], .data[["Model"]]))) +
+    p <- ggplot(frame,
+                aes(group = interaction(.data[[group_var]], .data[["Model"]]))) +
         scale_y_continuous(trans = ytrans, breaks = ybreaks,
                            labels = prettyNum, name = ylab,
                            limits = ylim) +
@@ -307,6 +334,8 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                       linetype = .data[["Model"]])) +
         scale_colour_manual(values = linecolour) +
         scale_linetype_discrete(drop = FALSE)
+    make_mizer_plot(p, mizer_tooltip_vars(frame, group_var, x_var, y_var,
+                                          legend_var, extra = "Model"))
 }
 
 plotRelativeDataFrame <- function(frame1, frame2, params,
@@ -347,7 +376,7 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
     xbreaks <- waiver()
     if (xtrans == "log10") xbreaks <- log_breaks()
 
-    ggplot(frame, aes(group = .data[[group_var]])) +
+    p <- ggplot(frame, aes(group = .data[[group_var]])) +
         scale_y_continuous(name = "Relative difference", limits = ylim) +
         scale_x_continuous(trans = xtrans, breaks = xbreaks, name = xlab,
                            limits = xlim) +
@@ -356,6 +385,8 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
         geom_line(aes(x = .data[[x_var]], y = .data[["rel_diff"]],
                       colour = .data[[legend_var]])) +
         scale_colour_manual(values = linecolour)
+    make_mizer_plot(p, mizer_tooltip_vars(frame, group_var, x_var, "rel_diff",
+                                          legend_var))
 }
 
 relative_difference <- function(first, second) {
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index cefda6250..d1cfb6597 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -123,6 +123,22 @@ test_that("plotly wrappers return plotly objects for spectra and rate plots", {
     expect_s3_class(plotlyDiet(params, species = species[[1]]), "plotly")
 })
 
+test_that("ggplotly(plot(...)) uses concise mizer tooltips", {
+    p <- plot(getEncounter(NS_params), species = "Cod")
+    expect_s3_class(p, "mizer_plot")
+    gp <- ggplotly(p)
+    first_tip <- gp$x$data[[1]]$text[[1]]
+    expect_true(grepl("Species: Cod", first_tip, fixed = TRUE))
+    expect_true(grepl("w:", first_tip, fixed = TRUE))
+    expect_true(grepl("value:", first_tip, fixed = TRUE))
+    legend_matches <- gregexpr("Legend:", first_tip, fixed = TRUE)[[1]]
+    expect_lte(sum(legend_matches > 0), 1)
+
+    ggplot2::set_last_plot(p)
+    gp_last <- ggplotly()
+    expect_identical(gp_last$x$data[[1]]$text[[1]], first_tip)
+})
+
 test_that("plotSpectra2 compares spectra from params and sims", {
     p_params <- plotSpectra2(params, params, name1 = "Original",
                              name2 = "Changed", species = species,