1. Stakeholder Communication Tool
- The charts turn abstract scheduling policies into something tangible:
- SLA lines + tenant latencies → investors & customers can see QoS differentiation.
- Queue depth vs. GPU occupancy → shows why overprovisioning is wasteful and how your scheduler avoids it.
- Job completion distributions → fairness & predictability made concrete.
Instead of raw math, you’ve now got a demo deck that can go in front of:
- A CTO at a cloud company (“here’s how we guarantee tenant A’s LLM inference stays <30ms while tenant C’s batch analytics don’t starve”).
- A DoD program office (“this ensures mission-critical signals intelligence workloads preempt non-critical batch jobs”).
2. Experimentation Sandbox
- Seeds → variability exploration. You already ran seed sweeps and showed variance in SLA violations and job sharing.
- Interference matrices → you can add knobs (e.g., RF vs. NLP jobs interfere at X% more) and run comparative visualizations.
- Scaling → move from 3 tenants to 30 with synthetic traffic models (Poisson arrivals, bursty loads, periodic sweeps).
- Stress test → crank total arrival rate beyond GPU capacity and measure which tenants degrade, and how.
This lets you prove robustness before deploying into messy production workloads.
3. Productization Hook
Your harness is now a demo-able feature:
- Run it live in front of an audience → tweak SLAs interactively and watch graphs react.
- Wrap it in a simple web dashboard (Flask + React) → prospective customers can spin up tenants and see “their” jobs behave under your scheduler.
- Position it as a QoS enforcement simulator:
- Cloud GPU tenancy: Hugging Face, Lambda Labs, RunPod.
- Enterprise AI ops: inference on shared GPU clusters.
- Gov/defense: ISR workloads with strict timelines.
4. Competitive Differentiator
Most schedulers in the wild (e.g. Kubernetes + NVIDIA MPS) don’t show transparent tenant-level metrics.
Your harness makes that visibility a feature:
- “Not only do we schedule fairly, but here’s how you can prove it.”
- That becomes both a sales story and a compliance artifact (think NIST/FedRAMP).
5. Bridge to Real-World Signals
The beauty is: this isn’t limited to synthetic jobs.
Next step = swap MockScheduler with real:
- Torch jobs (
torch.cudaload). - RF sweeps (SoapySDR/BladeRF workloads).
- NLP inference (transformer forward passes).
Then your visualizations become a control tower for real GPU fleets.
cd /home/bgilbert/editable_files/SignalIntelligence && chmod +x test_scheduler.py && ./test_scheduler.py
=== GPU Possession Scheduler Test Harness ===
Target latency SLA: 50.0ms
=== Simulating default QoS (latency factor: 1.0) ===
[Step 00] Latency p95=16.7ms Occ=0.00 QPS=239.8 Batch=4 → Suggest=8
[Step 01] Latency p95=17.4ms Occ=0.11 QPS=348.6 Batch=8 → Suggest=16
[Step 02] Latency p95=30.2ms Occ=0.15 QPS=401.3 Batch=16 → Suggest=32
[Step 03] Latency p95=46.4ms Occ=0.25 QPS=464.3 Batch=32 → Suggest=32
[Step 04] Latency p95=49.3ms Occ=0.31 QPS=501.1 Batch=32 → Suggest=32
[Step 05] Latency p95=49.3ms Occ=0.35 QPS=528.1 Batch=32 → Suggest=32
[Step 06] Latency p95=50.5ms Occ=0.38 QPS=542.3 Batch=32 → Suggest=16
[Step 07] Latency p95=50.4ms Occ=0.37 QPS=542.4 Batch=16 → Suggest=8
[Step 08] Latency p95=50.3ms Occ=0.34 QPS=522.2 Batch=8 → Suggest=4
[Step 09] Latency p95=50.3ms Occ=0.31 QPS=491.3 Batch=4 → Suggest=2
[Step 10] Latency p95=50.2ms Occ=0.29 QPS=461.3 Batch=2 → Suggest=1
[Step 11] Latency p95=50.1ms Occ=0.26 QPS=428.3 Batch=1 → Suggest=1
[Step 12] Latency p95=50.0ms Occ=0.25 QPS=405.1 Batch=1 → Suggest=1
[Step 13] Latency p95=49.9ms Occ=0.23 QPS=382.8 Batch=1 → Suggest=1
[Step 14] Latency p95=49.9ms Occ=0.21 QPS=361.7 Batch=1 → Suggest=1
[Step 15] Latency p95=49.8ms Occ=0.20 QPS=347.0 Batch=1 → Suggest=1
[Step 16] Latency p95=49.7ms Occ=0.19 QPS=333.9 Batch=1 → Suggest=1
[Step 17] Latency p95=49.6ms Occ=0.18 QPS=319.8 Batch=1 → Suggest=1
[Step 18] Latency p95=49.5ms Occ=0.18 QPS=307.6 Batch=1 → Suggest=1
[Step 19] Latency p95=49.4ms Occ=0.17 QPS=299.0 Batch=1 → Suggest=1
[Step 20] Latency p95=49.4ms Occ=0.16 QPS=287.9 Batch=1 → Suggest=1
[Step 21] Latency p95=49.3ms Occ=0.16 QPS=280.2 Batch=1 → Suggest=1
[Step 22] Latency p95=49.3ms Occ=0.15 QPS=273.9 Batch=1 → Suggest=1
[Step 23] Latency p95=49.3ms Occ=0.14 QPS=266.4 Batch=1 → Suggest=1
[Step 24] Latency p95=49.3ms Occ=0.14 QPS=260.2 Batch=1 → Suggest=1
[Step 25] Latency p95=49.3ms Occ=0.13 QPS=253.8 Batch=1 → Suggest=1
[Step 26] Latency p95=49.2ms Occ=0.13 QPS=247.0 Batch=1 → Suggest=1
[Step 27] Latency p95=49.2ms Occ=0.13 QPS=241.1 Batch=1 → Suggest=1
[Step 28] Latency p95=49.2ms Occ=0.12 QPS=235.3 Batch=1 → Suggest=1
[Step 29] Latency p95=49.2ms Occ=0.12 QPS=231.1 Batch=1 → Suggest=1
=== Simulating high QoS (latency factor: 0.8) ===
[Step 00] Latency p95=7.5ms Occ=0.08 QPS=266.6 Batch=2 → Suggest=4
[Step 01] Latency p95=16.5ms Occ=0.05 QPS=250.9 Batch=4 → Suggest=8
[Step 02] Latency p95=19.5ms Occ=0.06 QPS=301.8 Batch=8 → Suggest=8
[Step 03] Latency p95=19.4ms Occ=0.08 QPS=378.5 Batch=8 → Suggest=8
[Step 04] Latency p95=19.3ms Occ=0.08 QPS=402.9 Batch=8 → Suggest=8
[Step 05] Latency p95=19.1ms Occ=0.09 QPS=419.0 Batch=8 → Suggest=8
[Step 06] Latency p95=21.8ms Occ=0.08 QPS=409.8 Batch=8 → Suggest=8
[Step 07] Latency p95=21.6ms Occ=0.08 QPS=421.2 Batch=8 → Suggest=8
[Step 08] Latency p95=21.5ms Occ=0.09 QPS=431.0 Batch=8 → Suggest=8
[Step 09] Latency p95=21.3ms Occ=0.08 QPS=444.9 Batch=8 → Suggest=8
[Step 10] Latency p95=21.2ms Occ=0.09 QPS=461.0 Batch=8 → Suggest=8
[Step 11] Latency p95=21.1ms Occ=0.08 QPS=468.0 Batch=8 → Suggest=8
[Step 12] Latency p95=21.5ms Occ=0.08 QPS=461.6 Batch=8 → Suggest=8
[Step 13] Latency p95=21.4ms Occ=0.08 QPS=463.5 Batch=8 → Suggest=8
[Step 14] Latency p95=21.3ms Occ=0.09 QPS=462.7 Batch=8 → Suggest=8
[Step 15] Latency p95=21.2ms Occ=0.09 QPS=461.5 Batch=8 → Suggest=8
[Step 16] Latency p95=21.1ms Occ=0.09 QPS=460.4 Batch=8 → Suggest=8
[Step 17] Latency p95=21.1ms Occ=0.08 QPS=466.9 Batch=8 → Suggest=8
[Step 18] Latency p95=21.0ms Occ=0.08 QPS=468.5 Batch=8 → Suggest=8
[Step 19] Latency p95=20.9ms Occ=0.09 QPS=470.3 Batch=8 → Suggest=8
=== Simulating low QoS (latency factor: 1.2) ===
[Step 00] Latency p95=25.4ms Occ=0.08 QPS=315.2 Batch=8 → Suggest=16
[Step 01] Latency p95=29.3ms Occ=0.18 QPS=428.4 Batch=16 → Suggest=32
[Step 02] Latency p95=51.3ms Occ=0.28 QPS=484.3 Batch=32 → Suggest=16
[Step 03] Latency p95=50.6ms Occ=0.30 QPS=483.6 Batch=16 → Suggest=8
[Step 04] Latency p95=49.6ms Occ=0.28 QPS=454.9 Batch=8 → Suggest=8
[Step 05] Latency p95=48.6ms Occ=0.25 QPS=441.9 Batch=8 → Suggest=8
[Step 06] Latency p95=47.6ms Occ=0.23 QPS=446.9 Batch=8 → Suggest=8
[Step 07] Latency p95=46.5ms Occ=0.22 QPS=433.6 Batch=8 → Suggest=8
[Step 08] Latency p95=45.5ms Occ=0.22 QPS=419.3 Batch=8 → Suggest=8
[Step 09] Latency p95=44.5ms Occ=0.22 QPS=409.1 Batch=8 → Suggest=8
[Step 10] Latency p95=43.5ms Occ=0.21 QPS=400.9 Batch=8 → Suggest=8
[Step 11] Latency p95=42.4ms Occ=0.20 QPS=393.4 Batch=8 → Suggest=8
[Step 12] Latency p95=41.4ms Occ=0.19 QPS=389.7 Batch=8 → Suggest=8
[Step 13] Latency p95=40.4ms Occ=0.19 QPS=393.2 Batch=8 → Suggest=8
[Step 14] Latency p95=39.4ms Occ=0.19 QPS=391.7 Batch=8 → Suggest=16
[Step 15] Latency p95=39.5ms Occ=0.19 QPS=396.0 Batch=16 → Suggest=32
[Step 16] Latency p95=54.4ms Occ=0.20 QPS=405.6 Batch=32 → Suggest=16
[Step 17] Latency p95=54.2ms Occ=0.21 QPS=409.0 Batch=16 → Suggest=8
[Step 18] Latency p95=54.0ms Occ=0.20 QPS=406.6 Batch=8 → Suggest=4
[Step 19] Latency p95=53.9ms Occ=0.20 QPS=399.3 Batch=4 → Suggest=2
=== Simulating load spike for default QoS ===
[Step 00] Latency p95=49.2ms Occ=0.12 QPS=230.3 Batch=2 → Suggest=2
[Step 01] Latency p95=49.2ms Occ=0.12 QPS=229.1 Batch=2 → Suggest=2
[Step 02] Latency p95=49.1ms Occ=0.11 QPS=229.6 Batch=2 → Suggest=2
[Step 03] Latency p95=49.1ms Occ=0.11 QPS=227.0 Batch=2 → Suggest=2
[Step 04] Latency p95=49.1ms Occ=0.11 QPS=227.3 Batch=2 → Suggest=2
[Step 05] Latency p95=49.1ms Occ=0.11 QPS=226.0 Batch=2 → Suggest=2
[Step 06] Latency p95=49.1ms Occ=0.10 QPS=225.0 Batch=2 → Suggest=2
[Step 07] Latency p95=49.0ms Occ=0.10 QPS=223.2 Batch=2 → Suggest=2
[Step 08] Latency p95=49.0ms Occ=0.10 QPS=221.3 Batch=2 → Suggest=2
[Step 09] Latency p95=49.0ms Occ=0.10 QPS=219.0 Batch=2 → Suggest=2
[Step 10] Latency p95=49.0ms Occ=0.10 QPS=216.8 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 11] Latency p95=48.9ms Occ=0.11 QPS=215.3 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 12] Latency p95=48.9ms Occ=0.11 QPS=213.0 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 13] Latency p95=48.9ms Occ=0.12 QPS=211.2 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 14] Latency p95=48.8ms Occ=0.12 QPS=208.8 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 15] Latency p95=48.8ms Occ=0.13 QPS=206.6 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 16] Latency p95=48.8ms Occ=0.13 QPS=204.7 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 17] Latency p95=48.7ms Occ=0.13 QPS=203.2 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 18] Latency p95=48.7ms Occ=0.14 QPS=201.1 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 19] Latency p95=48.7ms Occ=0.14 QPS=199.2 Batch=2 → Suggest=2 🔥 SPIKE 🔥
[Step 20] Latency p95=48.7ms Occ=0.14 QPS=196.8 Batch=2 → Suggest=2
[Step 21] Latency p95=48.7ms Occ=0.14 QPS=192.5 Batch=2 → Suggest=2
[Step 22] Latency p95=48.7ms Occ=0.13 QPS=186.3 Batch=2 → Suggest=2
[Step 23] Latency p95=39.8ms Occ=0.12 QPS=178.5 Batch=2 → Suggest=4
[Step 24] Latency p95=26.2ms Occ=0.11 QPS=172.0 Batch=4 → Suggest=8
[Step 25] Latency p95=21.1ms Occ=0.11 QPS=166.8 Batch=8 → Suggest=16
[Step 26] Latency p95=21.1ms Occ=0.10 QPS=164.0 Batch=16 → Suggest=32
[Step 27] Latency p95=21.1ms Occ=0.11 QPS=166.4 Batch=32 → Suggest=32
[Step 28] Latency p95=27.0ms Occ=0.11 QPS=172.9 Batch=32 → Suggest=32
[Step 29] Latency p95=40.5ms Occ=0.12 QPS=181.9 Batch=32 → Suggest=32
=== Simulating multi-QoS interaction ===
[Step 00] Shared occupancy: 0.33
HIGH: Latency p95=20.8ms Batch=4 QPS=452.5
DEFAULT: Latency p95=40.5ms Batch=4 QPS=181.7
LOW: Latency p95=53.7ms Batch=4 QPS=391.1
[Step 01] Shared occupancy: 0.55
HIGH: Latency p95=22.5ms Batch=8 QPS=439.9
DEFAULT: Latency p95=40.5ms Batch=4 QPS=183.8
LOW: Latency p95=52.8ms Batch=2 QPS=379.6
[Step 02] Shared occupancy: 0.69
HIGH: Latency p95=22.8ms Batch=8 QPS=433.4
DEFAULT: Latency p95=40.5ms Batch=4 QPS=185.3
LOW: Latency p95=51.8ms Batch=1 QPS=367.3
[Step 03] Shared occupancy: 0.60
HIGH: Latency p95=26.8ms Batch=8 QPS=427.3
DEFAULT: Latency p95=40.5ms Batch=4 QPS=187.5
LOW: Latency p95=50.9ms Batch=1 QPS=354.6
[Step 04] Shared occupancy: 0.61
HIGH: Latency p95=27.8ms Batch=8 QPS=421.3
DEFAULT: Latency p95=40.5ms Batch=4 QPS=190.5
LOW: Latency p95=50.0ms Batch=1 QPS=342.8
[Step 05] Shared occupancy: 0.44
HIGH: Latency p95=27.8ms Batch=8 QPS=418.6
DEFAULT: Latency p95=40.5ms Batch=4 QPS=191.8
LOW: Latency p95=49.1ms Batch=1 QPS=331.5
[Step 06] Shared occupancy: 0.56
HIGH: Latency p95=27.8ms Batch=8 QPS=414.6
DEFAULT: Latency p95=40.5ms Batch=4 QPS=193.7
LOW: Latency p95=48.1ms Batch=1 QPS=321.2
[Step 07] Shared occupancy: 0.65
HIGH: Latency p95=28.5ms Batch=8 QPS=409.2
DEFAULT: Latency p95=40.5ms Batch=4 QPS=195.4
LOW: Latency p95=47.2ms Batch=1 QPS=311.4
[Step 08] Shared occupancy: 0.56
HIGH: Latency p95=28.5ms Batch=8 QPS=407.2
DEFAULT: Latency p95=40.5ms Batch=4 QPS=197.5
LOW: Latency p95=46.3ms Batch=1 QPS=302.2
[Step 09] Shared occupancy: 0.69
HIGH: Latency p95=28.9ms Batch=8 QPS=402.8
DEFAULT: Latency p95=40.5ms Batch=4 QPS=197.8
LOW: Latency p95=45.4ms Batch=1 QPS=294.1
[Step 10] Shared occupancy: 0.64
HIGH: Latency p95=28.9ms Batch=8 QPS=398.9
DEFAULT: Latency p95=40.5ms Batch=4 QPS=200.3
LOW: Latency p95=44.4ms Batch=1 QPS=286.2
[Step 11] Shared occupancy: 0.62
HIGH: Latency p95=28.9ms Batch=8 QPS=396.5
DEFAULT: Latency p95=40.5ms Batch=4 QPS=202.4
LOW: Latency p95=43.5ms Batch=1 QPS=278.7
[Step 12] Shared occupancy: 0.69
HIGH: Latency p95=28.9ms Batch=8 QPS=393.7
DEFAULT: Latency p95=40.5ms Batch=4 QPS=202.8
LOW: Latency p95=42.6ms Batch=1 QPS=272.2
[Step 13] Shared occupancy: 0.54
HIGH: Latency p95=28.9ms Batch=8 QPS=390.7
DEFAULT: Latency p95=40.5ms Batch=4 QPS=204.5
LOW: Latency p95=41.7ms Batch=1 QPS=265.6
[Step 14] Shared occupancy: 0.69
HIGH: Latency p95=29.3ms Batch=8 QPS=387.1
DEFAULT: Latency p95=40.5ms Batch=4 QPS=206.3
LOW: Latency p95=40.7ms Batch=1 QPS=259.3
[Step 15] Shared occupancy: 0.45
HIGH: Latency p95=29.3ms Batch=8 QPS=384.5
DEFAULT: Latency p95=40.5ms Batch=4 QPS=208.1
LOW: Latency p95=39.8ms Batch=2 QPS=253.9
[Step 16] Shared occupancy: 0.31
HIGH: Latency p95=29.2ms Batch=8 QPS=383.1
DEFAULT: Latency p95=40.5ms Batch=4 QPS=210.2
LOW: Latency p95=38.9ms Batch=4 QPS=249.8
[Step 17] Shared occupancy: 0.56
HIGH: Latency p95=29.1ms Batch=8 QPS=380.6
DEFAULT: Latency p95=40.5ms Batch=4 QPS=211.4
LOW: Latency p95=37.9ms Batch=8 QPS=248.2
[Step 18] Shared occupancy: 0.55
HIGH: Latency p95=30.3ms Batch=8 QPS=376.9
DEFAULT: Latency p95=40.5ms Batch=4 QPS=213.3
LOW: Latency p95=37.7ms Batch=16 QPS=247.5
[Step 19] Shared occupancy: 0.65
HIGH: Latency p95=30.6ms Batch=8 QPS=372.9
DEFAULT: Latency p95=41.5ms Batch=4 QPS=213.5
LOW: Latency p95=46.2ms Batch=16 QPS=250.1
[Step 20] Shared occupancy: 0.65
HIGH: Latency p95=34.2ms Batch=8 QPS=369.0
DEFAULT: Latency p95=42.5ms Batch=4 QPS=211.5
LOW: Latency p95=50.8ms Batch=8 QPS=251.7
[Step 21] Shared occupancy: 0.55
HIGH: Latency p95=34.0ms Batch=8 QPS=367.3
DEFAULT: Latency p95=42.5ms Batch=4 QPS=210.3
LOW: Latency p95=50.6ms Batch=4 QPS=251.6
[Step 22] Shared occupancy: 0.54
HIGH: Latency p95=33.8ms Batch=8 QPS=365.6
DEFAULT: Latency p95=42.5ms Batch=4 QPS=208.8
LOW: Latency p95=50.3ms Batch=2 QPS=250.2
[Step 23] Shared occupancy: 0.51
HIGH: Latency p95=33.6ms Batch=8 QPS=364.9
DEFAULT: Latency p95=42.5ms Batch=4 QPS=209.8
LOW: Latency p95=50.1ms Batch=1 QPS=246.7
[Step 24] Shared occupancy: 0.44
HIGH: Latency p95=33.4ms Batch=8 QPS=363.9
DEFAULT: Latency p95=42.5ms Batch=4 QPS=209.2
LOW: Latency p95=49.8ms Batch=1 QPS=242.2
[Step 25] Shared occupancy: 0.38
HIGH: Latency p95=33.2ms Batch=8 QPS=363.4
DEFAULT: Latency p95=42.5ms Batch=4 QPS=209.5
LOW: Latency p95=49.6ms Batch=1 QPS=238.1
[Step 26] Shared occupancy: 0.42
HIGH: Latency p95=33.0ms Batch=8 QPS=363.1
DEFAULT: Latency p95=42.5ms Batch=4 QPS=208.9
LOW: Latency p95=49.3ms Batch=1 QPS=234.0
[Step 27] Shared occupancy: 0.54
HIGH: Latency p95=32.9ms Batch=8 QPS=363.0
DEFAULT: Latency p95=42.5ms Batch=4 QPS=209.5
LOW: Latency p95=49.1ms Batch=1 QPS=230.0
[Step 28] Shared occupancy: 0.41
HIGH: Latency p95=32.7ms Batch=8 QPS=361.2
DEFAULT: Latency p95=42.5ms Batch=4 QPS=209.6
LOW: Latency p95=48.8ms Batch=1 QPS=226.1
[Step 29] Shared occupancy: 0.68
HIGH: Latency p95=32.5ms Batch=8 QPS=361.6
DEFAULT: Latency p95=42.5ms Batch=4 QPS=210.2
LOW: Latency p95=48.6ms Batch=1 QPS=222.9
=== Simulation Complete ===
Summary of final recommendations:
HIGH: Suggested batch=8, P95 Latency=32.5ms, QPS=361.6
DEFAULT: Suggested batch=4, P95 Latency=42.5ms, QPS=210.2
LOW: Suggested batch=1, P95 Latency=48.6ms, QPS=222.9
2
Semi Related:
# We'll generate an updated trajectory inference module with TDoA support
# and a tiny demo that synthesizes TDoA measurements and runs inference.
from textwrap import dedent
from pathlib import Path
module_code = dedent(r"""
# rf_sequence_recovery_tdoa.py
# Trajectory inference from sparse RF sightings with full TDoA residuals.
# Drop-in: replace/augment your previous rf_sequence_recovery.py with this file.
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Dict, List, Tuple, Optional
import numpy as np
import math
import heapq
# --------------------------
# Constants
# --------------------------
C_MPS = 299_792_458.0 # speed of light (m/s)
# --------------------------
# Data containers
# --------------------------
@dataclass
class Sensor:
sensor_id: str
xy: np.ndarray # shape (2,), meters in a local ENU
clock_bias_s: float = 0.0 # optional fixed clock bias (sec); use 0 if unknown
@dataclass
class RFSighting:
t: float # seconds
xy: Optional[np.ndarray] = None # rough geofix (m); None if not available
sigma_pos: Optional[float] = None # 1-sigma of geofix (m)
sensor_id: Optional[str] = None
aoa_rad: Optional[float] = None # bearing at sensor (rad), 0 = +x
sigma_aoa: Optional[float] = None # rad
snr_db: Optional[float] = None
@dataclass
class TDoAMeasurement:
"""A single TDoA observation between two sensors (i, j) near time t.
tdoa_s = t_i - t_j (time of arrival difference in seconds).
sigma_s = standard deviation (seconds).
Optionally set weight to down/up-weight this pair (e.g., by cross-corr. peak quality).
"""
t: float
sensor_i: str
sensor_j: str
tdoa_s: float
sigma_s: float
weight: float = 1.0
# --------------------------
# Mobility graph
# --------------------------
class MobilityGraph:
def __init__(self, nodes_xy: np.ndarray, adj: List[List[int]],
edge_dist: np.ndarray, edge_vmax: np.ndarray,
default_vmax: float = 15.0):
self.nodes_xy = np.asarray(nodes_xy, dtype=float) # (N,2)
self.adj = adj # list of neighbor indices
self.edge_dist = edge_dist # (E,) meters, matched with adj pairs
self.edge_vmax = edge_vmax # (E,) m/s
self.default_vmax = float(default_vmax)
@staticmethod
def grid(x_min: float, x_max: float, y_min: float, y_max: float, step: float,
default_vmax: float = 15.0, eight_connected: bool = True) -> "MobilityGraph":
xs = np.arange(x_min, x_max + 1e-9, step, dtype=float)
ys = np.arange(y_min, y_max + 1e-9, step, dtype=float)
X, Y = np.meshgrid(xs, ys)
nodes = np.stack([X.ravel(), Y.ravel()], axis=1)
nx, ny = len(xs), len(ys)
def idx(ix, iy): return iy*nx + ix
adj: List[List[int]] = [[] for _ in range(nodes.shape[0])]
edges = []
vmaxs = []
for iy in range(ny):
for ix in range(nx):
u = idx(ix, iy)
for dx, dy in [(1,0),(0,1)] + ([(1,1),(-1,1)] if eight_connected else []):
jx, jy = ix+dx, iy+dy
if 0 <= jx < nx and 0 <= jy < ny:
v = idx(jx, jy)
adj[u].append(v)
adj[v].append(u)
dist = np.linalg.norm(nodes[v]-nodes[u])
edges.append(dist)
vmaxs.append(default_vmax)
return MobilityGraph(nodes, adj, np.array(edges, dtype=float),
np.array(vmaxs, dtype=float), default_vmax=default_vmax)
# --------------------------
# Loss helpers
# --------------------------
def wrap_angle(a: float) -> float:
"""Wrap angle to [-pi, pi]."""
a = (a + np.pi) % (2*np.pi) - np.pi
return a
def huber_loss(r: float, delta: float = 2.5) -> float:
"""Robustify extreme outliers (unitless residual)."""
r = float(r)
ad = abs(r)
if ad <= delta:
return 0.5 * r * r
return delta * (ad - 0.5*delta)
# --------------------------
# Trajectory inference (beam-Viterbi on graph with observation model)
# --------------------------
class TrajectoryInferrer:
def __init__(self, graph: MobilityGraph, sensors: Dict[str, Sensor],
dt: float = 5.0, beam_width: int = 256,
vmax_hard: float = 40.0, stay_bias: float = 0.25,
huber_delta: float = 3.0, c_mps: float = C_MPS):
self.G = graph
self.S = sensors
self.dt = float(dt)
self.beam_width = int(beam_width)
self.vmax_hard = float(vmax_hard)
self.stay_bias = float(stay_bias)
self.huber_delta = float(huber_delta)
self.c_mps = float(c_mps)
# ----------------------
# Public API
# ----------------------
def infer(self,
sightings: List[RFSighting],
tdoas: Optional[List[TDoAMeasurement]] = None,
t0: float = 0.0, t1: float = 300.0,
window: float = 10.0) -> Dict[str, np.ndarray]:
"""Return MAP path, per-step entropy proxy, and times."""
t_axis = np.arange(t0, t1 + 1e-6, self.dt, dtype=float)
# Pre-index measurements by time for efficiency
s_by_t = self._bucketize(sightings, t_axis, window)
d_by_t = self._bucketize(tdoas or [], t_axis, window)
# Initialize beam uniformly (or from first measurement peak if desired)
N = self.G.nodes_xy.shape[0]
# Simple heuristic: seed beam by top-k nodes closest to first-pos sighting if available
init_scores = np.zeros(N, dtype=float)
if len(sightings) > 0 and any(s.xy is not None for s in sightings):
xy0 = min([s for s in sightings if s.xy is not None], key=lambda s: s.t).xy
d = np.linalg.norm(self.G.nodes_xy - xy0[None,:], axis=1)
init_scores = 0.5*(d/ (np.median(d)+1e-6))**2
init_beam = [(init_scores[i], i, -1) for i in np.argpartition(init_scores, self.beam_width)[:self.beam_width]]
prev = {i: (cost, parent) for cost, i, parent in init_beam}
backptr: List[Dict[int, Tuple[float,int]]] = []
for step, t in enumerate(t_axis):
# Expand
cand: Dict[int, Tuple[float,int]] = {}
for i, (cost_i, parent_i) in prev.items():
# Stay option
self._accumulate_candidate(cand, i, cost_i + self._trans_cost(i, i), i)
# Move options
for j in self.G.adj[i]:
tc = self._trans_cost(i, j)
if tc == np.inf:
continue
self._accumulate_candidate(cand, j, cost_i + tc, i)
# Observation cost at this time
obs = s_by_t.get(step, [])
dmeas = d_by_t.get(step, [])
for j in list(cand.keys()):
oc = self._obs_cost(self.G.nodes_xy[j], t, obs, dmeas)
cost_prev, p = cand[j]
cand[j] = (cost_prev + oc, p)
# Prune to beam
if len(cand) == 0:
# Fall back—keep previous
cand = prev.copy()
top = sorted(cand.items(), key=lambda kv: kv[1][0])[:self.beam_width]
backptr.append({node:(cost,parent) for node,(cost,parent) in top})
prev = {node:(cost,parent) for node,(cost,parent) in top}
# Backtrack best
last_nodes = min(prev.items(), key=lambda kv: kv[1][0])[0]
path_idx = [last_nodes]
for bp in reversed(backptr):
node = path_idx[-1]
cost, parent = bp[node]
if parent < 0:
break
path_idx.append(parent)
path_idx = path_idx[::-1]
path_xy = self.G.nodes_xy[np.array(path_idx[:len(t_axis)])]
# Entropy proxy: log-sum-exp of beam vs best at each step
ent = np.zeros(len(backptr), dtype=float)
for k, bp in enumerate(backptr):
costs = np.array([c for (_, (c,_)) in bp.items()])
m = np.min(costs)
ent[k] = m + np.log(np.exp(-(costs - m)).sum())
ent = ent - np.min(ent)
return {
"times": t_axis,
"path_xy": path_xy,
"entropy": ent,
"path_idx": np.array(path_idx[:len(t_axis)]),
}
# ----------------------
# Internals
# ----------------------
def _bucketize(self, items, t_axis, window):
buckets = {}
if items is None:
return buckets
half = window/2.0
for it in items:
# find nearest time index
idx = int(np.clip(round((it.t - t_axis[0]) / self.dt), 0, len(t_axis)-1))
# only assign if within window
t_center = t_axis[idx]
if abs(it.t - t_center) <= half:
buckets.setdefault(idx, []).append(it)
return buckets
def _trans_cost(self, i: int, j: int) -> float:
if i == j:
return self.stay_bias
dx = self.G.nodes_xy[j] - self.G.nodes_xy[i]
dist = float(np.linalg.norm(dx))
v = dist / max(self.dt, 1e-6)
if v > self.vmax_hard:
return np.inf
# soft penalty above default vmax
pen = max(0.0, v - self.G.default_vmax)
return 0.5 * (pen / (0.25*self.G.default_vmax + 1e-6))**2 + dist*1e-3 # small distance prior
def _accumulate_candidate(self, cand: Dict[int, Tuple[float,int]], node: int,
cost: float, parent: int):
cur = cand.get(node, (np.inf, -1))
if cost < cur[0]:
cand[node] = (cost, parent)
# Observation model
def _obs_cost(self, xy_node: np.ndarray, t: float,
sightings: List[RFSighting],
tdoas: List[TDoAMeasurement]) -> float:
cost = 0.0
# Rough fixes
for s in sightings:
if s.xy is not None and (s.sigma_pos is not None) and s.sigma_pos > 0:
r = np.linalg.norm(xy_node - s.xy) / s.sigma_pos
cost += huber_loss(r, self.huber_delta)
# AoA residual if available
if s.aoa_rad is not None and s.sensor_id is not None and s.sigma_aoa:
sens = self.S.get(s.sensor_id)
if sens is not None:
vec = xy_node - sens.xy
pred = math.atan2(vec[1], vec[0])
da = wrap_angle(pred - s.aoa_rad) / max(s.sigma_aoa, 1e-6)
w = 1.0
if s.snr_db is not None:
# map SNR to 0.5..2.0 weight
w = float(np.clip(0.5 + 0.05*(s.snr_db - 10.0), 0.5, 2.0))
cost += w * huber_loss(da, self.huber_delta)
# TDoA residuals
for d in tdoas:
si = self.S.get(d.sensor_i)
sj = self.S.get(d.sensor_j)
if (si is None) or (sj is None):
continue
# geometric time difference for node XY
ri = float(np.linalg.norm(xy_node - si.xy)) / self.c_mps + si.clock_bias_s
rj = float(np.linalg.norm(xy_node - sj.xy)) / self.c_mps + sj.clock_bias_s
tau_pred = ri - rj # seconds
# residual (pred - measured)
res = (tau_pred - d.tdoa_s) / max(d.sigma_s, 1e-9)
cost += float(d.weight) * huber_loss(res, self.huber_delta)
return cost
""")
demo_code = dedent(r"""
# demo_rf_sequence_recovery_tdoa.py
# Tiny synthetic demo: emit a path, generate AoA + TDoA, and recover.
import numpy as np
import matplotlib.pyplot as plt
from rf_sequence_recovery_tdoa import (
MobilityGraph, Sensor, RFSighting, TDoAMeasurement, TrajectoryInferrer, C_MPS
)
rng = np.random.default_rng(7)
# 1) Build grid graph
G = MobilityGraph.grid(-2000, 2000, -2000, 2000, step=100.0, default_vmax=20.0)
# 2) Sensors
S = {
"A": Sensor("A", np.array([-1500.0, 0.0])),
"B": Sensor("B", np.array([ 1500.0, 200.0])),
"C": Sensor("C", np.array([ 600.0,-1400.0])),
}
# 3) True path (smooth curve)
T = 160
dt = 2.0
t_axis = np.arange(0, T*dt, dt)
px = -1200 + 1800*np.sin(np.linspace(0, 2.5, len(t_axis)))
py = -800 + 1200*np.sin(np.linspace(0, 1.25, len(t_axis)) + 0.5)
true_xy = np.stack([px, py], axis=1)
# helper: nearest grid node
def snap(xy):
d = np.linalg.norm(G.nodes_xy - xy[None,:], axis=1)
return G.nodes_xy[np.argmin(d)]
# 4) Generate sightings (every ~8 steps)
sightings = []
for k, t in enumerate(t_axis):
if k % 8 == 0:
# rough fix with 150m sigma
xy = true_xy[k] + rng.normal(0, 150.0, size=2)
s = RFSighting(t=float(t), xy=xy, sigma_pos=150.0)
sightings.append(s)
if k % 5 == 0:
# AoA at a random sensor with 6deg sigma
sid = rng.choice(list(S.keys()))
sens = S[sid]
vec = true_xy[k] - sens.xy
aoa = np.arctan2(vec[1], vec[0]) + rng.normal(0, np.deg2rad(6.0))
sightings.append(RFSighting(t=float(t), sensor_id=sid, aoa_rad=float(aoa),
sigma_aoa=np.deg2rad(6.0), snr_db=15.0 + 10*rng.random()))
# 5) Generate TDoA pairs (A-B, A-C) every ~4 steps
tdoas = []
sigma_t = 30e-9 # 30 ns (~9 m range-diff sigma)
for k, t in enumerate(t_axis):
if k % 4 == 0:
for (i,j) in [("A","B"),("A","C")]:
ri = np.linalg.norm(true_xy[k] - S[i].xy)/C_MPS + S[i].clock_bias_s
rj = np.linalg.norm(true_xy[k] - S[j].xy)/C_MPS + S[j].clock_bias_s
tau = ri - rj + rng.normal(0, sigma_t)
tdoas.append(TDoAMeasurement(t=float(t), sensor_i=i, sensor_j=j,
tdoa_s=float(tau), sigma_s=float(sigma_t), weight=1.0))
# 6) Infer
inf = TrajectoryInferrer(G, S, dt=dt, beam_width=256, vmax_hard=45.0)
result = inf.infer(sightings, tdoas, t0=float(t_axis[0]), t1=float(t_axis[-1]), window=6.0)
# 7) Plot
plt.figure(figsize=(8,7))
plt.scatter([S[k].xy[0] for k in S], [S[k].xy[1] for k in S], c="k", marker="^", label="Sensors")
for sid in S:
plt.text(S[sid].xy[0]+20, S[sid].xy[1]+20, sid, fontsize=9)
plt.plot(true_xy[:,0], true_xy[:,1], "g-", lw=2, label="True path")
plt.plot(result["path_xy"][:,0], result["path_xy"][:,1], "b--", lw=2, label="Inferred")
# rough fix points
rx = [s.xy[0] for s in sightings if s.xy is not None]
ry = [s.xy[1] for s in sightings if s.xy is not None]
plt.scatter(rx, ry, c="orange", s=20, alpha=0.6, label="Rough fixes")
plt.axis("equal")
plt.grid(True, ls="--", alpha=0.4)
plt.legend()
plt.title("Trajectory inference with AoA + TDoA residuals")
plt.tight_layout()
plt.savefig("rf_sequence_tdoa_demo.png")
# Uncertainty plot
import matplotlib.ticker as mt
plt.figure(figsize=(8,3))
plt.plot(result["times"][:len(result["entropy"])], result["entropy"], "m-")
plt.gca().yaxis.set_major_locator(mt.MaxNLocator(4))
plt.grid(True, ls="--", alpha=0.4)
plt.xlabel("Time (s)")
plt.ylabel("Entropy (arb.)")
plt.title("Entropy proxy over time (lower = more confident)")
plt.tight_layout()
plt.savefig("rf_sequence_tdoa_entropy.png")
print("Saved: rf_sequence_tdoa_demo.png, rf_sequence_tdoa_entropy.png")
""")
Path("/mnt/data/rf_sequence_recovery_tdoa.py").write_text(module_code)
Path("/mnt/data/demo_rf_sequence_recovery_tdoa.py").write_text(demo_code)
print("Wrote /mnt/data/rf_sequence_recovery_tdoa.py and demo_rf_sequence_recovery_tdoa.py")