seismo-relay/tests/test_idf_binary_codec.py

"""Per-sample verification of the Thor / Micromate (series-4) IDF binary codec.

Ground truth is Thor's own CSV export, written next to each binary by the
Thor desktop application.  For waveforms the export carries a per-sample
block of four columns (Tran, Vert, Long, Mic) in in/s and psi -- the
series-4 equivalent of Blastware's ``_ASCII.TXT`` exports.

The full-corpus harness is ``scratch/verify_thor_against_csv.py``; these
tests pin the two constants that harness established so they cannot
regress silently.
"""
from __future__ import annotations

import csv
import os
import sys
from pathlib import Path

import pytest

sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))

from micromate.idf_file import (
    _GEO_LSB_IPS,
    geo_count_to_ips,
    read_idf_file,
)

FIXTURES = Path(__file__).parent / "fixtures" / "thor-idf"
IDFW = FIXTURES / "UM11719_20231219162723.IDFW"
IDFH = FIXTURES / "UM11719_20231219162648.IDFH"

GEO_CHANNELS = ("Tran", "Vert", "Long")

# tests/fixtures/ is gitignored, so a fresh checkout has no sample data.
# Skip rather than fail, matching test_idf_ascii_report.py.  To populate:
#
#   B="<thor-watcher>/example-data/THORDATA_example/THORDATA_example/UPMC Presby"
#   mkdir -p tests/fixtures/thor-idf
#   for f in UM11719/UM11719_20231219162723.IDFW \
#            UM11719/UM11719_20231219162648.IDFH \
#            UM13981/UM13981_20220207084555.IDFW \
#            UM13981/UM13981_20220207183102.IDFH \
#            UM13981/UM13981_20221202063059.IDFH; do
#     cp "$B/$f" tests/fixtures/thor-idf/
#     cp "$B/$(dirname $f)/CSV/$(basename $f).csv" tests/fixtures/thor-idf/
#   done
pytestmark = pytest.mark.skipif(
    not FIXTURES.is_dir() or not any(FIXTURES.glob("*.IDFW")),
    reason=f"Thor IDF fixtures not present under {FIXTURES}",
)


def _parse_export(path: Path):
    """Split a Thor CSV export into (header dict, per-sample rows)."""
    header, rows = {}, []
    with path.open(newline="", encoding="utf-8", errors="replace") as fh:
        for rec in csv.reader(fh):
            if len(rec) == 2:
                header[rec[0].strip()] = rec[1].strip()
            elif len(rec) >= 3:
                try:
                    rows.append([float(x) for x in rec])
                except ValueError:
                    pass
    return header, rows


def _header_float(header, key):
    return float(header[key].split()[0])


@pytest.fixture(scope="module")
def idfw_export():
    return _parse_export(IDFW.with_suffix(".IDFW.csv"))


# ─── The geo scale constant ────────────────────────────────────────────────


def test_geo_lsb_matches_thor_quantisation():
    """Thor's own export quantises geo samples to this LSB.

    Derived by maximising exact-match count over 1,046,016 paired samples
    (454 channel-events, 2 units); independently corroborated on 8
    production units via their device-reported PPV.  The historical value
    0.0003 read every series-4 geophone sample 3.3% low.
    """
    assert _GEO_LSB_IPS == pytest.approx(0.000310308, rel=1e-6)


def test_geo_lsb_is_not_the_legacy_value():
    # Guards against a revert to the truncated 0.0003 constant.
    assert abs(_GEO_LSB_IPS - 0.0003) > 1e-6


# ─── Per-sample fidelity ───────────────────────────────────────────────────


def test_waveform_channel_lengths_match_export(idfw_export):
    _header, rows = idfw_export
    result = read_idf_file(IDFW)
    for channel in GEO_CHANNELS:
        assert len(result.samples[channel]) == len(rows), (
            f"{channel} truncated: decoded {len(result.samples[channel])} "
            f"samples, export has {len(rows)}"
        )


def test_waveform_samples_match_export_exactly(idfw_export):
    """Every geo sample must reproduce Thor's exported value to 4 dp."""
    _header, rows = idfw_export
    result = read_idf_file(IDFW)
    for index, channel in enumerate(GEO_CHANNELS):
        decoded = result.samples[channel]
        expected = [row[index] for row in rows]
        mismatches = [
            (i, geo_count_to_ips(c), v)
            for i, (c, v) in enumerate(zip(decoded, expected))
            if abs(geo_count_to_ips(c) - v) >= 5e-5
        ]
        assert not mismatches, (
            f"{channel}: {len(mismatches)} of {len(expected)} samples differ; "
            f"first three {mismatches[:3]}"
        )


