rm3100/tests/test_calibrate.py

"""Calibration: the timing fit, the gain factor, and the file it is carried in.

Everything is checked against a planted answer. A calibration is a claim about a
particular part, so the tests that matter are the ones showing a value goes in
and the same value comes back out -- and that a file which cannot be fully
understood is refused rather than half-read into a plausible-looking correction.
"""

import json
import math

import numpy as np
import pytest

import calibrate
import capture
import rm3100
from synthetic import drifting_times, field_counts, write_capture

FIELD = (-40000.0, 16000.0, 12000.0)


def planted_period(cycle_count, counts_per_second, overhead):
    return rm3100.AXES * (cycle_count / counts_per_second + overhead)


def make(path, cycle_count=100, counts_per_second=88546.0, overhead=40.61e-6,
         rows=2000, mean=FIELD, noise_nt=5.0, seed=0):
    counts = field_counts(rows, mean, cycle_count, noise_nt, seed)
    return capture.load(write_capture(
        path, rows=rows, cycle_count=cycle_count, counts=counts,
        dt=planted_period(cycle_count, counts_per_second, overhead)))


def written(tmp_path, **fields):
    body = {"rm3100_calibration": calibrate.CALIBRATION_VERSION,
            "counts_per_second": 88546.0, "axis_overhead_s": 40.61e-6}
    body.update(fields)
    path = tmp_path / "cal.json"
    path.write_text(json.dumps(body))
    return str(path)


# --------------------------------------------------------------------------
# fit_timing
# --------------------------------------------------------------------------

def test_fit_timing_recovers_a_planted_pair_exactly():
    counts_per_second, overhead = 88546.0, 40.61e-6
    ccs = [100, 400]
    periods = [planted_period(cc, counts_per_second, overhead) for cc in ccs]
    fit = calibrate.fit_timing(ccs, periods, period_sd=1e-4)
    assert fit.counts_per_second == pytest.approx(counts_per_second, rel=1e-9)
    assert fit.axis_overhead_s == pytest.approx(overhead, rel=1e-9)
    # Two points fit both terms exactly, so there is nothing left over -- and
    # that is precisely why the residual cannot be used as a check.
    assert fit.residual_ppm == 0.0


def test_fit_timing_least_squares_over_more_points():
    counts_per_second, overhead = 92889.0, 38.11e-6
    ccs = [50, 100, 200, 400, 800]
    periods = [planted_period(cc, counts_per_second, overhead) for cc in ccs]
    fit = calibrate.fit_timing(ccs, periods)
    assert fit.counts_per_second == pytest.approx(counts_per_second, rel=1e-6)
    assert fit.axis_overhead_s == pytest.approx(overhead, rel=1e-6)


def test_fit_timing_needs_two_distinct_cycle_counts():
    """Two unknowns need two points; one would be a guess dressed as a fit."""
    with pytest.raises(calibrate.CalibrationError,
                       match="two distinct cycle counts"):
        calibrate.fit_timing([100, 100, 100], [3.5e-3] * 3, period_sd=1e-4)


def test_fit_timing_rejects_data_that_does_not_follow_the_model():
    """A higher cycle count that samples faster is not this chip."""
    with pytest.raises(calibrate.CalibrationError, match="not positive"):
        calibrate.fit_timing([100, 400], [13.0e-3, 3.5e-3], period_sd=1e-4)


def test_fit_timing_rejects_mismatched_lengths():
    with pytest.raises(calibrate.CalibrationError, match="against"):
        calibrate.fit_timing([100, 400], [3.5e-3], period_sd=1e-4)


def test_holding_the_divisor_at_spec_cannot_absorb_the_misfit():
    """The reason both terms are calibrated, not just the overhead.

    Solving for the overhead alone with the spec divisor gives a different
    answer at every cycle count, which is what makes the nominal model's error
    cycle-count dependent.
    """
    counts_per_second, overhead = 88546.0, 40.61e-6
    solved = [planted_period(cc, counts_per_second, overhead) / rm3100.AXES
              - cc / rm3100.COUNTS_PER_SECOND for cc in (100, 400)]
    assert solved[0] == pytest.approx(58.9e-6, abs=1e-6)
    assert solved[1] == pytest.approx(113.6e-6, abs=1e-6)


