From 73987eb179746cce1c288aa3d73da0d1e3b41541 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Thu, 16 Jul 2026 19:08:24 +0000
Subject: [PATCH 1/4] Add plot(GCIMSDataset): several samples on a shared color
 scale

Adds a new plot() method for GCIMSDataset that renders several samples
side by side (via cowplot::plot_grid()), all on the same intensity
color scale, so they're visually comparable.

The shared scale is controlled by a new intensity_range argument:
"global" (default), "ranged", a fixed c(min, max), or a list/vector
whose min/max are independently a number, "global" or "ranged". The
resolver (resolve_intensity_range()) only evaluates whichever of
"global"/"ranged" is actually referenced, and at most once.

"global" is backed by a per-sample intensity range cached on the
dataset (ds$IntensityRange, internal only), computed for free: it
rides along inside .extract_RIC_and_TIS_fun_extract, which already
loads each sample's full intensity matrix to compute RIC/TIS on every
realize() -- adding range(intmat) costs one extra range() call on data
already in memory, no new pass, no new DelayedOperation. Using the
sample's full uncropped range (rather than scoping to dt_range/
rt_range/sample selection) is deliberate: it's always a safe bound and
keeps the scale stable across differently-cropped or -subset calls.
It doesn't apply with remove_baseline = TRUE (the cache holds raw
intensities), which errors clearly toward intensity_range = "ranged"
instead.

plot(GCIMSSample) gained a fill_range parameter (default NULL, current
auto-scaling behavior unchanged) to make one panel's color scale
overridable, threaded through to mat_to_gplot() and mat_to_nativeRaster()
(which already had an unused rangex override hook).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019V3CHRGSzcEu3k6tpUFv56
---
 R/aaa-class-GCIMSDataset.R              |   6 +
 R/getTIS_getRIC-GCIMSDataset.R          |  21 +++-
 R/plot-GCIMSDataset.R                   | 139 ++++++++++++++++++++++++
 R/plot-GCIMSSample.R                    |  15 ++-
 man/GCIMSDataset.Rd                     |   5 +
 man/plot-GCIMSDataset-ANY-method.Rd     |  62 +++++++++++
 man/plot-GCIMSSample-ANY-method.Rd      |   8 +-
 tests/testthat/test-plot-GCIMSDataset.R | 124 +++++++++++++++++++++
 tests/testthat/test-plot-GCIMSSample.R  |  11 ++
 9 files changed, 384 insertions(+), 7 deletions(-)
 create mode 100644 R/plot-GCIMSDataset.R
 create mode 100644 man/plot-GCIMSDataset-ANY-method.Rd
 create mode 100644 tests/testthat/test-plot-GCIMSDataset.R

