rm3100/calibrate.py

#!/usr/bin/env python3
"""Convert a raw capture to nanotesla, optionally applying a calibration.

    ./.venv/bin/python calibrate.py raw.csv -o field.csv
    ./.venv/bin/python calibrate.py raw.csv -o field.csv --calibration bench.json
    ./.venv/bin/python calibrate.py raw.csv -o field.csv --trim 30

Without `--calibration` this only changes the format: counts become nanotesla
using the `tesla_per_count` the capture already carries, and nothing is
corrected. That is the default deliberately -- see "What is not corrected".

**Why a calibration is per-unit.** The chip times everything from one oscillator,
and `rm3100`'s model of it, `period = 3 x (cycle_count / C + overhead)`, uses a
specified C and an overhead measured on one rig. Neither fits every part. Solving
the model against one unit's two cycle counts gave C = 88,546 with a 40.6 us
overhead at 3.006 V, against the model's 90,000 and 68.7 us -- and C moved to
92,889 when the supply moved to 3.353 V, because the oscillator does.

That matters beyond timing. Gain moves with the oscillator, so a capture taken at
a different supply reads a different field for the same magnet. A calibration
records the timing pair and a reference oscillator frequency, and the correction
is how far this capture's oscillator sits from that reference.

Once the timing pair is right, **one capture is enough** to recover the
oscillator: rearranging the model gives `C = cycle_count / (period/3 - overhead)`,
and a 1 us error in the overhead moves C by only 0.09% at cycle count 100. With
the model's 28 us error it moves by 2.5%, and by a *different* amount at each
cycle count, which is what makes the uncalibrated figures unusable for this.

**Why the exponent is 1.** The measurement interval is clocked entirely by that
oscillator: the per-axis overhead measures 3.60 counts at 3.006 V and 3.54 at
3.353 V while its *duration* moves 6%, so the same number of counts is always
collected over a window whose length goes as 1/f. Gain then follows integration
time, giving `gain ~ (cc + n) / f` -- exponent exactly 1. Measured rates agree:
the supply rate ratio is the same at both cycle counts to 0.04%, which a
fixed-time overhead would not produce. It is favoured rather than proven, at 0.3
sigma against 1.2 for a fixed-time overhead; thermal drift is what limits it.

**Why there is a gain offset.** Table 3-1's `0.3671*cc + 1.5` LSB/uT implies the
gain goes as `cc + 4.086` in count units. Solving the same-supply cycle-count
contrast on one unit gives 0.90 and -0.09 -- near zero, so the overhead counts
cost time without integrating field. Leaving it at the datasheet value accounts
for only half of what the correction can remove.

**What is not corrected, and why nothing is by default.** All of it was measured
on four captures, two per supply, with a sensor that rotated between runs by 2.8
to 23 degrees, and every capture drifts 186 to 1241 ppm through its own thermal
warm-up. So it is applied only when a calibration is named, every factor arrives
with the standard deviation it earned, and the output records what was done.

The output is a different format from a raw capture and `capture.py` will refuse
to read it, which is intended: raw captures stay the one source of truth, and a
derived file that drifts from its source should be detectable rather than
plausible. The header carries the source's SHA-256 for that reason.
"""

import argparse
import csv
import hashlib
import json
import math
import os
import sys
from collections import namedtuple

import numpy as np

import capture
import characterize as ch
import rm3100

# Bumped when either file format changes incompatibly, so a reader can refuse
# rather than misinterpret.
CALIBRATION_VERSION = 2
CALIBRATED_VERSION = 1

# Read in blocks when hashing: a capture can be tens of megabytes and there is
# no reason to hold one in memory twice.
HASH_BLOCK = 1 << 20

# Sub-windows used to measure how much the oscillator moved during a capture.
# Eight over a five-minute run is ~40 s each, long enough to fit a rate against
# and short enough to see a thermal transient.
STABILITY_WINDOWS = 8


class CalibrationError(Exception):
    pass


class Estimate(namedtuple("Estimate", "value sd")):
    """A number and its standard deviation, carried together.

    Every corrected field here rests on a chain of measurements, and a residual
    only means something against the uncertainty of the thing that was supposed
    to remove it. Returning the two together is what makes that comparison
    possible at the point of use rather than in a comment.
    """

    __slots__ = ()

    @property
    def relative(self):
        return self.sd / abs(self.value) if self.value else float("inf")

    def __format__(self, spec):
        return f"{self.value:{spec or '.6g'}} +/- {self.sd:{spec or '.2g'}}"