# --------------------------------------------------------------------------
# oscillator_hz
# --------------------------------------------------------------------------

@pytest.mark.parametrize("cycle_count", [50, 100, 200, 400, 800])
def test_oscillator_is_recovered_from_a_single_capture(tmp_path, cycle_count):
    """One capture is enough once the overhead is known -- the whole point."""
    counts_per_second, overhead = 88546.0, 40.61e-6
    cap = make(tmp_path / "c.csv", cycle_count=cycle_count,
               counts_per_second=counts_per_second, overhead=overhead)
    cal = calibrate.calibration(counts_per_second, overhead)
    assert calibrate.oscillator_hz(cap, cal).value == pytest.approx(
        counts_per_second, rel=1e-5)


def test_overhead_error_costs_less_than_a_tenth_of_a_percent(tmp_path):
    """The sensitivity the single-capture path depends on.

    A 1 us error in the overhead must stay under 0.1% of the oscillator at the
    worst cycle count in use, or a calibration from one supply could not be
    applied to a capture from another.
    """
    counts_per_second, overhead = 88546.0, 40.61e-6
    cap = make(tmp_path / "c.csv", cycle_count=100,
               counts_per_second=counts_per_second, overhead=overhead)
    exact = calibrate.oscillator_hz(
        cap, calibrate.calibration(counts_per_second, overhead)).value
    off_by_1us = calibrate.oscillator_hz(
        cap, calibrate.calibration(counts_per_second, overhead + 1e-6)).value
    assert abs(off_by_1us / exact - 1) < 0.001


def test_the_nominal_overhead_costs_far_more_and_varies_with_cycle_count(tmp_path):
    """Why the nominal model cannot be used as a gain reference."""
    counts_per_second, overhead = 88546.0, 40.61e-6
    errors = []
    for cycle_count in (100, 400):
        cap = make(tmp_path / f"{cycle_count}.csv", cycle_count=cycle_count,
                   counts_per_second=counts_per_second, overhead=overhead)
        nominal = calibrate.calibration(counts_per_second,
                                        rm3100.AXIS_OVERHEAD_S)
        errors.append(calibrate.oscillator_hz(cap, nominal).value
                      / counts_per_second)
    # Both wrong, and wrong by different amounts -- so it does not cancel.
    assert abs(errors[0] - 1) > 0.02
    assert abs(errors[1] - 1) < 0.01


def test_oscillator_refuses_a_calibration_that_cannot_describe_the_capture(tmp_path):
    """An overhead longer than the whole per-axis time is not merely wrong."""
    cap = make(tmp_path / "c.csv", cycle_count=100)
    absurd = calibrate.calibration(88546.0, axis_overhead_s=1.0)
    with pytest.raises(calibrate.CalibrationError, match="does not describe"):
        calibrate.oscillator_hz(cap, absurd)


# --------------------------------------------------------------------------
# gain_factor
# --------------------------------------------------------------------------

def test_a_capture_at_its_own_reference_is_not_corrected(tmp_path):
    counts_per_second, overhead = 88546.0, 40.61e-6
    cap = make(tmp_path / "c.csv", counts_per_second=counts_per_second,
               overhead=overhead)
    cal = calibrate.calibration(counts_per_second, overhead)
    assert calibrate.gain_factor(cap, cal).value == pytest.approx(1.0, rel=1e-5)


def test_gain_factor_is_the_oscillator_ratio_at_exponent_one(tmp_path):
    overhead = 40.61e-6
    cap = make(tmp_path / "c.csv", counts_per_second=92889.0, overhead=overhead)
    cal = calibrate.calibration(88546.0, overhead, reference_oscillator_hz=88546.0)
    assert calibrate.gain_factor(cap, cal).value == pytest.approx(
        92889.0 / 88546.0, rel=1e-4)


def test_a_zero_exponent_disables_the_oscillator_term(tmp_path):
    overhead = 40.61e-6
    cap = make(tmp_path / "c.csv", counts_per_second=92889.0, overhead=overhead)
    cal = calibrate.calibration(88546.0, overhead, gain_exponent=0.0)
    assert calibrate.gain_factor(cap, cal).value == 1.0