diff --git a/R/aaa-class-GCIMSDataset.R b/R/aaa-class-GCIMSDataset.R
index d367742d..c95ecefb 100644
--- a/R/aaa-class-GCIMSDataset.R
+++ b/R/aaa-class-GCIMSDataset.R
@@ -79,6 +79,11 @@ GCIMSDataset <- R6::R6Class("GCIMSDataset",
     dt_ref = NULL, # numeric or NULL
     #' @field rt_ref A numeric retention time of reference
     rt_ref = NULL, # numeric or NULL
+    #' @field IntensityRange A matrix of n_samples x c(min, max) with each
+    #' sample's raw intensity range, extracted alongside TIS/RIC. Internal
+    #' only: no S4 generic exposes this, it only backs `plot(GCIMSDataset)`'s
+    #' shared color scale across samples.
+    IntensityRange = NULL, # matrix or NULL
     #' @field userData A list to store arbitrary data in the dataset
     userData = list(), # list
     # Methods:
@@ -221,6 +226,7 @@ GCIMSDataset <- R6::R6Class("GCIMSDataset",
       self$RIC <- NULL
       self$dt_ref <- NULL
       self$rt_ref <- NULL
+      self$IntensityRange <- NULL
       private$delayed_dataset$subset(samples)
       self
     },
diff --git a/R/getTIS_getRIC-GCIMSDataset.R b/R/getTIS_getRIC-GCIMSDataset.R
index 1cfa79b6..5ba253a2 100644
--- a/R/getTIS_getRIC-GCIMSDataset.R
+++ b/R/getTIS_getRIC-GCIMSDataset.R
@@ -114,7 +114,11 @@ setMethod(
   ric <- intmat[ric_pos, ]
   ric <- max(ric) - ric
   ric <- ric/sum(ric)
-  list(ric = ric, tis = tis, rt = rt, dt = dt)
+  # intmat is already loaded here, so this rides along for free -- it backs
+  # plot(GCIMSDataset)'s shared intensity scale across samples (see
+  # dataset_intensity_range()), without a dedicated extraction pass.
+  intensity_range <- range(intmat)
+  list(ric = ric, tis = tis, rt = rt, dt = dt, intensity_range = intensity_range)
 }
 
 .extract_RIC_and_TIS_fun_aggregate <- function(ds, objs) {
@@ -123,6 +127,7 @@ setMethod(
   tiss <- purrr::map(objs, "tis")
   dtimes <- purrr::map(objs, "dt")
   rtimes <- purrr::map(objs, "rt")
+  intensity_ranges <- purrr::map(objs, "intensity_range")
 
   dt_ref <- ds$dt_ref
   rt_ref <- ds$rt_ref
@@ -144,9 +149,23 @@ setMethod(
   }
   stopifnot(nrow(ds$TIS) == num_samples)
   stopifnot(nrow(ds$RIC) == num_samples)
+  ds$IntensityRange <- do.call(rbind, intensity_ranges)
+  colnames(ds$IntensityRange) <- c("min", "max")
   ds
 }
 
+#' The dataset's overall raw intensity range, realizing/refreshing the cache if needed
+#' @param object A [GCIMSDataset] object
+#' @return A length-2 numeric vector `c(min, max)`
+#' @noRd
+dataset_intensity_range <- function(object) {
+  if (object$hasDelayedOps() || is.null(object$IntensityRange)) {
+    object$extract_RIC_and_TIS()
+    object$realize()
+  }
+  range(object$IntensityRange)
+}
+
 #' Extract the Reverse Ion Chromatogram and Total Ion Spectrum from the samples
 #'
 #' @param object A GCIMSDataset object
diff --git a/R/plot-GCIMSDataset.R b/R/plot-GCIMSDataset.R
new file mode 100644
index 00000000..effd6b02
--- /dev/null
+++ b/R/plot-GCIMSDataset.R
@@ -0,0 +1,139 @@
+#' Resolve the `intensity_range` argument of `plot(GCIMSDataset)` into `c(min, max)`
+#'
+#' @param intensity_range One of `"global"`, `"ranged"`, a numeric vector of
+#' length 2, or a length-2 vector/list whose elements are independently a
+#' number, `"global"` or `"ranged"`.
+#' @param global_range A zero-argument function returning the cached,
+#' whole-dataset `c(min, max)`. Only called if actually referenced.
+#' @param ranged_range A zero-argument function returning the `c(min, max)`
+#' of exactly what will be plotted. Only called if actually referenced.
+#' @return A numeric vector `c(min, max)`
+#' @noRd
+resolve_intensity_range <- function(intensity_range, global_range, ranged_range) {
+  resolve_endpoint <- function(value, position) {
+    if (is.numeric(value) && length(value) == 1) {
+      return(value)
+    }
+    if (identical(value, "global")) {
+      return(global_range()[position])
+    }
+    if (identical(value, "ranged")) {
+      return(ranged_range()[position])
+    }
+    cli_abort('Each element of intensity_range should be a number, "global" or "ranged"')
+  }
+
+  if (is.character(intensity_range) && length(intensity_range) == 1) {
+    if (intensity_range == "global") {
+      return(global_range())
+    } else if (intensity_range == "ranged") {
+      return(ranged_range())
+    }
+    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
+  }
+
+  if (length(intensity_range) != 2) {
+    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
+  }
+
+  c(
+    resolve_endpoint(intensity_range[[1]], 1L),
+    resolve_endpoint(intensity_range[[2]], 2L)
+  )
+}
+
+#' Topographical plot of several samples of a GCIMSDataset, on a shared color scale
+#'
+#' Several samples are plotted side by side, on the same intensity color
+#' scale, so they are visually comparable.
+#'
+#' @param x A [GCIMSDataset] object
+#' @param sample A number, a string, or a vector of numbers/strings with the
+#' sample index(es) or name(s) to plot. If `NULL` (the default), all samples
+#' are plotted.
+#' @inheritParams dt_rt_range_normalization
+#' @param ... Ignored
+#' @param remove_baseline Set to `TRUE` to subtract each sample's estimated
+#' baseline first
+#' @param trans The transformation to the intensity values, see `plot,GCIMSSample-method`
+#' @param intensity_range Controls the shared color scale across all plotted
+#' samples. One of:
+#' - `"global"` (the default): each sample's full, uncropped intensity range,
+#'   cached on the dataset (computed once, alongside TIS/RIC, regardless of
+#'   `sample`/`dt_range`/`rt_range`). Free, and keeps the scale stable across
+#'   different calls with different `sample`/`dt_range`/`rt_range` selections.
+#'   Not valid together with `remove_baseline = TRUE`, since the cache holds
+#'   raw intensities.
+#' - `"ranged"`: the range of exactly what's plotted (respecting `dt_range`,
+#'   `rt_range` and `remove_baseline`). Requires loading every selected
+#'   sample's data twice (once to compute the range, once to plot).
+#' - A numeric vector of length 2, `c(min, max)`: fixed limits.
+#' - A length-2 vector/list whose elements are independently a number,
+#'   `"global"` or `"ranged"`, e.g. `list(min = 0, max = "global")`.
+#' @param ncol Number of columns for [cowplot::plot_grid()]. `NULL` lets `cowplot` decide.
+#' @return A combined plot with one panel per sample
+#' @export
+setMethod(
+  "plot",
+  "GCIMSDataset",
+  function(x, sample = NULL, dt_range = NULL, rt_range = NULL, ...,
+           remove_baseline = FALSE, trans = "cubic_root",
+           intensity_range = "global", ncol = NULL) {
+    require_pkgs("cowplot")
+    sample_names_all <- sampleNames(x)
+    if (is.null(sample)) {
+      sample <- sample_names_all
+    }
+    selected_names <- sample_name_or_number_to_both(sample, sample_names_all)$name
+
+    samples <- stats::setNames(
+      purrr::map(selected_names, function(nm) x$getSample(nm)),
+      selected_names
+    )
+
+    cropped_intensity <- function(s) {
+      dt <- dtime(s)
+      rt <- rtime(s)
+      idx <- dt_rt_range_normalization(dt, rt, dt_range, rt_range)
+      intmat <- intensity(s, idx)
+      if (isTRUE(remove_baseline)) {
+        basel <- baseline(s)[idx[["dt_logical"]], idx[["rt_logical"]]]
+        intmat <- intmat - basel
+      }
+      intmat
+    }
+
+    ranged_range_cache <- NULL
+    get_ranged_range <- function() {
+      if (is.null(ranged_range_cache)) {
+        ranges <- purrr::map(samples, function(s) range(cropped_intensity(s)))
+        ranged_range_cache <<- range(unlist(ranges))
+      }
+      ranged_range_cache
+    }
+    get_global_range <- function() {
+      if (isTRUE(remove_baseline)) {
+        cli_abort(
+          c(
+            "intensity_range can't use {.val global} together with {.code remove_baseline = TRUE}",
+            "i" = "The cached global range holds raw intensities. Use {.code intensity_range = \"ranged\"} instead"
+          )
+        )
+      }
+      dataset_intensity_range(x)
+    }
+
+    limits <- resolve_intensity_range(intensity_range, get_global_range, get_ranged_range)
+
+    panels <- purrr::map(selected_names, function(nm) {
+      plot(
+        samples[[nm]],
+        dt_range = dt_range, rt_range = rt_range,
+        remove_baseline = remove_baseline, trans = trans,
+        fill_range = limits
+      ) + ggplot2::labs(title = nm)
+    })
+
+    cowplot::plot_grid(plotlist = panels, ncol = ncol)
+  }
+)
diff --git a/R/plot-GCIMSSample.R b/R/plot-GCIMSSample.R
index a3c9e3eb..2224e008 100644
--- a/R/plot-GCIMSSample.R
+++ b/R/plot-GCIMSSample.R
@@ -7,6 +7,10 @@
 #' @param remove_baseline Set to `TRUE` to subtract the estimated baseline first
 #' @param trans The transformation to the intensity values. "cubic_root" is the default. "intensity" is also valid.
 #' See the `trans` argument in [ggplot2::continuous_scale()] for other possibilities.
+#' @param fill_range A numeric vector of length 2 with the intensity range the
+#' color scale should span, or `NULL` (the default) to use this sample's own
+#' `range()`. Set it explicitly to make several plots comparable on the same
+#' color scale.
 #' @return A plot of the GCIMSSample
 #' @examples
 #' dummy_obj <-GCIMSSample(
@@ -22,7 +26,7 @@
 setMethod(
   "plot",
   "GCIMSSample",
-  function(x, dt_range = NULL, rt_range = NULL, ..., remove_baseline = FALSE, trans = "cubic_root") {
+  function(x, dt_range = NULL, rt_range = NULL, ..., remove_baseline = FALSE, trans = "cubic_root", fill_range = NULL) {
   dt <- dtime(x)
   rt <- rtime(x)
   idx <- dt_rt_range_normalization(dt, rt, dt_range, rt_range)
@@ -38,12 +42,13 @@ setMethod(
     dt_max = idx$dt_ms_max,
     rt_min = idx$rt_s_min,
     rt_max = idx$rt_s_max,
-    trans = trans
+    trans = trans,
+    fill_range = fill_range
   )
 })
 
 
-mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt_max = NULL, trans = "cubic_root") {
+mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt_max = NULL, trans = "cubic_root", fill_range = NULL) {
   require_pkgs(c("farver", "viridisLite"))
   if (is.null(dt_min)) {
     dt_min <- as.numeric(rownames(intmat)[1L])
@@ -57,7 +62,7 @@ mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt
   if (is.null(rt_max)) {
     rt_max <- as.numeric(colnames(intmat)[ncol(intmat)])
   }
-  minmax <- range(intmat)
+  minmax <- if (is.null(fill_range)) range(intmat) else fill_range
 
   if (is.character(trans)) {
     trans_func <- paste0(trans, "_trans")
@@ -76,7 +81,7 @@ mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt
   colormap <- farver::encode_native(
     viridisLite::viridis(256L, direction = -1, option = "G")
   )
-  nr <- mat_to_nativeRaster(intmat_trans, colormap)
+  nr <- mat_to_nativeRaster(intmat_trans, colormap, rangex = trans$transform(minmax))
 
   # The geom_rect is fake and it is only used to force the fill legend to appear
   # The geom_rect  limits are used to help set the plot limits
diff --git a/man/GCIMSDataset.Rd b/man/GCIMSDataset.Rd
index bbe2408d..0e674dd3 100644
--- a/man/GCIMSDataset.Rd
+++ b/man/GCIMSDataset.Rd
@@ -96,6 +96,11 @@ dummy_dataset <- GCIMSDataset$new(
 
 \item{\code{rt_ref}}{A numeric retention time of reference}
 
+\item{\code{IntensityRange}}{A matrix of n_samples x c(min, max) with each
+sample's raw intensity range, extracted alongside TIS/RIC. Internal
+only: no S4 generic exposes this, it only backs \code{plot(GCIMSDataset)}'s
+shared color scale across samples.}
+
 \item{\code{userData}}{A list to store arbitrary data in the dataset}
 }
 \if{html}{\out{</div>}}
diff --git a/man/plot-GCIMSDataset-ANY-method.Rd b/man/plot-GCIMSDataset-ANY-method.Rd
new file mode 100644
index 00000000..9318d8e7
--- /dev/null
+++ b/man/plot-GCIMSDataset-ANY-method.Rd
@@ -0,0 +1,62 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot-GCIMSDataset.R
+\name{plot,GCIMSDataset,ANY-method}
+\alias{plot,GCIMSDataset,ANY-method}
+\title{Topographical plot of several samples of a GCIMSDataset, on a shared color scale}
+\usage{
+\S4method{plot}{GCIMSDataset,ANY}(
+  x,
+  sample = NULL,
+  dt_range = NULL,
+  rt_range = NULL,
+  ...,
+  remove_baseline = FALSE,
+  trans = "cubic_root",
+  intensity_range = "global",
+  ncol = NULL
+)
+}
+\arguments{
+\item{x}{A \link{GCIMSDataset} object}
+
+\item{sample}{A number, a string, or a vector of numbers/strings with the
+sample index(es) or name(s) to plot. If \code{NULL} (the default), all samples
+are plotted.}
+
+\item{dt_range}{The minimum and maximum drift times to extract (length 2 vector)}
+
+\item{rt_range}{The minimum and maximum retention times to extract (length 2 vector)}
+
+\item{...}{Ignored}
+
+\item{remove_baseline}{Set to \code{TRUE} to subtract each sample's estimated
+baseline first}
+
+\item{trans}{The transformation to the intensity values, see \verb{plot,GCIMSSample-method}}
+
+\item{intensity_range}{Controls the shared color scale across all plotted
+samples. One of:
+\itemize{
+\item \code{"global"} (the default): each sample's full, uncropped intensity range,
+cached on the dataset (computed once, alongside TIS/RIC, regardless of
+\code{sample}/\code{dt_range}/\code{rt_range}). Free, and keeps the scale stable across
+different calls with different \code{sample}/\code{dt_range}/\code{rt_range} selections.
+Not valid together with \code{remove_baseline = TRUE}, since the cache holds
+raw intensities.
+\item \code{"ranged"}: the range of exactly what's plotted (respecting \code{dt_range},
+\code{rt_range} and \code{remove_baseline}). Requires loading every selected
+sample's data twice (once to compute the range, once to plot).
+\item A numeric vector of length 2, \code{c(min, max)}: fixed limits.
+\item A length-2 vector/list whose elements are independently a number,
+\code{"global"} or \code{"ranged"}, e.g. \code{list(min = 0, max = "global")}.
+}}
+
+\item{ncol}{Number of columns for \code{\link[cowplot:plot_grid]{cowplot::plot_grid()}}. \code{NULL} lets \code{cowplot} decide.}
+}
+\value{
+A combined plot with one panel per sample
+}
+\description{
+Several samples are plotted side by side, on the same intensity color
+scale, so they are visually comparable.
+}
diff --git a/man/plot-GCIMSSample-ANY-method.Rd b/man/plot-GCIMSSample-ANY-method.Rd
index 49194da6..766e7e37 100644
--- a/man/plot-GCIMSSample-ANY-method.Rd
+++ b/man/plot-GCIMSSample-ANY-method.Rd
@@ -10,7 +10,8 @@
   rt_range = NULL,
   ...,
   remove_baseline = FALSE,
-  trans = "cubic_root"
+  trans = "cubic_root",
+  fill_range = NULL
 )
 }
 \arguments{
@@ -26,6 +27,11 @@
 
 \item{trans}{The transformation to the intensity values. "cubic_root" is the default. "intensity" is also valid.
 See the \code{trans} argument in \code{\link[ggplot2:continuous_scale]{ggplot2::continuous_scale()}} for other possibilities.}
+
+\item{fill_range}{A numeric vector of length 2 with the intensity range the
+color scale should span, or \code{NULL} (the default) to use this sample's own
+\code{range()}. Set it explicitly to make several plots comparable on the same
+color scale.}
 }
 \value{
 A plot of the GCIMSSample
diff --git a/tests/testthat/test-plot-GCIMSDataset.R b/tests/testthat/test-plot-GCIMSDataset.R
new file mode 100644
index 00000000..dae40053
--- /dev/null
+++ b/tests/testthat/test-plot-GCIMSDataset.R
@@ -0,0 +1,124 @@
+# --- resolve_intensity_range() -----------------------------------------
+
+make_counting_range <- function(value) {
+  n_calls <- 0L
+  fn <- function() {
+    n_calls <<- n_calls + 1L
+    value
+  }
+  attr(fn, "n_calls") <- function() n_calls
+  fn
+}
+
+test_that("resolve_intensity_range('global'/'ranged') calls only the matching closure, once", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_equal(resolve_intensity_range("global", global_range, ranged_range), c(0, 10))
+  expect_equal(attr(global_range, "n_calls")(), 1L)
+  expect_equal(attr(ranged_range, "n_calls")(), 0L)
+})
+
+test_that("resolve_intensity_range() with a fixed numeric vector never calls either closure", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_equal(resolve_intensity_range(c(1, 9), global_range, ranged_range), c(1, 9))
+  expect_equal(attr(global_range, "n_calls")(), 0L)
+  expect_equal(attr(ranged_range, "n_calls")(), 0L)
+})
+
+test_that("resolve_intensity_range() resolves mixed list(min=, max=) elements independently", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_equal(resolve_intensity_range(list(min = 1, max = "global"), global_range, ranged_range), c(1, 10))
+  expect_equal(attr(global_range, "n_calls")(), 1L)
+  expect_equal(attr(ranged_range, "n_calls")(), 0L)
+
+  expect_equal(resolve_intensity_range(list("ranged", 99), global_range, ranged_range), c(2, 99))
+  expect_equal(attr(ranged_range, "n_calls")(), 1L)
+})
+
+test_that("resolve_intensity_range() rejects invalid specs", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_error(resolve_intensity_range("nope", global_range, ranged_range), "intensity_range should be")
+  expect_error(resolve_intensity_range(c(1, 2, 3), global_range, ranged_range), "intensity_range should be")
+  expect_error(resolve_intensity_range(list("nope", 1), global_range, ranged_range), "should be a number")
+})
+
+# --- IntensityRange caching ----------------------------------------------
+
+make_range_dataset <- function() {
+  s1 <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
+  s2 <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(100:124, nrow = 5))
+  GCIMSDataset$new_from_list(samples = list(s1 = s1, s2 = s2), on_ram = TRUE, scratch_dir = NULL)
+}
+
+test_that("realizing a dataset populates IntensityRange alongside TIS/RIC", {
+  old_bpparam <- BiocParallel::bpparam()
+  BiocParallel::register(BiocParallel::SerialParam())
+  on.exit(BiocParallel::register(old_bpparam))
+
+  ds <- make_range_dataset()
+  ds$realize()
+
+  expect_equal(unname(ds$IntensityRange["s1", ]), c(1, 25))
+  expect_equal(unname(ds$IntensityRange["s2", ]), c(100, 124))
+})
+
+test_that("subsetting a dataset resets the cached IntensityRange", {
+  old_bpparam <- BiocParallel::bpparam()
+  BiocParallel::register(BiocParallel::SerialParam())
+  on.exit(BiocParallel::register(old_bpparam))
+
+  ds <- make_range_dataset()
+  ds$realize()
+  expect_false(is.null(ds$IntensityRange))
+
+  ds$subset("s1", inplace = TRUE)
+
+  expect_null(ds$IntensityRange)
+})
+
+# --- plot(GCIMSDataset) ----------------------------------------------------
+
+test_that("plot(GCIMSDataset) returns a combined plot using the cached global range by default", {
+  ds <- make_range_dataset()
+
+  p <- plot(ds)
+
+  expect_s3_class(p, "ggplot")
+})
+
+test_that("plot(GCIMSDataset, sample =) restricts which samples are plotted", {
+  ds <- make_range_dataset()
+
+  expect_no_error(plot(ds, sample = "s1"))
+  expect_no_error(plot(ds, sample = 1))
+})
+
+test_that("plot(GCIMSDataset, intensity_range = 'global') errors clearly with remove_baseline = TRUE", {
+  ds <- make_range_dataset()
+
+  expect_error(
+    plot(ds, remove_baseline = TRUE),
+    "global.*remove_baseline"
+  )
+})
+
+test_that("plot(GCIMSDataset, intensity_range = 'ranged') works with remove_baseline = TRUE", {
+  ds <- make_range_dataset()
+  ds <- estimateBaseline(ds, dt_peak_fwhm_ms = 1, dt_region_multiplier = 2, rt_length_s = 2)
+  ds$realize()
+
+  expect_no_error(plot(ds, remove_baseline = TRUE, intensity_range = "ranged"))
+})
+
+test_that("plot(GCIMSDataset, intensity_range = c(min, max)) accepts fixed limits without touching the dataset", {
+  ds <- make_range_dataset()
+
+  expect_no_error(plot(ds, intensity_range = c(0, 200)))
+})
diff --git a/tests/testthat/test-plot-GCIMSSample.R b/tests/testthat/test-plot-GCIMSSample.R
index 4687f3f0..72dac140 100644
--- a/tests/testthat/test-plot-GCIMSSample.R
+++ b/tests/testthat/test-plot-GCIMSSample.R
@@ -81,6 +81,17 @@ test_that("plot(remove_baseline = TRUE) on a GCIMSSample does not error once a b
   expect_s3_class(p, "ggplot")
 })
 
