From 447ab491707c310c6cb15c515d02b9ec48512e6d Mon Sep 17 00:00:00 2001 From: Remis Date: Wed, 2 Sep 2026 10:38:06 +0100 Subject: [PATCH] ekscn: read through dcio M2 of the record-layer plan. EKDIST carried a second implementation of the SCN binary decoding, and two copies of a format parser is two places for a record to be read differently. It was: the interval array is stored as 4-byte floats, and under NEP 50 the millisecond-to-second multiply stayed in single precision here, discarding most of the digits the file carries. read_data now delegates to dcio.formats.scn.read. Verified identical on both shipped records before the change -- intervals equal bit for bit, raw amplitudes and flags equal, and the same trailing-gap trim. The contract is unchanged: amplitudes come back as raw stored integers, not pA, because that is what this has always returned and what load_SCN_file scales by calfac2. dcio applies the calibration on read, so it is divided back out. The round trip is exercised on 181003S8.SCN, whose calfac2 is 0.0025 rather than 1, with a test asserting that it stays non-unit so the check cannot become vacuous. read_header stays EKDIST's own. It parses about 130 header fields -- the full SCAN header including display and fitting parameters -- where dcio names 18. Shimming it would lose information dcio does not currently model, so the duplication removed here is the data decoding, which is where the numerical disagreement lived. Co-Authored-By: Claude Opus 5 --- ekdist/ekscn.py | 83 +++++++++++++++++++--------------------- tests/test_ekscn_shim.py | 74 +++++++++++++++++++++++++++++++++++ 2 files changed, 113 insertions(+), 44 deletions(-) create mode 100644 tests/test_ekscn_shim.py diff --git a/ekdist/ekscn.py b/ekdist/ekscn.py index 8f7ca66..36cc3b2 100644 --- a/ekdist/ekscn.py +++ b/ekdist/ekscn.py @@ -603,52 +603,47 @@ def read_header (fname, verbose=False): return header def read_data(fname, header): + """Read idealised intervals, amplitudes and flags from an SCN file. + + Delegates to :func:`dcio.formats.scn.read`, which is the one SCN reader in + the DCPROGS stack. EKDIST had a second implementation of the same binary + decoding, and two copies of a format parser is two places for a record to + be read differently -- which is what happened: the interval array is + stored as 4-byte floats, and under NEP 50 the millisecond-to-second + multiply stayed in single precision here, discarding most of the digits + the file actually carries. + + The contract is unchanged. Amplitudes come back as the raw stored integers, + not in pA, because that is what this has always returned and what + ``SingleChannelRecord.load_SCN_file`` scales by ``calfac2``. dcio applies + the calibration on read, so it is divided back out; the values are small + integers and the round trip is exact. + + Parameters + ---------- + fname : str + Path to the SCN file. + header : dict + As returned by :func:`read_header`. Only ``calfac2`` is used; dcio + reads the file's own header for everything else. + + Returns + ------- + intervals : ndarray of float64 + Durations in seconds. + amplitudes : ndarray + Raw stored amplitudes, before calibration. + flags : ndarray + Property flags; bit 3 set means the duration is unusable. """ - Read idealised data- intervals, amplitudes, flags- rom SCN file. - Data= - real*4 tint(1...nint) 4nint bytes - integer*2 iampl(1..nint) 2nint bytes - integer*1 iprops(1..nint) nint bytes - Total storage needed = 7 * nint bytes - integer*1 iprops(i) holds properties of ith duration and amplitude - (integer*1 has range -128 to +127 (bit 7 set gives -128; can use bits 0-6) - 0 = all OK; - 1 = amplitude dubious = bit 0; - 2 = amplitude fixed = bit 1; - 4 = amplitude of opening constrained (see fixamp) = bit 2; - 8 = duration unusable = bit 3; etc - and keep sum of values of more than one property is true. - """ + from dcio.formats.scn import read as _dcio_read + + record = _dcio_read(fname) + calfac = record.header.calfac2 or 1.0 + return (np.asarray(record.intervals, dtype=np.float64), + np.asarray(record.amplitudes) / calfac, + np.asarray(record.flags)) - tint = array ('f') # 4 byte float - iampl = array ('h') # 2 byte integer - iprops = array('b') # 1 byte integer - - f=open(fname, 'rb') - f.seek(header['ioffset']-1) - tint.fromfile(f, header['nint']) - iampl.fromfile(f, header['nint']) - iprops.fromfile(f, header['nint']) - f.close() - - if header['iscanver'] > 0: - gapnotfound = True - while gapnotfound: - if iampl[-1] == 0: - gapnotfound = False - iprops[-1] = 8 - else: - tint.pop() - iampl.pop() - iprops.pop() - - # float64 before scaling, not after. tint is array('f'), and under NEP 50 - # a float32 array times a Python float stays float32, so the ms-to-s - # multiply was done in single precision and threw away most of the - # digits the stored value carries. Over A-10 that is 2e-5 s in a 435 s - # record -- 4.6e-8 relative, and it reached the fitted likelihood. - return (np.array(tint, dtype=np.float64) * 1e-3, - np.array(iampl), np.array(iprops)) def write(intervals, amplitudes, flags, calfac=1.0, ffilt=-1.0, rms=0.0, treso=0.0, tresg=0.0, Emem=0.0, diff --git a/tests/test_ekscn_shim.py b/tests/test_ekscn_shim.py new file mode 100644 index 0000000..91b7d77 --- /dev/null +++ b/tests/test_ekscn_shim.py @@ -0,0 +1,74 @@ +"""ekscn reads through dcio. + +M2 of the record-layer plan. EKDIST carried a second implementation of the SCN +binary decoding, and two copies of a format parser is two places for a record +to be read differently. It was: the interval array is stored as 4-byte floats, +and under NEP 50 the millisecond-to-second multiply stayed in single precision +here, discarding most of the digits the file carries. +""" + +import numpy as np +import pytest + +from dcio.formats import scn as dcio_scn +from ekdist import ekrecord, ekscn + +FILES = ["./tests/AChsim.scn", "./tests/181003S8.SCN"] + + +@pytest.mark.parametrize("path", FILES) +class TestAgreesWithDcio: + + def test_intervals_are_identical(self, path): + """Not close -- identical. Both read the same bytes; any difference + would be a decoding difference, not rounding.""" + header = ekscn.read_header(path) + intervals, _, _ = ekscn.read_data(path, header) + assert np.array_equal(intervals, dcio_scn.read(path).intervals) + + def test_flags_are_identical(self, path): + header = ekscn.read_header(path) + _, _, flags = ekscn.read_data(path, header) + record = dcio_scn.read(path) + assert np.array_equal(flags.astype(int), record.flags.astype(int)) + + def test_amplitudes_round_trip_through_the_calibration(self, path): + """read_data returns raw stored amplitudes, as it always has, so + load_SCN_file's multiply by calfac2 stays correct. dcio applies the + calibration on read, so it is divided back out.""" + header = ekscn.read_header(path) + _, raw, _ = ekscn.read_data(path, header) + record = dcio_scn.read(path) + np.testing.assert_allclose(raw * record.header.calfac2, + record.amplitudes, rtol=1e-12) + + def test_trailing_partial_interval_trimmed_the_same_way(self, path): + """An experimental record is truncated back to a final shutting. Both + readers do it, and they must agree on how many intervals remain.""" + header = ekscn.read_header(path) + intervals, _, _ = ekscn.read_data(path, header) + assert len(intervals) == len(dcio_scn.read(path).intervals) + + +class TestRecordStillLoads: + + @pytest.mark.parametrize("path", FILES) + def test_load_scn_file_end_to_end(self, path): + rec = ekrecord.SingleChannelRecord() + rec.load_SCN_file(path) + assert len(rec.itint) > 0 + assert len(rec.itint) == len(rec.iampl) == len(rec.iprop) + + def test_amplitudes_are_calibrated_exactly_once(self): + """The file with a non-unit calfac2 is the one that would expose a + double or missing application.""" + path = "./tests/181003S8.SCN" + rec = ekrecord.SingleChannelRecord() + rec.load_SCN_file(path) + np.testing.assert_allclose(rec.iampl, dcio_scn.read(path).amplitudes, + rtol=1e-12) + + def test_calfac2_is_not_unity_in_that_file(self): + """Guards the test above: if this ever became 1.0, the check would + pass while proving nothing.""" + assert dcio_scn.read("./tests/181003S8.SCN").header.calfac2 != 1.0