Made of little things.

#!/usr/bin/env python3
"""Compare an original-game log with a port replay, frame by frame.
    python3 tools/sim/compare.py original.txt port.txt [every]
"""
import math, sys

def load(p):
    rows = []
    for line in open(p, newline=''):
        for part in line.replace('\r', '\n').split('\n'):
            r = part.split()
            if len(r) >= 33:
                rows.append([float(x) for x in r[:33]])
    return rows


def pair_by_sim(orig, port, tol=0.02):
    """Pair each original row with the port row nearest in simTime (within tol)."""
    out, j = [], 0
    for o in orig:
        while j + 1 < len(port) and abs(port[j + 1][1] - o[1]) <= abs(port[j][1] - o[1]):
            j += 1
        if port and abs(port[j][1] - o[1]) <= tol:
            out.append((o, port[j]))
    return out

orig, port = load(sys.argv[1]), load(sys.argv[2])
pairs = pair_by_sim(orig, port)
every = int(sys.argv[3]) if len(sys.argv) > 3 else 10
print(f"{'sim':>6} | {'orig pos':>22} {'spd':>5} {'tilt':>5} {'spr':>5} | {'port pos':>22} {'spd':>5} {'tilt':>5} {'spr':>5} | {'dpos':>6} keys")
def tilt(r):  # angle of the body's up axis from world up, from axis-angle
    ax, ay, az, ang = r[5], r[6], r[7], math.radians(r[8])
    n = math.sqrt(ax*ax + ay*ay + az*az) or 1
    ax, ay, az = ax/n, ay/n, az/n
    c, s = math.cos(ang), math.sin(ang)
    upz = c + az*az*(1-c)
    return math.degrees(math.acos(max(-1, min(1, upz))))
for o, p in pairs[::every]:
    d = math.dist(o[2:5], p[2:5])
    so, sp = math.hypot(o[9], o[10]), math.hypot(p[9], p[10])
    keys = ''.join(str(int(k)) for k in o[27:32])
    print(f"{o[1]:6.2f} | {o[2]:7.0f}{o[3]:8.0f}{o[4]:7.1f} {so:5.0f} {tilt(o):5.1f} {o[16]:5.2f} | {p[2]:7.0f}{p[3]:8.0f}{p[4]:7.1f} {sp:5.0f} {tilt(p):5.1f} {p[16]:5.2f} | {d:6.1f} {keys}")

# Summary: position error and heading error over the run.
def heading(r):
    ax, ay, az, ang = r[5], r[6], r[7], math.radians(r[8])
    n = math.sqrt(ax*ax + ay*ay + az*az) or 1
    ax, ay, az = ax/n, ay/n, az/n
    c, s = math.cos(ang), math.sin(ang)
    fx = ax*ay*(1-c) - az*s      # rotated (0,1,0)
    fy = c + ay*ay*(1-c)
    return math.degrees(math.atan2(-fx, fy))
errs = [math.dist(o[2:5], p[2:5]) for o, p in pairs]
herr = [abs((heading(o) - heading(p) + 180) % 360 - 180) for o, p in pairs]
first_div = next((o[1] for (o, p), e in zip(pairs, errs) if e > 50), None)
print(f"SUMMARY frames={len(pairs)} meanPosErr={sum(errs)/len(errs):.1f} maxPosErr={max(errs):.1f} meanHdgErr={sum(herr)/len(herr):.1f} firstDiverge(>50u)@sim={first_div}")