def test_the_correction_is_exactly_one_at_the_reference_cycle_count(tmp_path):
    """The offset is only known up to a scale, so the reference pins it."""
    counts_per_second, overhead = 88546.0, 40.61e-6
    for reference in (50, 100, 400):
        cap = make(tmp_path / f"{reference}.csv", cycle_count=reference,
                   counts_per_second=counts_per_second, overhead=overhead)
        for offset in (-0.5, 0.0, 0.9, 4.086):
            cal = calibrate.calibration(
                counts_per_second, overhead, reference_cycle_count=reference,
                gain_offset_counts=offset)
            assert calibrate.gain_factor(cap, cal).value == pytest.approx(
                1.0, rel=1e-5)


def test_the_offset_corrects_the_datasheet_gain_shape(tmp_path):
    """The datasheet implies gain ~ (cc + 4.086); this unit measures ~0.

    Two captures whose counts embed gain = A(cc + n) must read the same field
    once corrected with that n, and must not with the datasheet's.
    """
    counts_per_second, overhead, n = 88546.0, 40.61e-6, 0.0
    field = 47000.0
    raw, good, bad = {}, {}, {}
    for cycle_count in (100, 400):
        # The counts this chip would produce if its gain really went as (cc+n).
        per_axis = field / math.sqrt(3) * 0.3671 * (cycle_count + n)
        counts = np.tile([per_axis, per_axis, per_axis], (2000, 1))
        cap = capture.load(write_capture(
            tmp_path / f"{cycle_count}.csv", rows=2000,
            cycle_count=cycle_count, counts=counts,
            dt=planted_period(cycle_count, counts_per_second, overhead)))
        raw[cycle_count] = cap.total.mean()
        for name, offset, into in (("good", n, good), ("bad", 4.086, bad)):
            cal = calibrate.calibration(counts_per_second, overhead,
                                        reference_cycle_count=100,
                                        gain_offset_counts=offset)
            into[cycle_count] = cap.total.mean() * calibrate.gain_factor(cap, cal).value

    # Uncorrected the two disagree, because gain_model has the wrong offset.
    assert abs(raw[400] / raw[100] - 1) > 0.02
    assert good[400] == pytest.approx(good[100], rel=1e-4)
    assert abs(bad[400] / bad[100] - 1) > 0.02


def test_a_gain_offset_that_zeroes_the_gain_is_refused(tmp_path):
    cap = make(tmp_path / "c.csv", cycle_count=100)
    cal = calibrate.calibration(88546.0, 40.61e-6, reference_cycle_count=100,
                                gain_offset_counts=-100.0)
    with pytest.raises(calibrate.CalibrationError, match="at or below zero"):
        calibrate.gain_factor(cap, cal)


# --------------------------------------------------------------------------
# Propagated uncertainty
# --------------------------------------------------------------------------

def test_zero_input_uncertainty_gives_zero_output_uncertainty(tmp_path):
    """A calibration claiming perfect knowledge says so, rather than guessing."""
    cap = make(tmp_path / "c.csv", cycle_count=400)
    cal = calibrate.calibration(88546.0, 40.61e-6)
    assert calibrate.gain_factor(cap, cal, period_sd=0.0).sd == 0.0
    assert calibrate.oscillator_hz(cap, cal, period_sd=0.0).sd == 0.0


def test_the_offset_uncertainty_vanishes_at_the_reference_cycle_count(tmp_path):
    """At the reference the shape term is 1 by construction, so it contributes
    nothing however badly the offset is known."""
    cap = make(tmp_path / "c.csv", cycle_count=100)
    cal = calibrate.calibration(88546.0, 40.61e-6, reference_cycle_count=100,
                                gain_offset_counts_sd=10.0)
    assert calibrate.gain_factor(cap, cal, period_sd=0.0).sd == pytest.approx(0.0)


