Simulates drift due to packet timing affecting moving physical object such as our bullet

bullet_interp_bench.py

#!/usr/bin/env python3
"""Bullet interpolation benchmark: Portland datacenter -> London, UK network model.

Models a 200 m/s projectile as the simulator would present it: 45 Hz physics ticks,
divergence-triggered terse updates (velocity-only prediction credit, matching the
observed update stream), U16-quantized payloads, packet bundling, and a
transatlantic backbone path with hop queuing, bufferbloat, packet loss, sim stalls
and viewer frame hitches. Updates are stamped at viewer pump time (no timestamp in
the protocol), exactly like ObjectUpdate traffic.

Viewer correction modes compared against ground truth:
  vanilla   - current viewer: snap to update at arrival time (llviewerobject.cpp:2527)
  ping_only - the disabled PingInterpolate: blind shift forward by estimated ping/2
  blend     - vanilla logical state, rendered position exponentially smoothed (TC=0.15s)
  timing    - timing inversion: solve update age from the along-track residual, clamp
              to a jitter window, back-date the anchor; fall back to vanilla snap when
              the residual is not timing-shaped (perpendicular gate)
"""

import math
import random
import statistics as st

DT_PHYS = 1.0 / 45.0
DT_FRAME = 1.0 / 60.0
FLIGHT_SECONDS = 12.0
GRAVITY = (0.0, 0.0, -9.81)

ROUTE_KM = 8100.0 * 1.35
FIBER_KMS_PER_MS = 200.0
PROPAGATION_MS = ROUTE_KM / FIBER_KMS_PER_MS
HOPS = 16
CONGESTED_HOPS = 2

PACKET_FLUSH_S = 0.020
MAX_BUNDLE = 6
LOSS_PCT = 0.005
BURST_START_P = 0.0005
BURST_MAX = 5

VIEWER_HITCH_P = 0.002
VIEWER_HITCH_MS = (40.0, 130.0)
PING_UPDATE_S = 1.0

QUANT_STEP_XY = 512.0 / 65535.0
QUANT_STEP_Z = 4196.0 / 65535.0
QUANT_STEP_VEL = 512.0 / 65535.0

SIM_POS_DIV_THRESHOLD_M = 0.05
SPEED_GATE_MPS = 5.0
PERP_GATE_M = 0.25
WINDOW_S = 0.200
PRIOR_BOUND_S = 0.250
BLEND_TC_S = 0.15

N_SEEDS = 8
BASE_SEED = 1000


def v_add(a, b, s=1.0):
    return (a[0] + b[0] * s, a[1] + b[1] * s, a[2] + b[2] * s)


def v_sub(a, b):
    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])


def v_scale(a, s):
    return (a[0] * s, a[1] * s, a[2] * s)


def v_dot(a, b):
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def v_len(a):
    return math.sqrt(v_dot(a, a))


def v_norm(a):
    l = v_len(a)
    return (a[0] / l, a[1] / l, a[2] / l) if l > 1e-9 else (0.0, 0.0, 0.0)


def quantize(val, step):
    return round(val / step) * step


def quantize_vec(p):
    return (
        quantize(p[0], QUANT_STEP_XY),
        quantize(p[1], QUANT_STEP_XY),
        quantize(p[2], QUANT_STEP_Z),
    )


def quantize_vel(v):
    return tuple(quantize(c, QUANT_STEP_VEL) for c in v)


class Sim:
    def __init__(self, scenario, rng):
        self.scenario = scenario
        self.rng = rng
        self.t = 0.0
        self.p = scenario["p0"]
        self.v = scenario["v0"]
        self.a = scenario["a"]
        self.pending = []
        self.packets = []
        self.last_sent = None
        self.next_flush = PACKET_FLUSH_S
        self.next_tick = DT_PHYS

    def truth(self, t):
        s = self.scenario
        dt = max(0.0, t)
        return v_add(v_add(s["p0"], s["v0"], dt), s["a"], 0.5 * dt * dt)

    def tick_until(self, now):
        while self.next_tick <= now:
            dt = DT_PHYS
            self.t = self.next_tick
            self.p = v_add(v_add(self.p, self.v, dt), self.a, 0.5 * dt * dt)
            self.v = v_add(self.v, self.a, dt)
            self.next_tick += dt
            if self.scenario["sim_stalls"] and self.rng.random() < 0.004:
                self.next_tick += self.rng.uniform(0.04, 0.2)
            if self.last_sent is None:
                self.emit()
                continue
            t_since = self.t - self.last_sent[0]
            pred = v_add(self.last_sent[1], self.last_sent[2], t_since)
            if v_len(v_sub(self.p, pred)) > SIM_POS_DIV_THRESHOLD_M:
                self.emit()

    def emit(self):
        self.last_sent = (self.t, self.p, self.v, self.a)
        v_avg = v_add(self.v, self.a, -0.5 * DT_PHYS)
        self.pending.append(
            {
                "t_send": self.t,
                "p": quantize_vec(self.p),
                "v": quantize_vel(v_avg),
                "a": quantize_vel(self.a),
            }
        )

    def flush_until(self, now):
        while self.next_flush <= now:
            if self.pending:
                bundle = self.pending[:MAX_BUNDLE]
                self.pending = self.pending[MAX_BUNDLE:]
                self.packets.append({"t_send": self.next_flush, "updates": bundle, "t_arrival": None})
            self.next_flush += PACKET_FLUSH_S


