"""R5 — Timely detection: probability of interruption (EASI-style model, after Garcia / Sandia).
Adversary path = ordered elements, each with detection probability p_i, delay t_i (s) and std of delay.
Detection at element i is 'timely' if the remaining path delay (from the detection point; we assume detection
at the start of the element, as in EASI's conservative 'B' location choice for barriers) exceeds response time RFT.
P_I = sum_i p_i * prod_{j<i}(1-p_j) * P(remaining_delay_i - RFT > 0), with normal approximations for both.
Scenarios are ILLUSTRATIVE: a remote substation with fence, yard and control house.
Outputs research/out/probability_of_interruption.json."""
import json, math
from scipy.stats import norm

def p_interrupt(path, rft_mean, rft_sd):
    total, miss = 0.0, 1.0
    for i, (name, p, t, sd) in enumerate(path):
        rem = sum(e[2] for e in path[i:]); rem_var = sum(e[3] ** 2 for e in path[i:])
        z = (rem - rft_mean) / math.sqrt(rem_var + rft_sd ** 2)
        total += miss * p * norm.cdf(z); miss *= (1 - p)
    return total

base = [  # name, p_detect, delay_s, sd_s
    ("approach 300 m open ground", 0.05, 150, 45),
    ("climb fence", 0.80, 20, 6),
    ("cross yard 60 m", 0.30, 40, 12),
    ("defeat control-house door", 0.50, 180, 54),
    ("disable equipment", 0.10, 300, 90),
]
outer = [("approach 300 m open ground", 0.85, 150, 45)] + base[1:]     # wide-area radar + perched EO/IR on approach
outer_verify_delay = [("approach 300 m open ground", 0.85, 150, 45), ("climb fence", 0.80, 20, 6), ("cross yard 60 m", 0.6, 40, 12)] + base[3:]
scen = {"fence-centric": base, "approach detection": outer, "approach + yard sentry": outer_verify_delay}
rfts = [180, 300, 480, 600, 900, 1200, 1800]
out = {"scenarios": {k: [{"rft_s": r, "PI": round(p_interrupt(v, r, 0.3 * r), 4)} for r in rfts] for k, v in scen.items()},
       "paths": {k: [list(e) for e in v] for k, v in scen.items()}, "total_delay_s": {k: sum(e[2] for e in v) for k, v in scen.items()}}
json.dump(out, open("out/probability_of_interruption.json", "w"), indent=1)
for k, v in out["scenarios"].items(): print(k, [(x["rft_s"], x["PI"]) for x in v])