Bucket A (1/2): refuse corrupt statistics: blocks, counts, ranges, tails, plane identity

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

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 |

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:

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(