def test_the_offset_uncertainty_scales_linearly_away_from_the_reference(tmp_path):
    cap = make(tmp_path / "c.csv", cycle_count=400)
    one = calibrate.calibration(88546.0, 40.61e-6, reference_cycle_count=100,
                                gain_offset_counts_sd=0.5)
    two = calibrate.calibration(88546.0, 40.61e-6, reference_cycle_count=100,
                                gain_offset_counts_sd=1.0)
    a = calibrate.gain_factor(cap, one, period_sd=0.0)
    b = calibrate.gain_factor(cap, two, period_sd=0.0)
    assert b.sd == pytest.approx(2 * a.sd, rel=1e-9)
    assert a.relative > 0.001          # and it is not negligible


def test_the_oscillator_term_scales_with_the_exponent(tmp_path):
    cap = make(tmp_path / "c.csv", cycle_count=100, counts_per_second=92889.0)
    one = calibrate.calibration(88546.0, 40.61e-6, axis_overhead_s_sd=1e-6,
                                gain_exponent=1.0)
    two = calibrate.calibration(88546.0, 40.61e-6, axis_overhead_s_sd=1e-6,
                                gain_exponent=2.0)
    a = calibrate.gain_factor(cap, one, period_sd=0.0)
    b = calibrate.gain_factor(cap, two, period_sd=0.0)
    assert b.relative == pytest.approx(2 * a.relative, rel=1e-9)


def test_rate_stability_is_zero_for_a_constant_period(tmp_path):
    """Not exactly zero: the capture format stores host time to a microsecond,
    which is ~0.05 ppm of scatter here. Real drifts run 186 to 1241 ppm, so a
    1 ppm bar separates "nothing" from anything worth reporting."""
    cap = make(tmp_path / "c.csv", rows=4000)
    assert calibrate.rate_stability(cap) < 1e-6


def test_rate_stability_recovers_a_planted_ramp(tmp_path):
    """The oscillator warms up and slows; this is what measures that."""
    rows, dt, ramp = 8000, 1 / 300.0, 0.002        # 2000 ppm across the run
    index = np.arange(rows)
    # Period growing linearly, so the instantaneous rate falls by `ramp`.
    times = np.cumsum(dt * (1 + ramp * index / rows))
    path = tmp_path / "c.csv"
    write_capture(path, rows=rows, dt=dt)
    text = path.read_text().splitlines()
    header = [l for l in text if l.startswith("#")] + [text[len(
        [l for l in text if l.startswith("#")])]]
    body = text[len(header):]
    rebuilt = header + [
        ",".join([row.split(",")[0], f"{1_700_000_000.0 + t:.6f}"]
                 + row.split(",")[2:])
        for row, t in zip(body, times)]
    path.write_text("\n".join(rebuilt) + "\n")
    cap = capture.load(path)
    assert calibrate.rate_stability(cap) == pytest.approx(ramp, rel=0.15)


def test_a_two_point_fit_refuses_to_invent_a_confidence():
    """It has no residual, so it cannot estimate its own uncertainty."""
    ccs = [100, 400]
    periods = [planted_period(cc, 88546.0, 40.61e-6) for cc in ccs]
    with pytest.raises(calibrate.CalibrationError, match="two-point fit"):
        calibrate.fit_timing(ccs, periods)


def test_three_points_carry_a_residual_that_can_catch_a_bad_one():
    clean = [planted_period(cc, 88546.0, 40.61e-6) for cc in (100, 200, 400)]
    good = calibrate.fit_timing([100, 200, 400], clean)
    assert good.residual_ppm < 1.0
    nudged = list(clean)
    nudged[1] *= 1.002
    bad = calibrate.fit_timing([100, 200, 400], nudged)
    assert bad.residual_ppm > 100.0
    assert bad.counts_per_second_sd > good.counts_per_second_sd


# --------------------------------------------------------------------------
# The calibration file
# --------------------------------------------------------------------------

def test_a_calibration_round_trips(tmp_path):
    cal = calibrate.calibration(
        88546.0, 40.61e-6, reference_oscillator_hz=88000.0,
        reference_cycle_count=200, gain_exponent=1.4, gain_offset_counts=0.9,
        counts_per_second_sd=70.0, axis_overhead_s_sd=1.1e-6,
        gain_exponent_sd=0.2, gain_offset_counts_sd=0.5,
        note="bench", created="2026-08-24")
    path = tmp_path / "cal.json"
    calibrate.save_calibration(path, cal)
    assert calibrate.load_calibration(path) == cal


