Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions AGENTS.md
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,8 @@ src/morphobricks_normalization/
issues.py Issue: warnings returned to the calling brick (channel for the checks)
tests/ pytest with simulated planes (tests/simulation.py); tests/data holds
a schema-version-2 fixture written by the Milestone 2 code
tools/dump_outputs.py bit-exact dump of the outputs for valid inputs and its SHA-256, to
show a change leaves them unchanged (run on base and change, compare)
```

The six modules `dtypes` to `normalization` were moved from morphobricks-preprocessing
Expand Down
9 changes: 9 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,15 @@ histogram, check = histogram_with_range_check(lambda: (plane,), plane.dtype, set
warnings = range_check_issues([("img.tif::z=0,t=0", 0, check)])
```

### Unreleased (next tag)

- **`PlaneHistogram`** no longer fills a missing `pixel_min` or `pixel_max` from
`binning_min` or `binning_max` (#47). A caller that builds a histogram of a plane with
valid pixels passes both extrema explicitly, e.g. `pixel_min=lo, pixel_max=hi` for a
min-max histogram. A `plane_histograms.csv` row with one blank intensity limit and valid
pixels is refused with a `CsvFormatError` naming the column. Valid inputs give
bit-identical outputs (`python tools/dump_outputs.py` prints the digest to compare).

### Migrating to v0.3.3

No valid statistics change: on the same inputs, histograms, `plane_histograms.csv` bytes,
Expand Down
6 changes: 4 additions & 2 deletions docs/contract.md
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,9 @@ filled in:
- 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.
beyond it. Each of the four is given: a missing `pixel_min` or `pixel_max` is refused,
never taken from `binning_min` or `binning_max`, even where a min-max histogram would
have them equal (#47). In the statistics file the message names the blank column.
- 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
Expand Down Expand Up @@ -296,7 +298,7 @@ The names that change in Milestone 5:

| Module | Milestone 5 |
|---|---|
| `histogram.py` | `PlaneHistogram` gains `pixel_min` and `pixel_max` as stored fields (default `binning_min`, `binning_max`), `underflow` and `overflow` (`HistogramTail(count, sum, sum_of_squares)`, empty by default) and `has_tails`; `is_exact` and `max_representative_error` include the tails. `compute_plane_histogram`, `compute_block_histogram` and `estimate_histogram_memory` take `settings`; new `histogram_with_range_check`. |
| `histogram.py` | `PlaneHistogram` gains `pixel_min` and `pixel_max` as stored fields (defaulting to `binning_min`, `binning_max` up to `v0.3.3`; from #47 they must be given whenever the plane has valid pixels), `underflow` and `overflow` (`HistogramTail(count, sum, sum_of_squares)`, empty by default) and `has_tails`; `is_exact` and `max_representative_error` include the tails. `compute_plane_histogram`, `compute_block_histogram` and `estimate_histogram_memory` take `settings`; new `histogram_with_range_check`. |
| `ranges.py` (new) | `HistogramSettings`, `RangeCheck`, `range_check_issues`. |
| `csv_io.py` | `SCHEMA_VERSION = "4"`; `CHANNEL_FEATURES` gains the six tail features (§6.2). |
| `pooling.py` | `PooledHistogram.values`/`counts` include one pair per tail; `bin_values`, `bin_counts` and `tails` give the bins alone and the tails' sums; `WeightedHistogram` the same, weighted. |
Expand Down
24 changes: 21 additions & 3 deletions src/morphobricks_normalization/csv_io.py
Original file line number Diff line number Diff line change
Expand Up @@ -294,7 +294,8 @@ def _record(
for channel in range(n_channels):
values = row[n_meta + channel * width : n_meta + (channel + 1) * width]
try:
histograms.append(_histogram(values, dtype, n_bins, meta["total_pixel_count"]))
total = meta["total_pixel_count"]
histograms.append(_histogram(values, dtype, n_bins, total, channel))
except ValueError as error:
raise CsvFormatError(
f"{path}, line {line}: the bin counts of channel {channel} are inconsistent: "
Expand Down Expand Up @@ -331,10 +332,27 @@ 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."""
def _histogram(
values: list[str], dtype: np.dtype, n_bins: int, total: str, channel: int
) -> PlaneHistogram:
"""Read one channel's columns; the histogram checks its own invariants.

A blank intensity limit of a channel with valid pixels is refused here, so
that the message names its column (#47).
"""
number = int if is_integer(dtype) else float
pixel_min, pixel_max, valid, lo, hi = values[:5]
blank = [
channel_column(feature, channel)
for feature, value in zip(CHANNEL_FEATURES[:5], values[:5], strict=True)
if value == "" and feature != "Histogram_ValidPixelCount"
]
if blank and int(valid) > 0:
raise ValueError(
f"{', '.join(blank)} {'is' if len(blank) == 1 else 'are'} blank, but the channel has "
f"{int(valid)} valid pixel(s): a plane with valid pixels needs its four finite "
"intensity limits, and none is filled in from another."
)
try:
counts = np.array(values[len(CHANNEL_FEATURES) :], dtype=np.int64)
except OverflowError:
Expand Down
16 changes: 8 additions & 8 deletions src/morphobricks_normalization/histogram.py
Original file line number Diff line number Diff line change
Expand Up @@ -123,8 +123,10 @@ class PlaneHistogram:
valid_pixel_count
Finite pixels: those counted in ``counts`` and in the two tails.
pixel_min, pixel_max
Exact finite minimum and maximum of the plane (None if it has no finite
pixels). Default to ``binning_min`` and ``binning_max``.
Exact finite minimum and maximum of the plane; None only if it has no
finite pixels. They are never taken from ``binning_min`` and
``binning_max``: a plane with finite pixels must give both, even when
they equal the binning limits (min-max range).
underflow, overflow
The pixels below ``binning_min`` and above ``binning_max``; empty unless
the robust range was used. All are checked by
Expand All @@ -144,10 +146,7 @@ class PlaneHistogram:
overflow: HistogramTail = EMPTY_TAIL

def __post_init__(self) -> None:
if self.pixel_min is None:
object.__setattr__(self, "pixel_min", self.binning_min)
if self.pixel_max is None:
object.__setattr__(self, "pixel_max", self.binning_max)
# No limit is filled in from another: a missing extremum is refused (#47).
check_histogram(self)

@property
Expand Down Expand Up @@ -472,7 +471,8 @@ def min_max() -> PlaneHistogram:
counts = _counts_from_per_value(
table[first : last + 1], _integer_bin_starts(last + 1 - first, n), n
)
return PlaneHistogram(dtype, n, smallest + first, smallest + last, counts, total, total)
lo, hi = smallest + first, smallest + last
return PlaneHistogram(dtype, n, lo, hi, counts, total, total, pixel_min=lo, pixel_max=hi)

if last + 1 - first <= n or not (settings.robust or check):
return min_max(), None
Expand Down Expand Up @@ -589,7 +589,7 @@ def _bin_pass(
else:
counts = _float_bins(blocks, dtype, n, lo, hi, valid, extrema, selection, under, over)
if extrema is None:
return PlaneHistogram(dtype, n, lo, hi, counts, total, valid)
return PlaneHistogram(dtype, n, lo, hi, counts, total, valid, pixel_min=lo, pixel_max=hi)
return PlaneHistogram(
dtype,
n,
Expand Down
19 changes: 13 additions & 6 deletions src/morphobricks_normalization/integrity.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,8 @@
# 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
# The four intensity limits of a histogram, in the order they must be in.
LIMIT_NAMES = ("pixel_min", "binning_min", "binning_max", "pixel_max")


def checked_blocks(blocks: Iterable[np.ndarray], dtype: np.dtype) -> Iterator[np.ndarray]:
Expand Down Expand Up @@ -138,8 +140,9 @@ def check_histogram(histogram: PlaneHistogram) -> None:
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.
binning_max <= pixel_max``, each given (none is filled in from another), a
limit inside its extremum exactly where a tail lies beyond it, and tails
whose sums fit their range.

Raises
------
Expand Down Expand Up @@ -183,12 +186,16 @@ 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
):
bad = [
name
for name, limit in zip(LIMIT_NAMES, limits, strict=True)
if limit is None or not math.isfinite(float(limit))
]
if bad:
raise ValueError(
"A plane with valid pixels needs its four finite intensity limits (minimum, binning "
f"minimum, binning maximum, maximum); got {list(limits)}."
f"minimum, binning maximum, maximum); {', '.join(bad)} "
f"{'is' if len(bad) == 1 else 'are'} missing or not finite, got {list(limits)}."
)
pixel_min, binning_min, binning_max, pixel_max = limits
if not pixel_min <= binning_min <= binning_max <= pixel_max:
Expand Down
50 changes: 48 additions & 2 deletions tests/test_integrity.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
write_plane_histograms,
)
from morphobricks_normalization.csv_io import CsvFormatError
from morphobricks_normalization.integrity import LIMIT_NAMES

ROBUST = HistogramSettings(histogram_range="robust")
rng = np.random.default_rng(44)
Expand Down Expand Up @@ -112,7 +113,9 @@ def test_empty_blocks_of_any_type_are_skipped():

def singleton(count, value=0):
counts = np.array([count], dtype=np.int64)
return PlaneHistogram(np.dtype(np.uint8), 1, value, value, counts, count, count)
return PlaneHistogram(
np.dtype(np.uint8), 1, value, value, counts, count, count, pixel_min=value, pixel_max=value
)


def test_pooled_counts_stay_exact_above_2_to_the_53():
Expand Down Expand Up @@ -142,6 +145,8 @@ def histogram(**changes):
counts=np.array([1, 1], np.int64),
total_pixel_count=2,
valid_pixel_count=2,
pixel_min=0,
pixel_max=1,
)
return PlaneHistogram(**{**fields, **changes})

Expand Down Expand Up @@ -170,6 +175,28 @@ def test_impossible_histograms_are_refused(changes, message):
histogram(**changes)


@pytest.mark.parametrize("field", LIMIT_NAMES)
def test_a_single_missing_limit_is_refused_not_filled_in(field):
# Before #47 was fixed, a missing pixel_min or pixel_max was taken from the binning limit.
with pytest.raises(ValueError, match=f"four finite intensity limits .*; {field} is missing"):
histogram(**{field: None})


@pytest.mark.parametrize("field", ["pixel_min", "pixel_max"])
def test_a_missing_extremum_of_a_robust_histogram_is_refused(field):
good = robust_plane_histogram()
fields = {name: getattr(good, name) for name in good.__dataclass_fields__}
with pytest.raises(ValueError, match=f"{field} is missing"):
PlaneHistogram(**{**fields, field: None})


def test_computed_min_max_histograms_give_their_extrema():
for plane in (np.array([3, 9, 4], np.uint8), np.array([1.5, -2.0, np.nan], np.float32)):
hist = compute_plane_histogram(plane, 4)
assert (hist.pixel_min, hist.pixel_max) == (hist.binning_min, hist.binning_max)
assert (hist.pixel_min, hist.pixel_max) == (np.nanmin(plane), np.nanmax(plane))


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)
Expand All @@ -182,7 +209,8 @@ def test_wrapping_int64_sum_is_refused():

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
no_limits = dict.fromkeys(LIMIT_NAMES)
assert histogram(**empty, **no_limits).is_empty
with pytest.raises(ValueError, match="without valid pixels can have no intensity range"):
histogram(**empty)

Expand Down Expand Up @@ -299,6 +327,24 @@ def test_impossible_counts_and_limits_in_the_file_are_refused(tmp_path, changes,
read_plane_histograms(path)


@pytest.mark.parametrize("column", LIMITS)
def test_a_single_blank_limit_in_the_file_is_refused(tmp_path, column):
# The default record's limits are min-max, so a blank extremum equals its binning limit
# and was silently filled in from it before #47 was fixed.
path = edited(tmp_path, [record()], {column: ""})
message = f"line 2: the bin counts of channel 0 .*{column} is blank, but the channel has 4"
with pytest.raises(CsvFormatError, match=message):
read_plane_histograms(path)


def test_blank_limits_of_a_channel_without_valid_pixels_are_read(tmp_path):
empty = compute_plane_histogram(np.full((2, 2), np.nan, np.float32), 2)
path = write_plane_histograms([record(histograms=[empty])], tmp_path / "plane_histograms.csv")
(back,) = read_plane_histograms(path)
assert back.histograms[0].is_empty
assert [getattr(back.histograms[0], name) for name in LIMIT_NAMES] == [None] * 4


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)}
Expand Down
2 changes: 2 additions & 0 deletions tests/test_schema_v2_equivalence.py
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,8 @@ def version_2_records(name):
counts=np.array([row[f"histogram_bin_{i:05d}"] for i in range(bin_count)], np.int64),
total_pixel_count=int(row["total_pixel_count"]),
valid_pixel_count=int(row["valid_pixel_count"]),
pixel_min=number(row["pixel_min"]) if row["pixel_min"] else None,
pixel_max=number(row["pixel_max"]) if row["pixel_max"] else None,
)
labels = (row["condition"], row["biological_repeat"], row["technical_repeat"])
key = (row["source_relative_path"], int(row["t"]), int(row["c"]))
Expand Down
Loading
Loading