+test_that("plot(fill_range = ) overrides the auto-computed color scale limits", {
+  s <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
+
+  p_auto <- plot(s)
+  p_fixed <- plot(s, fill_range = c(0, 200))
+
+  tr <- cubic_root_trans()
+  expect_equal(p_auto$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 25)))
+  expect_equal(p_fixed$scales$get_scales("fill")$get_limits(), tr$transform(c(0, 200)))
+})
+
 test_that("plot() on a GCIMSSample accepts trans as a string, or as a transform object directly", {
   s <- make_sample()
 

From 5172084e39bde1c03f3daef46f3ada2934a10a6d Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Thu, 16 Jul 2026 19:21:41 +0000
Subject: [PATCH 2/4] Unify plot(GCIMSSample)'s fill_range into intensity_range

plot(GCIMSDataset) already resolved intensity_range via "global"/
"ranged"/fixed/list and passed the result straight into
plot(GCIMSSample, fill_range = ...) -- same knob, two names depending
which layer you were at. Renamed plot(GCIMSSample)'s parameter to
intensity_range and gave it the same resolve_intensity_range()
vocabulary (moved to utils-plot.R now that both methods share it).

This isn't just a rename: "global" vs "ranged" is a real distinction
even for one sample once dt_range/rt_range crop the view. "ranged"
(the new default, matching prior behavior exactly) auto-scales to just
what's shown; "global" scales against the sample's own full range
instead, e.g. to see how strong a cropped region is relative to the
whole sample. GCIMSSample has no cache (deliberately, given the
staleness risk discussed earlier for per-sample derived data), so
"global" is just range(intensity(object)) computed on demand -- cheap
enough for a user-invoked plot call. Errors clearly, same as at the
dataset level, if combined with remove_baseline = TRUE.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019V3CHRGSzcEu3k6tpUFv56
---
 R/plot-GCIMSDataset.R                  | 46 +-------------------------
 R/plot-GCIMSSample.R                   | 39 +++++++++++++++++-----
 R/utils-plot.R                         | 45 +++++++++++++++++++++++++
 man/plot-GCIMSSample-ANY-method.Rd     | 19 ++++++++---
 tests/testthat/test-plot-GCIMSSample.R | 41 +++++++++++++++++++++--
 5 files changed, 130 insertions(+), 60 deletions(-)

diff --git a/R/plot-GCIMSDataset.R b/R/plot-GCIMSDataset.R
index effd6b02..2a93a6b8 100644
--- a/R/plot-GCIMSDataset.R
+++ b/R/plot-GCIMSDataset.R
@@ -1,47 +1,3 @@
-#' Resolve the `intensity_range` argument of `plot(GCIMSDataset)` into `c(min, max)`
-#'
-#' @param intensity_range One of `"global"`, `"ranged"`, a numeric vector of
-#' length 2, or a length-2 vector/list whose elements are independently a
-#' number, `"global"` or `"ranged"`.
-#' @param global_range A zero-argument function returning the cached,
-#' whole-dataset `c(min, max)`. Only called if actually referenced.
-#' @param ranged_range A zero-argument function returning the `c(min, max)`
-#' of exactly what will be plotted. Only called if actually referenced.
-#' @return A numeric vector `c(min, max)`
-#' @noRd
-resolve_intensity_range <- function(intensity_range, global_range, ranged_range) {
-  resolve_endpoint <- function(value, position) {
-    if (is.numeric(value) && length(value) == 1) {
-      return(value)
-    }
-    if (identical(value, "global")) {
-      return(global_range()[position])
-    }
-    if (identical(value, "ranged")) {
-      return(ranged_range()[position])
-    }
-    cli_abort('Each element of intensity_range should be a number, "global" or "ranged"')
-  }
-
-  if (is.character(intensity_range) && length(intensity_range) == 1) {
-    if (intensity_range == "global") {
-      return(global_range())
-    } else if (intensity_range == "ranged") {
-      return(ranged_range())
-    }
-    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
-  }
-
-  if (length(intensity_range) != 2) {
-    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
-  }
-
-  c(
-    resolve_endpoint(intensity_range[[1]], 1L),
-    resolve_endpoint(intensity_range[[2]], 2L)
-  )
-}
-
 #' Topographical plot of several samples of a GCIMSDataset, on a shared color scale
 #'
 #' Several samples are plotted side by side, on the same intensity color
@@ -130,7 +86,7 @@ setMethod(
         samples[[nm]],
         dt_range = dt_range, rt_range = rt_range,
         remove_baseline = remove_baseline, trans = trans,
-        fill_range = limits
+        intensity_range = limits
       ) + ggplot2::labs(title = nm)
     })
 
