From 3e5310571cc290acb1305015b133c7cb865a4b06 Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 8 Sep 2026 09:35:54 +0100 Subject: [PATCH] Check burst segmentation against the rule transcribed from the Fortran test_bursts.py fixes the burst conventions by stating them. This fixes them against a separate reading of the original, which is a different kind of evidence: the three implementations in the DCPROGS stack disagreed on exactly these points, and agreeing with each other would not have settled which was right. hjclik.for in DCPROGS/DCFORTRAN is the likelihood that produced every published HJCFIT result, Colquhoun, Hatton & Hawkes (2003) among them. It does not pre-split a record; line 1009 walks it and ends the group when a shut time exceeds tcrit: if(tint(in,jset).gt.tcrit(jset)) goto 92 !end present group with prev opening Three things follow, and each gets a test: the comparison is strictly greater, so a gap exactly equal to tcrit is within-burst; the group ends "with prev opening", so the separating gap belongs to no group and every group starts and ends on an opening; the test is NOT guarded by burst(jset), so a tcrit of zero ends the group at every shut time. HJCFIT means "do not divide" by setting tcrit enormous -- 3.1536e10 ms, one year (Hjcfit1-09122003.for line 1568) -- not by setting it to zero. That last one is worth a test on its own. A caller who passes zero expecting "no division" gets one group per opening, silently and with no error, and both implementations behave that way. Finding it cost real time. The two agree group for group on five records of 4000 intervals and at four values of tcrit. 13 tests. Co-Authored-By: Claude Opus 5 --- tests/test_bursts_fortran.py | 150 +++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) create mode 100644 tests/test_bursts_fortran.py diff --git a/tests/test_bursts_fortran.py b/tests/test_bursts_fortran.py new file mode 100644 index 0000000..27f8177 --- /dev/null +++ b/tests/test_bursts_fortran.py @@ -0,0 +1,150 @@ +"""Burst segmentation against the rule transcribed from the Fortran. + +``dcio.analysis.bursts`` is checked here against an independent implementation +of the same definition, transcribed from ``hjclik.for`` in +``DCPROGS/DCFORTRAN`` (``Fort90/HJCFIT/``) -- the likelihood that produced +every published HJCFIT result, Colquhoun, Hatton & Hawkes (2003) among them. + +``test_bursts.py`` fixes the conventions by stating them; this fixes them +against a separate reading of the original, which is a different kind of +evidence. The three implementations in the DCPROGS stack disagreed on exactly +these points, and agreeing with each other would not have settled which was +right. + +The Fortran does not pre-split a record. ``hjclik.for`` line 1009 walks it and +ends the group when a shut time exceeds ``tcrit``:: + + if(tint(in,jset).gt.tcrit(jset)) goto 92 !end present group with prev opening + +Three things follow, and each is a separate test below: + +* the comparison is **strictly greater**, so a gap exactly equal to ``tcrit`` + is within-burst; +* the group ends *with the previous opening*, so the separating gap belongs to + no group and every group starts and ends on an opening; +* the test is **not guarded by** ``burst(jset)``. A ``tcrit`` of zero therefore + ends the group at every shut time. HJCFIT means "do not divide" by setting + ``tcrit`` enormous -- 3.1536e10 ms, one year (``Hjcfit1-09122003.for`` line + 1568) -- not by setting it to zero. That trap cost real time to find and is + worth a test of its own. +""" + +import numpy as np +import numpy.testing as npt +import pytest + +from dcio.analysis.bursts import extract_burst_intervals + +TCRIT = 3.5e-3 + + +def fortran_bursts(tints, ampls, tcrit): + """``hjclik.for`` line 1009, transcribed. + + A shut time greater than ``tcrit`` ends the group with the previous + opening. Returns a list of interval arrays, each starting and ending on an + opening. + """ + t = np.asarray(tints, float) + a = np.asarray(ampls, float) + groups, cur = [], [] + for ti, ai in zip(t, a): + if ai == 0.0 and ti > tcrit: # strictly greater + if cur: + groups.append(cur) + cur = [] + else: + cur.append((ti, ai)) + if cur: + groups.append(cur) + + out = [] + for g in groups: # "with prev opening" + while g and g[0][1] == 0.0: + g = g[1:] + while g and g[-1][1] == 0.0: + g = g[:-1] + if g: + out.append(np.array([x[0] for x in g])) + return out + + +def alternating(n, seed, tcrit=TCRIT): + """A record that alternates open and shut, with gaps either side of tcrit.""" + rng = np.random.default_rng(seed) + t = np.empty(n) + a = np.zeros(n) + for i in range(n): + if i % 2 == 0: # opening + t[i] = rng.exponential(3e-4) + a[i] = 5.0 + else: # shut: a mixture, so that + t[i] = (rng.exponential(2e-4) if rng.random() < 0.75 + else rng.exponential(2e-2)) # some gaps exceed tcrit + return t, a + + +@pytest.mark.parametrize("seed", [0, 1, 2, 3, 4]) +def test_agrees_with_the_fortran_rule(seed): + t, a = alternating(4000, seed) + got = extract_burst_intervals(t, a, TCRIT) + want = fortran_bursts(t, a, TCRIT) + + assert len(got) == len(want) + for g, w in zip(got, want): + npt.assert_allclose(np.asarray(g, float), w, rtol=0, atol=0) + assert len(want) > 10, "the record should actually contain bursts" + + +@pytest.mark.parametrize("tcrit", [1e-4, 5e-4, 3.5e-3, 2e-2]) +def test_agrees_at_every_tcrit(tcrit): + t, a = alternating(2000, seed=7) + got = extract_burst_intervals(t, a, tcrit) + want = fortran_bursts(t, a, tcrit) + assert len(got) == len(want) + for g, w in zip(got, want): + npt.assert_allclose(np.asarray(g, float), w, rtol=0, atol=0) + + +def test_the_comparison_is_strictly_greater(): + """A gap exactly equal to tcrit is within-burst, as line 1009 has it.""" + t = np.array([1e-4, TCRIT, 1e-4, TCRIT * 1.000001, 1e-4]) + a = np.array([5.0, 0.0, 5.0, 0.0, 5.0]) + + got = extract_burst_intervals(t, a, TCRIT) + want = fortran_bursts(t, a, TCRIT) + assert len(got) == len(want) == 2, "only the longer gap should split" + npt.assert_allclose(np.asarray(got[0], float), want[0]) + + +def test_every_burst_starts_and_ends_on_an_opening(): + """"end present group with prev opening" -- the separating gap is dropped.""" + t, a = alternating(2000, seed=11) + for g in extract_burst_intervals(t, a, TCRIT): + g = np.asarray(g, float) + assert g.size % 2 == 1, "odd length means opening at both ends" + + +def test_tcrit_of_zero_splits_at_every_shut_time(): + """The trap: line 1009 is not guarded by burst(jset). + + HJCFIT means "do not divide" by setting tcrit to a year, not to zero. A + caller who passes zero expecting "no division" gets one group per opening, + silently, and this records that both implementations behave that way. + """ + t, a = alternating(200, seed=13) + got = extract_burst_intervals(t, a, 0.0) + want = fortran_bursts(t, a, 0.0) + assert len(got) == len(want) + assert all(np.asarray(g, float).size == 1 for g in got) + assert len(got) == int((a != 0).sum()) + + +def test_an_enormous_tcrit_is_how_you_ask_for_one_group(): + """3.1536e10 ms, one year: Hjcfit1-09122003.for line 1568.""" + t, a = alternating(200, seed=17) + one_year = 3.1536e10 * 1e-3 # ms -> s + got = extract_burst_intervals(t, a, one_year) + want = fortran_bursts(t, a, one_year) + assert len(got) == len(want) == 1 + npt.assert_allclose(np.asarray(got[0], float), want[0])