class Network:
    def __init__(self, scenario, rng):
        self.scenario = scenario
        self.rng = rng
        self.bloat = False
        self.burst_left = 0

    def sample_latency_ms(self):
        rng = self.rng
        if self.scenario["bufferbloat"]:
            if self.bloat:
                if rng.random() < 0.05:
                    self.bloat = False
            elif rng.random() < 0.002:
                self.bloat = True
        lat = PROPAGATION_MS
        for _ in range(HOPS):
            lat += rng.expovariate(1.0)
        for _ in range(CONGESTED_HOPS):
            lat += rng.expovariate(0.25)
        if self.bloat:
            lat += rng.uniform(15.0, 50.0)
        if self.scenario["spikes"] and rng.random() < 0.004:
            lat += rng.uniform(30.0, 120.0)
        return lat

    def lost(self):
        rng = self.rng
        if not self.scenario["loss"]:
            return False
        if self.burst_left > 0:
            self.burst_left -= 1
            return True
        if rng.random() < BURST_START_P:
            self.burst_left = rng.randint(1, BURST_MAX)
            return True
        return rng.random() < LOSS_PCT


class Viewer:
    def __init__(self, scenario, mode, rng):
        self.scenario = scenario
        self.mode = mode
        self.rng = rng
        self.t = 0.0
        self.anchor = None
        self.ping_ms = 2.0 * (PROPAGATION_MS + HOPS * 1.0 + CONGESTED_HOPS * 4.0)
        self.l_prior_s = self.ping_ms / 2000.0
        self.taus = []
        self.accepted = 0
        self.rejected = 0
        self.tau_center = None
        self.prev_v = None
        self.last_arrival = None
        self.smooth_p = None
        self.last_t = 0.0

    def observe_rtt(self, mean_one_way_ms):
        self.ping_ms = self.ping_ms * 0.7 + (mean_one_way_ms * 2.0) * 0.3 + self.rng.gauss(0.0, 4.0)

    def expected(self, t):
        p, v, a, t0 = self.anchor
        dt = t - t0
        return v_add(v_add(p, v, dt), a, 0.5 * dt * dt)

    def render(self, t, dt_frame):
        if self.mode == "blend" and self.anchor is not None:
            target = self.expected(t)
            if self.smooth_p is None:
                self.smooth_p = target
            k = 1.0 - math.exp(-dt_frame / BLEND_TC_S)
            self.smooth_p = v_add(self.smooth_p, v_sub(target, self.smooth_p), k)
            return self.smooth_p
        return self.expected(t) if self.anchor is not None else self.scenario["p0"]

    def apply_update(self, upd, t_arrival):
        p_recv, v_recv, a_recv = upd["p"], upd["v"], upd["a"]
        if self.anchor is None:
            self.anchor = (p_recv, v_recv, a_recv, t_arrival)
            self.smooth_p = p_recv
            return

        if self.mode in ("vanilla", "blend"):
            self.anchor = (p_recv, v_recv, a_recv, t_arrival)
            return

        if self.mode == "ping_only":
            l0 = self.ping_ms / 2000.0
            shifted = v_add(v_add(p_recv, v_recv, l0), a_recv, 0.5 * l0 * l0)
            self.anchor = (shifted, v_recv, a_recv, t_arrival)
            return

        if self.mode == "timing":
            l0 = self.ping_ms / 2000.0
            p_exp = self.expected(t_arrival - l0)
            r = v_sub(p_recv, p_exp)
            speed = v_len(v_recv)
            if speed < SPEED_GATE_MPS:
                self.anchor = (p_recv, v_recv, a_recv, t_arrival)
                return
            v_hat = v_norm(v_recv)
            along = v_dot(r, v_hat)
            perp = v_sub(r, v_scale(v_hat, along))
            tau = along / speed
            self.taus.append(tau)
            center = self.tau_center if self.tau_center is not None else tau
            dt_upd = t_arrival - self.last_arrival if self.last_arrival is not None else 0.0
            dv = v_len(v_sub(v_recv, self.prev_v)) if self.prev_v is not None else 0.0
            real_change = dt_upd > 0.0 and dv > max(1.0, 3.0 * v_len(a_recv) * dt_upd)
            if abs(tau - center) <= WINDOW_S and v_len(perp) < PERP_GATE_M and not real_change:
                self.accepted += 1
                l_hat = min(max(l0 - tau, l0 - PRIOR_BOUND_S), l0 + PRIOR_BOUND_S)
                self.anchor = (p_recv, v_recv, a_recv, t_arrival - l_hat)
            else:
                self.rejected += 1
                self.anchor = (p_recv, v_recv, a_recv, t_arrival)
            self.tau_center = center * 0.8 + tau * 0.2
            self.prev_v = v_recv
            self.last_arrival = t_arrival