Calibration = namedtuple("Calibration", [
    "counts_per_second",            # C in the timing model, this unit's own
    "axis_overhead_s",              # the model's fixed per-axis term
    "reference_oscillator_hz",      # C at which the gain is taken as correct
    "reference_cycle_count",        # cycle count at which it is taken as correct
    "gain_exponent",                # gain ~ (cc + offset) / f ** exponent
    "gain_offset_counts",           # the "+ n" -- see gain_factor()
    "counts_per_second_sd",
    "axis_overhead_s_sd",
    "gain_exponent_sd",
    "gain_offset_counts_sd",
    "note",
    "created",
])


def calibration(counts_per_second, axis_overhead_s,
                reference_oscillator_hz=None, reference_cycle_count=100,
                gain_exponent=1.0, gain_offset_counts=0.0,
                counts_per_second_sd=0.0, axis_overhead_s_sd=0.0,
                gain_exponent_sd=0.0, gain_offset_counts_sd=0.0,
                note="", created=""):
    """Build a Calibration, defaulting the reference to the measured value.

    A calibration whose reference is its own oscillator frequency corrects
    nothing on the unit and supply it was taken at, and corrects other captures
    relative to it. That is the useful default: the reference is a choice of
    where "right" is, and the only defensible choice without an absolute field
    standard is the condition the gain was characterised at.

    The uncertainties all default to zero, so a calibration that makes no claim
    about its own accuracy still works -- and every field it corrects then
    reports a standard deviation of zero, which is a visible claim rather than a
    silent one.
    """
    return Calibration(
        counts_per_second=float(counts_per_second),
        axis_overhead_s=float(axis_overhead_s),
        reference_oscillator_hz=float(reference_oscillator_hz
                                      if reference_oscillator_hz is not None
                                      else counts_per_second),
        reference_cycle_count=int(reference_cycle_count),
        gain_exponent=float(gain_exponent),
        gain_offset_counts=float(gain_offset_counts),
        counts_per_second_sd=float(counts_per_second_sd),
        axis_overhead_s_sd=float(axis_overhead_s_sd),
        gain_exponent_sd=float(gain_exponent_sd),
        gain_offset_counts_sd=float(gain_offset_counts_sd),
        note=str(note), created=str(created),
    )


def _finite(value, field):
    """A field that must parse as a real number. Booleans are not numbers here.

    JSON's `true` reaches Python as a bool, which floats to 1.0 without
    complaint -- a silent way for a malformed calibration to become a plausible
    one.
    """
    if isinstance(value, bool) or not isinstance(value, (int, float, str)):
        raise CalibrationError(f"{field}: {value!r} is not a number")
    try:
        number = float(value)
    except ValueError:
        raise CalibrationError(f"{field}: {value!r} is not a number")
    if not math.isfinite(number):
        raise CalibrationError(f"{field}: {value!r} is not finite")
    return number


def _positive(value, field):
    number = _finite(value, field)
    if number <= 0:
        raise CalibrationError(f"{field}: {number} must be positive")
    return number


def _non_negative(value, field):
    number = _finite(value, field)
    if number < 0:
        raise CalibrationError(f"{field}: {number} must not be negative")
    return number


