diff --git a/CHANGELOG.md b/CHANGELOG.md index 6b36f3f..43f0592 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,8 @@ new refusal or a numerical-range fix (see the pull requests for the exact messag different pixel types are refused before pooling, `Parameters.select` of a strict subset of R, G and B is refused, and `apply_rgb` uses one float32 copy (#78). - `range_check_issues` refuses malformed channels and mixed settings (#79). +- New `iter_image_summaries(records)`: streaming, bounded-memory per-image summaries for rows + grouped by image (contiguity required, otherwise refused); `image_summaries` is unchanged (#58). ## Release notes moved from the README diff --git a/docs/contract.md b/docs/contract.md index dca7379..0513c57 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -601,6 +601,16 @@ library writes no file (N50). The columns are the shared columns `source_relativ | `Intensity_SaturatedFraction_C` | Pixels equal to the type's maximum over valid pixels, for integer types whose histogram holds the maximum value in a bin of its own; blank otherwise (float images, or inexact top bins) | | `Histogram_Exact_C` | Whether every pooled histogram of the channel was exact; blank (`None`) for a channel without finite pixels, which has no statistics to be exact | +`iter_image_summaries(records, percentile=(1.0, 99.8))` is the streaming form for large tables +(issue #58): a generator that holds one image's planes at a time and yields the same +`ImageSummary` that `image_summaries` gives for that image. It requires every Z and T row of an +image to be contiguous in `records`; it never sorts or stages rows, and refuses (`ValueError`) +an image whose rows are not together, a source file under two sets of labels, repeated +`plane_id`s and images with different numbers of channels. Summaries come in the order the +images first appear (the eager order for rows sorted by condition, repeat and path). Channel +counts are only compared as images are read, so an error can follow summaries already yielded: +a run that raised is incomplete and its output must be discarded. `image_summaries` is unchanged. + The checks below call these values `median`, `p_low`, `p_high` and `saturated_fraction`, the attributes of `ImageSummary` that hold one value per channel. A channel without finite pixels has blank values and is not checked; `unchecked_channels(summaries)` returns one diff --git a/src/morphobricks_normalization/__init__.py b/src/morphobricks_normalization/__init__.py index 81f5c1e..bf6a650 100644 --- a/src/morphobricks_normalization/__init__.py +++ b/src/morphobricks_normalization/__init__.py @@ -48,6 +48,7 @@ ImageSummary, IncompatibleSummariesError, image_summaries, + iter_image_summaries, summary_columns, ) @@ -85,6 +86,7 @@ "group_records", "histogram_with_range_check", "image_summaries", + "iter_image_summaries", "intensity_checks", "make_plane_id", "min_max", diff --git a/src/morphobricks_normalization/summaries.py b/src/morphobricks_normalization/summaries.py index 064ff77..78f1da1 100644 --- a/src/morphobricks_normalization/summaries.py +++ b/src/morphobricks_normalization/summaries.py @@ -19,8 +19,9 @@ from __future__ import annotations -from collections.abc import Iterable +from collections.abc import Iterable, Iterator from dataclasses import dataclass +from itertools import groupby from typing import Any import numpy as np @@ -187,9 +188,7 @@ def image_summaries( IncompatibleSummariesError For an image whose Z or T planes have different pixel types (issue #54). """ - low, high = (float(q) for q in percentile) - if not 0 <= low < high <= 100: - raise ValueError(f"Expected 0 <= low < high <= 100; got low={low}, high={high}.") + low, high = _check_percentile(percentile) records = list(records) duplicate_planes(records, "These rows") if len({r.channel_count for r in records}) > 1: @@ -204,6 +203,102 @@ def image_summaries( return [_summary(planes, (low, high)) for _, planes in sorted(images.items())] +def _check_percentile(percentile: tuple[float, float]) -> tuple[float, float]: + low, high = (float(q) for q in percentile) + if not 0 <= low < high <= 100: + raise ValueError(f"Expected 0 <= low < high <= 100; got low={low}, high={high}.") + return low, high + + +def iter_image_summaries( + records: Iterable[PlaneRecord], percentile: tuple[float, float] = DEFAULT_PERCENTILE +) -> Iterator[ImageSummary]: + """Summarize images one at a time from rows grouped by image, in bounded memory (issue #58). + + The streaming counterpart of :func:`image_summaries`: it holds only the planes of the + image being read, emits its :class:`ImageSummary` when the next image starts, and drops + its histograms. Each summary equals the one :func:`image_summaries` gives for that image. + + Precondition: all Z and T rows of an image are contiguous (next to each other in + ``records``). Rows are never sorted or staged here, so input that breaks the precondition is + refused, not summarized from part of an image. Summaries come in the order the images first + appear, which is the order of :func:`image_summaries` when the rows are already sorted by + condition, biological repeat, technical repeat and path. + + Parameters + ---------- + records + Rows of the statistics file (one XYC plane each), grouped by image. + percentile + ``(low, high)`` in percent, ``0 <= low < high <= 100``; default ``(1.0, 99.8)``. + + Yields + ------ + ImageSummary + One per image. + + Raises + ------ + ValueError + At the call, for invalid percentiles. While iterating: for an image whose rows are not + contiguous, a source file listed under two sets of labels, rows that repeat a + ``plane_id``, or images with different numbers of channels (N43). + IncompatibleSummariesError + For an image whose Z or T planes have different pixel types (issue #54). + + Notes + ----- + Whether all images have the same number of channels is only known as they are read, so an + error can come after earlier summaries were emitted. Those summaries are valid on their + own, but a run that ended with an error is incomplete: discard its output; do not publish + it as the result for the experiment. :func:`image_summaries` validates every row before + returning anything. + """ + low, high = _check_percentile(percentile) + return _stream_summaries(records, (low, high)) + + +def _stream_summaries( + records: Iterable[PlaneRecord], percentile: tuple[float, float] +) -> Iterator[ImageSummary]: + """Group rows by adjacent image, check each group and yield its summary.""" + seen_keys: set[tuple[str, ...]] = set() + labels_of_path: dict[str, tuple[str, ...]] = {} + first_channels: tuple[str, int] | None = None + for key, group in groupby(records, key=lambda r: r.group_key("image")): + planes = list(group) + path = key[-1] + if key in seen_keys: + raise ValueError( + f"The rows of image {path} are not together: other images come between its " + "planes. The streaming summary reads one image at a time and never reorders " + "rows, so it cannot pool them. Sort the rows by condition, biological repeat, " + "technical repeat and path, or use image_summaries()." + ) + seen_keys.add(key) + if labels_of_path.setdefault(path, key[:-1]) != key[:-1]: + raise ValueError( + f"The source file {path} appears under different condition or repeat labels, " + "so which image its planes belong to is not defined. Correct the labels, or run " + "the preprocessing brick again on this experiment." + ) + duplicate_planes(planes, "These rows") + counts = {r.channel_count for r in planes} + if first_channels is None: + first_channels = (path, planes[0].channel_count) + counts.add(first_channels[1]) + if len(counts) > 1: + raise ValueError( + different_channel_counts( + [(first_channels[0], first_channels[1])] + + [(r.source_relative_path, r.channel_count) for r in planes], + "These images", + ) + ) + yield _summary(planes, percentile) + del planes + + def _check_one_dtype(planes: list[PlaneRecord]) -> None: """Refuse an image whose Z/T planes have different pixel types (issue #54).""" by_dtype: dict[str, list[str]] = {} diff --git a/tests/test_summaries.py b/tests/test_summaries.py index 80ff6e8..01a79ad 100644 --- a/tests/test_summaries.py +++ b/tests/test_summaries.py @@ -1,6 +1,8 @@ """Per-image intensity summaries from plane_histograms.csv rows (contract ยง7.1, N50).""" +import gc import re +import weakref import numpy as np import pytest @@ -11,6 +13,7 @@ PlaneRecord, compute_plane_histogram, image_summaries, + iter_image_summaries, read_plane_histograms, summary_columns, write_plane_histograms, @@ -245,3 +248,85 @@ def test_different_images_may_differ_in_pixel_type(): record("b.tif", [np.full((8, 8), 5, np.uint16)]), ] assert [s.dtype for s in image_summaries(records)] == [np.uint8, np.uint16] + + +# --- streaming API (issue #58) ------------------------------------------------------------- + + +def grouped_records(n_images=4, planes_per_image=3, channels=1, seed=20): + """Rows grouped by image and sorted as image_summaries sorts its output.""" + rng = np.random.default_rng(seed) + rows = [] + for i in range(n_images): + labels = ("control", "rep1", SINGLE) + for z in range(planes_per_image): + planes = [simulate_plane(rng, np.uint16, (10, 10)) for _ in range(channels)] + rows.append(record(f"img{i}.tif", planes, z=z, labels=labels)) + return rows + + +def test_streaming_summaries_equal_the_eager_ones(): + rows = grouped_records(channels=2) + assert list(iter_image_summaries(rows)) == image_summaries(rows) + assert list(iter_image_summaries(rows, (5, 95))) == image_summaries(rows, (5, 95)) + + +def test_streaming_emits_one_image_at_a_time_and_drops_its_histograms(): + rows = grouped_records(n_images=3) + pulled = [] + first_hist = weakref.ref(rows[0].histograms[0]) + + def source(): + for row in rows: + pulled.append(row.plane_id) + yield row + + stream = iter_image_summaries(source()) + next(stream) + # The first summary is out once the second image has started, not after all rows. + assert len(pulled) == 4 < len(rows) + del rows[:3] # the caller's own list drops image 0; the stream must not keep it + gc.collect() + next(stream) + gc.collect() + assert first_hist() is None + + +def test_streaming_refuses_rows_that_are_not_together(): + rows = grouped_records(n_images=2) + mixed = [rows[0], rows[3], rows[1], rows[4], rows[2], rows[5]] + with pytest.raises(ValueError, match="not together"): + list(iter_image_summaries(mixed)) + # The eager API still accepts any order. + assert len(image_summaries(mixed)) == 2 + + +def test_streaming_refuses_a_file_under_two_label_sets(): + plane = simulate_plane(np.random.default_rng(1), np.uint16, (8, 8)) + rows = [ + record("a.tif", [plane], z=0, labels=("control", "rep1", SINGLE)), + record("a.tif", [plane], z=1, labels=("treated", "rep1", SINGLE)), + ] + with pytest.raises(ValueError, match="different condition or repeat labels"): + list(iter_image_summaries(rows)) + + +def test_streaming_refuses_repeated_planes_and_different_channel_counts(): + rows = grouped_records(n_images=1) + with pytest.raises(ValueError, match="more than once"): + list(iter_image_summaries([*rows, rows[0]])) + one = grouped_records(n_images=1, channels=1) + two = [ + r for r in grouped_records(n_images=2, channels=2) if r.source_relative_path == "img1.tif" + ] + with pytest.raises(ValueError, match="different numbers of channels"): + list(iter_image_summaries([*one, *two])) + + +def test_streaming_validates_percentiles_at_the_call(): + with pytest.raises(ValueError, match="0 <= low < high <= 100"): + iter_image_summaries([], percentile=(99.8, 1.0)) + + +def test_streaming_no_records_no_summaries(): + assert list(iter_image_summaries([])) == []