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