def test_waveform_ppv_matches_export(idfw_export):
    header, _rows = idfw_export
    result = read_idf_file(IDFW)
    for channel, attr in (
        ("Tran", "transverse_ips"),
        ("Vert", "vertical_ips"),
        ("Long", "longitudinal_ips"),
    ):
        decoded = getattr(result.event.peaks, attr)
        assert decoded == pytest.approx(
            _header_float(header, f"{channel}PPV"), abs=5e-5
        ), f"{channel} PPV disagrees with Thor's export"


# ─── Histogram path shares the same scale ──────────────────────────────────


def test_histogram_peaks_match_export():
    header, _rows = _parse_export(IDFH.with_suffix(".IDFH.csv"))
    result = read_idf_file(IDFH)
    assert result.intervals, "IDFH decoded no intervals"
    for channel, attr in (
        ("Tran", "transverse_ips"),
        ("Vert", "vertical_ips"),
        ("Long", "longitudinal_ips"),
    ):
        decoded = getattr(result.event.peaks, attr)
        expected = _header_float(header, f"{channel}PPV")
        # Histogram peaks are stored per-interval, so the export's PPV is
        # reproduced within one quantisation step rather than exactly.
        assert decoded == pytest.approx(expected, abs=2 * _GEO_LSB_IPS), (
            f"{channel} histogram peak {decoded} vs export {expected}"
        )


# ─── Regressions found 2026-09-10 ──────────────────────────────────────────

IDFH_LONG = FIXTURES / "UM13981_20220207183102.IDFH"      # 719 intervals
IDFH_SENTINEL = FIXTURES / "UM13981_20221202063059.IDFH"  # holds an unwritten slot
IDFW_RAW16 = FIXTURES / "UM13981_20220207084555.IDFW"     # segment 0 is MODE_RAW16


def test_histogram_decodes_past_250_intervals():
    """The segment validator must not require a zero counter high byte.

    The interval counter is a uint16 cumulative index.  Requiring its high
    byte to be zero rejected every segment past interval 255, capping each
    histogram at 250 intervals and truncating any run longer than ~4 hours —
    frequently discarding the part that held the peak.
    """
    result = read_idf_file(IDFH_LONG)
    header, _rows = _parse_export(IDFH_LONG.with_suffix(".IDFH.csv"))
    expected = float(header["NumberOfIntervals"])
    assert len(result.intervals) == 719
    assert len(result.intervals) == pytest.approx(expected, abs=1.0)


def test_histogram_ignores_unwritten_interval_slot():
    """A never-written interval keeps its ±full-scale seed and must be dropped.

    Counting it fabricates a 10.0 in/s peak on every channel, which then wins
    the max-over-intervals and poisons the whole file's PPV.
    """
    header, _rows = _parse_export(IDFH_SENTINEL.with_suffix(".IDFH.csv"))
    result = read_idf_file(IDFH_SENTINEL)
    for channel, attr in (
        ("Tran", "transverse_ips"),
        ("Vert", "vertical_ips"),
        ("Long", "longitudinal_ips"),
    ):
        decoded = getattr(result.event.peaks, attr)
        assert decoded < 1.0, f"{channel} peak {decoded} looks like the ±FS seed"
        assert decoded == pytest.approx(
            _header_float(header, f"{channel}PPV"), abs=2 * _GEO_LSB_IPS
        )


def test_waveform_raw16_segment_zero_is_decoded():
    """Segment-0 records can be raw int16 (MODE_RAW16, 10-byte header).

    That mode was absent from the dispatch, so the record fell through
    unhandled and the channel silently lost its first 512 samples.
    """
    rows = _parse_export(IDFW_RAW16.with_suffix(".IDFW.csv"))[1]
    result = read_idf_file(IDFW_RAW16)
    for index, channel in enumerate(GEO_CHANNELS):
        decoded = result.samples[channel]
        assert len(decoded) == len(rows), f"{channel} lost segment 0"
        expected = [row[index] for row in rows]
        bad = sum(
            1 for c, v in zip(decoded, expected)
            if abs(geo_count_to_ips(c) - v) >= 5e-5
        )
        assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"


