"""R6 — Cost of a 24/7 guard post vs an autonomous post: Monte Carlo + sensitivity.
Wage anchor: BLS OOH median hourly wage for security guards, May 2024 = $18.46 (verified source).
All other distributions are ASSUMPTIONS stated in the article (triangular: low, mode, high).
FTE factor = 8760 / productive hours per FTE. 100,000 draws, seed=5. Outputs research/out/guard_cost.json."""
import json, numpy as np
rng = np.random.default_rng(5)
n = 100_000
tri = lambda lo, mo, hi: rng.triangular(lo, mo, hi, n)
wage = tri(16.0, 18.46, 28.0)                 # $/h paid to the guard
burden = tri(1.18, 1.28, 1.45)                # payroll taxes, insurance, benefits multiplier
markup = tri(1.20, 1.35, 1.60)                # contractor overhead + margin
prod_h = tri(1800, 1890, 2000)                # productive hours per FTE-year (after PTO, holidays, sick, training)
ot_share = tri(0.00, 0.05, 0.15)              # share of hours paid at 1.5x
supervision = tri(0.05, 0.08, 0.12)           # supervisor overhead as share of labor
bill = wage * burden * markup * (1 + 0.5 * ot_share) * (1 + supervision)
post_cost = bill * 8760
fte = 8760 / prod_h
# autonomous alternative per post replaced
units_per_post = tri(0.5, 0.67, 1.0)
rate = tri(8.5, 9.80, 11.0)                   # $/sentry-hour (Tencer list price as mode)
hive_share = tri(0.06, 0.07, 0.10)            # Roost cost as share of sentry cost
monitor = tri(0.0, 2.0, 6.0)                  # $/h incremental monitoring per post replaced
auto_cost = 8760 * (units_per_post * rate * (1 + hive_share) + monitor)
saving = 1 - auto_cost / post_cost
pct = lambda a: {p: round(float(np.percentile(a, p)), 3) for p in (5, 25, 50, 75, 95)}
out = {"n": n, "bill_rate": pct(bill), "post_cost": pct(post_cost), "fte": pct(fte), "auto_cost": pct(auto_cost), "saving": pct(saving),
       "p_saving_gt_50": round(float((saving > 0.5).mean()), 4), "p_saving_lt_0": round(float((saving < 0).mean()), 4),
       "fte_naive": round(8760 / 2080, 3)}
# tornado: swing each input between its 5th and 95th pct with others at median, on saving
inputs = {"wage": wage, "burden": burden, "markup": markup, "ot_share": ot_share, "supervision": supervision,
          "units_per_post": units_per_post, "rate": rate, "hive_share": hive_share, "monitor": monitor}
med = {k: float(np.median(v)) for k, v in inputs.items()}
def sav(d):
    b = d["wage"] * d["burden"] * d["markup"] * (1 + 0.5 * d["ot_share"]) * (1 + d["supervision"])
    a = d["units_per_post"] * d["rate"] * (1 + d["hive_share"]) + d["monitor"]
    return 1 - a / b
torn = []
for k, v in inputs.items():
    lo, hi = np.percentile(v, 5), np.percentile(v, 95)
    torn.append({"input": k, "lo": round(float(sav({**med, k: lo})), 4), "hi": round(float(sav({**med, k: hi})), 4)})
torn.sort(key=lambda r: -abs(r["hi"] - r["lo"]))
out["tornado"] = torn; out["saving_at_median"] = round(float(sav(med)), 4)
# histogram of saving
h, edges = np.histogram(saving, bins=30, range=(0.3, 0.95))
out["hist"] = {"counts": h.tolist(), "edges": [round(float(e), 3) for e in edges]}
json.dump(out, open("out/guard_cost.json", "w"), indent=1)
print({k: out[k] for k in ("bill_rate", "post_cost", "fte", "auto_cost", "saving", "p_saving_gt_50", "p_saving_lt_0", "saving_at_median")}); print(torn[:4])