def load_calibration(path):
    """Read a calibration file, refusing anything it cannot fully understand."""
    try:
        with open(path) as handle:
            raw = json.load(handle)
    except json.JSONDecodeError as exc:
        raise CalibrationError(f"{path}: not valid JSON ({exc})")
    if not isinstance(raw, dict):
        raise CalibrationError(f"{path}: expected an object, got "
                               f"{type(raw).__name__}")
    version = raw.get("rm3100_calibration")
    if version != CALIBRATION_VERSION:
        raise CalibrationError(
            f"{path}: rm3100_calibration is {version!r}, this reads "
            f"{CALIBRATION_VERSION}. Version 1 carried a per-cycle-count gain "
            "table; version 2 replaces it with gain_offset_counts, which is "
            "not derivable from it -- re-measure rather than convert.")
    missing = {"counts_per_second", "axis_overhead_s"} - raw.keys()
    if missing:
        raise CalibrationError(f"{path}: missing {', '.join(sorted(missing))}")

    cycle_count = raw.get("reference_cycle_count", 100)
    if isinstance(cycle_count, bool) or not isinstance(cycle_count, (int, str)):
        raise CalibrationError(
            f"{path}: reference_cycle_count {cycle_count!r} is not an integer")
    try:
        cycle_count = int(cycle_count)
    except ValueError:
        raise CalibrationError(
            f"{path}: reference_cycle_count {cycle_count!r} is not an integer")
    if cycle_count < 1:
        raise CalibrationError(
            f"{path}: reference_cycle_count {cycle_count} must be at least 1")

    return Calibration(
        counts_per_second=_positive(raw["counts_per_second"],
                                    "counts_per_second"),
        # The overhead is a duration and may legitimately be zero, but not
        # negative: a measurement cannot take less time than its own counts.
        axis_overhead_s=_non_negative(raw["axis_overhead_s"],
                                      "axis_overhead_s"),
        reference_oscillator_hz=_positive(
            raw.get("reference_oscillator_hz", raw["counts_per_second"]),
            "reference_oscillator_hz"),
        reference_cycle_count=cycle_count,
        gain_exponent=_finite(raw.get("gain_exponent", 1.0), "gain_exponent"),
        gain_offset_counts=_finite(raw.get("gain_offset_counts", 0.0),
                                   "gain_offset_counts"),
        counts_per_second_sd=_non_negative(
            raw.get("counts_per_second_sd", 0.0), "counts_per_second_sd"),
        axis_overhead_s_sd=_non_negative(
            raw.get("axis_overhead_s_sd", 0.0), "axis_overhead_s_sd"),
        gain_exponent_sd=_non_negative(
            raw.get("gain_exponent_sd", 0.0), "gain_exponent_sd"),
        gain_offset_counts_sd=_non_negative(
            raw.get("gain_offset_counts_sd", 0.0), "gain_offset_counts_sd"),
        note=str(raw.get("note", "")),
        created=str(raw.get("created", "")),
    )


def save_calibration(path, cal):
    """Write a calibration file."""
    if cal.axis_overhead_s < 0:
        raise CalibrationError(
            f"axis_overhead_s {cal.axis_overhead_s} is negative")
    body = {"rm3100_calibration": CALIBRATION_VERSION,
            "note": cal.note, "created": cal.created}
    body.update({field: getattr(cal, field) for field in Calibration._fields
                 if field not in ("note", "created")})
    with open(path, "w") as handle:
        json.dump(body, handle, indent=2, sort_keys=False)
        handle.write("\n")
    return path


