ekscn: read through dcio (M2)

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,

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