"""How much a simulated crowd looks like the real one. Both files: frame, id, x, y (0.4 s steps, metres).
Three distributions, each compared by the 1-Wasserstein distance between simulated and real:
  speed        every step's speed (m/s)
  nearest      distance from each person to the nearest other person present in the same frame (m)
  straightness each person's straight-line distance from first to last point over path length (0-1)
Also reported: share of frames in which some pair is closer than 0.3 m (contact).

  python3 peds_realism.py REAL.txt SIM.txt [SIM2.txt ...]   prints one JSON line per simulation
"""
import json, sys
import numpy as np


def load(path):
    d = np.loadtxt(path)
    frames = np.unique(d[:, 0]); step = np.min(np.diff(frames)) if len(frames) > 1 else 1
    d[:, 0] = np.round(d[:, 0] / step)
    return d


def w1(a, b):
    a, b = np.sort(np.asarray(a, float)), np.sort(np.asarray(b, float))
    if len(a) == 0 or len(b) == 0: return float("nan")
    q = np.linspace(0, 1, 2001)
    return float(np.mean(np.abs(np.quantile(a, q) - np.quantile(b, q))))


def stats(d):
    speed, nearest, straight, contact_frames, frames = [], [], [], 0, 0
    for pid in np.unique(d[:, 1]):
        t = d[d[:, 1] == pid]; t = t[np.argsort(t[:, 0])]
        if len(t) < 2: continue
        step = np.diff(t[:, 0]); seg = np.linalg.norm(np.diff(t[:, 2:4], axis=0), axis=1)
        ok = step == 1; speed += list(seg[ok] / 0.4)
        L = seg.sum(); chord = np.linalg.norm(t[-1, 2:4] - t[0, 2:4])
        if L > 1.0: straight.append(chord / L)
    for f in np.unique(d[:, 0]):
        p = d[d[:, 0] == f][:, 2:4]
        if len(p) < 2: continue
        frames += 1
        dist = np.linalg.norm(p[:, None, :] - p[None, :, :], axis=2) + np.eye(len(p)) * 1e9
        m = dist.min(axis=1); nearest += list(m)
        contact_frames += int(m.min() < 0.3)
    return {"speed": speed, "nearest": nearest, "straight": straight, "contact": contact_frames / max(1, frames)}


if __name__ == "__main__":
    real = stats(load(sys.argv[1]))
    for sim_path in sys.argv[2:]:
        s = stats(load(sim_path))
        out = {k: round(w1(real[k], s[k]), 4) for k in ("speed", "nearest", "straight")}
        out.update({"contact_real": round(real["contact"], 3), "contact_sim": round(s["contact"], 3), "sim": sim_path})
        print(json.dumps(out))