From 77a95770d6661b9243642bc008df554f4ded47a4 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 3 Nov 2026 14:02:06 +0000 Subject: [PATCH] Refuse corrupt statistics: blocks, counts, ranges, tails, plane identity Bucket-A integrity fixes for the statistics side of the brick. Valid inputs give bit-identical outputs; invalid ones are refused with a named error instead of being silently reinterpreted, wrapped or double-counted. - Blocks: every non-empty block must hold the plane's pixel type; a block in the other byte order is converted by value, so big-endian 8/16-bit planes are no longer read as native-endian bytes (#25, #26). - New integrity module: PlaneHistogram and HistogramTail check their invariants when made (int64 non-negative counts, exact unwrapped totals, valid <= total, all-or-none ordered finite limits, feasible tail sums with float64 overflow accepted only where the range allows it) (#45, #46, #47). - plane_histograms.csv: a header-only table, overflowing counts and repeated plane_id rows are refused; write_plane_histograms, scope_groups, fit and image_summaries refuse repeated planes too (#43, #45, #48). - PlaneRecord: z and t must be whole numbers 0 or above (#44). - PooledHistogram adds counts as integers (exact above 2**53) and refuses a group above 2**63 - 1 pixels (#31). - scope_groups refuses empty or non-text repeat-group labels and two biological repeats of one condition in one repeat group (#42). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G3qp9AxjpGPw9xfWAUJNqc --- AGENTS.md | 2 + docs/contract.md | 44 ++- src/morphobricks_normalization/csv_io.py | 125 +++--- src/morphobricks_normalization/fitting.py | 7 +- src/morphobricks_normalization/histogram.py | 49 +-- src/morphobricks_normalization/integrity.py | 236 +++++++++++ src/morphobricks_normalization/pooling.py | 41 +- src/morphobricks_normalization/records.py | 52 ++- src/morphobricks_normalization/scopes.py | 41 +- src/morphobricks_normalization/summaries.py | 6 +- tests/test_csv_io.py | 3 +- tests/test_integrity.py | 408 ++++++++++++++++++++ tests/test_scopes.py | 21 + 13 files changed, 924 insertions(+), 111 deletions(-) create mode 100644 src/morphobricks_normalization/integrity.py create mode 100644 tests/test_integrity.py diff --git a/AGENTS.md b/AGENTS.md index 8a0df27..87f5432 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -31,6 +31,8 @@ src/morphobricks_normalization/ robust range with HistogramTail underflow/overflow and histogram_with_range_check, quantiles by radix selection (Milestone 5) ranges.py HistogramSettings, RangeCheck, range_check_issues (Milestone 5) + integrity.py invariants every PlaneHistogram and HistogramTail checks when made + (counts, limits, tail sums) and the pixel type of every block records.py PlaneRecord (one plane_histograms.csv row: an XYC plane, one histogram per channel), plane_id, group_records csv_io.py write/read plane_histograms.csv schema version 4, channels in diff --git a/docs/contract.md b/docs/contract.md index 8e6e761..4bf4056 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -145,6 +145,35 @@ plane, so a row never depends on whether, or in which blocks, its plane was read robust range (§1.4) this holds for the limits, the counts and the tail counts too; only the tails' float sums may differ in the last bits, because they are added block by block. +Every non-empty block must hold the plane's pixel type: a block of another type (another +kind, width or signedness) is refused with `ValueError` naming the block and both types, +because reading its bytes as the plane's type would change its values and its number of +pixels. A block stored in the other byte order (for example a big-endian TIFF read on a +little-endian computer) holds the same type and is converted by value. + +##### Invariants of a histogram + +`PlaneHistogram` and `HistogramTail` check these rules when they are made, whether +computed, built by a calling brick or read from the statistics file (module +`integrity.py`); a histogram that breaks one is refused with `ValueError`, never clipped or +filled in: + +- `counts` is a one-dimensional int64 array of `bin_count` non-negative values; + `total_pixel_count` and `valid_pixel_count` are whole numbers from 0 to 2**63 − 1, with + `valid ≤ total`; the bins and the two tails hold exactly the valid pixels (an exact + integer sum, never one that wraps). +- Without valid pixels, the four limits are blank and the tails empty. With valid pixels, + the four limits are finite and ordered `pixel_min ≤ binning_min ≤ binning_max ≤ + pixel_max`, and a binning limit differs from its extremum exactly when a tail lies + beyond it. +- An empty tail has zero sums. A populated tail has a sum that is a number and a sum of + squares of 0 or more, with `count · sum_of_squares ≥ sum²`, a mean within its range and + a sum of squares possible for values in it (relative rounding slack 1e-9). A sum of + squares (or sum) that overflows float64 to infinity is accepted only when the tail's + range makes the overflow possible (values beyond about 1e154). +- Pooling adds the counts as integers, so a group's pixel count is exact; a group of more + than 2**63 − 1 finite pixels is refused. + #### 1.2.3 Group statistics Groups pool their planes' (representative value, count) pairs; per-plane statistics are @@ -407,7 +436,9 @@ y_c = (x_c − offset) / scale for c = R, G, B - **Pairing:** `repeat_groups` maps `(condition, biological_repeat)` to a label chosen in the preprocessing brick (`experiment_design.csv`, follow-up P3). Every `(condition, biological_repeat)` of the statistics must be in the mapping when one is given; a missing - one is an error (§9.1). + one is an error (§9.1). Its labels are non-empty text, and a repeat group holds at most + one biological repeat of each condition: two repeats of one condition mapped to the same + label are an error (§9.1), never merged. Entries for pairs without rows are ignored. - **`single`:** a condition without biological-repeat folders has the label `single`, which forms a group with the other conditions' `single` rows, or with the repeat it is paired with. @@ -503,10 +534,10 @@ they are zero. |---|---|---| | `position_id` | text or blank | The XY position within its technical repeat (MorphoBricks decision D02): from reliable metadata, or the image file name without extension when one file holds exactly one position, as morphobricks-preprocessing writes it. Blank only in tables written before it was filled. | | `image_file` | text | File name. | -| `z`, `t` | integer | Zero-based plane coordinates, always written; 0 for an axis the image does not have (preprocessing contract §2). | +| `z`, `t` | integer | Zero-based plane coordinates, always written; 0 for an axis the image does not have (preprocessing contract §2). Whole numbers 0 or above; `PlaneRecord` refuses booleans, fractions, text and negative numbers. | | `condition`, `biological_repeat`, `technical_repeat` | text | Folder labels; `single` when absent. Always read as text. | | `source_relative_path` | text | Path relative to the experiment root, `/`-separated. | -| `plane_id` | text | `source_relative_path::z=…,t=…`, always with both coordinates (`control/img.tif::z=0,t=0` for a 2D image). Row key. | +| `plane_id` | text | `source_relative_path::z=…,t=…`, always with both coordinates (`control/img.tif::z=0,t=0` for a 2D image). Row key: unique in a table, and in the rows given to `scope_groups`, `fit` and `image_summaries`. | | `dtype` | text | NumPy dtype name, e.g. `uint16`; the same for every channel. | | `bit_depth` | integer | Storage bits (8, 16, 32, 64). | | `channel_count` | integer | Number of channels; the same in every row. | @@ -661,6 +692,13 @@ Raised as exceptions before any pixel is normalized; the calling brick blocks th | Percentiles outside `0 ≤ low < high ≤ 100`; unknown `pooling` or `image_type` | `ValueError` | Existing message for percentiles; the allowed values otherwise | | Statistics file with unexpected columns, schema version or inconsistent counts | `CsvFormatError` | Existing messages | | Statistics file row whose `plane_id`, `image_file` or `bit_depth` differs from what its `source_relative_path`, `z`, `t` and `dtype` give | `CsvFormatError` | Names the file, line, column, stored value and implied value, and asks to run the preprocessing brick again | +| Statistics file with the columns but no rows | `CsvFormatError` | Names the file and asks to run the preprocessing brick again | +| Two rows with the same `plane_id`: in the statistics file, or given to `write_plane_histograms`, `scope_groups`, `fit` or `image_summaries` | `CsvFormatError` (file), `ValueError` | Names the plane and both lines (file) or every repeated plane, and asks to remove the repeats or run the preprocessing brick again | +| A histogram or tail that breaks the invariants of §1.2.2 (negative or wrapping counts, valid above total, missing or disordered limits, impossible tail sums) | `ValueError`; `CsvFormatError` in the file | Names the rule and the values; in the file also the line and channel | +| A plane coordinate `z` or `t` that is not a whole number 0 or above | `ValueError`; `CsvFormatError` in the file | Names the image, coordinate and value | +| A block of another pixel type than the plane's | `ValueError` | Names the block and both pixel types | +| `repeat_groups` with an empty or non-text label, or two biological repeats of one condition in one repeat group | `ValueError` | Names the pairs and the label, and asks to correct the experiment design in the preprocessing brick | +| A group of more than 2**63 − 1 finite pixels | `ValueError` | Asks to pool fewer planes | | Statistics file of schema version 1 or 2 (one row per channel) | `CsvFormatError` | "This plane_histograms.csv was written with schema version {version}, one row per channel. Run the preprocessing brick again on this experiment: it rewrites the table with the channels in columns (schema version 4)." | | Statistics file of schema version 3 (N56) | `CsvFormatError` | "This plane_histograms.csv was written with schema version 3, without the underflow and overflow columns. Run the preprocessing brick again on this experiment: it rewrites the table with these columns (schema version 4)." | | `HistogramSettings` with an unknown range, invalid quantiles or `min_central_bins` below 1 | `ValueError` | Names the value and the allowed ones | diff --git a/src/morphobricks_normalization/csv_io.py b/src/morphobricks_normalization/csv_io.py index 150e98a..504c64d 100644 --- a/src/morphobricks_normalization/csv_io.py +++ b/src/morphobricks_normalization/csv_io.py @@ -29,7 +29,7 @@ from .dtypes import bit_depth, check_supported, is_integer from .histogram import EMPTY_TAIL, HistogramTail, PlaneHistogram -from .records import PlaneRecord, different_channel_counts +from .records import PlaneRecord, different_channel_counts, duplicate_planes # Version 4: the underflow and overflow columns of the robust histogram range # (N56). Version 3 put one row per XYC plane, the channels in columns (N41); @@ -118,7 +118,8 @@ def write_plane_histograms(records: Iterable[PlaneRecord], path: str | Path) -> Raises ------ ValueError - If there are no records, or their channel counts or bin counts differ. + If there are no records, two records have the same ``plane_id``, or their + channel counts or bin counts differ. """ records = list(records) if not records: @@ -135,6 +136,7 @@ def write_plane_histograms(records: Iterable[PlaneRecord], path: str | Path) -> raise ValueError( f"All planes in one table must use the same bin count; found {sorted(bin_counts)}." ) + duplicate_planes(records, "The planes to write") path = Path(path) n_channels, n_bins = channel_counts.pop(), bin_counts.pop() header = [*METADATA_COLUMNS] @@ -162,9 +164,11 @@ def read_plane_histograms(path: str | Path) -> list[PlaneRecord]: CsvFormatError If the header, schema version or a row does not match the format, including a row whose ``plane_id``, ``image_file`` or ``bit_depth`` differs from what its - path, coordinates and pixel type give; if the table was written with an - earlier schema version (1 and 2 had one row per channel, 3 no underflow and - overflow columns); or if its rows have different numbers of channels. + path, coordinates and pixel type give, or whose statistics break the rules + of :mod:`~morphobricks_normalization.integrity`; if the table was written with + an earlier schema version (1 and 2 had one row per channel, 3 no underflow and + overflow columns); if it has no rows, or two rows with the same ``plane_id``; + or if its rows have different numbers of channels. """ with open(path, newline="", encoding="utf-8") as fh: reader = csv.reader(fh) @@ -187,16 +191,24 @@ def read_plane_histograms(path: str | Path) -> list[PlaneRecord]: if layout is None: raise CsvFormatError(f"{path} has unexpected per-channel column names.") n_channels, n_bins = layout - records, other_counts = [], [] + records, other_counts, lines = [], [], {} for line, row in enumerate(reader, 2): record = _record(row, n_channels, n_bins, path, line) - if isinstance(record, PlaneRecord): - records.append(record) - else: + if not isinstance(record, PlaneRecord): other_counts.append(record) + continue + if record.plane_id in lines: + raise CsvFormatError(_duplicate_message(path, record.plane_id, lines, line)) + lines[record.plane_id] = line + records.append(record) if other_counts: listed = [(r.source_relative_path, r.channel_count) for r in records] + other_counts raise CsvFormatError(f"{path}: " + different_channel_counts(listed, "Its images")) + if not records: + raise CsvFormatError( + f"{path} has the plane-histogram columns but no rows, so it holds no planes. The " + "table is incomplete; run the preprocessing brick again on this experiment." + ) return records @@ -256,7 +268,7 @@ def _row(record: PlaneRecord) -> list[str]: def _record( - row: list[str], n_channels: int, n_bins: int, path, line: int + row: list[str], n_channels: int, n_bins: int, path: str | Path, line: int ) -> PlaneRecord | tuple[str, int]: """Read one row; a row with another channel count is returned as ``(path, count)``.""" n_meta = len(METADATA_COLUMNS) @@ -276,35 +288,19 @@ def _record( if channel_count != n_channels: return meta["source_relative_path"], channel_count dtype = check_supported(meta["dtype"]) - number = int if is_integer(dtype) else float + if int(meta["bin_count"]) != n_bins: + raise ValueError(f"its bin_count is {meta['bin_count']}, but it has {n_bins} bins.") histograms = [] for channel in range(n_channels): values = row[n_meta + channel * width : n_meta + (channel + 1) * width] - pixel_min, pixel_max, valid, lo, hi = values[:5] - underflow, overflow = _tail(values[5:8]), _tail(values[8:11]) - histogram = PlaneHistogram( - dtype=dtype, - bin_count=int(meta["bin_count"]), - binning_min=_value(lo, number), - binning_max=_value(hi, number), - counts=np.array(values[len(CHANNEL_FEATURES) :], dtype=np.int64), - total_pixel_count=int(meta["total_pixel_count"]), - valid_pixel_count=int(valid), - pixel_min=_value(pixel_min, number), - pixel_max=_value(pixel_max, number), - underflow=underflow, - overflow=overflow, - ) - in_bins = int(histogram.counts.sum()) - if ( - histogram.bin_count != n_bins - or in_bins + underflow.count + overflow.count != histogram.valid_pixel_count - or not _consistent_range(histogram) - ): + try: + histograms.append(_histogram(values, dtype, n_bins, meta["total_pixel_count"])) + except ValueError as error: raise CsvFormatError( - f"{path}, line {line}: the bin counts of channel {channel} are inconsistent." - ) - histograms.append(histogram) + f"{path}, line {line}: the bin counts of channel {channel} are inconsistent: " + f"{error} The table was changed after it was written; run the preprocessing " + "brick again on this experiment to rewrite it." + ) from error record = PlaneRecord( source_relative_path=meta["source_relative_path"], condition=meta["condition"], @@ -335,37 +331,44 @@ def _record( raise CsvFormatError(f"{path}, line {line}: {error}") from error +def _histogram(values: list[str], dtype: np.dtype, n_bins: int, total: str) -> PlaneHistogram: + """Read one channel's columns; the histogram checks its own invariants.""" + number = int if is_integer(dtype) else float + pixel_min, pixel_max, valid, lo, hi = values[:5] + try: + counts = np.array(values[len(CHANNEL_FEATURES) :], dtype=np.int64) + except OverflowError: + raise ValueError("a bin count is larger than 2**63 - 1.") from None + return PlaneHistogram( + dtype=dtype, + bin_count=n_bins, + binning_min=_value(lo, number), + binning_max=_value(hi, number), + counts=counts, + total_pixel_count=int(total), + valid_pixel_count=int(valid), + pixel_min=_value(pixel_min, number), + pixel_max=_value(pixel_max, number), + underflow=_tail(values[5:8]), + overflow=_tail(values[8:11]), + ) + + +def _duplicate_message(path: str | Path, plane_id: str, lines: dict[str, int], line: int) -> str: + return ( + f"{path}, line {line}: plane {plane_id} is already in line {lines[plane_id]}. Every " + "plane (an image at one Z and T position) has one row; a repeated row would count its " + "pixels twice. Run the preprocessing brick again on this experiment to rewrite the table." + ) + + def _tail_text(tail: HistogramTail) -> tuple[str, str, str]: return str(tail.count), repr(float(tail.sum)), repr(float(tail.sum_of_squares)) def _tail(texts: list[str]) -> HistogramTail: - count, total, squares = int(texts[0]), float(texts[1]), float(texts[2]) - if count < 0: - raise ValueError(f"a pixel count cannot be negative; got {count}.") - if count == 0: - if total != 0 or squares != 0: - raise ValueError("an empty underflow or overflow must have zero sums.") - return EMPTY_TAIL - return HistogramTail(count, total, squares) - - -def _consistent_range(histogram: PlaneHistogram) -> bool: - """Check that a binning limit differs from the extremum only where a tail lies beyond it.""" - if histogram.binning_min is None or histogram.pixel_min is None: - return ( - histogram.binning_min is None - and histogram.binning_max is None - and histogram.pixel_min is None - and histogram.pixel_max is None - and not histogram.has_tails - ) - low = (histogram.pixel_min == histogram.binning_min) == histogram.underflow.is_empty - high = (histogram.pixel_max == histogram.binning_max) == histogram.overflow.is_empty - ordered = ( - histogram.pixel_min <= histogram.binning_min <= histogram.binning_max <= histogram.pixel_max - ) - return low and high and ordered + tail = HistogramTail(int(texts[0]), float(texts[1]), float(texts[2])) + return EMPTY_TAIL if tail.is_empty else tail def _earlier_version_message(path, version: str) -> str: diff --git a/src/morphobricks_normalization/fitting.py b/src/morphobricks_normalization/fitting.py index a96b486..4944b9c 100644 --- a/src/morphobricks_normalization/fitting.py +++ b/src/morphobricks_normalization/fitting.py @@ -38,7 +38,7 @@ z_score, ) from .pooling import NoValidPixelsError, PooledHistogram, WeightedHistogram -from .records import PlaneRecord, different_channel_counts +from .records import PlaneRecord, different_channel_counts, duplicate_planes from .rgb import RGB_CHANNEL_COUNT, apply_rgb, check_rgb, rgb_percentile from .settings import NormalizationSettings, check_settings @@ -351,8 +351,8 @@ def fit( If a channel of the group has no finite pixel. ValueError If a fluorescence channel's scale is zero or not finite, the rows have - different numbers of channels (N43), or RGB images do not have exactly - three channels of an integer type (contract §5). + different numbers of channels (N43) or repeat a ``plane_id``, or RGB images + do not have exactly three channels of an integer type (contract §5). """ from_pixels = isinstance(statistics, np.ndarray) check_settings(settings, statistics_file=not from_pixels) @@ -423,6 +423,7 @@ def _records(statistics: Iterable[PlaneRecord], group: str | None) -> list[Plane "fit expects the plane_histograms.csv rows of one group (PlaneRecord objects) " "or a pixel array." ) + duplicate_planes(records, f"The rows of group '{group}'" if group else "The rows of this group") if len({record.channel_count for record in records}) > 1: subject = f"The images of group '{group}'" if group else "The images of this group" raise ValueError( diff --git a/src/morphobricks_normalization/histogram.py b/src/morphobricks_normalization/histogram.py index 2f06cc7..2341cd4 100644 --- a/src/morphobricks_normalization/histogram.py +++ b/src/morphobricks_normalization/histogram.py @@ -51,6 +51,7 @@ import numpy as np from .dtypes import check_supported, default_bin_count, is_integer, validate_bin_count +from .integrity import check_histogram, check_tail, checked_blocks from .ranges import CENTRAL_PERCENTILES, HistogramSettings, RangeCheck # Above this many distinct values, integer planes are binned with a search over @@ -75,13 +76,19 @@ class HistogramTail: count Number of pixels. sum, sum_of_squares - Sum of their values and of their squared values, in float64. + Sum of their values and of their squared values, in float64; checked + by :func:`~morphobricks_normalization.integrity.check_tail`. """ count: int = 0 sum: float = 0.0 sum_of_squares: float = 0.0 + def __post_init__(self) -> None: + checked = check_tail(self.count, self.sum, self.sum_of_squares) + for name, value in zip(("count", "sum", "sum_of_squares"), checked, strict=True): + object.__setattr__(self, name, value) + @property def is_empty(self) -> bool: """True when no pixel lies in this tail.""" @@ -120,7 +127,8 @@ class PlaneHistogram: pixels). Default to ``binning_min`` and ``binning_max``. underflow, overflow The pixels below ``binning_min`` and above ``binning_max``; empty unless - the robust range was used. + the robust range was used. All are checked by + :func:`~morphobricks_normalization.integrity.check_histogram`. """ dtype: np.dtype @@ -140,6 +148,7 @@ def __post_init__(self) -> None: object.__setattr__(self, "pixel_min", self.binning_min) if self.pixel_max is None: object.__setattr__(self, "pixel_max", self.binning_max) + check_histogram(self) @property def is_empty(self) -> bool: @@ -444,8 +453,7 @@ def _small_integer_histogram( bits = dtype.itemsize * 8 table = np.zeros(2**bits, dtype=np.int64) total = 0 - for block in blocks: - values = np.asarray(block).ravel() + for values in checked_blocks(blocks, dtype): if values.size == 0: continue total += int(values.size) @@ -459,14 +467,15 @@ def _small_integer_histogram( present = np.flatnonzero(table) first, last = int(present[0]), int(present[-1]) smallest = int(np.iinfo(dtype).min) - exact = last + 1 - first <= n - if exact or not (settings.robust or check): - starts = _integer_bin_starts(last + 1 - first, n) - counts = _counts_from_per_value(table[first : last + 1], starts, n) - histogram = PlaneHistogram( - dtype, n, smallest + first, smallest + last, counts, total, total + + def min_max() -> PlaneHistogram: + counts = _counts_from_per_value( + table[first : last + 1], _integer_bin_starts(last + 1 - first, n), n ) - return histogram, None + return PlaneHistogram(dtype, n, smallest + first, smallest + last, counts, total, total) + + if last + 1 - first <= n or not (settings.robust or check): + return min_max(), None # The per-value counts give every order statistic exactly, without another pass. cumulative = np.cumsum(table) @@ -495,11 +504,7 @@ def order(q: float) -> _Order: overflow=_table_tail(table[high + 1 : last + 1], smallest + high + 1), ) else: - starts = _integer_bin_starts(last + 1 - first, n) - counts = _counts_from_per_value(table[first : last + 1], starts, n) - histogram = PlaneHistogram( - dtype, n, smallest + first, smallest + last, counts, total, total - ) + histogram = min_max() return histogram, _range_check(histogram, orders, settings) if check else None @@ -521,8 +526,7 @@ def _wide_histogram( selection = _RadixSelection(dtype.itemsize * 8) if settings.robust or check else None total = valid = 0 lo = hi = None - for block in blocks(): - values = np.asarray(block).ravel() + for values in checked_blocks(blocks(), dtype): total += int(values.size) finite = values if integer else values[np.isfinite(values)] if finite.size == 0: @@ -544,8 +548,7 @@ def _wide_histogram( # last round of the selection shares the binning pass. shared_round = 0 if settings.robust else 1 while selection.rounds_left > shared_round: - for block in blocks(): - values = np.asarray(block).ravel() + for values in checked_blocks(blocks(), dtype): finite = values if integer else values[np.isfinite(values)] if finite.size: selection.feed(_ordered_keys(finite, dtype)) @@ -619,8 +622,7 @@ def _integer_bins( per_value = np.zeros(width, dtype=np.int64) if width <= _PER_VALUE_COUNT_LIMIT else None bin_starts = np.array(starts[:-1], dtype=np.uint64) counts = np.zeros(n, dtype=np.int64) - for block in blocks(): - values = np.asarray(block).ravel() + for values in checked_blocks(blocks(), dtype): if values.size == 0: continue if selection is not None: @@ -660,8 +662,7 @@ def _float_bins( return counts # Halving keeps hi - lo finite even for ranges near the float64 limits. half_range = hi * 0.5 - lo * 0.5 - for block in blocks(): - values = np.asarray(block).ravel() + for values in checked_blocks(blocks(), dtype): finite = values[np.isfinite(values)] if finite.size == 0: continue diff --git a/src/morphobricks_normalization/integrity.py b/src/morphobricks_normalization/integrity.py new file mode 100644 index 0000000..6ef325e --- /dev/null +++ b/src/morphobricks_normalization/integrity.py @@ -0,0 +1,236 @@ +"""Invariants of the statistics: pixel blocks, bin counts, ranges and tails (contract §1.2, §6.2). + +:class:`~morphobricks_normalization.histogram.PlaneHistogram` and +:class:`~morphobricks_normalization.histogram.HistogramTail` check these rules +when they are made, so a histogram computed from pixels, built by a calling +brick or read from ``plane_histograms.csv`` obeys the same ones. A statistic +that breaks a rule is refused with a ``ValueError``; nothing is clipped, +rounded or filled in. + +The sums of a tail are float64 accumulations of the tail's pixel values, so +the checks that compare them allow a relative rounding slack of +``MOMENT_TOLERANCE``. A float64 plane whose values exceed about ``1e154`` can +overflow its sum of squares to infinity (and its sum beyond ``1.8e308``); such +an overflow is accepted only where the tail's range makes it possible. +""" + +from __future__ import annotations + +import math +import numbers +import sys +from collections.abc import Iterable, Iterator +from typing import TYPE_CHECKING + +import numpy as np + +if TYPE_CHECKING: # the histogram module imports this one + from .histogram import HistogramTail, PlaneHistogram + +INT64_MAX = 2**63 - 1 +# Relative rounding slack of the float64 tail sums: summing N values in float64 +# loses at most about N * 2.2e-16 relative, far below this for any plane. +MOMENT_TOLERANCE = 1e-9 +# Below this float64 estimate, the int64 sum of non-negative counts cannot have +# wrapped (its rounding is far smaller than the 2**62 left before 2**63). +_SAFE_FLOAT_TOTAL = 2.0**62 + + +def checked_blocks(blocks: Iterable[np.ndarray], dtype: np.dtype) -> Iterator[np.ndarray]: + """Yield every block flattened, as values of ``dtype`` in native byte order. + + Parameters + ---------- + blocks + The blocks of one plane. + dtype + The plane's pixel type, in native byte order. + + Yields + ------ + numpy.ndarray + The block's pixels, one dimension. A block stored in the other byte + order is converted by value; empty blocks are passed on unchecked. + + Raises + ------ + ValueError + If a non-empty block holds another pixel type: reading its bytes as + ``dtype`` would change its values and its number of pixels. + """ + for index, block in enumerate(blocks): + values = np.asarray(block).ravel() + if values.size and values.dtype.newbyteorder("=") != dtype: + raise ValueError( + f"Block {index} of the plane holds {values.dtype.name} pixels, but the plane's " + f"pixel type is {dtype.name}. Give every block in the plane's pixel type; " + "converting it would change the pixel values." + ) + if not values.dtype.isnative: + values = values.astype(dtype) + yield values + + +def pixel_count(value: object, name: str) -> int: + """Return ``value`` as a Python int if it is a whole number of pixels. + + Raises + ------ + ValueError + For booleans, fractions, text, negative numbers and numbers above + ``2**63 - 1`` (the int64 counts of the statistics). + """ + if isinstance(value, bool | np.bool_) or not isinstance(value, numbers.Integral): + raise ValueError(f"The {name} must be a whole number of pixels; got {value!r}.") + number = int(value) + if not 0 <= number <= INT64_MAX: + raise ValueError(f"The {name} must be between 0 and 2**63 - 1 pixels; got {number}.") + return number + + +def exact_total(counts: np.ndarray) -> int: + """Return the exact sum of non-negative int64 counts, as a Python int (never wrapped).""" + if float(np.sum(counts, dtype=np.float64)) < _SAFE_FLOAT_TOTAL: + return int(counts.sum()) + return sum(counts.tolist()) + + +def check_tail(count: object, total: object, squares: object) -> tuple[int, float, float]: + """Check one tail's count and sums on their own; return them as int and floats. + + Raises + ------ + ValueError + For a count that is not a whole number of pixels, sums of an empty tail + that are not zero, a NaN sum, a negative sum of squares, an infinite sum + with a finite sum of squares, or sums no set of values can have + (``count * sum_of_squares < sum**2``). + """ + count = pixel_count(count, "pixel count of an underflow or overflow") + total, squares = float(total), float(squares) + if count == 0: + if total != 0 or squares != 0: + raise ValueError("an empty underflow or overflow must have zero sums.") + return count, total, squares + if math.isnan(total) or not squares >= 0: + raise ValueError( + f"an underflow or overflow of {count} pixel(s) has sum {total!r} and sum of squares " + f"{squares!r}; the sum must be a number and the sum of squares 0 or more." + ) + if math.isinf(total) and math.isfinite(squares): + raise ValueError( + f"an underflow or overflow has an infinite sum {total!r} but a finite sum of " + f"squares {squares!r}." + ) + if math.isfinite(squares) and count * squares < total * total * (1 - MOMENT_TOLERANCE): + raise ValueError( + f"an underflow or overflow of {count} pixel(s) has sum {total!r} and sum of squares " + f"{squares!r}, which no {count} values can have (their variance would be negative)." + ) + return count, total, squares + + +def check_histogram(histogram: PlaneHistogram) -> None: + """Check a :class:`~morphobricks_normalization.histogram.PlaneHistogram`'s invariants. + + The counts are a one-dimensional int64 array of ``bin_count`` non-negative + values; the total and valid pixel counts are whole numbers with ``valid <= + total``; the bins and tails hold exactly the valid pixels. A histogram + without valid pixels has no limits, empty bins and empty tails; one with + valid pixels has four finite limits ``pixel_min <= binning_min <= + binning_max <= pixel_max``, a limit inside its extremum exactly where a tail + lies beyond it, and tails whose sums fit their range. + + Raises + ------ + ValueError + Naming the rule that is broken and the values that break it. + """ + counts = histogram.counts + if not isinstance(counts, np.ndarray) or counts.dtype != np.int64 or counts.ndim != 1: + raise ValueError("The bin counts must be a one-dimensional int64 array.") + if counts.size != histogram.bin_count: + raise ValueError( + f"There are {counts.size} bin counts for a bin count of {histogram.bin_count}." + ) + if counts.size and int(counts.min()) < 0: + raise ValueError(f"A bin count cannot be negative; got {int(counts.min())}.") + total = pixel_count(histogram.total_pixel_count, "total pixel count") + valid = pixel_count(histogram.valid_pixel_count, "valid pixel count") + if valid > total: + raise ValueError(f"The plane has {valid} valid pixels but only {total} pixels in all.") + under, over = histogram.underflow, histogram.overflow + in_bins = exact_total(counts) + if in_bins + under.count + over.count != valid: + raise ValueError( + f"The bins hold {in_bins} pixel(s) and the tails {under.count + over.count}, but " + f"the plane has {valid} valid pixel(s)." + ) + limits = ( + histogram.pixel_min, + histogram.binning_min, + histogram.binning_max, + histogram.pixel_max, + ) + if valid == 0: + if any(limit is not None for limit in limits) or not (under.is_empty and over.is_empty): + raise ValueError("A plane without valid pixels can have no intensity range or tails.") + return + _check_range(limits, under, over) + + +def _check_range( + limits: tuple[int | float | None, ...], under: HistogramTail, over: HistogramTail +) -> None: + """Check the four limits of a histogram with valid pixels, and its tails against them.""" + if any(limit is None for limit in limits) or not all( + math.isfinite(float(limit)) for limit in limits + ): + raise ValueError( + "A plane with valid pixels needs its four finite intensity limits (minimum, binning " + f"minimum, binning maximum, maximum); got {list(limits)}." + ) + pixel_min, binning_min, binning_max, pixel_max = limits + if not pixel_min <= binning_min <= binning_max <= pixel_max: + raise ValueError( + "The intensity limits must be ordered minimum <= binning minimum <= binning maximum " + f"<= maximum; got {list(limits)}." + ) + for name, tail, extremum, limit in ( + ("underflow", under, pixel_min, binning_min), + ("overflow", over, pixel_max, binning_max), + ): + if (extremum == limit) != tail.is_empty: + raise ValueError( + f"The binning limit {limit!r} and the extremum {extremum!r} differ exactly when " + f"pixels lie beyond the limit, but the {name} holds {tail.count} pixel(s)." + ) + if not tail.is_empty: + _check_tail_range(name, tail, *sorted((float(extremum), float(limit)))) + + +def _check_tail_range(name: str, tail: HistogramTail, low: float, high: float) -> None: + """Check that a tail's mean and sum of squares are possible for values in ``[low, high]``.""" + slack = MOMENT_TOLERANCE * max(abs(low), abs(high)) + largest = max(low * low, high * high) + smallest = 0.0 if low <= 0 <= high else min(low * low, high * high) + if math.isinf(tail.sum): + possible = ( + tail.count * max(abs(low), abs(high)) * (1 + MOMENT_TOLERANCE) >= sys.float_info.max + ) + else: + possible = low - slack <= tail.sum / tail.count <= high + slack + if math.isinf(tail.sum_of_squares): + possible = possible and tail.count * largest * (1 + MOMENT_TOLERANCE) >= sys.float_info.max + else: + possible = possible and ( + tail.count * smallest * (1 - MOMENT_TOLERANCE) + <= tail.sum_of_squares + <= tail.count * largest * (1 + MOMENT_TOLERANCE) + ) + if not possible: + raise ValueError( + f"The {name} of {tail.count} pixel(s) between {low!r} and {high!r} has sum " + f"{tail.sum!r} and sum of squares {tail.sum_of_squares!r}, which values in that " + "range cannot have." + ) diff --git a/src/morphobricks_normalization/pooling.py b/src/morphobricks_normalization/pooling.py index b279a27..abff086 100644 --- a/src/morphobricks_normalization/pooling.py +++ b/src/morphobricks_normalization/pooling.py @@ -39,6 +39,7 @@ import numpy as np from .histogram import PlaneHistogram +from .integrity import INT64_MAX from .records import PlaneRecord DEFAULT_PERCENTILES = (0.1, 1.0, 5.0, 25.0, 50.0, 75.0, 95.0, 99.0, 99.9) @@ -88,12 +89,18 @@ class PooledHistogram: If none of the planes has a finite pixel. ValueError If a record has more than one channel and ``channel`` is not given, or - has no channel ``channel``. + has no channel ``channel``; or if the planes hold more than ``2**63 - 1`` + finite pixels together, more than the int64 counts can hold. + + Notes + ----- + The counts are added as integers, so ``pixel_count`` is exactly the sum of + the planes' valid pixel counts, however large. """ def __init__( self, histograms: Iterable[PlaneHistogram | PlaneRecord], channel: int | None = None - ): + ) -> None: values, weights = [], [] bin_values, bin_weights = [], [] # (count, sum, sum of squares) of every non-empty tail. @@ -129,20 +136,19 @@ def __init__( "None of the selected planes contains finite pixel values, so no statistics " "can be calculated for this group." ) + if sum(int(w.sum()) for w in weights) > INT64_MAX: + raise ValueError( + "These planes hold more than 2**63 - 1 finite pixels together, more than the " + "pixel counts can hold. Pool fewer planes per group." + ) # Planes binned the same way share representative values; merging them # keeps the list short (a few thousand entries for typical integer data). - merged, inverse = np.unique(np.concatenate(values), return_inverse=True) - self.values = merged - self.counts = np.bincount(inverse, weights=np.concatenate(weights)).astype(np.int64) + self.values, self.counts = _merge(values, weights) self.pixel_count = int(self.counts.sum()) self._cumulative = np.cumsum(self.counts) if self.tails: # The bins alone, for the mean and standard deviation (the tails use their sums). - merged, inverse = np.unique(np.concatenate(bin_values), return_inverse=True) - self.bin_values = merged - self.bin_counts = np.bincount( - inverse, weights=np.concatenate(bin_weights), minlength=merged.size - ).astype(np.int64) + self.bin_values, self.bin_counts = _merge(bin_values, bin_weights) else: self.bin_values, self.bin_counts = self.values, self.counts @@ -211,6 +217,15 @@ def summary(self, percentiles: Sequence[float] = DEFAULT_PERCENTILES) -> PooledS ) +def _merge(values: list[np.ndarray], counts: list[np.ndarray]) -> tuple[np.ndarray, np.ndarray]: + """Merge (value, int64 count) pairs: sorted distinct values and their exact count sums.""" + merged, inverse = np.unique(np.concatenate(values), return_inverse=True) + total = np.zeros(merged.size, dtype=np.int64) + # Integer addition, never float64 weights: counts above 2**53 stay exact. + np.add.at(total, inverse.ravel(), np.concatenate(counts)) + return merged, total + + class WeightedHistogram: """Several pooled histograms added with a weight each (contract §4.3). @@ -227,6 +242,12 @@ class WeightedHistogram: channel. weights One positive, finite weight per part. + + Notes + ----- + The weighted counts are float64, as the weights are fractions in general, + so they carry float64 rounding (exact for whole counts below ``2**53``); + ``pixel_count`` adds the parts' exact integer counts. """ def __init__(self, parts: Sequence[PooledHistogram], weights: Sequence[float]): diff --git a/src/morphobricks_normalization/records.py b/src/morphobricks_normalization/records.py index 1688626..6d6d1b8 100644 --- a/src/morphobricks_normalization/records.py +++ b/src/morphobricks_normalization/records.py @@ -2,9 +2,12 @@ from __future__ import annotations -from collections.abc import Sequence +import numbers +from collections.abc import Iterable, Sequence from dataclasses import dataclass +import numpy as np + from .histogram import PlaneHistogram # Label used for hierarchy levels a branch of the folder tree does not have. @@ -40,7 +43,8 @@ class PlaneRecord: count and size. A single histogram is accepted for a one-channel plane. z, t Zero-based plane coordinates. An axis the image does not have, or has - with length 1, has coordinate 0; None is stored as 0. + with length 1, has coordinate 0; None is stored as 0. Booleans, + fractions, text and negative numbers are refused with ``ValueError``. position_id The XY position within its technical repeat (MorphoBricks decision D02): from reliable metadata, or the file stem when one file is one position; None @@ -60,8 +64,7 @@ def __post_init__(self) -> None: # Every plane has both coordinates (contract §6.2), so a plane is # identified the same way whether or not its file has the axis. for name in ("z", "t"): - if getattr(self, name) is None: - object.__setattr__(self, name, 0) + object.__setattr__(self, name, _coordinate(getattr(self, name), name, self)) histograms = self.histograms if isinstance(histograms, PlaneHistogram): histograms = (histograms,) @@ -133,6 +136,47 @@ def _check_channels(self) -> None: ) +def _coordinate(value: object, name: str, record: PlaneRecord) -> int: + """Return a plane coordinate as an int: None is 0, other values must be whole and >= 0.""" + if value is None: + return 0 + if isinstance(value, bool | np.bool_) or not isinstance(value, numbers.Integral) or value < 0: + raise ValueError( + f"Plane of {record.source_relative_path}: its {name} coordinate must be a whole number " + f"0 or above (coordinates count from 0); got {value!r}." + ) + return int(value) + + +def duplicate_planes(records: Iterable[PlaneRecord], subject: str) -> None: + """Refuse records that share a ``plane_id``: each plane's pixels must count once. + + Parameters + ---------- + records + The records of one table or group. + subject + Starts the message, e.g. ``"The planes to write"``. + + Raises + ------ + ValueError + Naming every repeated ``plane_id`` and how often it appears. + """ + seen: dict[str, int] = {} + for record in records: + seen[record.plane_id] = seen.get(record.plane_id, 0) + 1 + repeated = [f"{plane_id} ({n} times)" for plane_id, n in seen.items() if n > 1] + if repeated: + shown = ", ".join(repeated[:_MAX_LISTED]) + more = f" and {len(repeated) - _MAX_LISTED} more" if len(repeated) > _MAX_LISTED else "" + raise ValueError( + f"{subject} hold some planes more than once: {shown}{more}. Every plane (an image " + "at one Z and T position) must appear once, or its pixels count twice; remove the " + "repeated rows, or run the preprocessing brick again on this experiment." + ) + + def make_plane_id(source_relative_path: str, *, z: int | None = 0, t: int | None = 0) -> str: """Build a plane identifier such as ``"cond/rep1/img.tif::z=3,t=0"``. diff --git a/src/morphobricks_normalization/scopes.py b/src/morphobricks_normalization/scopes.py index 1055058..c8a201f 100644 --- a/src/morphobricks_normalization/scopes.py +++ b/src/morphobricks_normalization/scopes.py @@ -24,7 +24,7 @@ from collections.abc import Iterable, Mapping from .issues import REPEAT_MISSING_CONDITIONS, WARNING, Issue -from .records import PlaneRecord +from .records import PlaneRecord, duplicate_planes from .settings import NormalizationSettings RepeatGroups = Mapping[tuple[str, str], str] @@ -110,13 +110,19 @@ def scope_groups( ------ ValueError As :func:`group_key`; all missing ``(condition, biological_repeat)`` - pairs are named in one message. + pairs are named in one message. Also when two rows share a + ``plane_id``, or ``repeat_groups`` gives a present pair a label that is + not text or is empty, or puts two biological repeats of one condition in + the same repeat group (contract §4.1: one repeat per condition and group). """ records = list(records) + duplicate_planes(records, "The rows of plane_histograms.csv") if settings.scope == "biological_repeat" and repeat_groups is not None: - missing = {(r.condition, r.biological_repeat) for r in records} - set(repeat_groups) + pairs = {(r.condition, r.biological_repeat) for r in records} + missing = pairs - set(repeat_groups) if missing: raise _missing_pairs_error(missing) + _check_pairing(pairs, repeat_groups) groups: dict[tuple, list[PlaneRecord]] = {} for record in records: groups.setdefault(group_key(record, settings, repeat_groups), []).append(record) @@ -135,6 +141,35 @@ def _repeat_group(record: PlaneRecord, repeat_groups: RepeatGroups | None) -> st return repeat_groups[pair] +def _check_pairing(pairs: set[tuple[str, str]], repeat_groups: RepeatGroups) -> None: + """Refuse labels that are not text, and two repeats of one condition in one repeat group.""" + unlabelled = sorted( + p for p in pairs if not isinstance(repeat_groups[p], str) or not repeat_groups[p] + ) + if unlabelled: + names = ", ".join(f"{c}/{r} ({repeat_groups[(c, r)]!r})" for c, r in unlabelled) + raise ValueError( + f"The experiment design gives these biological repeats no repeat group name: {names}. " + "Name the repeat group of every biological repeat in the preprocessing brick's " + "experiment design, then run again." + ) + members: dict[tuple[str, str], set[str]] = {} + for condition, repeat in pairs: + members.setdefault((repeat_groups[(condition, repeat)], condition), set()).add(repeat) + conflicts = [ + f"{label}: {condition}/{' and '.join(sorted(repeats))}" + for (label, condition), repeats in sorted(members.items()) + if len(repeats) > 1 + ] + if conflicts: + raise ValueError( + "A repeat group holds one biological repeat of each condition, but the experiment " + "design puts several repeats of one condition in the same group " + f"({'; '.join(conflicts)}). Pair each biological repeat with one repeat of every " + "other condition in the preprocessing brick's experiment design, then run again." + ) + + def _missing_pairs_error(pairs: set[tuple[str, str]]) -> ValueError: names = ", ".join(f"{condition}/{repeat}" for condition, repeat in sorted(pairs)) return ValueError( diff --git a/src/morphobricks_normalization/summaries.py b/src/morphobricks_normalization/summaries.py index 925b2d8..be95831 100644 --- a/src/morphobricks_normalization/summaries.py +++ b/src/morphobricks_normalization/summaries.py @@ -28,7 +28,7 @@ from .csv_io import channel_column from .dtypes import is_integer from .pooling import NoValidPixelsError, PooledHistogram -from .records import PlaneRecord, different_channel_counts +from .records import PlaneRecord, different_channel_counts, duplicate_planes from .settings import DEFAULT_PERCENTILE SHARED_COLUMNS = ( @@ -157,12 +157,14 @@ def image_summaries( Raises ------ ValueError - For invalid percentiles, or images with different numbers of channels (N43). + For invalid percentiles, rows that repeat a ``plane_id``, or images with + different numbers of channels (N43). """ 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}.") records = list(records) + duplicate_planes(records, "These rows") if len({r.channel_count for r in records}) > 1: raise ValueError( different_channel_counts( diff --git a/tests/test_csv_io.py b/tests/test_csv_io.py index 8d59db8..64319cd 100644 --- a/tests/test_csv_io.py +++ b/tests/test_csv_io.py @@ -207,7 +207,8 @@ def failing_row(rec): original_row = csv_io._row monkeypatch.setattr(csv_io, "_row", failing_row) with pytest.raises(RuntimeError): - write_plane_histograms([record(), record()], path) + # Two different planes: one plane twice is refused before writing (issue #43). + write_plane_histograms([record("a/1.tif"), record("a/2.tif")], path) assert path.read_bytes() == before assert [p.name for p in tmp_path.iterdir()] == ["t.csv"] # no temporary files left diff --git a/tests/test_integrity.py b/tests/test_integrity.py new file mode 100644 index 0000000..c480a2b --- /dev/null +++ b/tests/test_integrity.py @@ -0,0 +1,408 @@ +"""Regression tests for the integrity of the statistics (issues #25, #26, #31, #43-#48).""" + +import csv + +import numpy as np +import pytest + +from morphobricks_normalization import ( + HistogramSettings, + HistogramTail, + NormalizationSettings, + PlaneHistogram, + PlaneRecord, + PooledHistogram, + compute_block_histogram, + compute_plane_histogram, + fit, + histogram_with_range_check, + image_summaries, + read_plane_histograms, + scope_groups, + write_plane_histograms, +) +from morphobricks_normalization.csv_io import CsvFormatError + +ROBUST = HistogramSettings(histogram_range="robust") +rng = np.random.default_rng(44) + + +def assert_same(a, b): + assert a.dtype == b.dtype and a.bin_count == b.bin_count + assert (a.pixel_min, a.binning_min, a.binning_max, a.pixel_max) == ( + b.pixel_min, + b.binning_min, + b.binning_max, + b.pixel_max, + ) + assert (a.total_pixel_count, a.valid_pixel_count) == (b.total_pixel_count, b.valid_pixel_count) + np.testing.assert_array_equal(a.counts, b.counts) + assert (a.underflow, a.overflow) == (b.underflow, b.overflow) + + +def plane_of(dtype): + if np.dtype(dtype).kind == "f": + plane = rng.normal(-3, 50, (40, 30)).astype(dtype) + plane.ravel()[:3] = [np.nan, -1e6, 1e6] + return plane + info = np.iinfo(dtype) + plane = rng.integers(max(info.min, -3000), min(info.max, 3000), (40, 30)).astype(dtype) + plane.ravel()[:2] = [info.min, info.max] + return plane + + +def blocks_of(plane): + return lambda: (plane[:7], plane[7:20, ::2], plane[7:20, 1::2], plane[20:]) + + +# --- Pixel blocks (#25, #26) ------------------------------------------------------------------- + + +@pytest.mark.parametrize("dtype", ["u1", "i1", "u2", "i2", "u4", "i4", "u8", "i8", "f4", "f8"]) +@pytest.mark.parametrize("settings", [None, ROBUST]) +def test_other_byte_order_gives_the_histogram_of_the_values(dtype, settings): + native = plane_of(np.dtype(dtype)) + for order in "<>": + swapped = native.astype(np.dtype(dtype).newbyteorder(order)) + np.testing.assert_array_equal(swapped, native) # same values, other bytes + expected = compute_plane_histogram(native, 50, settings=settings) + assert_same(compute_plane_histogram(swapped, 50, settings=settings), expected) + blockwise = compute_block_histogram( + blocks_of(swapped), swapped.dtype, 50, settings=settings + ) + assert_same(blockwise, expected) + checked, check = histogram_with_range_check(blocks_of(swapped), swapped.dtype, 50) + assert_same(checked, compute_plane_histogram(native, 50)) + assert check == histogram_with_range_check(blocks_of(native), native.dtype, 50)[1] + + +def test_big_endian_uint16_values_are_not_reinterpreted(): + histogram = compute_plane_histogram(np.array([1, 2, 3], dtype=">u2")) + assert (histogram.pixel_min, histogram.pixel_max) == (1, 3) + assert histogram.counts[:3].tolist() == [1, 1, 1] + + +@pytest.mark.parametrize( + "blocks, declared, observed", + [ + ((np.array([1, 2], np.uint8),), np.uint16, "uint8"), + ((np.array([1.5], np.float64),), np.float32, "float64"), + ((np.array([-1, 2], np.int16),), np.uint16, "int16"), + ((np.array([1], np.uint32),), np.int32, "uint32"), + ((np.array([1, 2], np.uint16), np.array([3], np.int64)), np.uint16, "int64"), + ], +) +def test_blocks_of_another_pixel_type_are_refused(blocks, declared, observed): + index = len(blocks) - 1 + message = f"Block {index} of the plane holds {observed} pixels, but the plane's pixel type" + with pytest.raises(ValueError, match=message): + compute_block_histogram(lambda: blocks, declared) + with pytest.raises(ValueError, match=message): + histogram_with_range_check(lambda: blocks, declared, settings=ROBUST) + + +def test_empty_blocks_of_any_type_are_skipped(): + blocks = (np.array([], np.float64), np.array([1, 2], np.uint16)) + histogram = compute_block_histogram(lambda: blocks, np.uint16) + assert histogram.valid_pixel_count == int(histogram.counts.sum()) == 2 + + +# --- Pooled counts (#31) ----------------------------------------------------------------------- + + +def singleton(count, value=0): + counts = np.array([count], dtype=np.int64) + return PlaneHistogram(np.dtype(np.uint8), 1, value, value, counts, count, count) + + +def test_pooled_counts_stay_exact_above_2_to_the_53(): + big = 2**53 + 1 + pooled = PooledHistogram([singleton(big)]) + assert pooled.pixel_count == big and pooled.counts.tolist() == [big] + overlapping = PooledHistogram([singleton(big), singleton(big), singleton(3, 7)]) + assert overlapping.counts.tolist() == [2 * big, 3] + assert overlapping.pixel_count == 2 * big + 3 + assert overlapping.percentile(100) == 7.0 and overlapping.percentile(0) == 0.0 + + +def test_pooled_counts_beyond_int64_are_refused(): + with pytest.raises(ValueError, match="more than 2\\*\\*63 - 1 finite pixels"): + PooledHistogram([singleton(2**62), singleton(2**62)]) + + +# --- Histogram invariants (#45, #47) ----------------------------------------------------------- + + +def histogram(**changes): + fields = dict( + dtype=np.dtype(np.uint8), + bin_count=2, + binning_min=0, + binning_max=1, + counts=np.array([1, 1], np.int64), + total_pixel_count=2, + valid_pixel_count=2, + ) + return PlaneHistogram(**{**fields, **changes}) + + +@pytest.mark.parametrize( + "changes, message", + [ + ({"counts": np.array([3, -1], np.int64)}, "cannot be negative"), + ({"counts": np.array([1.0, 1.0])}, "int64 array"), + ({"counts": np.array([[1, 1]], np.int64)}, "int64 array"), + ({"counts": np.array([2], np.int64)}, "1 bin counts for a bin count of 2"), + ({"valid_pixel_count": 3}, "3 valid pixels but only 2"), + ({"valid_pixel_count": -1}, "between 0 and 2"), + ({"total_pixel_count": 2.0}, "whole number"), + ({"valid_pixel_count": True}, "whole number"), + ({"total_pixel_count": 3, "valid_pixel_count": 3}, "hold 2 pixel"), + ({"binning_min": None, "binning_max": None}, "four finite intensity limits"), + ({"binning_max": None, "pixel_max": None}, "four finite intensity limits"), + ({"binning_min": 1, "binning_max": 0}, "must be ordered"), + ({"binning_max": float("nan")}, "four finite intensity limits"), + ({"pixel_min": -1}, "differ exactly when"), + ], +) +def test_impossible_histograms_are_refused(changes, message): + with pytest.raises(ValueError, match=message): + histogram(**changes) + + +def test_wrapping_int64_sum_is_refused(): + # The int64 sum of these counts wraps to 3 (issue #45). + counts = np.array([2**62] * 4 + [3], np.int64) + assert int(counts.sum()) == 3 + with pytest.raises(ValueError, match=f"bins hold {2**64 + 3} pixel"): + histogram( + bin_count=5, binning_max=4, counts=counts, total_pixel_count=3, valid_pixel_count=3 + ) + + +def test_a_plane_without_valid_pixels_has_no_limits(): + empty = dict(counts=np.zeros(2, np.int64), valid_pixel_count=0) + assert histogram(**empty, binning_min=None, binning_max=None).is_empty + with pytest.raises(ValueError, match="without valid pixels can have no intensity range"): + histogram(**empty) + + +# --- Tails (#46) ------------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "tail, message", + [ + ((2, float("nan"), 1.0), "the sum must be a number"), + ((2, 1.0, -1.0), "sum of squares 0 or more"), + ((2, 1.0, float("nan")), "sum of squares 0 or more"), + ((2, float("inf"), 4.0), "infinite sum"), + ((2, 10.0, 1.0), "variance would be negative"), + ((-1, 0.0, 0.0), "between 0 and 2"), + ((0, 1.0, 0.0), "empty underflow or overflow must have zero sums"), + ], +) +def test_impossible_tail_sums_are_refused(tail, message): + with pytest.raises(ValueError, match=message): + HistogramTail(*tail) + + +def test_signed_and_constant_tails_are_accepted(): + assert HistogramTail(2, -6.0, 18.0).mean == -3.0 + third = 0.1 + 0.2 # rounding of the sums must not refuse a constant tail + assert HistogramTail(3, 3 * third, 3 * third * third).count == 3 + + +def robust_plane_histogram(dtype=np.float64, outliers=(-1e6, 1e6)): + # Ten distinct outliers at each end, so the robust limits leave some below and above. + plane = rng.normal(0, 1, 40_000).astype(dtype) + plane[:10] = np.linspace(outliers[0], outliers[0] / 10, 10) + plane[10:20] = np.linspace(outliers[1] / 10, outliers[1], 10) + return compute_plane_histogram(plane, 100, settings=ROBUST) + + +def test_tails_must_fit_their_range(): + good = robust_plane_histogram() + assert not good.underflow.is_empty and not good.overflow.is_empty + lo, above = good.underflow, good.binning_min + 1.0 + for bad in ( + HistogramTail(lo.count, lo.count * above, lo.count * above * above), # above the limit + HistogramTail(lo.count, lo.sum, lo.sum_of_squares * 10), # values beyond the minimum + HistogramTail(lo.count, lo.sum, float("inf")), # cannot overflow in this range + ): + with pytest.raises(ValueError, match="which values in that range cannot have"): + PlaneHistogram( + good.dtype, + good.bin_count, + good.binning_min, + good.binning_max, + good.counts, + good.total_pixel_count, + good.valid_pixel_count, + good.pixel_min, + good.pixel_max, + bad, + good.overflow, + ) + + +@pytest.mark.filterwarnings("ignore:overflow encountered:RuntimeWarning") # issue #33 +def test_a_sum_of_squares_that_overflows_float64_is_accepted_where_possible(): + # Values near 1e200 overflow their squares (issue #33): legitimate, not corruption. + result = robust_plane_histogram(outliers=(-1e200, 1e200)) + assert result.underflow.sum_of_squares == np.inf and result.overflow.count >= 1 + + +# --- The statistics file (#43, #44, #45, #46, #47, #48) ---------------------------------------- + + +def record(path="c/r/a.tif", histograms=None, **coords): + if histograms is None: + histograms = [compute_plane_histogram(np.array([[0, 1], [1, 255]], np.uint8), 2)] + return PlaneRecord(path, "c", "r", "t", histograms, **coords) + + +def edited(tmp_path, records, changes, row=0): + path = write_plane_histograms(records, tmp_path / "plane_histograms.csv") + with open(path, newline="", encoding="utf-8") as fh: + header, *rows = list(csv.reader(fh)) + for column, value in changes.items(): + rows[row][header.index(column)] = value + with open(path, "w", newline="", encoding="utf-8") as fh: + csv.writer(fh).writerows([header, *rows]) + return path + + +LIMITS = ( + "Intensity_MinIntensity_C0", + "Intensity_MaxIntensity_C0", + "Histogram_BinningMin_C0", + "Histogram_BinningMax_C0", +) + + +@pytest.mark.parametrize( + "changes, message", + [ + ({"Histogram_Count_C0_00000": "4", "Histogram_Count_C0_00001": "-1"}, "negative"), + ({"Histogram_Count_C0_00000": str(2**63)}, "larger than 2\\*\\*63 - 1"), + ({"Histogram_Count_C0_00000": "1.5"}, "invalid literal"), + ({"Histogram_ValidPixelCount_C0": "-3", "total_pixel_count": "-3"}, "between 0"), + ({"Histogram_ValidPixelCount_C0": "5"}, "5 valid pixels but only 4"), + (dict.fromkeys(LIMITS, ""), "four finite intensity limits"), + ({"Intensity_MaxIntensity_C0": "", "Histogram_BinningMax_C0": ""}, "four finite"), + ], +) +def test_impossible_counts_and_limits_in_the_file_are_refused(tmp_path, changes, message): + path = edited(tmp_path, [record()], changes) + with pytest.raises(CsvFormatError, match=f"line 2: the bin counts of channel 0 .*{message}"): + read_plane_histograms(path) + + +def test_wrapping_counts_in_the_file_are_refused(tmp_path): + histograms = [compute_plane_histogram(np.arange(5, dtype=np.uint8), 5)] + changes = {f"Histogram_Count_C0_{i:05d}": str(2**62) for i in range(4)} + changes |= {"Histogram_Count_C0_00004": "3", "Histogram_ValidPixelCount_C0": "3"} + path = edited(tmp_path, [record(histograms=histograms)], changes | {"total_pixel_count": "3"}) + with pytest.raises(CsvFormatError, match=f"bins hold {2**64 + 3} pixel"): + read_plane_histograms(path) + + +@pytest.mark.parametrize( + "column, value", + [ + ("Histogram_UnderflowSum_C0", "nan"), + ("Histogram_UnderflowSumOfSquares_C0", "-1"), + ("Histogram_UnderflowSumOfSquares_C0", "inf"), + ("Histogram_OverflowSum_C0", "inf"), + ("Histogram_OverflowSum_C0", "0.0"), + ], +) +def test_impossible_sums_of_a_populated_tail_are_refused(tmp_path, column, value): + path = edited(tmp_path, [record(histograms=[robust_plane_histogram()])], {column: value}) + with pytest.raises(CsvFormatError, match="line 2: the bin counts of channel 0"): + read_plane_histograms(path) + + +@pytest.mark.filterwarnings("ignore:overflow encountered:RuntimeWarning") # issue #33 +def test_populated_tails_and_overflowed_squares_round_trip(tmp_path): + records = [ + record("c/r/a.tif", [robust_plane_histogram()]), + record("c/r/b.tif", [robust_plane_histogram(outliers=(-1e200, 1e200))]), + ] + back = read_plane_histograms(write_plane_histograms(records, tmp_path / "t.csv")) + for a, b in zip(records, back, strict=True): + assert_same(a.histograms[0], b.histograms[0]) + + +def test_a_table_without_rows_is_refused(tmp_path): + path = write_plane_histograms([record()], tmp_path / "plane_histograms.csv") + lines = path.read_text(encoding="utf-8").splitlines(keepends=True) + path.write_text(lines[0], encoding="utf-8") + with pytest.raises(CsvFormatError, match="has the plane-histogram columns but no rows"): + read_plane_histograms(path) + + +@pytest.mark.parametrize("conflicting", [False, True]) +def test_a_repeated_plane_row_is_refused(tmp_path, conflicting): + a, b = record("c/r/a.tif"), record("c/r/b.tif") + path = write_plane_histograms([a, b], tmp_path / "plane_histograms.csv") + lines = path.read_text(encoding="utf-8").splitlines(keepends=True) + repeated = lines[2] + if conflicting: # same plane, other counts: still one plane twice + other = [compute_plane_histogram(np.array([[9, 9], [9, 200]], np.uint8), 2)] + table = write_plane_histograms([record("c/r/b.tif", other)], tmp_path / "other.csv") + repeated = table.read_text(encoding="utf-8").splitlines(keepends=True)[1] + assert repeated != lines[2] + path.write_text("".join([*lines, repeated]), encoding="utf-8") + with pytest.raises( + CsvFormatError, match="line 4: plane c/r/b.tif::z=0,t=0 is already in line 3" + ): + read_plane_histograms(path) + + +def test_planes_that_differ_in_z_t_or_folder_are_not_repeats(tmp_path): + records = [ + record("c/r/a.tif"), + record("c/r/a.tif", t=1), + record("c/r/a.tif", z=2, t=1), + record("c/s/a.tif"), + ] + back = read_plane_histograms(write_plane_histograms(records, tmp_path / "t.csv")) + assert [r.plane_id for r in back] == [r.plane_id for r in records] + + +def test_repeated_planes_are_refused_by_every_entry_point(tmp_path): + a, b = record("c/r/a.tif"), record("c/r/b.tif") + message = "more than once: c/r/b.tif::z=0,t=0 \\(2 times\\)" + with pytest.raises(ValueError, match=message): + write_plane_histograms([a, b, b], tmp_path / "t.csv") + assert not (tmp_path / "t.csv").exists() + with pytest.raises(ValueError, match=message): + scope_groups([a, b, b], NormalizationSettings()) + with pytest.raises(ValueError, match=message): + fit([a, b, b], NormalizationSettings(method="minmax")) + with pytest.raises(ValueError, match=message): + image_summaries([a, b, b]) + + +@pytest.mark.parametrize("value", [-1, True, 1.5, 1.0, "1", np.float64(2.0), np.bool_(True)]) +@pytest.mark.parametrize("name", ["z", "t"]) +def test_plane_coordinates_must_be_whole_numbers_from_zero(name, value): + with pytest.raises( + ValueError, match=f"its {name} coordinate must be a whole number 0 or above" + ): + record(**{name: value}) + + +def test_valid_plane_coordinates(): + plane = record(z=np.int64(3), t=None) + assert (plane.z, plane.t) == (3, 0) and type(plane.z) is int + assert record(z=2**40, t=7).plane_id == f"c/r/a.tif::z={2**40},t=7" + + +@pytest.mark.parametrize("column", ["z", "t"]) +def test_negative_coordinates_in_the_file_are_refused(tmp_path, column): + path = edited(tmp_path, [record()], {column: "-1", "plane_id": "c/r/a.tif::z=0,t=0"}) + with pytest.raises(CsvFormatError, match=f"line 2: .*its {column} coordinate must be"): + read_plane_histograms(path) diff --git a/tests/test_scopes.py b/tests/test_scopes.py index dc04f45..6b6702d 100644 --- a/tests/test_scopes.py +++ b/tests/test_scopes.py @@ -86,6 +86,27 @@ def test_every_pair_must_be_in_the_pairing(): group_key(records[1], REPEAT, {("control", "r1"): "A"}) +def test_two_repeats_of_one_condition_in_one_repeat_group_are_refused(): + # Issue #42: a repeat group holds one biological repeat per condition. + records = [record("control", "r1"), record("control", "r2"), record("treated", "r3")] + pairing = {("control", "r1"): "A", ("control", "r2"): "A", ("treated", "r3"): "A"} + with pytest.raises(ValueError, match=r"same group \(A: control/r1 and r2\)"): + scope_groups(records, REPEAT, pairing) + # Several rows of one repeat (other images, time points) are fine, so are unused entries. + records = [record("control", "r1"), record("control", "r1", image="b.tif", t=1)] + records.append(record("treated", "r3")) + unused = {("control", "r1"): "A", ("treated", "r3"): "A", ("drug", "r9"): "A"} + groups, issues = scope_groups(records, REPEAT, unused) + assert list(groups) == [("A",)] and len(groups[("A",)]) == 3 and issues == [] + + +@pytest.mark.parametrize("label", ["", None, 1]) +def test_repeat_group_labels_must_be_text(label): + records = [record("control", "r1"), record("treated", "r1")] + with pytest.raises(ValueError, match="no repeat group name: control/r1"): + scope_groups(records, REPEAT, {("control", "r1"): label, ("treated", "r1"): "A"}) + + def test_a_repeat_missing_from_some_conditions_is_used_with_a_warning(): design = {"control": ("r1", "r2"), "treated": ("r1",), "drug": ("r1",)} records = make_design_records(simulate_design(seed=3, design=design, n_channels=1))