"""R4: Nuisance alarms, coincidence fusion and operator load.
A perimeter of Z zones, each with two independent sensors (fence vibration F, radar R) and an optional
verification stage (sentry EO/IR classification). Nuisance alarms are Poisson; rates scale with wind.
Analytical: coincidence of two Poisson streams within window w has rate ~ 2*l1*l2*w (for l*w << 1).
Simulation: 10,000 hours per configuration checks the analytic rates. Seeded (seed=3).
All rates are ASSUMED illustrative inputs, not field measurements. Outputs research/out/false_alarm_fusion.json."""
import json, numpy as np
rng = np.random.default_rng(3)
Z = 40
Pd = {"F": 0.90, "R": 0.95}
lam = {"calm": {"F": 0.20, "R": 0.10}, "windy": {"F": 1.50, "R": 0.40}}   # nuisance alarms / zone / hour
Pv_detect, Pv_pass = 0.97, 0.03     # verification: P(confirm real intruder), P(nuisance passes as real)

def sim_coincidence(l1, l2, w_s, hours=10000):
    n1 = rng.poisson(l1 * hours); n2 = rng.poisson(l2 * hours)
    a = np.sort(rng.uniform(0, hours * 3600, n1)); b = np.sort(rng.uniform(0, hours * 3600, n2))
    if n1 == 0 or n2 == 0: return 0.0
    idx = np.searchsorted(b, a)
    near = np.zeros(n1, bool)
    for off in (-1, 0):
        j = np.clip(idx + off, 0, n2 - 1); near |= np.abs(b[j] - a) <= w_s
    return near.sum() / hours

out = {"assumptions": {"zones": Z, "Pd": Pd, "nuisance_rates": lam, "verify": [Pv_detect, Pv_pass]}, "rows": []}
for cond, L in lam.items():
    for w in [10, 30, 60, 120]:
        lf, lr = L["F"], L["R"]
        rows = {
            "OR": (1 - (1 - Pd["F"]) * (1 - Pd["R"]), (lf + lr) * Z),
            "AND": (Pd["F"] * Pd["R"], 2 * lf * lr * (w / 3600) * Z),
            "OR+verify": ((1 - (1 - Pd["F"]) * (1 - Pd["R"])) * Pv_detect, (lf + lr) * Z * Pv_pass),
        }
        sim_and = sim_coincidence(lf, lr, w) * Z
        for name, (pd, fa) in rows.items():
            out["rows"].append({"cond": cond, "window_s": w, "logic": name, "Pd": round(pd, 4), "alarms_per_hr": round(fa, 3),
                                "alarms_per_shift_12h": round(fa * 12, 1), "sim_check": round(sim_and, 3) if name == "AND" else None})
# operator saturation: alerts/hour an operator can assess at ~90 s each with 50% utilisation ceiling
out["operator_capacity_per_hr"] = round(3600 * 0.5 / 90, 1)
json.dump(out, open("out/false_alarm_fusion.json", "w"), indent=1)
for r in out["rows"]:
    if r["window_s"] in (30,): print(r)