def run_flight(scenario, mode, seed):
    rng = random.Random(seed)
    sim = Sim(scenario, rng)
    net = Network(scenario, rng)
    viewer = Viewer(scenario, mode, rng)

    frames = []
    updates = 0
    lost = 0
    delivered_latencies = []
    next_ping = PING_UPDATE_S

    while viewer.t < FLIGHT_SECONDS:
        dt_frame = max(0.004, rng.gauss(DT_FRAME, 0.0015))
        if scenario["viewer_hitches"] and rng.random() < VIEWER_HITCH_P:
            dt_frame += rng.uniform(*VIEWER_HITCH_MS) / 1000.0
        viewer.t += dt_frame
        now = viewer.t

        sim.tick_until(now)
        sim.flush_until(now)

        kept = []
        for pkt in sim.packets:
            if pkt["t_arrival"] is None:
                if net.lost():
                    lost += len(pkt["updates"])
                    continue
                pkt["t_arrival"] = pkt["t_send"] + net.sample_latency_ms() / 1000.0
            if pkt["t_arrival"] <= now:
                delivered_latencies.append((pkt["t_arrival"] - pkt["t_send"]) * 1000.0)
                for upd in pkt["updates"]:
                    updates += 1
                    viewer.apply_update(upd, now)
            else:
                kept.append(pkt)
        sim.packets = kept

        if now >= next_ping and delivered_latencies:
            viewer.observe_rtt(st.mean(delivered_latencies[-20:]))
            next_ping += PING_UPDATE_S

        truth = sim.truth(now)
        if viewer.anchor is not None:
            rendered = viewer.render(now, dt_frame)
            frames.append((now, rendered, truth))

    return {"frames": frames, "updates": updates, "lost": lost, "taus": viewer.taus}


def summarize(scenario, seeds, modes):
    results = {}
    for mode in modes:
        rows = []
        for seed in seeds:
            m = run_flight(scenario, mode, seed)
            offs = [v_sub(rp, tp) for _, rp, tp in m["frames"]]
            mean_off = tuple(st.mean(o[i] for o in offs) for i in range(3))
            dyn = [v_sub(o, mean_off) for o in offs]
            wobble = sorted(v_len(d) for d in dyn)
            v_hat = v_norm(scenario["v0"])
            lateral = [abs(v_dot(d, v_hat)) for d in dyn]
            frames = m["frames"]
            snaps = 0
            snap_max = 0.0
            for i in range(1, len(frames)):
                dtf = frames[i][0] - frames[i - 1][0]
                dr = v_len(v_sub(frames[i][1], frames[i - 1][1]))
                dt = v_len(v_sub(frames[i][2], frames[i - 1][2]))
                anomaly = abs(dr - dt)
                snap_max = max(snap_max, anomaly)
                if anomaly > 1.0:
                    snaps += 1
            speeds = [
                v_len(v_sub(frames[i][1], frames[i - 1][1])) / (frames[i][0] - frames[i - 1][0])
                for i in range(1, len(frames))
            ]
            tspeeds = [
                v_len(v_sub(frames[i][2], frames[i - 1][2])) / (frames[i][0] - frames[i - 1][0])
                for i in range(1, len(frames))
            ]
            jit = [a - b for a, b in zip(speeds, tspeeds)]
            jit_rms = math.sqrt(st.mean(j * j for j in jit)) if jit else 0.0
            rows.append(
                {
                    "wobble_rms": st.mean(wobble),
                    "wobble_p95": wobble[int(0.95 * len(wobble))],
                    "wobble_max": wobble[-1],
                    "lateral_rms": st.mean(lateral),
                    "speed_cv": (st.pstdev(speeds) / st.mean(speeds)) if speeds else 0.0,
                    "jit_rms": jit_rms,
                    "snap_frames": snaps,
                    "snap_max": snap_max,
                    "updates": m["updates"],
                    "lost": m["lost"],
                }
            )
        results[mode] = rows
    return results


