From 5154ef3c2e8e4e4ebb652b025d1424126d6244ad Mon Sep 17 00:00:00 2001 From: Remis Date: Wed, 2 Sep 2026 07:43:06 +0100 Subject: [PATCH] Add golden-value regression tests for the record layer Every package in the stack reads its records through this code, so a change here reaches SCALCS, EKDIST and HJCFIT at once. Pinning it here rather than in HJCFIT is the same work with three times the coverage. Five cases over the two shipped records, each fixing resolved-interval counts, period counts, burst counts and openings, mean burst length, mean Popen within bursts, and the histogram binning. Resolutions are chosen so the dead time actually bites. The simulated record's shortest interval is 125 us, so a golden case at 25 or 100 us would merge nothing and exercise the concatenation logic not at all -- one of the invariants now asserts that, so a future case cannot quietly test nothing. Alongside the exact numbers, invariants that hold for any record at any resolution: periods alternate strictly, every burst has an odd interval count, burst openings account for every open period, no within-burst gap exceeds tcrit, the histogram counts every interval, the sign of tcrit is ignored. Those keep meaning something if the example files are ever replaced. Two of them were wrong when first written, and the failures were worth following rather than relaxing: Total time is not conserved outright. impose_resolution starts at the first interval that is both resolvable and usable, because anything earlier has a start that cannot be located, and the time before it is dropped -- 75.8 us of a 102.7 s record on the experimental example at 100 us. The invariant now accounts for that head explicitly rather than allowing a fudge factor that would have hidden it. Resolved intervals are not all at least tres long. Two documented exceptions carry FLAG_UNUSABLE: the negative sentinel for an unfinished final opening, and the short shut emitted when the record ends on a brief shutting that terminates an open run. The invariant is now about usable intervals, which is what a consumer relies on. 83 tests; 345 passing overall. Co-Authored-By: Claude Opus 5 --- tests/conftest.py | 6 ++ tests/test_golden.py | 250 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 256 insertions(+) create mode 100644 tests/test_golden.py diff --git a/tests/conftest.py b/tests/conftest.py index 22c1251..689297c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -45,3 +45,9 @@ def legacy_experimental_scn() -> Path: if not p.exists(): pytest.skip("glyr_experimental.scn not found.") return p + + +@pytest.fixture() +def examples_dir() -> Path: + """Directory holding the example SCN records.""" + return EXAMPLES_SCN diff --git a/tests/test_golden.py b/tests/test_golden.py new file mode 100644 index 0000000..61ce4c2 --- /dev/null +++ b/tests/test_golden.py @@ -0,0 +1,250 @@ +"""Golden values for the record layer. + +Every package in the DCPROGS stack reads its records through this code, so a +change here reaches SCALCS, EKDIST and HJCFIT at once. These are the numbers +that pin it: fixed outputs for the shipped example records at several +resolutions and critical times. + +Two kinds of assertion, deliberately mixed. + +The *golden* ones are exact and will move if any of the arithmetic changes -- +resolution imposition, period grouping, the burst convention, the binning. +That is what they are for. When one moves, the question to answer is whether +the change was intended, not whether the number should be updated. + +The *invariant* ones hold for any record at any resolution and say what the +layer guarantees rather than what it happened to produce. They are the ones +that keep meaning something if the example files are ever replaced. + +Resolutions are chosen so the dead time actually bites. The simulated record's +shortest interval is 125 us, so anything at or below 100 us merges nothing and +would exercise the concatenation logic not at all. +""" + +from dataclasses import dataclass + +import numpy as np +import pytest + +from dcio.analysis import bursts_from_record, from_scn, log_bin_histogram +from dcio.formats.scn import FLAG_UNUSABLE, read + + +@dataclass(frozen=True) +class Case: + filename: str + tres: float + tcrit: float + n_resolved: int + n_open: int + n_shut: int + n_bursts: int + n_openings: int + mean_burst_ms: float + mean_popen: float + open_nbdec: int + open_bins: int + shut_nbdec: int + shut_bins: int + + def __str__(self): + return f"{self.filename.split('_')[1][:4]}-{self.tres * 1e6:.0f}us-{self.tcrit * 1e3:.0f}ms" + + +CASES = [ + Case("glyr_experimental.scn", 2.5e-05, 0.004, + n_resolved=3232, n_open=731, n_shut=730, + n_bursts=63, n_openings=731, + mean_burst_ms=98.629089, mean_popen=0.991578, + open_nbdec=8, open_bins=29, shut_nbdec=8, shut_bins=53), + Case("glyr_experimental.scn", 0.0001, 0.004, + n_resolved=2408, n_open=90, n_shut=89, + n_bursts=62, n_openings=90, + mean_burst_ms=100.030402, mean_popen=0.997952, + open_nbdec=5, open_bins=20, shut_nbdec=5, shut_bins=30), + Case("glyr_experimental.scn", 0.0005, 0.01, + n_resolved=1812, n_open=67, n_shut=66, + n_bursts=56, n_openings=67, + mean_burst_ms=111.190894, mean_popen=0.995910, + open_nbdec=5, open_bins=17, shut_nbdec=5, shut_bins=27), + Case("glyr_simulated.scn", 0.0005, 0.004, + n_resolved=33377, n_open=16688, n_shut=16687, + n_bursts=9270, n_openings=16688, + mean_burst_ms=111.686132, mean_popen=0.993102, + open_nbdec=12, open_bins=40, shut_nbdec=12, shut_bins=52), + Case("glyr_simulated.scn", 0.001, 0.01, + n_resolved=23095, n_open=11547, n_shut=11546, + n_bursts=7021, n_openings=11547, + mean_burst_ms=149.705263, mean_popen=0.983202, + open_nbdec=12, open_bins=37, shut_nbdec=12, shut_bins=49), +] + + +@pytest.fixture(params=CASES, ids=str) +def case(request, examples_dir): + c = request.param + path = examples_dir / c.filename + if not path.is_file(): + pytest.skip(f"{c.filename} not found") + return c, from_scn(read(path), tres=c.tres) + + +# -------------------------------------------------------------------------- +# Golden: exact outputs +# -------------------------------------------------------------------------- + +class TestGoldenValues: + + def test_resolved_interval_count(self, case): + c, rec = case + assert rec.n_resolved == c.n_resolved + + def test_period_counts(self, case): + c, rec = case + assert (rec.periods.n_open, rec.periods.n_shut) == (c.n_open, c.n_shut) + + def test_burst_count_and_openings(self, case): + c, rec = case + lengths, nops = bursts_from_record(rec, tcrit=c.tcrit) + assert len(lengths) == c.n_bursts + assert int(nops.sum()) == c.n_openings + + def test_mean_burst_length(self, case): + c, rec = case + lengths, _ = bursts_from_record(rec, tcrit=c.tcrit) + assert lengths.mean() * 1e3 == pytest.approx(c.mean_burst_ms, abs=1e-6) + + def test_mean_popen_within_bursts(self, case): + c, rec = case + segs = bursts_from_record(rec, tcrit=c.tcrit, intervals_only=True) + popen = np.mean([s[0::2].sum() / s.sum() for s in segs]) + assert popen == pytest.approx(c.mean_popen, abs=1e-6) + + def test_histogram_binning(self, case): + c, rec = case + for kind, X, nbdec, nbins in ( + ("open", rec.periods.open_intervals, c.open_nbdec, c.open_bins), + ("shut", rec.periods.shut_intervals, c.shut_nbdec, c.shut_bins), + ): + counts, _, got = log_bin_histogram(X, c.tres) + assert got == nbdec, kind + assert len(counts) == nbins, kind + + +# -------------------------------------------------------------------------- +# Invariant: what the layer guarantees, for any record at any resolution +# -------------------------------------------------------------------------- + +class TestInvariants: + + def test_resolution_actually_bites(self, case): + """A golden case that merges nothing tests the concatenation logic not + at all. The simulated record's shortest interval is 125 us, so this + would silently pass for every tres below that.""" + c, rec = case + assert rec.n_resolved < len(rec.intervals) + + def test_resolved_intervals_are_resolvable_unless_flagged(self, case): + """Every resolved interval is at least one dead time long, unless it is + flagged unusable. + + Two kinds of exception exist and both are deliberate: a negative + sentinel for an unfinished final opening, and the short shut that + impose_resolution emits when the last interval of the record is a brief + shutting that terminates an open run. Both carry FLAG_UNUSABLE, so the + invariant is about the usable intervals -- which is what a consumer + actually relies on.""" + c, rec = case + usable = (rec.resolved_flags & FLAG_UNUSABLE) == 0 + finite = rec.resolved_intervals[usable & (rec.resolved_intervals > 0)] + assert (finite >= c.tres - 1e-15).all() + + def test_total_time_is_conserved_after_the_first_resolvable_interval(self, case): + """Concatenation moves time between intervals; it must not create or + destroy any. + + The head of the record is the one legitimate loss. impose_resolution + starts at the first interval that is both resolvable and usable, since + anything earlier has a start that cannot be located, and the time + before it is dropped. On the experimental example at 100 us that is + 75.8 us of a 102.7 s record -- small, but not rounding, and an + invariant that ignored it would be false rather than approximate.""" + c, rec = case + flags = rec.flags.astype(np.int32).copy() + flags[rec.intervals < 0] |= FLAG_UNUSABLE + resolvable = np.where((rec.intervals > c.tres) & (flags < FLAG_UNUSABLE))[0] + head = rec.intervals[:int(resolvable[0])] + discarded = head[head > 0].sum() + + raw = rec.intervals[rec.intervals > 0].sum() + resolved = rec.resolved_intervals[rec.resolved_intervals > 0].sum() + assert resolved + discarded == pytest.approx(raw, rel=1e-12) + + def test_periods_strictly_alternate(self, case): + """set_periods merges runs of the same conductance, so the period list + alternates open/shut with no two of a kind adjacent. Burst + segmentation depends on this.""" + c, rec = case + is_open = rec.periods.amplitudes != 0.0 + assert not (is_open[:-1] & is_open[1:]).any() + assert not (~is_open[:-1] & ~is_open[1:]).any() + + def test_periods_start_and_end_open(self, case): + c, rec = case + assert rec.periods.n_open == rec.periods.n_shut + 1 + + def test_every_burst_has_an_odd_interval_count(self, case): + """Each burst starts and ends on an opening. The missed-events + likelihood is a product of matrices alternating A->F and F->A; an + even-length burst would end on a shut.""" + c, rec = case + for seq in bursts_from_record(rec, tcrit=c.tcrit, intervals_only=True): + assert len(seq) % 2 == 1 + + def test_burst_openings_equal_period_openings(self, case): + """Every open period belongs to exactly one burst, because the record + begins and ends on an opening and nothing between is discarded.""" + c, rec = case + _, nops = bursts_from_record(rec, tcrit=c.tcrit) + assert int(nops.sum()) == rec.periods.n_open + + def test_no_burst_contains_a_gap_longer_than_tcrit(self, case): + c, rec = case + for seq in bursts_from_record(rec, tcrit=c.tcrit, intervals_only=True): + if len(seq) > 1: + assert (seq[1::2] <= c.tcrit).all() + + def test_histogram_counts_every_interval(self, case): + c, rec = case + for X in (rec.periods.open_intervals, rec.periods.shut_intervals): + counts, edges, _ = log_bin_histogram(X, c.tres) + assert counts.sum() == len(X) + assert X.max() <= edges[-1] + + def test_tcrit_sign_is_ignored(self, case): + c, rec = case + pos, _ = bursts_from_record(rec, tcrit=c.tcrit) + neg, _ = bursts_from_record(rec, tcrit=-c.tcrit) + np.testing.assert_allclose(pos, neg) + + +class TestGoldenTableItself: + """The table is only worth having if it covers what it claims to.""" + + def test_both_example_records_are_covered(self): + assert {c.filename for c in CASES} == { + "glyr_experimental.scn", "glyr_simulated.scn"} + + def test_experimental_record_carries_unusable_intervals(self, examples_dir): + """The experimental record is the one that exercises flag handling; + if it ever stops containing unusable intervals, the burst convention + for them is no longer covered by any golden case.""" + path = examples_dir / "glyr_experimental.scn" + if not path.is_file(): + pytest.skip("glyr_experimental.scn not found") + rec = read(path) + assert int((rec.flags & FLAG_UNUSABLE != 0).sum()) > 0 + + def test_cases_span_more_than_one_resolution_per_record(self): + for name in {c.filename for c in CASES}: + assert len({c.tres for c in CASES if c.filename == name}) > 1