def test_the_reference_defaults_to_the_measured_oscillator():
    """A calibration corrects nothing at the condition it was taken at."""
    cal = calibrate.calibration(88546.0, 40.61e-6)
    assert cal.reference_oscillator_hz == 88546.0


def test_a_missing_version_is_refused(tmp_path):
    path = tmp_path / "cal.json"
    path.write_text(json.dumps({"counts_per_second": 1.0,
                                "axis_overhead_s": 0.0}))
    with pytest.raises(calibrate.CalibrationError, match="rm3100_calibration"):
        calibrate.load_calibration(path)


def test_a_future_version_is_refused_rather_than_guessed(tmp_path):
    with pytest.raises(calibrate.CalibrationError, match="rm3100_calibration"):
        calibrate.load_calibration(
            written(tmp_path, rm3100_calibration=calibrate.CALIBRATION_VERSION + 1))


def test_truncated_json_is_refused(tmp_path):
    path = tmp_path / "cal.json"
    path.write_text('{"rm3100_calibration": 1, "counts_per_second":')
    with pytest.raises(calibrate.CalibrationError, match="not valid JSON"):
        calibrate.load_calibration(path)


def test_a_json_list_is_refused(tmp_path):
    path = tmp_path / "cal.json"
    path.write_text("[1, 2, 3]")
    with pytest.raises(calibrate.CalibrationError, match="expected an object"):
        calibrate.load_calibration(path)


@pytest.mark.parametrize("field", ["counts_per_second", "axis_overhead_s"])
def test_a_missing_required_field_names_itself(tmp_path, field):
    body = {"rm3100_calibration": calibrate.CALIBRATION_VERSION,
            "counts_per_second": 88546.0, "axis_overhead_s": 40.61e-6}
    del body[field]
    path = tmp_path / "cal.json"
    path.write_text(json.dumps(body))
    with pytest.raises(calibrate.CalibrationError, match=field):
        calibrate.load_calibration(path)


@pytest.mark.parametrize("bad", [0, -1, "abc", None, True, [1], float("nan")])
def test_a_bad_count_rate_is_refused(tmp_path, bad):
    path = tmp_path / "cal.json"
    # json.dump cannot write nan as valid JSON, so write it literally.
    body = ('{"rm3100_calibration": 2, "axis_overhead_s": 4e-05, '
            f'"counts_per_second": {json.dumps(bad) if bad == bad else "NaN"}}}')
    path.write_text(body)
    with pytest.raises(calibrate.CalibrationError, match="counts_per_second"):
        calibrate.load_calibration(path)


@pytest.mark.parametrize("bad", ["abc", None, True])
def test_a_bad_exponent_is_refused(tmp_path, bad):
    with pytest.raises(calibrate.CalibrationError, match="gain_exponent"):
        calibrate.load_calibration(written(tmp_path, gain_exponent=bad))


def test_a_non_integer_reference_cycle_count_is_refused(tmp_path):
    with pytest.raises(calibrate.CalibrationError,
                       match="reference_cycle_count"):
        calibrate.load_calibration(
            written(tmp_path, reference_cycle_count="fast"))


def test_a_negative_uncertainty_is_refused(tmp_path):
    with pytest.raises(calibrate.CalibrationError,
                       match="gain_offset_counts_sd"):
        calibrate.load_calibration(
            written(tmp_path, gain_offset_counts_sd=-1.0))


def test_a_version_1_file_is_refused_with_an_explanation(tmp_path):
    """Version 1's per-cycle-count table has no version 2 equivalent."""
    path = tmp_path / "old.json"
    path.write_text(json.dumps({"rm3100_calibration": 1,
                                "counts_per_second": 88546.0,
                                "axis_overhead_s": 40.61e-6,
                                "gain_scale_by_cycle_count": {"400": 0.977}}))
    with pytest.raises(calibrate.CalibrationError, match="re-measure"):
        calibrate.load_calibration(path)