SCENARIOS = {
    "A vertical 200m/s gravity, harsh transatlantic": dict(
        p0=(128.0, 128.0, 40.0), v0=(0.0, 0.0, 200.0), a=(0.0, 0.0, -9.81),
        sim_stalls=True, viewer_hitches=True, bufferbloat=True, spikes=True, loss=True,
    ),
    "B vertical 200m/s gravity, clean link": dict(
        p0=(128.0, 128.0, 40.0), v0=(0.0, 0.0, 200.0), a=(0.0, 0.0, -9.81),
        sim_stalls=False, viewer_hitches=False, bufferbloat=False, spikes=False, loss=False,
    ),
    "C horizontal 200m/s buoyant, sim goes quiet": dict(
        p0=(128.0, 128.0, 40.0), v0=(200.0, 0.0, 0.0), a=(0.0, 0.0, 0.0),
        sim_stalls=False, viewer_hitches=False, bufferbloat=False, spikes=False, loss=False,
    ),
}

MODES = ["vanilla", "ping_only", "blend", "timing"]


def main():
    global WINDOW_S
    print(f"PDX->LHR: geodesic 8100 km x1.35 route, fiber {FIBER_KMS_PER_MS:.0f} km/ms -> propagation {PROPAGATION_MS:.1f} ms one-way")
    print(f"backbone: {HOPS} hops (exp 1 ms) + {CONGESTED_HOPS} congested (exp 4 ms); bufferbloat +15-50 ms episodes; spikes +30-120 ms (0.4%); loss {LOSS_PCT*100:.1f}% + bursts")
    print(f"sim: {1/DT_PHYS:.0f} Hz ticks with stalls, divergence threshold {SIM_POS_DIV_THRESHOLD_M} m (velocity-only prediction), flush {PACKET_FLUSH_S*1000:.0f} ms bundling <= {MAX_BUNDLE}")
    print(f"viewer: 60 fps (sigma 1.5 ms) + hitches, pump-time stamping, U16 quantization; ping estimated from delivered-packet mean every {PING_UPDATE_S:.0f} s")
    print(f"timing mode: window {WINDOW_S*1000:.0f} ms (centered), perp gate {PERP_GATE_M} m, speed gate {SPEED_GATE_MPS} m/s; {N_SEEDS} seeds x {FLIGHT_SECONDS:.0f} s\n")

    header = f"{'scenario':<48} {'mode':<10} {'wob rms':>8} {'wob p95':>8} {'wob max':>8} {'spd err':>8} {'snap>1m':>8} {'upd':>5}"
    print(header)
    print("-" * len(header))
    for name, sc in SCENARIOS.items():
        results = summarize(sc, range(BASE_SEED, BASE_SEED + N_SEEDS), MODES)
        for mode in MODES:
            rows = results[mode]
            agg = lambda k: st.mean(r[k] for r in rows)
            print(
                f"{name:<48} {mode:<10} {agg('wobble_rms'):>8.3f} {agg('wobble_p95'):>8.3f} "
                f"{agg('wobble_max'):>8.3f} {agg('jit_rms'):>8.3f} "
                f"{agg('snap_frames'):>8.1f} {agg('updates'):>5.0f}"
            )
        print()

    seed_runs = {m: run_flight(SCENARIOS[next(iter(SCENARIOS))], m, BASE_SEED) for m in MODES}
    print("tau stats (timing mode, scenario A, first seed): "
          f"n={len(seed_runs['timing']['taus'])}, "
          f"mean {st.mean(seed_runs['timing']['taus'])*1000:.1f} ms, "
          f"std {st.pstdev(seed_runs['timing']['taus'])*1000:.1f} ms, "
          f"max|tau| {max(abs(t) for t in seed_runs['timing']['taus'])*1000:.1f} ms" if seed_runs['timing']['taus'] else "no taus")

    print("\nwindow sensitivity (scenario A, timing mode, mean wobble rms over seeds):")
    for w_ms in (60.0, 120.0, 200.0, 300.0, 400.0):
        WINDOW_S = w_ms / 1000.0
        rows = summarize(SCENARIOS[next(iter(SCENARIOS))], range(BASE_SEED, BASE_SEED + N_SEEDS), ["timing"])["timing"]
        print(f"  window {w_ms:>5.0f} ms: wobble rms {st.mean(r['wobble_rms'] for r in rows):.3f}  "
              f"p95 {st.mean(r['wobble_p95'] for r in rows):.3f}  max {st.mean(r['wobble_max'] for r in rows):.3f}")


if __name__ == "__main__":
    main()

There might be more to the source of the issue, this was just a one-shot attempt at a specific angle to the problem. I put this implementation into practice but the bullets still wobbled.