def test_body_offset_search_is_not_quadratic():
    """The body scan must stay cheap enough for bulk ingest.

    MODE_RAW16 is (0x00, 0x00), so scanning for candidate *preambles* treats
    every run of three zero bytes as a body start and trial-decodes each one
    (~0.5 s/file measured).  The search anchors on record headers instead.
    """
    import time

    start = time.perf_counter()
    for _ in range(3):
        read_idf_file(IDFW_RAW16)
    elapsed = (time.perf_counter() - start) / 3
    assert elapsed < 0.15, f"body-offset search took {elapsed*1000:.0f} ms/file"


# ─── Mic-disabled (3-channel) units, found 2026-09-10 ──────────────────────

IDFW_3CH = FIXTURES / "UM20147_20250531135901.IDFW"   # body head below old floor
IDFH_3CH = FIXTURES / "UM20147_20250330070110.IDFH"   # 56-byte interval records


def test_three_channel_waveform_decodes_all_geo_channels():
    """A mic-disabled unit's shorter header moves the record chain head.

    Its head sits at 0x0dba, below the old ``_BODY_SCAN_FLOOR`` of 0x0E00, so
    the scan could not see it and fell through to the *Vert* segment-0 record
    — decoding a body shifted one position around the channel rotation, which
    surfaced as Vert being exactly 512 samples short.
    """
    rows = _parse_export(IDFW_3CH.with_suffix(".IDFW.csv"))[1]
    result = read_idf_file(IDFW_3CH)
    for index, channel in enumerate(GEO_CHANNELS):
        decoded = result.samples[channel]
        assert len(decoded) == len(rows), (
            f"{channel}: {len(decoded)} samples, export has {len(rows)}"
        )
        expected = [row[index] for row in rows]
        bad = sum(
            1 for c, v in zip(decoded, expected)
            if abs(geo_count_to_ips(c) - v) >= 5e-5
        )
        assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"
    # Mic is genuinely absent on these units, not merely undecoded.
    assert not result.samples.get("MicL")


def test_three_channel_histogram_uses_56_byte_intervals():
    """Interval stride is 16 bytes per channel + an 8-byte tail, not a constant.

    A mic-disabled unit packs 56-byte records, so assuming 72 read 7 intervals
    out of every 10-interval segment and then walked off alignment into
    garbage, which decoded as ~10 in/s peaks.  The true count comes from the
    segment's cumulative interval counter.
    """
    header, _rows = _parse_export(IDFH_3CH.with_suffix(".IDFH.csv"))
    result = read_idf_file(IDFH_3CH)
    expected_intervals = float(header["NumberOfIntervals"])
    assert len(result.intervals) == pytest.approx(expected_intervals, abs=1.0)
    assert {iv.n_channels for iv in result.intervals} == {3}
    for channel, attr in (
        ("Tran", "transverse_ips"),
        ("Vert", "vertical_ips"),
        ("Long", "longitudinal_ips"),
    ):
        decoded = getattr(result.event.peaks, attr)
        assert decoded < 1.0, f"{channel} peak {decoded} looks like walked-off garbage"
        assert decoded == pytest.approx(
            _header_float(header, f"{channel}PPV"), rel=0.02
        )


# ─── `40 NN` blocks with NN > 8, verified 2026-09-11 ───────────────────────

IDFW_WIDE40 = FIXTURES / "UM12947_20250806134504.IDFW"


def test_wide_forty_nn_block_does_not_truncate_channels():
    """Loud events use `40 NN` blocks with NN well above the old cap of 8.

    ``data_block_len()`` rejected NN > 0x08, which halted the block walk
    part-way through a record.  The walker stops at the first unrecognised
    tag instead of raising, so this surfaced as silently short channels —
    here Tran 1812 / Vert 2132 / Long 2324 where the export has 2324 for all
    three.  The affected files use NN of 12, 16, 20 ... up to 196.
    """
    rows = _parse_export(IDFW_WIDE40.with_suffix(".IDFW.csv"))[1]
    result = read_idf_file(IDFW_WIDE40)
    for index, channel in enumerate(GEO_CHANNELS):
        decoded = result.samples[channel]
        assert len(decoded) == len(rows), (
            f"{channel}: {len(decoded)} samples, export has {len(rows)}"
        )
        expected = [row[index] for row in rows]
        bad = sum(
            1 for c, v in zip(decoded, expected)
            if abs(geo_count_to_ips(c) - v) >= 5e-5
        )
        assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"