From abbda7a05e6f91f93b3761f454b26afef9682899 Mon Sep 17 00:00:00 2001 From: Josh Poole Date: Tue, 15 Sep 2026 08:37:16 +0100 Subject: [PATCH] 20260915 - Fit R and Q out of the innovations, along the prediction interval The records exist; nothing reads them. This reads them, and it exists now rather than after the capture because a day's analysis written while the data is sitting there is a day the calibration is not landing. R and Q cancel, so neither can be fitted alone and NIS moves the wrong way for both. What separates them is that only one grows with how long the filter had been predicting: var(nu | L) = R + a * L^p, with R the intercept. Fixing p makes the rest linear in R and a, so a grid over p reduces the whole thing to a 2x2 solve per candidate and it stays stdlib-only, which is what lets it be copied onto a node like live_score. p is the finding, not just a nuisance parameter. Live data said the model error grows as L^1.96 where white jerk, which is what Q implements, would be L^5. An exponent near 2 says no jerk value fits and the adaptive scale cannot converge, and that is a different repair from raising range_jerk. Scatter is also reported per axis against the S the filter claimed. NIS sums the two axes into one number, which is precisely how an oversized R and an undersized Q have been hiding in it. The innovation record now carries the occupancy map's verdict per detection, so the interference that is three quarters of the tracks at both nodes is dropped on evidence rather than on a track's Doppler holding still. The old heuristic stays for records written before the map landed, and it is worth less than it looks: it costs the genuine constant-Doppler aircraft too. Verified by synthesising records from a known R and a known velocity error and asserting both come back. From the plan's own measured figures it recovers 0.141 and 0.181 cells against 0.15 and 0.18, and an exponent of 2.05 against 1.96. Co-Authored-By: Claude Opus 5 (1M context) --- retina_tracker/calibrate.py | 333 ++++++++++++++++++++++++++++++++++++ retina_tracker/output.py | 5 + tests/test_calibrate.py | 211 +++++++++++++++++++++++ tests/test_innovations.py | 1 + 4 files changed, 550 insertions(+) create mode 100644 retina_tracker/calibrate.py create mode 100644 tests/test_calibrate.py diff --git a/retina_tracker/calibrate.py b/retina_tracker/calibrate.py new file mode 100644 index 0000000..bb25051 --- /dev/null +++ b/retina_tracker/calibrate.py @@ -0,0 +1,333 @@ +"""Fit R and Q from recorded innovations, separating two constants that cancel. + +R and Q are both wrong today, in opposite directions. R is around 45x too large +in delay, while the model error it is standing in for grows as the square of the +prediction interval rather than the fifth power white jerk implies. The sum comes +out plausible, which is why nothing has complained, and it is why neither can be +fitted alone: a single scalar like NIS moves one way for an oversized R and the +other for an undersized Q. + +What separates them is the prediction interval. The innovation covariance is +HPH' + R, so scatter measured against how long the filter had been predicting is +a constant plus a growing term: + + var(nu | L) = R + a * L^p + +R is the intercept, and (a, p) are the shape and size of the model error. A grid +over p makes the rest a two-parameter linear fit, which is all this needs. + +Deliberately stdlib-only, like live_score: it is meant to be copied onto a node +and run under whatever python3 is there. + + python -m retina_tracker.calibrate innovations.jsonl --blah2-config config.yml + python -m retina_tracker.calibrate jn1.jsonl ffb.jsonl --fc 177000000 +""" + +import argparse +import json +import math +import statistics +import sys +from collections import defaultdict + +SPEED_OF_LIGHT = 299792458.0 + +# Median of chi-squared with two degrees of freedom. A matched filter's NIS +# should sit here; the plan measured 1.083 live, which is the signature of the +# oversized R holding the whole thing down. +CHI2_2DOF_MEDIAN = 2.0 * math.log(2.0) + +# 1.4826 * MAD estimates sigma for a Gaussian, and unlike a standard deviation +# it does not let one mis-association set the answer. Calibration data comes +# from a live node, so it has mis-associations in it. +MAD_TO_SIGMA = 1.4826 + +# Prediction intervals are quantised by the frame rate anyway; rounding groups +# the ones that differ only by clock jitter. +INTERVAL_QUANTUM_S = 0.05 + +MIN_PER_INTERVAL = 30 +MIN_INTERVALS = 3 + +# A track whose Doppler never leaves one bin is the fallback test for +# interference, used only on records written before the map carried its verdict. +FIXED_DOPPLER_BINS = 1.0 +MIN_RECORDS_TO_JUDGE_A_TRACK = 8 + + +def wavelength_km(fc_hz): + return SPEED_OF_LIGHT / fc_hz / 1000.0 + + +def delay_cell_km(fs_hz): + return SPEED_OF_LIGHT / fs_hz / 1000.0 + + +def load_blah2_capture(path): + """Read fc, fs and cpi out of a blah2 config without a YAML dependency.""" + found = {} + section = [] + with open(path) as handle: + for line in handle: + if not line.strip() or line.lstrip().startswith("#"): + continue + indent = len(line) - len(line.lstrip()) + stripped = line.strip() + depth = indent // 2 + section = section[:depth] + if stripped.endswith(":"): + section.append(stripped[:-1]) + continue + key, _, value = stripped.partition(":") + value = value.strip() + if section[:1] == ["capture"] and key in ("fc", "fs"): + found[key] = float(value) + elif section[:2] == ["process", "data"] and key == "cpi": + found["cpi"] = float(value) + return found + + +def load_records(path): + records = [] + with open(path) as handle: + for line in handle: + line = line.strip() + if line: + records.append(json.loads(line)) + return records + + +def interfering_tracks(records): + """Tracks whose Doppler never left one bin, for records with no verdict. + + Only a fallback. It costs the genuine constant-Doppler aircraft along with + the tone, which biases what is left toward manoeuvring targets, so it is + worth something for R and worth less for Q. Records written by a tracker + carrying the occupancy map say so per detection instead. + """ + by_track = defaultdict(list) + for record in records: + doppler = record.get("doppler") + if doppler is not None: + by_track[record.get("birth")].append(doppler) + return { + birth + for birth, dopplers in by_track.items() + if len(dopplers) >= MIN_RECORDS_TO_JUDGE_A_TRACK and max(dopplers) - min(dopplers) < FIXED_DOPPLER_BINS + } + + +def usable(records): + """Records that describe a real measurement of a real target. + + Drops the interference, which is most of the record at an interfered node, + and the degenerate updates, whose innovation is NaN because a singular + covariance meant the measurement was skipped. + """ + fallback = interfering_tracks(records) if not any("interfering" in r for r in records) else set() + kept = [] + for record in records: + if record.get("interfering") or record.get("birth") in fallback: + continue + innovation = record.get("innovation") or [] + if len(innovation) != 2 or any(x is None or not math.isfinite(x) for x in innovation): + continue + if not record.get("dt"): + continue + kept.append(record) + return kept + + +def prediction_interval(record): + """How long the filter had been predicting when this innovation was taken. + + dt is the interval of the last predict and n_missed the frames coasted + before this association landed, so this assumes a steady frame rate across + a coast. It is the quantity R and Q are separated along, and getting it + from the record is the whole reason the tracker emits n_missed. + """ + return record["dt"] * (1 + record.get("n_missed", 0)) + + +def mad_sigma(values): + centre = statistics.median(values) + return MAD_TO_SIGMA * statistics.median([abs(v - centre) for v in values]) + + +def group_by_interval(records, component): + """Scatter of one innovation component against prediction interval.""" + buckets = defaultdict(list) + for record in records: + interval = round(prediction_interval(record) / INTERVAL_QUANTUM_S) * INTERVAL_QUANTUM_S + buckets[round(interval, 3)].append(record["innovation"][component]) + return sorted( + (interval, len(values), mad_sigma(values) ** 2) + for interval, values in buckets.items() + if interval > 0 and len(values) >= MIN_PER_INTERVAL + ) + + +def fit_noise_model(points, exponents=None): + """Least squares for var = R + a * L^p, weighted by sample count. + + Linear in R and a once p is fixed, so a grid over p reduces the whole thing + to a 2x2 solve per candidate. p is what says whether Q has the right shape: + white jerk, which is what the filter implements, would be 5. + """ + if len(points) < MIN_INTERVALS: + return None + if exponents is None: + exponents = [0.5 + 0.05 * i for i in range(91)] + + best = None + for p in exponents: + s_w = s_x = s_xx = s_y = s_xy = 0.0 + for interval, count, variance in points: + x = interval**p + w = float(count) + s_w += w + s_x += w * x + s_xx += w * x * x + s_y += w * variance + s_xy += w * x * variance + determinant = s_w * s_xx - s_x * s_x + if abs(determinant) < 1e-30: + continue + r = (s_y * s_xx - s_x * s_xy) / determinant + a = (s_w * s_xy - s_x * s_y) / determinant + if r < 0 or a < 0: + continue + cost = sum(count * (variance - (r + a * interval**p)) ** 2 for interval, count, variance in points) + if best is None or cost < best[0]: + best = (cost, r, a, p) + + if best is None: + return None + _cost, r, a, p = best + return {"R": r, "a": a, "exponent": p} + + +def normalised_scatter(records, component): + """Actual scatter over what the filter claimed, per axis. + + NIS collapses both axes into one number, which is how an oversized R and an + undersized Q have been hiding in it. Above one means the filter is more + confident than it has earned on that axis; below one means the reverse. + """ + ratios = [ + record["innovation"][component] / math.sqrt(record["s_diag"][component]) + for record in records + if record.get("s_diag") and record["s_diag"][component] > 0 + ] + return mad_sigma(ratios) if len(ratios) >= MIN_PER_INTERVAL else None + + +def calibrate(records, fc_hz, fs_hz, cpi_s): + kept = usable(records) + report = { + "n_records": len(records), + "n_used": len(kept), + "n_interfering": sum(1 for r in records if r.get("interfering")), + "n_tracks": len({r.get("birth") for r in kept}), + "delay_cell_km": delay_cell_km(fs_hz), + "doppler_bin_hz": 1.0 / cpi_s, + "wavelength_km": wavelength_km(fc_hz), + "nis_median": statistics.median([r["nis"] for r in kept if r.get("nis") is not None]) if kept else None, + "axes": {}, + } + if not kept: + return report + + for name, component in (("delay", 0), ("doppler", 1)): + fit = fit_noise_model(group_by_interval(kept, component)) + axis = {"scatter_vs_claimed": normalised_scatter(kept, component)} + if fit: + sigma = math.sqrt(fit["R"]) + if name == "delay": + axis["sigma"] = sigma + axis["unit"] = "km" + axis["cells"] = sigma / report["delay_cell_km"] + else: + # The innovation is range rate, which is what R and S are in. + sigma_hz = sigma / report["wavelength_km"] + axis["sigma"] = sigma_hz + axis["unit"] = "Hz" + axis["cells"] = sigma_hz / report["doppler_bin_hz"] + axis["exponent"] = fit["exponent"] + # At p = 2 the growing term is a velocity error held for the whole + # interval, which is the shape the live data showed. + axis["error_rate"] = math.sqrt(fit["a"]) + report["axes"][name] = axis + + return report + + +def format_report(name, report): + lines = [f"{name}"] + lines.append( + f" records {report['n_used']} used of {report['n_records']} " + f"({report['n_interfering']} interfering), {report['n_tracks']} tracks" + ) + lines.append(f" cell {report['delay_cell_km']:.4f} km x {report['doppler_bin_hz']:.2f} Hz") + if report["nis_median"] is not None: + lines.append(f" NIS median {report['nis_median']:.3f}, matched {CHI2_2DOF_MEDIAN:.3f}") + + for axis_name in ("delay", "doppler"): + axis = report["axes"].get(axis_name) + if not axis: + lines.append(f" {axis_name:<8}not enough spread in prediction interval to separate R from Q") + continue + if "sigma" in axis: + lines.append( + f" {axis_name:<8}R = {axis['sigma']:.4f} {axis['unit']} = {axis['cells']:.3f} cells, " + f"model error grows as L^{axis['exponent']:.2f} at {axis['error_rate']:.4g}" + ) + if axis["scatter_vs_claimed"] is not None: + lines.append(f" {'':<8}scatter is {axis['scatter_vs_claimed']:.2f}x what the filter claimed") + return "\n".join(lines) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + parser.add_argument("files", nargs="+", help="innovations.jsonl written by --innovations") + parser.add_argument("--blah2-config", help="blah2 config.yml to read fc, fs and cpi from") + parser.add_argument("--fc", type=float, help="Centre frequency in Hz") + parser.add_argument("--fs", type=float, default=2000000.0, help="Sample rate in Hz (default: 2e6)") + parser.add_argument("--cpi", type=float, default=0.5, help="Coherent processing interval in s (default: 0.5)") + parser.add_argument("--json", action="store_true", help="Emit the report as JSON") + args = parser.parse_args(argv) + + capture = load_blah2_capture(args.blah2_config) if args.blah2_config else {} + fc = args.fc or capture.get("fc") + if not fc: + parser.error("centre frequency unknown: pass --fc or --blah2-config") + fs = capture.get("fs", args.fs) + cpi = capture.get("cpi", args.cpi) + + reports = {} + for path in args.files: + reports[path] = calibrate(load_records(path), fc, fs, cpi) + + if args.json: + print(json.dumps(reports, indent=2)) + return 0 + + for path, report in reports.items(): + print(format_report(path, report)) + print() + + cells = [ + report["axes"][axis]["cells"] + for report in reports.values() + for axis in ("delay", "doppler") + if report["axes"].get(axis, {}).get("cells") is not None + ] + if len(cells) > 1: + spread = max(cells) / min(cells) + print(f"k across {len(cells)} axis-site pairs: {min(cells):.3f} to {max(cells):.3f} cells, {spread:.2f}x apart") + print("One k is only shippable if that spread is small. Two axes agreeing is time-bandwidth duality holding.") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/retina_tracker/output.py b/retina_tracker/output.py index 662ce33..1e0438c 100644 --- a/retina_tracker/output.py +++ b/retina_tracker/output.py @@ -172,6 +172,11 @@ def write_residual(self, track_id, timestamp, track, detection): # events file, which a track that never confirmed is absent from. "delay": detection.get("delay"), "doppler": detection.get("doppler"), + # The occupancy map's verdict on that measurement, so calibration + # can drop the interference exactly rather than inferring it from a + # track's Doppler holding still. Present whether or not the node is + # suppressing on it. + "interfering": bool(detection.get("interfering")), "n_missed": track.last_n_missed, "q_scale": track.last_q_scale, "innovation": [float(x) for x in residual.innovation], diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py new file mode 100644 index 0000000..a672803 --- /dev/null +++ b/tests/test_calibrate.py @@ -0,0 +1,211 @@ +"""Fitting R and Q back out of recorded innovations. + +The two cancel: R is around 45x too large in delay while the model error it +stands in for grows as the square of the prediction interval, not the fifth +power white jerk implies, so the sum looks plausible and NIS alone cannot tell +which is wrong. What separates them is that only one of the two grows with how +long the filter had been predicting. + +These tests synthesise records from a known R and a known velocity error and +assert the fit recovers both, which is the only way to know the tool is +measuring rather than asserting before any of it is pointed at a node. +""" + +import json +import random + +import pytest + +from retina_tracker.calibrate import ( + calibrate, + fit_noise_model, + group_by_interval, + load_blah2_capture, + prediction_interval, + usable, + wavelength_km, +) +from retina_tracker.output import InnovationWriter +from retina_tracker.tracker import Tracker + +FC = 177000000.0 +FS = 2000000.0 +CPI = 0.5 + +TRUE_DELAY_SIGMA_KM = 0.022 +TRUE_DOPPLER_SIGMA_HZ = 0.356 +TRUE_VELOCITY_ERROR = 0.015 + +INTERVALS = (0.5, 1.0, 1.5, 2.0, 3.0) +PER_INTERVAL = 1200 + + +def record(interval, delay_innovation, rate_innovation, doppler=-120.0, birth=0, **overrides): + record = { + "track_id": "250915-000001", + "birth": birth, + "timestamp": 1718747745000, + "dt": 0.5, + "n_missed": round(interval / 0.5) - 1, + "snr": 16.0, + "delay": 20.0, + "doppler": doppler, + "interfering": False, + "q_scale": 1.0, + "innovation": [delay_innovation, rate_innovation], + "s_diag": [1.0, 5e-6], + "nis": 1.0, + } + record.update(overrides) + return record + + +def synthesise(seed=1, delay_sigma=TRUE_DELAY_SIGMA_KM, velocity_error=TRUE_VELOCITY_ERROR, exponent=2.0): + """Innovations from a filter whose R is `delay_sigma` and whose model error + is a velocity held for the whole prediction interval.""" + rng = random.Random(seed) + rate_sigma = TRUE_DOPPLER_SIGMA_HZ * wavelength_km(FC) + records = [] + for interval in INTERVALS: + growth = velocity_error * interval ** (exponent / 2.0) + delay_spread = (delay_sigma**2 + growth**2) ** 0.5 + rate_spread = (rate_sigma**2 + (growth / 50.0) ** 2) ** 0.5 + for i in range(PER_INTERVAL): + records.append( + record( + interval, + rng.gauss(0.0, delay_spread), + rng.gauss(0.0, rate_spread), + doppler=-120.0 + i % 40, + ) + ) + return records + + +class TestTheFitRecoversWhatWasPutIn: + def test_r_comes_back_from_the_intercept(self): + report = calibrate(synthesise(), FC, FS, CPI) + assert report["axes"]["delay"]["sigma"] == pytest.approx(TRUE_DELAY_SIGMA_KM, rel=0.2) + + def test_r_comes_back_in_resolution_cells(self): + """The number that has to agree across sites, and the one a shipped + constant is expressed in.""" + report = calibrate(synthesise(), FC, FS, CPI) + assert report["axes"]["delay"]["cells"] == pytest.approx(0.147, rel=0.25) + + def test_the_doppler_axis_comes_back_in_hertz(self): + """The innovation is a range rate. Reported in Hz it is comparable with + the configured value and with the other axis in cells.""" + report = calibrate(synthesise(), FC, FS, CPI) + assert report["axes"]["doppler"]["sigma"] == pytest.approx(TRUE_DOPPLER_SIGMA_HZ, rel=0.25) + + def test_the_shape_of_the_model_error_comes_back(self): + """Two, not five. An exponent near five would mean white jerk really is + what the filter is missing and range_jerk merely needs raising.""" + report = calibrate(synthesise(), FC, FS, CPI) + assert report["axes"]["delay"]["exponent"] == pytest.approx(2.0, abs=0.4) + + def test_the_size_of_the_model_error_comes_back(self): + report = calibrate(synthesise(), FC, FS, CPI) + assert report["axes"]["delay"]["error_rate"] == pytest.approx(TRUE_VELOCITY_ERROR, rel=0.3) + + def test_a_larger_r_is_reported_larger(self): + """The oversized R the tracker ships with today, against the measured + one, on data whose model error is identical.""" + loose = calibrate(synthesise(seed=2, delay_sigma=0.2), FC, FS, CPI) + tight = calibrate(synthesise(seed=2), FC, FS, CPI) + assert loose["axes"]["delay"]["sigma"] > 5 * tight["axes"]["delay"]["sigma"] + + def test_white_jerk_is_told_apart_from_a_velocity_error(self): + report = calibrate(synthesise(seed=3, exponent=5.0), FC, FS, CPI) + assert report["axes"]["delay"]["exponent"] > 3.0 + + +class TestSeparatingTheTwo: + def test_a_constant_scatter_reads_as_all_r_and_no_q(self): + report = calibrate(synthesise(seed=4, velocity_error=0.0), FC, FS, CPI) + assert report["axes"]["delay"]["sigma"] == pytest.approx(TRUE_DELAY_SIGMA_KM, rel=0.15) + assert report["axes"]["delay"]["error_rate"] < TRUE_DELAY_SIGMA_KM + + def test_the_interval_is_what_carries_the_separation(self): + """Without a spread of prediction intervals there is one equation and + two unknowns, and the honest answer is to decline.""" + records = [r for r in synthesise(seed=5) if prediction_interval(r) == 0.5] + assert "sigma" not in calibrate(records, FC, FS, CPI)["axes"]["delay"] + + def test_a_coasted_update_is_a_longer_prediction(self): + assert prediction_interval({"dt": 0.5, "n_missed": 3}) == pytest.approx(2.0) + + def test_scatter_is_reported_against_what_the_filter_claimed(self): + """Per axis, because NIS sums the two and that is how an oversized R + and an undersized Q have been hiding inside it.""" + report = calibrate(synthesise(seed=6), FC, FS, CPI) + # Far below one is the signature being hunted: the filter is carrying + # an R so large that its innovations never come close to filling it. + assert report["axes"]["delay"]["scatter_vs_claimed"] < 0.1 + + +class TestWhatIsExcluded: + def test_interference_is_dropped_on_the_records_own_verdict(self): + records = synthesise(seed=7)[:200] + for r in records[:120]: + r["interfering"] = True + assert len(usable(records)) == 80 + + def test_a_degenerate_update_is_not_a_measurement(self): + """A singular covariance skipped the measurement, so its innovation is + NaN and there is nothing in it to fit.""" + records = [record(0.5, float("nan"), float("nan")), record(0.5, 0.01, 1e-5)] + assert len(usable(records)) == 1 + + def test_older_records_fall_back_to_a_tracks_doppler_holding_still(self): + """Written before the tracker carried the occupancy map's verdict.""" + tone = [record(0.5, 0.01, 1e-5, doppler=27.9, birth=1) for _ in range(12)] + aircraft = [record(0.5, 0.01, 1e-5, doppler=-120.0 + 4 * i, birth=2) for i in range(12)] + for r in tone + aircraft: + del r["interfering"] + kept = usable(tone + aircraft) + assert {r["birth"] for r in kept} == {2} + + def test_the_fallback_is_not_used_when_a_verdict_is_present(self): + """A real target can hold one Doppler bin. Where the map has spoken, + the heuristic that would convict it must not run.""" + holding = [record(0.5, 0.01, 1e-5, doppler=27.9, birth=1) for _ in range(12)] + assert len(usable(holding)) == 12 + + def test_a_thin_interval_is_not_fitted(self): + assert fit_noise_model([(0.5, 40, 1e-4), (1.0, 40, 2e-4)]) is None + + def test_a_sparse_interval_is_not_a_measurement(self): + records = [record(0.5, 0.01, 1e-5) for _ in range(5)] + assert group_by_interval(records, 0) == [] + + +class TestTheCaptureConfig: + def test_all_three_parameters_come_out_of_blah2s_own_file(self, tmp_path): + path = tmp_path / "config.yml" + path.write_text( + "capture:\n fs: 2000000\n fc: 177000000\n device:\n type: 'RspDuo'\n" + "process:\n data:\n cpi: 0.5\n buffer: 1.5\n" + ) + assert load_blah2_capture(str(path)) == {"fs": 2000000.0, "fc": 177000000.0, "cpi": 0.5} + + +class TestItReadsWhatTheTrackerWrites: + """The contract between the two is a set of field names in a file. A test + that synthesises its own records would never notice one being renamed.""" + + def test_a_real_recording_calibrates(self, tmp_path): + path = tmp_path / "innovations.jsonl" + writer = InnovationWriter(str(path), max_bytes=0) + tracker = Tracker(innovation_writer=writer) + for i in range(60): + tracker.process_frame([{"delay": 20.0 - 0.05 * i, "doppler": -120.0, "snr": 16.0}], 1718747745000 + i * 500) + writer.close() + + records = [json.loads(line) for line in path.read_text().splitlines() if line.strip()] + report = calibrate(records, FC, FS, CPI) + + assert report["n_records"] == len(records) + assert report["n_used"] > 0 + assert report["nis_median"] is not None diff --git a/tests/test_innovations.py b/tests/test_innovations.py index cd7780a..4efc228 100644 --- a/tests/test_innovations.py +++ b/tests/test_innovations.py @@ -93,6 +93,7 @@ def test_a_record_carries_what_calibration_needs(self, tmp_path): "snr", "delay", "doppler", + "interfering", "n_missed", "q_scale", "innovation",