def test_true_is_not_accepted_as_a_number(tmp_path):
    """JSON's true floats to 1.0 in Python, which would pass silently."""
    with pytest.raises(calibrate.CalibrationError, match="counts_per_second"):
        calibrate.load_calibration(written(tmp_path, counts_per_second=True))


# --------------------------------------------------------------------------
# Conversion and provenance
# --------------------------------------------------------------------------

def read_calibrated(path):
    """Parse a calibrated CSV back into (header dict, rows)."""
    meta, rows = {}, []
    with open(path) as handle:
        for line in handle:
            if line.startswith("#"):
                key, _, value = line[1:].partition(":")
                meta[key.strip()] = value.strip()
            else:
                rows.append(line.rstrip("\n").split(","))
    return meta, rows[0], rows[1:]


def test_uncorrected_output_matches_capture_exactly(tmp_path):
    source = tmp_path / "raw.csv"
    cap = make(source)
    out = tmp_path / "out.csv"
    calibrate.write_calibrated(cap, out, source)
    meta, header, rows = read_calibrated(out)
    assert meta["gain_factor"] == "1.0"
    assert meta["calibration"] == "none"
    assert header == ["sample_index", "elapsed_s",
                      "x_nT", "y_nT", "z_nT", "total_nT", "warning"]
    assert len(rows) == len(cap.sample_index)
    assert float(rows[0][2]) == pytest.approx(cap.x[0], abs=5e-4)
    assert float(rows[0][5]) == pytest.approx(cap.total[0], abs=5e-4)


def test_the_total_column_is_the_norm_of_the_three(tmp_path):
    source = tmp_path / "raw.csv"
    cap = make(source)
    out = tmp_path / "out.csv"
    calibrate.write_calibrated(cap, out, source, factor=1.049)
    _, _, rows = read_calibrated(out)
    for row in rows[:50]:
        x, y, z, total = (float(v) for v in row[2:6])
        assert total == pytest.approx(math.sqrt(x*x + y*y + z*z), abs=1e-2)


def test_the_factor_scales_the_field(tmp_path):
    source = tmp_path / "raw.csv"
    cap = make(source)
    out = tmp_path / "out.csv"
    calibrate.write_calibrated(cap, out, source, factor=2.0)
    _, _, rows = read_calibrated(out)
    assert float(rows[0][2]) == pytest.approx(cap.x[0] * 2.0, abs=1e-3)


def test_the_source_digest_matches_and_changes_with_the_source(tmp_path):
    """The answer to a derived file drifting from its source unnoticed."""
    source = tmp_path / "raw.csv"
    cap = make(source)
    out = tmp_path / "out.csv"
    calibrate.write_calibrated(cap, out, source)
    meta, _, _ = read_calibrated(out)
    assert meta["source_sha256"] == calibrate.source_digest(source)

    source.write_text(source.read_text() + "\n")
    assert meta["source_sha256"] != calibrate.source_digest(source)


def test_flags_survive_the_conversion(tmp_path):
    source = tmp_path / "raw.csv"
    write_capture(source, rows=200,
                  flags={40: capture.WARN_MISSED,
                         41: f"{capture.WARN_MISSED} {capture.WARN_AMBIGUOUS}",
                         42: capture.WARN_AMBIGUOUS})
    cap = capture.load(source)
    out = tmp_path / "out.csv"
    calibrate.write_calibrated(cap, out, source)
    _, _, rows = read_calibrated(out)
    assert rows[40][6] == capture.WARN_MISSED
    assert rows[41][6] == f"{capture.WARN_MISSED} {capture.WARN_AMBIGUOUS}"
    assert rows[42][6] == capture.WARN_AMBIGUOUS
    assert rows[39][6] == ""


def test_a_calibrated_file_is_not_a_capture(tmp_path):
    """Different format, and capture.py must say so rather than misread it."""
    source = tmp_path / "raw.csv"
    cap = make(source)
    out = tmp_path / "out.csv"
    calibrate.write_calibrated(cap, out, source)
    with pytest.raises(capture.CaptureError, match="no capture header"):
        capture.load(out)