def window_rates(cap, seconds):
    """Sample rate fitted independently on each consecutive window of a capture.

    Returns (centres, rates) -- seconds from the start of the capture, and Hz --
    with one point per window that held enough samples to fit.

    `Capture.dt_true` is a single straight line through a run whose rate is
    moving. This is the same fit repeated on pieces short enough that the
    movement within any one of them is small, which is what turns the drift from
    an error on the whole-capture number into something measurable in its own
    right.
    """
    seconds = float(seconds)
    start = float(cap.elapsed[0])
    if seconds <= 0 or cap.duration < 2 * seconds:
        return np.array([]), np.array([])
    centres, rates = [], []
    for w in range(int(cap.duration // seconds)):
        low = start + w * seconds
        try:
            rate = cap.restrict(low, low + seconds).true_rate_hz
        except capture.CaptureError:
            continue
        rates.append(rate)
        centres.append((w + 0.5) * seconds)
    return np.array(centres), np.array(rates)


def rate_stability(cap, windows=STABILITY_WINDOWS):
    """Fractional spread of the sample rate across sub-windows of a capture.

    The chip's RC oscillator warms up, and every capture here drifts downward
    through its own run -- 186 ppm on the quietest, 1241 ppm on the worst. A
    single fitted `dt_true` is the average over that transient, so this is the
    honest error bar on it, and the input to every uncertainty downstream.

    It is not the same thing as `Capture.drift_limited`, which asks whether one
    straight line describes the capture at all. This asks how far the rate moved
    while it did. Nor is it `rate_allan`, which asks how much of that movement
    is real: a spread this wide is only meaningful against the noise of the
    measurement that produced it.
    """
    if len(cap.sample_index) < 64 * windows:
        return 0.0
    _, rates = window_rates(cap, cap.duration / windows)
    if len(rates) < 2:
        return 0.0
    return float((rates.max() - rates.min()) / rates.mean())


def rate_allan(cap, taus):
    """Two-sample (Allan) deviation of the measured sample rate, fractional.

    For each averaging time tau the rate is fitted on consecutive tau-second
    windows and the deviation taken between neighbours:

        sigma_y(tau) = sqrt( mean( (y[k+1] - y[k])^2 ) / 2 )

    Differences between neighbours rather than a spread about the mean, for the
    same reason `characterize.white_sd` differences rather than taking a
    standard deviation: on a series that drifts, the spread measures the drift.
    Here both effects are wanted and the shape of the curve separates them.
    Falling with tau is measurement noise averaging down; rising with tau is
    drift the measurement has begun to resolve, and a rise proportional to tau
    specifically is a deterministic frequency ramp rather than a random walk.

    The minimum is therefore the best the rate can be known, and the tau where
    it falls is how long it is worth measuring for.

    Returns (taus, sigma_y) over the taus that yielded three or more windows.
    Two windows give a single difference, which is not an estimate of anything.
    """
    kept, devs = [], []
    for tau in taus:
        _, rates = window_rates(cap, tau)
        if len(rates) < 3:
            continue
        y = rates / rates.mean()
        kept.append(float(tau))
        devs.append(float(np.sqrt(np.mean(np.diff(y) ** 2) / 2)))
    return np.array(kept), np.array(devs)


TimingFit = namedtuple("TimingFit", "counts_per_second axis_overhead_s "
                                    "counts_per_second_sd axis_overhead_s_sd "
                                    "residual_ppm points")


def fit_timing(cycle_counts, periods, period_sd=None):
    """Solve `period = AXES * (cycle_count / C + overhead)` for both terms.

    `periods` are whole sample periods in seconds, not per-axis. `period_sd` is
    the fractional uncertainty on each -- `rate_stability()` is where it comes
    from.

    Both terms are needed. Holding C at the specified 90,000 and solving for the
    overhead alone gave 58.9 us at cycle count 100 and 113.6 us at 400 on the
    same unit -- no single value fits, because the misfit is in both.

    **Two points determine both terms exactly and leave no residual**, so a bad
    point cannot be detected: `residual_ppm` is 0 and means nothing. That case
    requires `period_sd`, so the uncertainty is at least stated rather than
    fabricated from a perfect fit. Three or more points carry their own.
    """
    counts = np.asarray(cycle_counts, dtype=float)
    per_axis = np.asarray(periods, dtype=float) / rm3100.AXES
    if counts.shape != per_axis.shape:
        raise CalibrationError(
            f"{counts.size} cycle counts against {per_axis.size} periods")
    distinct = len(np.unique(counts))
    if distinct < 2:
        raise CalibrationError(
            "need at least two distinct cycle counts to fit both the count "
            "rate and the per-axis overhead")

    if distinct == 2 and len(counts) == 2:
        if period_sd is None:
            raise CalibrationError(
                "a two-point fit has no residual and cannot detect a bad "
                "point, so it cannot estimate its own uncertainty. Pass "
                "period_sd (from rate_stability) or sweep a third cycle count.")
        slope, intercept = np.polyfit(counts, per_axis, 1)
        if slope <= 0:
            raise CalibrationError(
                "fitted count rate is not positive; these captures do not "
                "follow the timing model")
        count_rate = 1.0 / float(slope)
        # C leans on the *difference* of the two per-axis times, so a fractional
        # error on each is amplified by how close together they are.
        span = abs(per_axis[1] - per_axis[0])
        rate_sd = count_rate * period_sd * float(np.hypot(*per_axis)) / span
        overhead_sd = (per_axis.min() * period_sd
                       + counts.min() / count_rate * rate_sd / count_rate)
        return TimingFit(count_rate, float(intercept), rate_sd, overhead_sd,
                         0.0, len(counts))

    (slope, intercept), cov = np.polyfit(counts, per_axis, 1, cov=True)
    if slope <= 0:
        raise CalibrationError(
            "fitted count rate is not positive; these captures do not follow "
            "the timing model")
    count_rate = 1.0 / float(slope)
    slope_sd, overhead_sd = math.sqrt(cov[0, 0]), math.sqrt(cov[1, 1])
    residual = per_axis - (slope * counts + intercept)
    return TimingFit(
        count_rate, float(intercept),
        # d(1/slope) = dslope / slope^2
        float(slope_sd) / slope ** 2, float(overhead_sd),
        float(np.abs(residual / per_axis).max() * 1e6), len(counts))


def oscillator_hz(cap, cal, period_sd=None):
    """This capture's oscillator frequency, with its uncertainty.

    Inverts `period = AXES * (cycle_count / C + overhead)`. The overhead comes
    from the calibration because it cannot be recovered from a single cycle
    count -- one equation, two unknowns -- and it is the term the nominal model
    gets most wrong.

    `period_sd` is the fractional uncertainty on the measured period; it
    defaults to measuring it from the capture with `rate_stability()`.
    """
    per_axis = cap.dt_true / rm3100.AXES - cal.axis_overhead_s
    if per_axis <= 0:
        raise CalibrationError(
            f"{cap.path}: a {cap.dt_true * 1e3:.3f} ms period leaves no time "
            f"for {cap.cycle_count} counts once the calibration's "
            f"{cal.axis_overhead_s * 1e6:.1f} us overhead is removed. The "
            "calibration does not describe this capture.")
    if period_sd is None:
        period_sd = rate_stability(cap)
    value = cap.cycle_count / per_axis
    # Two independent contributions: how well the period was measured, and how
    # well the overhead subtracted from it is known. The second dominates at low
    # cycle counts, where the overhead is a larger share of the per-axis time.
    from_period = cap.dt_true * period_sd / rm3100.AXES
    relative = math.hypot(from_period, cal.axis_overhead_s_sd) / per_axis
    return Estimate(value, value * relative)


def gain_factor(cap, cal, period_sd=None):
    """Multiply measured nanotesla by this to put them on the reference gain.

    Three terms, each correcting something the capture's own conversion assumed:

      oscillator   how far this capture's clock sits from the reference, raised
                   to the fitted exponent. The measurement interval is clocked
                   entirely by that oscillator -- the per-axis overhead measures
                   3.57 counts at both supplies, constant, while its *duration*
                   moves 6% -- so the same number of counts is always collected
                   over a window whose length goes as 1/f. That predicts an
                   exponent of exactly 1.
      gain shape   Table 3-1's `0.3671*cc + 1.5` LSB/uT implies the gain goes as
                   `cc + 4.086` in count units. Measured on one unit it is 0 to
                   0.9, so the overhead counts cost time without integrating
                   field. `gain_offset_counts` is that "+ n".
      reference    the shape term is only known up to an overall scale, since
                   `gain = A(cc + n)` leaves A free. `reference_cycle_count`
                   pins it: the correction is exactly 1 there.
    """
    oscillator = oscillator_hz(cap, cal, period_sd)
    ratio = oscillator.value / cal.reference_oscillator_hz
    reference_cc, offset = cal.reference_cycle_count, cal.gain_offset_counts
    if cap.cycle_count + offset <= 0 or reference_cc + offset <= 0:
        raise CalibrationError(
            f"gain_offset_counts {offset} puts the gain at or below zero for "
            f"cycle count {min(cap.cycle_count, reference_cc)}")
    shape = ((rm3100.gain_lsb_per_tesla(cap.cycle_count)
              / rm3100.gain_lsb_per_tesla(reference_cc))
             * ((reference_cc + offset) / (cap.cycle_count + offset)))
    value = ratio ** cal.gain_exponent * shape

    # In logs, so the terms add in quadrature as independent contributions.
    from_oscillator = cal.gain_exponent * oscillator.relative
    from_exponent = cal.gain_exponent_sd * abs(math.log(ratio)) if ratio > 0 else 0.0
    from_offset = cal.gain_offset_counts_sd * abs(
        1.0 / (reference_cc + offset) - 1.0 / (cap.cycle_count + offset))
    return Estimate(value, value * math.sqrt(from_oscillator ** 2
                                             + from_exponent ** 2
                                             + from_offset ** 2))


def source_digest(path):
    """SHA-256 of a file, so a derived copy can be tied back to its source."""
    digest = hashlib.sha256()
    with open(path, "rb") as handle:
        for block in iter(lambda: handle.read(HASH_BLOCK), b""):
            digest.update(block)
    return digest.hexdigest()


def write_calibrated(cap, path, source, cal=None, factor=None,
                     calibration_path=None):
    """Write nanotesla and a header recording where they came from.

    `factor` is an Estimate; its standard deviation goes into the header beside
    it, so a reader can tell a real residual from calibration slop without
    having the calibration to hand.
    """
    if factor is None:
        factor = Estimate(1.0, 0.0)
    elif not isinstance(factor, Estimate):
        factor = Estimate(float(factor), 0.0)
    header = [
        ("rm3100_calibrated", CALIBRATED_VERSION),
        ("source", os.path.basename(source)),
        ("source_sha256", source_digest(source)),
        ("samples", len(cap.sample_index)),
        ("cycle_count", cap.cycle_count),
        ("tesla_per_count", repr(cap.tesla_per_count)),
        ("measured_rate_hz", repr(cap.true_rate_hz)),
        ("calibration", os.path.basename(calibration_path)
                        if calibration_path else "none"),
        ("gain_factor", repr(factor.value)),
        ("gain_factor_sd", repr(factor.sd)),
    ]
    if cal is not None:
        oscillator = oscillator_hz(cap, cal)
        header += [
            ("oscillator_hz", repr(oscillator.value)),
            ("oscillator_hz_sd", repr(oscillator.sd)),
            ("reference_oscillator_hz", repr(cal.reference_oscillator_hz)),
            ("reference_cycle_count", cal.reference_cycle_count),
            ("gain_exponent", repr(cal.gain_exponent)),
            ("gain_offset_counts", repr(cal.gain_offset_counts)),
        ]

    x, y, z = (cap.x * factor.value, cap.y * factor.value,
               cap.z * factor.value)
    total = np.sqrt(x**2 + y**2 + z**2)
    # The two flags are independent, so a row can carry both. Placeholder rows
    # keep their interpolated values here, exactly as capture.py produced them,
    # and the flag is what says the values were reconstructed.
    flags = [" ".join(
                 ([capture.WARN_MISSED] if missed else [])
                 + ([capture.WARN_AMBIGUOUS] if ambiguous else []))
             for missed, ambiguous in zip(cap.missed, cap.ambiguous)]

    with open(path, "w", newline="") as handle:
        for key, value in header:
            handle.write(f"# {key}: {value}\n")
        writer = csv.writer(handle)
        writer.writerow(["sample_index", "elapsed_s",
                         "x_nT", "y_nT", "z_nT", "total_nT", "warning"])
        writer.writerows(
            (index, f"{t:.6f}", f"{a:.3f}", f"{b:.3f}", f"{c:.3f}",
             f"{n:.3f}", flag)
            for index, t, a, b, c, n, flag
            in zip(cap.sample_index, cap.elapsed, x, y, z, total, flags))
    return path


class Calibrated:
    """A capture already converted to nanotesla, read back from disk.

    Deliberately not a `capture.Capture`: the raw counts and the host clock are
    gone, so the time base cannot be re-fitted and `restrict()` would be a lie.
    What survives is what a plot or a spectrum needs, plus the provenance of how
    it got that way.
    """

    def __init__(self, path, meta, sample_index, elapsed, x, y, z, total):
        self.path = path
        self.meta = meta
        self.sample_index = sample_index
        self.elapsed = elapsed
        self.x, self.y, self.z, self.total = x, y, z, total
        self.cycle_count = int(meta["cycle_count"])
        self.tesla_per_count = float(meta["tesla_per_count"])
        self.lsb_nt = self.tesla_per_count * rm3100.NT_PER_TESLA
        self.true_rate_hz = float(meta["measured_rate_hz"])
        self.dt_true = 1.0 / self.true_rate_hz
        self.gain_factor = float(meta["gain_factor"])
        self.gain_factor_sd = float(meta.get("gain_factor_sd", 0.0))
        self.calibration_name = meta.get("calibration", "none")
        self.source = meta.get("source", "")
        self.source_sha256 = meta.get("source_sha256", "")

    @property
    def duration(self):
        return float(self.elapsed[-1] - self.elapsed[0])

    def axes(self):
        return [("x", self.x), ("y", self.y), ("z", self.z),
                ("total", self.total)]

    def describe(self):
        if self.gain_factor == 1.0:
            return "converted to nT, uncorrected"
        return (f"gain x{self.gain_factor:.6f} +/- {self.gain_factor_sd:.6f} "
                f"from {self.calibration_name}")


def load_calibrated(path):
    """Read a file written by write_calibrated()."""
    meta, rows = {}, []
    with open(path, newline="") as handle:
        while True:
            position = handle.tell()
            line = handle.readline()
            if not line:
                break
            if not line.startswith("#"):
                handle.seek(position)
                break
            key, _, value = line[1:].partition(":")
            if value:
                meta[key.strip()] = value.strip()
        rows = list(csv.DictReader(handle))

    version = meta.get("rm3100_calibrated")
    if version is None:
        raise CalibrationError(
            f"{path}: no rm3100_calibrated header. A raw capture goes through "
            "capture.load() instead.")
    if int(version) != CALIBRATED_VERSION:
        raise CalibrationError(
            f"{path}: rm3100_calibrated is {version}, this reads "
            f"{CALIBRATED_VERSION}")
    missing = {"cycle_count", "tesla_per_count", "measured_rate_hz",
               "gain_factor"} - meta.keys()
    if missing:
        raise CalibrationError(f"{path}: header missing "
                               f"{', '.join(sorted(missing))}")
    if not rows:
        raise CalibrationError(f"{path}: no samples")

    column = lambda name: np.array([float(r[name]) for r in rows])
    return Calibrated(
        path, meta,
        np.array([int(r["sample_index"]) for r in rows], dtype=np.int64),
        column("elapsed_s"), column("x_nT"), column("y_nT"), column("z_nT"),
        column("total_nT"))


def load_any(path):
    """Load a raw capture or a calibrated one, whichever this is.

    Which format a file is, is a property of the file, not of the caller -- so
    tools take a path and get back something with axes() and a time base either
    way, rather than each having to know.
    """
    with open(path) as handle:
        head = handle.read(4096)
    if "rm3100_calibrated" in head:
        return load_calibrated(path)
    return capture.load(path)


def main():
    ap = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("csv", help="raw capture written by logger.py")
    ap.add_argument("-o", "--output", required=True, metavar="CSV",
                    help="where to write the nanotesla capture")
    ap.add_argument("--calibration", default=None, metavar="JSON",
                    help="apply this calibration; without it nothing is "
                         "corrected and only the units change")
    ap.add_argument("--trim", type=float, default=0.0, metavar="SECONDS",
                    help="drop this many seconds from BOTH ends before "
                         "converting (default: %(default)s)")
    args = ap.parse_args()

    # Check the destination before reading a capture that may be tens of
    # megabytes, so a mistyped directory fails in a second rather than a minute.
    parent = os.path.dirname(os.path.abspath(args.output)) or "."
    if not os.path.isdir(parent):
        sys.exit(f"{args.output}: {parent} is not a directory")
    if not os.access(parent, os.W_OK):
        sys.exit(f"{args.output}: {parent} is not writable")

    try:
        cal = load_calibration(args.calibration) if args.calibration else None
        cap = capture.load(args.csv)
        cap, note = ch.trimmed(cap, args.trim)
        factor = gain_factor(cap, cal) if cal else Estimate(1.0, 0.0)
    except (OSError, ValueError, CalibrationError,
            capture.CaptureError) as exc:
        sys.exit(str(exc))

    if note:
        print(f"  {note}")
    if cal is None:
        print("no calibration given: converting units only, correcting nothing")
    else:
        measured = oscillator_hz(cap, cal)
        stability = rate_stability(cap)
        print(f"calibration {args.calibration}"
              + (f" -- {cal.note}" if cal.note else ""))
        print(f"  rate        {cap.true_rate_hz:.4f} Hz, drifting "
              f"{stability * 1e6:.0f} ppm across the run")
        print(f"  oscillator  {measured.value:,.0f} +/- {measured.sd:,.0f} Hz "
              f"against a {cal.reference_oscillator_hz:,.0f} Hz reference "
              f"({measured.value / cal.reference_oscillator_hz - 1:+.2%})")
        if cap.cycle_count != cal.reference_cycle_count:
            print(f"  gain shape  cycle count {cap.cycle_count} against a "
                  f"{cal.reference_cycle_count} reference, offset "
                  f"{cal.gain_offset_counts:+.2f} counts")
        print(f"  gain factor {factor.value:.6f} +/- {factor.sd:.6f} "
              f"({factor.relative:.2%})   |B| {cap.total.mean():,.0f} -> "
              f"{cap.total.mean() * factor.value:,.0f} nT")

    try:
        out = write_calibrated(cap, args.output, args.csv, cal, factor,
                               args.calibration)
    except (OSError, CalibrationError) as exc:
        sys.exit(str(exc))
    print(f"-> {out}")
    return 0


if __name__ == "__main__":
    sys.exit(main())