"""Monte Carlo: P(FCC adopts final rules opening at least part of the 482 MHz D2D proposal by a deadline).
Clock starts at NPRM adoption, 2026-10-29. All durations in days."""
import numpy as np, datetime as dt, sys
rng=np.random.default_rng(42); N=400_000
START=dt.date(2026,10,29)
def ln(med,s,n): return med*np.exp(rng.normal(0,s,n))

def run(p_stall_a=2,p_stall_b=8, post_med=165, post_s=0.35, split_p=0.6, shutdown_p=0.30, verbose=False):
    fr = np.clip(ln(30,0.45,N),10,120)                       # adoption -> Federal Register publication
    cycle = np.full(N,60.0)                                   # comments +30, replies +60 after FR
    ext = rng.random(N)<0.40; cycle += ext*rng.uniform(15,30,N)   # extension request granted (contested item)
    p_stall = rng.beta(p_stall_a,p_stall_b,N)                 # parameter uncertainty on stall risk
    stall = rng.random(N)<p_stall
    post = np.where(stall, ln(540,0.6,N), ln(post_med,post_s,N))  # replies -> order: priority track vs languishing NPRM
    nosplit = rng.random(N)>split_p                            # FCC does NOT split off contested 2.5/3.45 GHz bands
    post += nosplit*rng.uniform(30,120,N)                     # extra interference analysis / negotiation
    sd = rng.random(N)<shutdown_p                             # lapse in appropriations during the window
    post += sd*np.clip(ln(20,0.6,N),3,70)
    snap = rng.uniform(0,30,N)                                # wait for next open meeting / circulation vote
    total = fr+cycle+post+snap
    d = np.array([ (START+dt.timedelta(days=0)).toordinal() ])[0] + total
    def p(by): return np.mean(d <= dt.date(*by).toordinal())
    out = dict(oct27=p((2027,10,31)), dec27=p((2027,12,31)), jun28=p((2028,6,30)),
               med=dt.date.fromordinal(int(np.median(d))), p10=dt.date.fromordinal(int(np.percentile(d,10))),
               p90=dt.date.fromordinal(int(np.percentile(d,90))), stall_share=stall.mean())
    return out

base=run()
print("BASE  by Oct 31 2027: %.1f%%  by Dec 31 2027: %.1f%%  by Jun 30 2028: %.1f%%" % (base['oct27']*100, base['dec27']*100, base['jun28']*100))
print("      median %s   10th pct %s   90th pct %s   stall share %.2f" % (base['med'],base['p10'],base['p90'],base['stall_share']))
print(f"correlated_mean={base['dec27']:.3f}")   # scripts/check_forecast_models.py reads this line
print("\nSENSITIVITY (P by Dec 31 2027)")
for name,kw in [("stall risk ~10%",dict(p_stall_a=1,p_stall_b=9)),("stall risk ~30%",dict(p_stall_a=3,p_stall_b=7)),
                ("post-reply median 4.5 mo",dict(post_med=135)),("post-reply median 7.5 mo",dict(post_med=225)),
                ("FCC always splits bands",dict(split_p=1.0)),("FCC never splits",dict(split_p=0.0)),
                ("no shutdown",dict(shutdown_p=0.0)),("shutdown 60%",dict(shutdown_p=0.6))]:
    r=run(**kw); print("  %-26s %.1f%%   (Oct 31: %.1f%%)" % (name, r['dec27']*100, r['oct27']*100))