def test_the_header_records_the_calibration_that_was_applied(tmp_path):
    source = tmp_path / "raw.csv"
    cap = make(source, counts_per_second=92889.0, overhead=40.61e-6)
    cal = calibrate.calibration(88546.0, 40.61e-6)
    cal_path = tmp_path / "bench.json"
    calibrate.save_calibration(cal_path, cal)
    out = tmp_path / "out.csv"
    calibrate.write_calibrated(cap, out, source, cal,
                               calibrate.gain_factor(cap, cal), cal_path)
    meta, _, _ = read_calibrated(out)
    assert meta["calibration"] == "bench.json"
    assert float(meta["gain_factor"]) == pytest.approx(92889.0 / 88546.0,
                                                       rel=1e-4)
    assert float(meta["oscillator_hz"]) == pytest.approx(92889.0, rel=1e-4)
    assert float(meta["reference_oscillator_hz"]) == 88546.0


# --------------------------------------------------------------------------
# End to end
# --------------------------------------------------------------------------

def run(monkeypatch, capsys, argv):
    monkeypatch.setattr("sys.argv", ["calibrate.py"] + argv)
    assert calibrate.main() == 0
    return capsys.readouterr()


def test_cli_converts_without_a_calibration(tmp_path, monkeypatch, capsys):
    source = tmp_path / "raw.csv"
    make(source)
    out = tmp_path / "out.csv"
    result = run(monkeypatch, capsys, [str(source), "-o", str(out)])
    assert "correcting nothing" in result.out
    assert out.exists()


def test_cli_applies_a_calibration(tmp_path, monkeypatch, capsys):
    source = tmp_path / "raw.csv"
    make(source, counts_per_second=92889.0, overhead=40.61e-6)
    cal_path = tmp_path / "bench.json"
    calibrate.save_calibration(cal_path,
                               calibrate.calibration(88546.0, 40.61e-6))
    out = tmp_path / "out.csv"
    result = run(monkeypatch, capsys,
                 [str(source), "-o", str(out), "--calibration", str(cal_path)])
    assert "gain factor 1.04" in result.out


def test_cli_fails_before_reading_when_the_output_is_unwritable(tmp_path):
    source = tmp_path / "raw.csv"
    make(source)
    import sys as _sys
    argv = [str(source), "-o", str(tmp_path / "nope" / "out.csv")]
    _sys.argv = ["calibrate.py"] + argv
    with pytest.raises(SystemExit, match="not a directory"):
        calibrate.main()


def test_cli_reports_a_bad_calibration_rather_than_writing_output(tmp_path):
    source = tmp_path / "raw.csv"
    make(source)
    bad = tmp_path / "bad.json"
    bad.write_text("{}")
    out = tmp_path / "out.csv"
    import sys as _sys
    _sys.argv = ["calibrate.py", str(source), "-o", str(out),
                 "--calibration", str(bad)]
    with pytest.raises(SystemExit, match="rm3100_calibration"):
        calibrate.main()
    assert not out.exists()


# --------------------------------------------------------------------------
# Measuring the oscillator's drift
# --------------------------------------------------------------------------

DRIFT_DT = 1.0 / 250.0
# 102 s: ten 10 s windows with room to spare. Landing exactly on a multiple of
# the window would make the window count depend on whether the planted drift
# happens to push the duration over it.
DRIFT_ROWS = 25_500


def drifting(path, ppm_per_second=0.0, jitter_s=0.0, rows=DRIFT_ROWS,
             dt=DRIFT_DT, seed=0):
    """A capture whose sample period ramps, with optional host-side jitter."""
    return capture.load(write_capture(
        path, rows=rows, dt=dt, counts=field_counts(rows, FIELD, 100, 5.0, seed),
        times=drifting_times(rows, dt, ppm_per_second, jitter_s, seed)))


def test_window_rates_recovers_a_planted_ramp(tmp_path):
    """The rate falls at the planted rate, in fractional terms per second."""
    ppm_per_second = 8.0
    cap = drifting(tmp_path / "d.csv", ppm_per_second=ppm_per_second)
    centres, rates = calibrate.window_rates(cap, 10.0)
    assert len(rates) == 10
    # A period ramping up as (1 + r*t) is a rate falling as (1 - r*t).
    slope = np.polyfit(centres, rates / rates.mean(), 1)[0]
    assert slope == pytest.approx(-ppm_per_second * 1e-6, rel=0.05)