diff --git a/R/plot-GCIMSSample.R b/R/plot-GCIMSSample.R
index 2224e008..019f4633 100644
--- a/R/plot-GCIMSSample.R
+++ b/R/plot-GCIMSSample.R
@@ -7,10 +7,17 @@
 #' @param remove_baseline Set to `TRUE` to subtract the estimated baseline first
 #' @param trans The transformation to the intensity values. "cubic_root" is the default. "intensity" is also valid.
 #' See the `trans` argument in [ggplot2::continuous_scale()] for other possibilities.
-#' @param fill_range A numeric vector of length 2 with the intensity range the
-#' color scale should span, or `NULL` (the default) to use this sample's own
-#' `range()`. Set it explicitly to make several plots comparable on the same
-#' color scale.
+#' @param intensity_range Controls the color scale limits. One of:
+#' - `"ranged"` (the default): the range of exactly what's plotted (respecting
+#'   `dt_range`, `rt_range` and `remove_baseline`).
+#' - `"global"`: this sample's own full, uncropped, raw intensity range,
+#'   regardless of `dt_range`/`rt_range` -- useful to see how strong a cropped
+#'   region is relative to the whole sample. Not valid together with
+#'   `remove_baseline = TRUE`.
+#' - A numeric vector of length 2, `c(min, max)`: fixed limits. Set this
+#'   explicitly to make several plots comparable on the same color scale.
+#' - A length-2 vector/list whose elements are independently a number,
+#'   `"global"` or `"ranged"`, e.g. `list(min = 0, max = "global")`.
 #' @return A plot of the GCIMSSample
 #' @examples
 #' dummy_obj <-GCIMSSample(