def test_window_rates_is_flat_without_drift(tmp_path):
    cap = drifting(tmp_path / "f.csv")
    _, rates = calibrate.window_rates(cap, 10.0)
    assert np.ptp(rates) / rates.mean() < 1e-9


def test_window_rates_centres_span_the_capture(tmp_path):
    cap = drifting(tmp_path / "c.csv")
    centres, rates = calibrate.window_rates(cap, 10.0)
    assert len(centres) == len(rates)
    assert centres[0] == pytest.approx(5.0)
    assert centres[-1] == pytest.approx(95.0)


@pytest.mark.parametrize("window", [0.0, -1.0, 60.0, 1e6])
def test_window_rates_declines_impossible_windows(tmp_path, window):
    """Non-positive, or too long to fit two windows in the capture."""
    cap = drifting(tmp_path / "s.csv")
    centres, rates = calibrate.window_rates(cap, window)
    assert len(centres) == 0 and len(rates) == 0


def test_rate_stability_recovers_a_planted_ramp(tmp_path):
    """Peak to peak across eight windows, against the ramp that produced it."""
    cap = drifting(tmp_path / "r.csv", ppm_per_second=8.0)
    # Eight windows of 12.5 s: centres 6.25 s and 93.75 s apart, so the spread
    # is the ramp over 87.5 s, not over the full 100 s.
    expected = 8.0e-6 * cap.duration * (1 - 1 / 8)
    assert calibrate.rate_stability(cap) == pytest.approx(expected, rel=0.05)


def test_rate_stability_is_zero_on_a_uniform_grid(tmp_path):
    assert calibrate.rate_stability(drifting(tmp_path / "u.csv")) < 1e-9


def test_rate_stability_declines_a_capture_too_short_to_window(tmp_path):
    cap = drifting(tmp_path / "t.csv", rows=200)
    assert calibrate.rate_stability(cap) == 0.0


def test_rate_allan_of_a_linear_drift_rises_with_tau(tmp_path):
    """sigma_y(tau) = |D| * tau / sqrt(2) for a deterministic frequency ramp.

    This is the shape the real captures show above a few seconds, and it is what
    separates a warming oscillator from a random walk.
    """
    ppm_per_second = 8.0
    cap = drifting(tmp_path / "a.csv", ppm_per_second=ppm_per_second)
    taus, devs = calibrate.rate_allan(cap, [4.0, 8.0, 16.0])
    assert list(taus) == [4.0, 8.0, 16.0]
    for tau, dev in zip(taus, devs):
        assert dev == pytest.approx(ppm_per_second * 1e-6 * tau / math.sqrt(2),
                                    rel=0.05)


def test_rate_allan_falls_with_tau_when_only_jitter_is_present(tmp_path):
    """Independent read jitter averages down; drift does not. Opposite slopes."""
    cap = drifting(tmp_path / "j.csv", jitter_s=2e-3)
    taus, devs = calibrate.rate_allan(cap, [4.0, 16.0])
    assert devs[0] > devs[1] * 2


def test_rate_allan_is_negligible_on_a_uniform_grid(tmp_path):
    cap = drifting(tmp_path / "z.csv")
    _, devs = calibrate.rate_allan(cap, [4.0, 8.0])
    assert np.all(devs < 1e-9)


def test_rate_allan_drops_taus_that_yield_too_few_windows(tmp_path):
    """Two windows give one difference, which is not an estimate of anything."""
    cap = drifting(tmp_path / "d2.csv")
    taus, devs = calibrate.rate_allan(cap, [10.0, 40.0, 200.0])
    assert list(taus) == [10.0]      # 40 s gives two windows, 200 s none
    assert len(devs) == 1


def test_rate_allan_accepts_an_empty_request(tmp_path):
    taus, devs = calibrate.rate_allan(drifting(tmp_path / "e.csv"), [])
    assert len(taus) == 0 and len(devs) == 0