@@ -26,7 +33,7 @@
 setMethod(
   "plot",
   "GCIMSSample",
-  function(x, dt_range = NULL, rt_range = NULL, ..., remove_baseline = FALSE, trans = "cubic_root", fill_range = NULL) {
+  function(x, dt_range = NULL, rt_range = NULL, ..., remove_baseline = FALSE, trans = "cubic_root", intensity_range = "ranged") {
   dt <- dtime(x)
   rt <- rtime(x)
   idx <- dt_rt_range_normalization(dt, rt, dt_range, rt_range)
@@ -36,6 +43,22 @@ setMethod(
     basel <- baseline(x)[idx$dt_logical, idx$rt_logical]
     intmat <- intmat - basel
   }
+
+  get_global_range <- function() {
+    if (isTRUE(remove_baseline)) {
+      cli_abort(
+        c(
+          "intensity_range can't use {.val global} together with {.code remove_baseline = TRUE}",
+          "i" = "The full-sample range is raw intensity. Use {.code intensity_range = \"ranged\"} instead"
+        )
+      )
+    }
+    range(intensity(x))
+  }
+  get_ranged_range <- function() range(intmat)
+
+  limits <- resolve_intensity_range(intensity_range, get_global_range, get_ranged_range)
+
   mat_to_gplot(
     intmat,
     dt_min = idx$dt_ms_min,
@@ -43,12 +66,12 @@ setMethod(
     rt_min = idx$rt_s_min,
     rt_max = idx$rt_s_max,
     trans = trans,
-    fill_range = fill_range
+    intensity_range = limits
   )
 })
 
 
-mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt_max = NULL, trans = "cubic_root", fill_range = NULL) {
+mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt_max = NULL, trans = "cubic_root", intensity_range = NULL) {
   require_pkgs(c("farver", "viridisLite"))
   if (is.null(dt_min)) {
     dt_min <- as.numeric(rownames(intmat)[1L])
@@ -62,7 +85,7 @@ mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt
   if (is.null(rt_max)) {
     rt_max <- as.numeric(colnames(intmat)[ncol(intmat)])
   }
-  minmax <- if (is.null(fill_range)) range(intmat) else fill_range
+  minmax <- if (is.null(intensity_range)) range(intmat) else intensity_range
 
   if (is.character(trans)) {
     trans_func <- paste0(trans, "_trans")
diff --git a/R/utils-plot.R b/R/utils-plot.R
index 41559e5e..3b850fb5 100644
--- a/R/utils-plot.R
+++ b/R/utils-plot.R
@@ -35,6 +35,51 @@ mat_to_nativeRaster <- function(x, colormap, rangex = NULL)  {
 }
 
 
+#' Resolve an `intensity_range` argument (shared by `plot,GCIMSSample-method`
+#' and `plot,GCIMSDataset-method`) into `c(min, max)`
+#'
+#' @param intensity_range One of `"global"`, `"ranged"`, a numeric vector of
+#' length 2, or a length-2 vector/list whose elements are independently a
+#' number, `"global"` or `"ranged"`.
+#' @param global_range A zero-argument function returning the full/uncropped
+#' `c(min, max)`. Only called if actually referenced.
+#' @param ranged_range A zero-argument function returning the `c(min, max)`
+#' of exactly what will be plotted. Only called if actually referenced.
+#' @return A numeric vector `c(min, max)`
+#' @noRd
+resolve_intensity_range <- function(intensity_range, global_range, ranged_range) {
+  resolve_endpoint <- function(value, position) {
+    if (is.numeric(value) && length(value) == 1) {
+      return(value)
+    }
+    if (identical(value, "global")) {
+      return(global_range()[position])
+    }
+    if (identical(value, "ranged")) {
+      return(ranged_range()[position])
+    }
+    cli_abort('Each element of intensity_range should be a number, "global" or "ranged"')
+  }
+
+  if (is.character(intensity_range) && length(intensity_range) == 1) {
+    if (intensity_range == "global") {
+      return(global_range())
+    } else if (intensity_range == "ranged") {
+      return(ranged_range())
+    }
+    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
+  }
+
+  if (length(intensity_range) != 2) {
+    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
+  }
+
+  c(
+    resolve_endpoint(intensity_range[[1]], 1L),
+    resolve_endpoint(intensity_range[[2]], 2L)
+  )
+}
+
 #' Make a plot interactive
 #'
 #' Wraps the `plt` with [plotly::ggplotly()] and sets the `xaxis` and `yaxis`
diff --git a/man/plot-GCIMSSample-ANY-method.Rd b/man/plot-GCIMSSample-ANY-method.Rd
index 766e7e37..5c94c8c3 100644
--- a/man/plot-GCIMSSample-ANY-method.Rd
+++ b/man/plot-GCIMSSample-ANY-method.Rd
@@ -11,7 +11,7 @@
   ...,
   remove_baseline = FALSE,
   trans = "cubic_root",
-  fill_range = NULL
+  intensity_range = "ranged"
 )
 }
 \arguments{
@@ -28,10 +28,19 @@
 \item{trans}{The transformation to the intensity values. "cubic_root" is the default. "intensity" is also valid.
 See the \code{trans} argument in \code{\link[ggplot2:continuous_scale]{ggplot2::continuous_scale()}} for other possibilities.}
 
-\item{fill_range}{A numeric vector of length 2 with the intensity range the
-color scale should span, or \code{NULL} (the default) to use this sample's own
-\code{range()}. Set it explicitly to make several plots comparable on the same
-color scale.}
+\item{intensity_range}{Controls the color scale limits. One of:
+\itemize{
+\item \code{"ranged"} (the default): the range of exactly what's plotted (respecting
+\code{dt_range}, \code{rt_range} and \code{remove_baseline}).
+\item \code{"global"}: this sample's own full, uncropped, raw intensity range,
+regardless of \code{dt_range}/\code{rt_range} -- useful to see how strong a cropped
+region is relative to the whole sample. Not valid together with
+\code{remove_baseline = TRUE}.
+\item A numeric vector of length 2, \code{c(min, max)}: fixed limits. Set this
+explicitly to make several plots comparable on the same color scale.
+\item A length-2 vector/list whose elements are independently a number,
+\code{"global"} or \code{"ranged"}, e.g. \code{list(min = 0, max = "global")}.
+}}
 }
 \value{
 A plot of the GCIMSSample
diff --git a/tests/testthat/test-plot-GCIMSSample.R b/tests/testthat/test-plot-GCIMSSample.R
index 72dac140..eb081ca9 100644
--- a/tests/testthat/test-plot-GCIMSSample.R
+++ b/tests/testthat/test-plot-GCIMSSample.R
@@ -81,17 +81,54 @@ test_that("plot(remove_baseline = TRUE) on a GCIMSSample does not error once a b
   expect_s3_class(p, "ggplot")
 })
 
-test_that("plot(fill_range = ) overrides the auto-computed color scale limits", {
+test_that("plot(intensity_range = ) overrides the auto-computed color scale limits", {
   s <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
 
   p_auto <- plot(s)
-  p_fixed <- plot(s, fill_range = c(0, 200))
+  p_fixed <- plot(s, intensity_range = c(0, 200))
 
   tr <- cubic_root_trans()
   expect_equal(p_auto$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 25)))
   expect_equal(p_fixed$scales$get_scales("fill")$get_limits(), tr$transform(c(0, 200)))
 })
 
+test_that("plot(intensity_range = 'ranged') is the default and scales to the cropped view", {
+  s <- GCIMSSample(drift_time = 1:10, retention_time = 1:10, data = matrix(1:100, nrow = 10))
+
+  p_cropped <- plot(s, dt_range = c(1, 1))
+
+  tr <- cubic_root_trans()
+  expect_equal(p_cropped$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 91)))
+})
+
+test_that("plot(intensity_range = 'global') scales to the sample's full range, ignoring dt_range/rt_range", {
+  s <- GCIMSSample(drift_time = 1:10, retention_time = 1:10, data = matrix(1:100, nrow = 10))
+
+  p_global <- plot(s, dt_range = c(1, 1), intensity_range = "global")
+
+  tr <- cubic_root_trans()
+  expect_equal(p_global$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 100)))
+})
+
+test_that("plot(intensity_range = 'global') errors clearly with remove_baseline = TRUE", {
+  s <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
+  s <- estimateBaseline(s, dt_peak_fwhm_ms = 1, dt_region_multiplier = 2, rt_length_s = 2)
+
+  expect_error(
+    plot(s, remove_baseline = TRUE, intensity_range = "global"),
+    "global.*remove_baseline"
+  )
+})
+
+test_that("plot(intensity_range = list(min=, max=)) resolves each endpoint independently", {
+  s <- GCIMSSample(drift_time = 1:10, retention_time = 1:10, data = matrix(1:100, nrow = 10))
+
+  p <- plot(s, dt_range = c(1, 1), intensity_range = list(min = "ranged", max = "global"))
+
+  tr <- cubic_root_trans()
+  expect_equal(p$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 100)))
+})
+
 test_that("plot() on a GCIMSSample accepts trans as a string, or as a transform object directly", {
   s <- make_sample()
 

From c93a4f4e05fb4b9ca26e25d51a2d053000b8a8ba Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Thu, 16 Jul 2026 20:01:24 +0000
Subject: [PATCH 3/4] Compute global+remove_baseline via two passes instead of
 erroring

Rename the internal cache field IntensityRange -> intensity_range, to
match the intensity_range argument it backs.

intensity_range = "global" combined with remove_baseline = TRUE
previously errored, since the dataset's cache (and the trivial
range(intensity(x)) shortcut on GCIMSSample) both hold raw intensity,
not baseline-removed. Instead of refusing the combination, compute it:
each selected sample's full, uncropped intensity minus its full
baseline, ranged across all of them. Same "global ignores dt_range/
rt_range" semantics as the raw case, just costing an extra pass
instead of being free. Memoized so a list(min=, max=) spec resolving
"global" for both endpoints only computes it once.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019V3CHRGSzcEu3k6tpUFv56
---
 R/aaa-class-GCIMSDataset.R              |  6 +++---
 R/getTIS_getRIC-GCIMSDataset.R          |  8 +++----
 R/plot-GCIMSDataset.R                   | 28 +++++++++++++------------
 R/plot-GCIMSSample.R                    | 18 +++++++---------
 man/GCIMSDataset.Rd                     |  2 +-
 man/plot-GCIMSDataset-ANY-method.Rd     | 11 +++++-----
 man/plot-GCIMSSample-ANY-method.Rd      |  8 +++----
 tests/testthat/test-plot-GCIMSDataset.R | 23 ++++++++++----------
 tests/testthat/test-plot-GCIMSSample.R  | 15 ++++++-------
 9 files changed, 59 insertions(+), 60 deletions(-)

diff --git a/R/aaa-class-GCIMSDataset.R b/R/aaa-class-GCIMSDataset.R
index c95ecefb..8cc58340 100644
--- a/R/aaa-class-GCIMSDataset.R
+++ b/R/aaa-class-GCIMSDataset.R
@@ -79,11 +79,11 @@ GCIMSDataset <- R6::R6Class("GCIMSDataset",
     dt_ref = NULL, # numeric or NULL
     #' @field rt_ref A numeric retention time of reference
     rt_ref = NULL, # numeric or NULL
-    #' @field IntensityRange A matrix of n_samples x c(min, max) with each
+    #' @field intensity_range A matrix of n_samples x c(min, max) with each
     #' sample's raw intensity range, extracted alongside TIS/RIC. Internal
     #' only: no S4 generic exposes this, it only backs `plot(GCIMSDataset)`'s
     #' shared color scale across samples.
-    IntensityRange = NULL, # matrix or NULL
+    intensity_range = NULL, # matrix or NULL
     #' @field userData A list to store arbitrary data in the dataset
     userData = list(), # list
     # Methods:
@@ -226,7 +226,7 @@ GCIMSDataset <- R6::R6Class("GCIMSDataset",
       self$RIC <- NULL
       self$dt_ref <- NULL
       self$rt_ref <- NULL
-      self$IntensityRange <- NULL
+      self$intensity_range <- NULL
       private$delayed_dataset$subset(samples)
       self
     },
diff --git a/R/getTIS_getRIC-GCIMSDataset.R b/R/getTIS_getRIC-GCIMSDataset.R
index 5ba253a2..0d7c3720 100644
--- a/R/getTIS_getRIC-GCIMSDataset.R
+++ b/R/getTIS_getRIC-GCIMSDataset.R
@@ -149,8 +149,8 @@ setMethod(
   }
   stopifnot(nrow(ds$TIS) == num_samples)
   stopifnot(nrow(ds$RIC) == num_samples)
-  ds$IntensityRange <- do.call(rbind, intensity_ranges)
-  colnames(ds$IntensityRange) <- c("min", "max")
+  ds$intensity_range <- do.call(rbind, intensity_ranges)
+  colnames(ds$intensity_range) <- c("min", "max")
   ds
 }
 
@@ -159,11 +159,11 @@ setMethod(
 #' @return A length-2 numeric vector `c(min, max)`
 #' @noRd
 dataset_intensity_range <- function(object) {
-  if (object$hasDelayedOps() || is.null(object$IntensityRange)) {
+  if (object$hasDelayedOps() || is.null(object$intensity_range)) {
     object$extract_RIC_and_TIS()
     object$realize()
   }
-  range(object$IntensityRange)
+  range(object$intensity_range)
 }
 
 #' Extract the Reverse Ion Chromatogram and Total Ion Spectrum from the samples
diff --git a/R/plot-GCIMSDataset.R b/R/plot-GCIMSDataset.R
index 2a93a6b8..cd0868a2 100644
--- a/R/plot-GCIMSDataset.R
+++ b/R/plot-GCIMSDataset.R
@@ -15,11 +15,12 @@
 #' @param intensity_range Controls the shared color scale across all plotted
 #' samples. One of:
 #' - `"global"` (the default): each sample's full, uncropped intensity range,
-#'   cached on the dataset (computed once, alongside TIS/RIC, regardless of
-#'   `sample`/`dt_range`/`rt_range`). Free, and keeps the scale stable across
-#'   different calls with different `sample`/`dt_range`/`rt_range` selections.
-#'   Not valid together with `remove_baseline = TRUE`, since the cache holds
-#'   raw intensities.
+#'   regardless of `sample`/`dt_range`/`rt_range`, keeping the scale stable
+#'   across calls with different selections. With `remove_baseline = FALSE`,
+#'   this comes from a cache computed once, alongside TIS/RIC, so it's free.
+#'   With `remove_baseline = TRUE`, the cache (raw intensities) doesn't apply,
+#'   so each selected sample's full, uncropped, baseline-removed range is
+#'   computed instead (loading every selected sample's data an extra time).
 #' - `"ranged"`: the range of exactly what's plotted (respecting `dt_range`,
 #'   `rt_range` and `remove_baseline`). Requires loading every selected
 #'   sample's data twice (once to compute the range, once to plot).
@@ -67,16 +68,17 @@ setMethod(
       }
       ranged_range_cache
     }
+    global_range_cache <- NULL
     get_global_range <- function() {
-      if (isTRUE(remove_baseline)) {
-        cli_abort(
-          c(
-            "intensity_range can't use {.val global} together with {.code remove_baseline = TRUE}",
-            "i" = "The cached global range holds raw intensities. Use {.code intensity_range = \"ranged\"} instead"
-          )
-        )
+      if (is.null(global_range_cache)) {
+        global_range_cache <<- if (isTRUE(remove_baseline)) {
+          ranges <- purrr::map(samples, function(s) range(intensity(s) - baseline(s)))
+          range(unlist(ranges))
+        } else {
+          dataset_intensity_range(x)
+        }
       }
-      dataset_intensity_range(x)
+      global_range_cache
     }
 
     limits <- resolve_intensity_range(intensity_range, get_global_range, get_ranged_range)
diff --git a/R/plot-GCIMSSample.R b/R/plot-GCIMSSample.R
index 019f4633..42c6ec5c 100644
--- a/R/plot-GCIMSSample.R
+++ b/R/plot-GCIMSSample.R
@@ -10,10 +10,10 @@
 #' @param intensity_range Controls the color scale limits. One of:
 #' - `"ranged"` (the default): the range of exactly what's plotted (respecting
 #'   `dt_range`, `rt_range` and `remove_baseline`).
-#' - `"global"`: this sample's own full, uncropped, raw intensity range,
-#'   regardless of `dt_range`/`rt_range` -- useful to see how strong a cropped
-#'   region is relative to the whole sample. Not valid together with
-#'   `remove_baseline = TRUE`.
+#' - `"global"`: this sample's own full, uncropped intensity range, regardless
+#'   of `dt_range`/`rt_range` -- useful to see how strong a cropped region is
+#'   relative to the whole sample. With `remove_baseline = TRUE`, this is the
+#'   full, uncropped, baseline-removed range.
 #' - A numeric vector of length 2, `c(min, max)`: fixed limits. Set this
 #'   explicitly to make several plots comparable on the same color scale.
 #' - A length-2 vector/list whose elements are independently a number,
@@ -46,14 +46,10 @@ setMethod(
 
   get_global_range <- function() {
     if (isTRUE(remove_baseline)) {
-      cli_abort(
-        c(
-          "intensity_range can't use {.val global} together with {.code remove_baseline = TRUE}",
-          "i" = "The full-sample range is raw intensity. Use {.code intensity_range = \"ranged\"} instead"
-        )
-      )
+      range(intensity(x) - baseline(x))
+    } else {
+      range(intensity(x))
     }
-    range(intensity(x))
   }
   get_ranged_range <- function() range(intmat)
 
diff --git a/man/GCIMSDataset.Rd b/man/GCIMSDataset.Rd
index 0e674dd3..5e2ca407 100644
--- a/man/GCIMSDataset.Rd
+++ b/man/GCIMSDataset.Rd
@@ -96,7 +96,7 @@ dummy_dataset <- GCIMSDataset$new(
 
 \item{\code{rt_ref}}{A numeric retention time of reference}
 
-\item{\code{IntensityRange}}{A matrix of n_samples x c(min, max) with each
+\item{\code{intensity_range}}{A matrix of n_samples x c(min, max) with each
 sample's raw intensity range, extracted alongside TIS/RIC. Internal
 only: no S4 generic exposes this, it only backs \code{plot(GCIMSDataset)}'s
 shared color scale across samples.}
diff --git a/man/plot-GCIMSDataset-ANY-method.Rd b/man/plot-GCIMSDataset-ANY-method.Rd
index 9318d8e7..da1ec7b2 100644
--- a/man/plot-GCIMSDataset-ANY-method.Rd
+++ b/man/plot-GCIMSDataset-ANY-method.Rd
@@ -38,11 +38,12 @@ baseline first}
 samples. One of:
 \itemize{
 \item \code{"global"} (the default): each sample's full, uncropped intensity range,
-cached on the dataset (computed once, alongside TIS/RIC, regardless of
-\code{sample}/\code{dt_range}/\code{rt_range}). Free, and keeps the scale stable across
-different calls with different \code{sample}/\code{dt_range}/\code{rt_range} selections.
-Not valid together with \code{remove_baseline = TRUE}, since the cache holds
-raw intensities.
+regardless of \code{sample}/\code{dt_range}/\code{rt_range}, keeping the scale stable
+across calls with different selections. With \code{remove_baseline = FALSE},
+this comes from a cache computed once, alongside TIS/RIC, so it's free.
+With \code{remove_baseline = TRUE}, the cache (raw intensities) doesn't apply,
+so each selected sample's full, uncropped, baseline-removed range is
+computed instead (loading every selected sample's data an extra time).
 \item \code{"ranged"}: the range of exactly what's plotted (respecting \code{dt_range},
 \code{rt_range} and \code{remove_baseline}). Requires loading every selected
 sample's data twice (once to compute the range, once to plot).
diff --git a/man/plot-GCIMSSample-ANY-method.Rd b/man/plot-GCIMSSample-ANY-method.Rd
index 5c94c8c3..c3fd3364 100644
--- a/man/plot-GCIMSSample-ANY-method.Rd
+++ b/man/plot-GCIMSSample-ANY-method.Rd
@@ -32,10 +32,10 @@ See the \code{trans} argument in \code{\link[ggplot2:continuous_scale]{ggplot2::
 \itemize{
 \item \code{"ranged"} (the default): the range of exactly what's plotted (respecting
 \code{dt_range}, \code{rt_range} and \code{remove_baseline}).
-\item \code{"global"}: this sample's own full, uncropped, raw intensity range,
-regardless of \code{dt_range}/\code{rt_range} -- useful to see how strong a cropped
-region is relative to the whole sample. Not valid together with
-\code{remove_baseline = TRUE}.
+\item \code{"global"}: this sample's own full, uncropped intensity range, regardless
+of \code{dt_range}/\code{rt_range} -- useful to see how strong a cropped region is
+relative to the whole sample. With \code{remove_baseline = TRUE}, this is the
+full, uncropped, baseline-removed range.
 \item A numeric vector of length 2, \code{c(min, max)}: fixed limits. Set this
 explicitly to make several plots comparable on the same color scale.
 \item A length-2 vector/list whose elements are independently a number,
diff --git a/tests/testthat/test-plot-GCIMSDataset.R b/tests/testthat/test-plot-GCIMSDataset.R
index dae40053..3ea18fff 100644
--- a/tests/testthat/test-plot-GCIMSDataset.R
+++ b/tests/testthat/test-plot-GCIMSDataset.R
@@ -49,7 +49,7 @@ test_that("resolve_intensity_range() rejects invalid specs", {
   expect_error(resolve_intensity_range(list("nope", 1), global_range, ranged_range), "should be a number")
 })
 
-# --- IntensityRange caching ----------------------------------------------
+# --- intensity_range caching ----------------------------------------------
 
 make_range_dataset <- function() {
   s1 <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
@@ -57,7 +57,7 @@ make_range_dataset <- function() {
   GCIMSDataset$new_from_list(samples = list(s1 = s1, s2 = s2), on_ram = TRUE, scratch_dir = NULL)
 }
 
-test_that("realizing a dataset populates IntensityRange alongside TIS/RIC", {
+test_that("realizing a dataset populates intensity_range alongside TIS/RIC", {
   old_bpparam <- BiocParallel::bpparam()
   BiocParallel::register(BiocParallel::SerialParam())
   on.exit(BiocParallel::register(old_bpparam))
@@ -65,22 +65,22 @@ test_that("realizing a dataset populates IntensityRange alongside TIS/RIC", {
   ds <- make_range_dataset()
   ds$realize()
 
-  expect_equal(unname(ds$IntensityRange["s1", ]), c(1, 25))
-  expect_equal(unname(ds$IntensityRange["s2", ]), c(100, 124))
+  expect_equal(unname(ds$intensity_range["s1", ]), c(1, 25))
+  expect_equal(unname(ds$intensity_range["s2", ]), c(100, 124))
 })
 
-test_that("subsetting a dataset resets the cached IntensityRange", {
+test_that("subsetting a dataset resets the cached intensity_range", {
   old_bpparam <- BiocParallel::bpparam()
   BiocParallel::register(BiocParallel::SerialParam())
   on.exit(BiocParallel::register(old_bpparam))
 
   ds <- make_range_dataset()
   ds$realize()
-  expect_false(is.null(ds$IntensityRange))
+  expect_false(is.null(ds$intensity_range))
 
   ds$subset("s1", inplace = TRUE)
 
-  expect_null(ds$IntensityRange)
+  expect_null(ds$intensity_range)
 })
 
 # --- plot(GCIMSDataset) ----------------------------------------------------
@@ -100,13 +100,12 @@ test_that("plot(GCIMSDataset, sample =) restricts which samples are plotted", {
   expect_no_error(plot(ds, sample = 1))
 })
 
-test_that("plot(GCIMSDataset, intensity_range = 'global') errors clearly with remove_baseline = TRUE", {
+test_that("plot(GCIMSDataset, intensity_range = 'global') computes a baseline-removed range with remove_baseline = TRUE", {
   ds <- make_range_dataset()
+  ds <- estimateBaseline(ds, dt_peak_fwhm_ms = 1, dt_region_multiplier = 2, rt_length_s = 2)
+  ds$realize()
 
-  expect_error(
-    plot(ds, remove_baseline = TRUE),
-    "global.*remove_baseline"
-  )
+  expect_no_error(plot(ds, remove_baseline = TRUE))
 })
 
 test_that("plot(GCIMSDataset, intensity_range = 'ranged') works with remove_baseline = TRUE", {
diff --git a/tests/testthat/test-plot-GCIMSSample.R b/tests/testthat/test-plot-GCIMSSample.R
index eb081ca9..fd3459f4 100644
--- a/tests/testthat/test-plot-GCIMSSample.R
+++ b/tests/testthat/test-plot-GCIMSSample.R
@@ -110,14 +110,15 @@ test_that("plot(intensity_range = 'global') scales to the sample's full range, i
   expect_equal(p_global$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 100)))
 })
 
-test_that("plot(intensity_range = 'global') errors clearly with remove_baseline = TRUE", {
-  s <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
-  s <- estimateBaseline(s, dt_peak_fwhm_ms = 1, dt_region_multiplier = 2, rt_length_s = 2)
+test_that("plot(intensity_range = 'global') computes the full, uncropped, baseline-removed range with remove_baseline = TRUE", {
+  s <- make_sample()
+  s <- estimateBaseline(s, dt_peak_fwhm_ms = 0.3, dt_region_multiplier = 6, rt_length_s = 10, remove = FALSE)
 
-  expect_error(
-    plot(s, remove_baseline = TRUE, intensity_range = "global"),
-    "global.*remove_baseline"
-  )
+  p <- plot(s, dt_range = c(0, 1), remove_baseline = TRUE, intensity_range = "global")
+
+  expected_range <- range(intensity(s) - baseline(s))
+  tr <- cubic_root_trans()
+  expect_equal(p$scales$get_scales("fill")$get_limits(), tr$transform(expected_range))
 })
 
 test_that("plot(intensity_range = list(min=, max=)) resolves each endpoint independently", {

From a08f1bb6e54db729e23576fe9753184246ca540c Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 17 Jul 2026 16:12:54 +0000
Subject: [PATCH 4/4] Add pagination to plot(GCIMSDataset): nrow, ncol, page

Replaces the unbounded "load and render every selected sample at once"
default with paginated rendering, per discussion:

- nrow/ncol define page capacity (nrow * ncol), replacing the earlier
  standalone ncol-only parameter -- matches ggforce::facet_wrap_paginate()'s
  convention rather than inventing a separate samples_per_page. When
  neither is given, both are picked from the sample count: an exact fit
  for 1-6 samples, otherwise capped at 3x3 (resolve_page_grid(), new in
  utils-plot.R). When exactly one is given, the other defaults to 3 if
  there are more than 9 samples, or just enough to fit everyone
  otherwise.
- page (default 1) selects which page to render; out-of-bounds pages
  error clearly instead of silently clamping or rendering nothing.
- sample= filters first (as before), pagination then slices that
  filtered set into pages -- no special-casing needed, just two
  sequential steps.
- Rendering now loads only the current page's samples, one at a time,
  discarding each raw sample right after building its (much smaller,
  native-raster-encoded) panel -- so peak memory during rendering is
  bounded by page size, not dataset size.
- intensity_range = "ranged" deliberately stays scoped to every
  selected sample across all pages, not just the current one: a
  page-local range would make pages incomparable to each other,
  defeating the point of paginating through a dataset to compare
  samples. This does mean it still costs a full pass over every
  selected sample's data (loaded and discarded one at a time, same
  memory bound as rendering) -- accepted as necessary for fair
  cross-page comparison. intensity_range = "global" is unaffected:
  still the cached, whole-dataset, page-independent range.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019V3CHRGSzcEu3k6tpUFv56
---
 R/plot-GCIMSDataset.R                   | 85 ++++++++++++++++++-------
 R/utils-plot.R                          | 36 +++++++++++
 man/plot-GCIMSDataset-ANY-method.Rd     | 40 ++++++++----
 tests/testthat/test-plot-GCIMSDataset.R | 79 +++++++++++++++++++++++
 4 files changed, 205 insertions(+), 35 deletions(-)

diff --git a/R/plot-GCIMSDataset.R b/R/plot-GCIMSDataset.R
index cd0868a2..9e8f0d30 100644
--- a/R/plot-GCIMSDataset.R
+++ b/R/plot-GCIMSDataset.R
@@ -1,52 +1,80 @@
 #' Topographical plot of several samples of a GCIMSDataset, on a shared color scale
 #'
 #' Several samples are plotted side by side, on the same intensity color
-#' scale, so they are visually comparable.
+#' scale, so they are visually comparable. Paginated: only one page's worth
+#' of samples is ever loaded into memory at a time for rendering.
 #'
 #' @param x A [GCIMSDataset] object
 #' @param sample A number, a string, or a vector of numbers/strings with the
 #' sample index(es) or name(s) to plot. If `NULL` (the default), all samples
-#' are plotted.
+#' are plotted (across as many pages as needed).
 #' @inheritParams dt_rt_range_normalization
 #' @param ... Ignored
 #' @param remove_baseline Set to `TRUE` to subtract each sample's estimated
 #' baseline first
 #' @param trans The transformation to the intensity values, see `plot,GCIMSSample-method`
 #' @param intensity_range Controls the shared color scale across all plotted
-#' samples. One of:
+#' samples (all pages, not just the current one). One of:
 #' - `"global"` (the default): each sample's full, uncropped intensity range,
-#'   regardless of `sample`/`dt_range`/`rt_range`, keeping the scale stable
-#'   across calls with different selections. With `remove_baseline = FALSE`,
-#'   this comes from a cache computed once, alongside TIS/RIC, so it's free.
-#'   With `remove_baseline = TRUE`, the cache (raw intensities) doesn't apply,
-#'   so each selected sample's full, uncropped, baseline-removed range is
-#'   computed instead (loading every selected sample's data an extra time).
+#'   regardless of `sample`/`dt_range`/`rt_range`. With `remove_baseline =
+#'   FALSE`, this comes from a cache computed once, alongside TIS/RIC, so
+#'   it's free. With `remove_baseline = TRUE`, the cache (raw intensities)
+#'   doesn't apply, so each selected sample's full, uncropped, baseline-removed
+#'   range is computed instead (loading every selected sample's data an extra
+#'   time).
 #' - `"ranged"`: the range of exactly what's plotted (respecting `dt_range`,
-#'   `rt_range` and `remove_baseline`). Requires loading every selected
+#'   `rt_range` and `remove_baseline`), across *every* selected sample, not
+#'   just the current page -- so the scale is still comparable page to page.
+#'   This is deliberately not page-scoped: a cheaper, page-local range would
+#'   make pages incomparable to each other, defeating the purpose of paging
+#'   through a dataset to compare samples. Costs loading every selected
 #'   sample's data twice (once to compute the range, once to plot).
 #' - A numeric vector of length 2, `c(min, max)`: fixed limits.
 #' - A length-2 vector/list whose elements are independently a number,
 #'   `"global"` or `"ranged"`, e.g. `list(min = 0, max = "global")`.
-#' @param ncol Number of columns for [cowplot::plot_grid()]. `NULL` lets `cowplot` decide.
-#' @return A combined plot with one panel per sample
+#' @param nrow,ncol Page grid shape, forwarded to [cowplot::plot_grid()]. Page
+#' capacity is `nrow * ncol`. If both are `NULL` (the default), they're picked
+#' from the number of selected samples: an exact fit for 1-6 samples (e.g. 5
+#' or 6 samples -> 2x3), otherwise 3x3 (so page capacity maxes out at 9, and
+#' more than 9 selected samples span multiple pages). If only one of
+#' `nrow`/`ncol` is given, the other becomes `3` when there are more than 9
+#' selected samples, or just enough to fit them all on one page otherwise.
+#' @param page Which page to plot, 1-based. Errors if out of bounds.
+#' @return A combined plot with one panel per sample, for the requested page
 #' @export
 setMethod(
   "plot",
   "GCIMSDataset",
   function(x, sample = NULL, dt_range = NULL, rt_range = NULL, ...,
            remove_baseline = FALSE, trans = "cubic_root",
-           intensity_range = "global", ncol = NULL) {
+           intensity_range = "global", nrow = NULL, ncol = NULL, page = 1) {
     require_pkgs("cowplot")
     sample_names_all <- sampleNames(x)
     if (is.null(sample)) {
       sample <- sample_names_all
     }
     selected_names <- sample_name_or_number_to_both(sample, sample_names_all)$name
+    num_samples <- length(selected_names)
+    if (num_samples == 0) {
+      cli_abort("No samples to plot")
+    }
 
-    samples <- stats::setNames(
-      purrr::map(selected_names, function(nm) x$getSample(nm)),
-      selected_names
-    )
+    grid_dims <- resolve_page_grid(nrow, ncol, num_samples)
+    nrow <- grid_dims[["nrow"]]
+    ncol <- grid_dims[["ncol"]]
+    page_size <- nrow * ncol
+    n_pages <- ceiling(num_samples / page_size)
+    if (page < 1 || page > n_pages) {
+      cli_abort(
+        c(
+          "{.arg page} = {page} is out of bounds",
+          "i" = "With {num_samples} sample(s) laid out on a {nrow}x{ncol} grid, there are {n_pages} page(s)"
+        )
+      )
+    }
+    page_start <- (page - 1) * page_size + 1
+    page_end <- min(page * page_size, num_samples)
+    page_names <- selected_names[page_start:page_end]
 
     cropped_intensity <- function(s) {
       dt <- dtime(s)
@@ -63,7 +91,13 @@ setMethod(
     ranged_range_cache <- NULL
     get_ranged_range <- function() {
       if (is.null(ranged_range_cache)) {
-        ranges <- purrr::map(samples, function(s) range(cropped_intensity(s)))
+        # Scoped to every selected sample (all pages), not just the current
+        # one, so the scale stays comparable across pages. Each sample is
+        # loaded, ranged and discarded one at a time to avoid holding them
+        # all in RAM at once.
+        ranges <- purrr::map(selected_names, function(nm) {
+          range(cropped_intensity(x$getSample(nm)))
+        })
         ranged_range_cache <<- range(unlist(ranges))
       }
       ranged_range_cache
@@ -72,7 +106,10 @@ setMethod(
     get_global_range <- function() {
       if (is.null(global_range_cache)) {
         global_range_cache <<- if (isTRUE(remove_baseline)) {
-          ranges <- purrr::map(samples, function(s) range(intensity(s) - baseline(s)))
+          ranges <- purrr::map(selected_names, function(nm) {
+            s <- x$getSample(nm)
+            range(intensity(s) - baseline(s))
+          })
           range(unlist(ranges))
         } else {
           dataset_intensity_range(x)
@@ -83,15 +120,19 @@ setMethod(
 
     limits <- resolve_intensity_range(intensity_range, get_global_range, get_ranged_range)
 
-    panels <- purrr::map(selected_names, function(nm) {
+    # Only the current page's samples are loaded for rendering, one at a
+    # time, so peak memory never holds more than one raw sample plus the
+    # (much smaller, native-raster-encoded) panels built so far.
+    panels <- purrr::map(page_names, function(nm) {
+      s <- x$getSample(nm)
       plot(
-        samples[[nm]],
+        s,
         dt_range = dt_range, rt_range = rt_range,
         remove_baseline = remove_baseline, trans = trans,
         intensity_range = limits
       ) + ggplot2::labs(title = nm)
     })
 
-    cowplot::plot_grid(plotlist = panels, ncol = ncol)
+    cowplot::plot_grid(plotlist = panels, nrow = nrow, ncol = ncol)
   }
 )
diff --git a/R/utils-plot.R b/R/utils-plot.R
index 3b850fb5..326b456a 100644
--- a/R/utils-plot.R
+++ b/R/utils-plot.R
@@ -35,6 +35,42 @@ mat_to_nativeRaster <- function(x, colormap, rangex = NULL)  {
 }
 
 
+#' Resolve the (`nrow`, `ncol`) page grid for `plot,GCIMSDataset-method`,
+#' auto-filling whichever of `nrow`/`ncol` isn't given
+#'
+#' @param nrow Number of rows, or `NULL`
+#' @param ncol Number of columns, or `NULL`
+#' @param num_samples Number of samples to lay out (before pagination)
+#' @return A list `list(nrow =, ncol =)`
+#' @noRd
+resolve_page_grid <- function(nrow, ncol, num_samples) {
+  if (is.null(nrow) && is.null(ncol)) {
+    dims <- if (num_samples <= 1) {
+      c(1L, 1L)
+    } else if (num_samples == 2) {
+      c(1L, 2L)
+    } else if (num_samples == 3) {
+      c(1L, 3L)
+    } else if (num_samples == 4) {
+      c(2L, 2L)
+    } else if (num_samples <= 6) {
+      c(2L, 3L)
+    } else {
+      c(3L, 3L)
+    }
+    return(list(nrow = dims[1], ncol = dims[2]))
+  }
+  if (is.null(nrow)) {
+    nrow <- if (num_samples > 9) 3L else ceiling(num_samples / ncol)
+    return(list(nrow = nrow, ncol = ncol))
+  }
+  if (is.null(ncol)) {
+    ncol <- if (num_samples > 9) 3L else ceiling(num_samples / nrow)
+    return(list(nrow = nrow, ncol = ncol))
+  }
+  list(nrow = nrow, ncol = ncol)
+}
+
 #' Resolve an `intensity_range` argument (shared by `plot,GCIMSSample-method`
 #' and `plot,GCIMSDataset-method`) into `c(min, max)`
 #'
diff --git a/man/plot-GCIMSDataset-ANY-method.Rd b/man/plot-GCIMSDataset-ANY-method.Rd
index da1ec7b2..6685ebea 100644
--- a/man/plot-GCIMSDataset-ANY-method.Rd
+++ b/man/plot-GCIMSDataset-ANY-method.Rd
@@ -13,7 +13,9 @@
   remove_baseline = FALSE,
   trans = "cubic_root",
   intensity_range = "global",
-  ncol = NULL
+  nrow = NULL,
+  ncol = NULL,
+  page = 1
 )
 }
 \arguments{
@@ -21,7 +23,7 @@
 
 \item{sample}{A number, a string, or a vector of numbers/strings with the
 sample index(es) or name(s) to plot. If \code{NULL} (the default), all samples
-are plotted.}
+are plotted (across as many pages as needed).}
 
 \item{dt_range}{The minimum and maximum drift times to extract (length 2 vector)}
 
@@ -35,29 +37,41 @@ baseline first}
 \item{trans}{The transformation to the intensity values, see \verb{plot,GCIMSSample-method}}
 
 \item{intensity_range}{Controls the shared color scale across all plotted
-samples. One of:
+samples (all pages, not just the current one). One of:
 \itemize{
 \item \code{"global"} (the default): each sample's full, uncropped intensity range,
-regardless of \code{sample}/\code{dt_range}/\code{rt_range}, keeping the scale stable
-across calls with different selections. With \code{remove_baseline = FALSE},
-this comes from a cache computed once, alongside TIS/RIC, so it's free.
-With \code{remove_baseline = TRUE}, the cache (raw intensities) doesn't apply,
-so each selected sample's full, uncropped, baseline-removed range is
-computed instead (loading every selected sample's data an extra time).
+regardless of \code{sample}/\code{dt_range}/\code{rt_range}. With \code{remove_baseline = FALSE}, this comes from a cache computed once, alongside TIS/RIC, so
+it's free. With \code{remove_baseline = TRUE}, the cache (raw intensities)
+doesn't apply, so each selected sample's full, uncropped, baseline-removed
+range is computed instead (loading every selected sample's data an extra
+time).
 \item \code{"ranged"}: the range of exactly what's plotted (respecting \code{dt_range},
-\code{rt_range} and \code{remove_baseline}). Requires loading every selected
+\code{rt_range} and \code{remove_baseline}), across \emph{every} selected sample, not
+just the current page -- so the scale is still comparable page to page.
+This is deliberately not page-scoped: a cheaper, page-local range would
+make pages incomparable to each other, defeating the purpose of paging
+through a dataset to compare samples. Costs loading every selected
 sample's data twice (once to compute the range, once to plot).
 \item A numeric vector of length 2, \code{c(min, max)}: fixed limits.
 \item A length-2 vector/list whose elements are independently a number,
 \code{"global"} or \code{"ranged"}, e.g. \code{list(min = 0, max = "global")}.
 }}
 
-\item{ncol}{Number of columns for \code{\link[cowplot:plot_grid]{cowplot::plot_grid()}}. \code{NULL} lets \code{cowplot} decide.}
+\item{nrow, ncol}{Page grid shape, forwarded to \code{\link[cowplot:plot_grid]{cowplot::plot_grid()}}. Page
+capacity is \code{nrow * ncol}. If both are \code{NULL} (the default), they're picked
+from the number of selected samples: an exact fit for 1-6 samples (e.g. 5
+or 6 samples -> 2x3), otherwise 3x3 (so page capacity maxes out at 9, and
+more than 9 selected samples span multiple pages). If only one of
+\code{nrow}/\code{ncol} is given, the other becomes \code{3} when there are more than 9
+selected samples, or just enough to fit them all on one page otherwise.}
+
+\item{page}{Which page to plot, 1-based. Errors if out of bounds.}
 }
 \value{
-A combined plot with one panel per sample
+A combined plot with one panel per sample, for the requested page
 }
 \description{
 Several samples are plotted side by side, on the same intensity color
-scale, so they are visually comparable.
+scale, so they are visually comparable. Paginated: only one page's worth
+of samples is ever loaded into memory at a time for rendering.
 }
diff --git a/tests/testthat/test-plot-GCIMSDataset.R b/tests/testthat/test-plot-GCIMSDataset.R
index 3ea18fff..5e16915b 100644
--- a/tests/testthat/test-plot-GCIMSDataset.R
+++ b/tests/testthat/test-plot-GCIMSDataset.R
@@ -49,6 +49,34 @@ test_that("resolve_intensity_range() rejects invalid specs", {
   expect_error(resolve_intensity_range(list("nope", 1), global_range, ranged_range), "should be a number")
 })
 
+# --- resolve_page_grid() --------------------------------------------------
+
+test_that("resolve_page_grid() with neither nrow nor ncol given follows the exact-fit-then-3x3-cap table", {
+  expected <- list(
+    `1` = c(1, 1), `2` = c(1, 2), `3` = c(1, 3), `4` = c(2, 2),
+    `5` = c(2, 3), `6` = c(2, 3), `7` = c(3, 3), `9` = c(3, 3), `25` = c(3, 3)
+  )
+  for (n in names(expected)) {
+    g <- resolve_page_grid(NULL, NULL, as.integer(n))
+    expect_equal(c(g$nrow, g$ncol), expected[[n]], info = n)
+  }
+})
+
+test_that("resolve_page_grid() with exactly one of nrow/ncol given fills the other", {
+  # num_samples <= 9: fit everyone on one page
+  expect_equal(resolve_page_grid(NULL, 2, 5)$nrow, 3)
+  expect_equal(resolve_page_grid(4, NULL, 9)$ncol, 3)
+
+  # num_samples > 9: missing dimension is fixed at 3, regardless of the given one
+  expect_equal(resolve_page_grid(NULL, 5, 20)$nrow, 3)
+  expect_equal(resolve_page_grid(7, NULL, 50)$ncol, 3)
+})
+
+test_that("resolve_page_grid() with both given uses them as-is", {
+  g <- resolve_page_grid(4, 5, 100)
+  expect_equal(c(g$nrow, g$ncol), c(4, 5))
+})
+
 # --- intensity_range caching ----------------------------------------------
 
 make_range_dataset <- function() {
@@ -121,3 +149,54 @@ test_that("plot(GCIMSDataset, intensity_range = c(min, max)) accepts fixed limit
 
   expect_no_error(plot(ds, intensity_range = c(0, 200)))
 })
+
+# --- pagination -------------------------------------------------------------
+
+make_paged_dataset <- function(n = 11) {
+  samples <- stats::setNames(
+    purrr::map(seq_len(n), function(i) {
+      GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(rep(i, 25), nrow = 5))
+    }),
+    paste0("s", seq_len(n))
+  )
+  GCIMSDataset$new_from_list(samples = samples, on_ram = TRUE, scratch_dir = NULL)
+}
+
+test_that("plot(GCIMSDataset) with more than 9 samples defaults to a 3x3 first page", {
+  ds <- make_paged_dataset(11)
+
+  expect_no_error(plot(ds))
+  expect_no_error(plot(ds, page = 1))
+  expect_no_error(plot(ds, page = 2))
+})
+
+test_that("plot(GCIMSDataset, page =) errors clearly when out of bounds", {
+  ds <- make_paged_dataset(11)
+
+  expect_error(plot(ds, page = 3), "page.*3.*out of bounds")
+  expect_error(plot(ds, page = 0), "page.*0.*out of bounds")
+})
+
+test_that("plot(GCIMSDataset, intensity_range = 'ranged') is computed over every selected sample, not just the current page", {
+  # 11 samples, values 1..11 -> page 1 (samples 1-9) would only see 1..9 if
+  # ranged were page-scoped, but sample 11 (value 11) is on page 2.
+  ds <- make_paged_dataset(11)
+
+  captured <- list()
+  testthat::local_mocked_bindings(
+    mat_to_gplot = function(intmat, ..., intensity_range = NULL) {
+      captured[[length(captured) + 1]] <<- intensity_range
+      ggplot2::ggplot()
+    },
+    .package = "GCIMS"
+  )
+
+  plot(ds, intensity_range = "ranged", page = 1)
+  range_page1 <- captured[[1]]
+  captured <<- list()
+  plot(ds, intensity_range = "ranged", page = 2)
+  range_page2 <- captured[[1]]
+
+  expect_equal(range_page1, c(1, 11))
+  expect_equal(range_page2, c(1, 11))
+})