"""R2 — Drone-in-a-box fleet sizing vs perch-and-hold, with weather.
Monte Carlo over one simulated year at 1-minute resolution.
Wind: hourly AR(1) process mapped to a Weibull(k=2, c=6 m/s) marginal (persistence phi=0.9/h).
Rain: two-state Markov chain, ~8% of hours wet, mean wet spell 3 h.
DiaB: 35 min sortie (4 min each way transit => 27 min on station), 55 min recharge; flies only if wind < 12 m/s and dry.
Perch: on watch whenever perched; must re-perch every 4 h (battery/rotation), transit needs wind < 12 m/s;
       if transit is weather-blocked it stays perched (tolerates 15.3 m/s sustained); above that it is blind.
Two DiaB policies: naive fixed-interval staggering (with 5-min weather retry) and relief scheduling.
Seeded (seed=7) so results are reproducible. Outputs research/out/dib_fleet.json."""
import json, math
import numpy as np
from scipy.stats import norm, weibull_min

rng = np.random.default_rng(7)
H = 24 * 365
phi = 0.9
z = np.zeros(H)
for t in range(1, H):
    z[t] = phi * z[t - 1] + math.sqrt(1 - phi ** 2) * rng.standard_normal()
wind = weibull_min.ppf(norm.cdf(z), 2, scale=6.0)            # m/s, hourly mean
wet = np.zeros(H, bool)
p_wd, p_dw = 1 / 3, (0.08 / 0.92) / 3                            # leave wet after ~3 h; stationary 8%
for t in range(1, H):
    wet[t] = (rng.random() > p_wd) if wet[t - 1] else (rng.random() < p_dw)
flyable_h = (wind < 12) & ~wet
perch_ok_h = wind < 15.3
M = H * 60
flyable = np.repeat(flyable_h, 60); perch_ok = np.repeat(perch_ok_h, 60)

SORTIE, TRANSIT, RECHARGE = 35, 4, 55
def dib_coverage(n):
    covered = np.zeros(M, bool)
    period = SORTIE + RECHARGE
    offsets = [round(i * period / n) for i in range(n)]
    for off in offsets:
        t = off
        while t < M:
            if flyable[t]:
                end = min(t + SORTIE, M)
                covered[t + TRANSIT:end - TRANSIT] = True
                t += period
            else:
                t += 5          # retry in 5 minutes
    return covered.mean()

def perch_coverage(n_units=1):
    covered = np.zeros(M, bool)
    t = 0
    while t < M:
        # perched for up to 240 min, then needs a flyable minute to relocate (≈2 min blind)
        end = min(t + 240, M)
        covered[t:end] = perch_ok[t:end]
        t = end
        while t < M and not flyable[t]:
            covered[t] = perch_ok[t]
            t += 1
        t += 2
    return covered.mean()

res = {"weather": {"flyable_frac": round(float(flyable_h.mean()), 4), "wet_frac": round(float(wet.mean()), 4),
                   "wind_mean": round(float(wind.mean()), 2), "frac_wind_over_12": round(float((wind >= 12).mean()), 4),
                   "frac_wind_over_15_3": round(float((wind >= 15.3).mean()), 4)}}
def dib_relief(n):
    """Relief scheduling: launch the next charged aircraft so it arrives on station as the current one departs."""
    covered = np.zeros(M, bool)
    ready = [0] * n                     # minute each dock's aircraft is charged
    station_until = -1                  # minute the current on-station aircraft leaves
    for t in range(M):
        if station_until - t <= TRANSIT and flyable[t]:
            i = min(range(n), key=lambda k: ready[k])
            if ready[i] <= t:
                start = max(t + TRANSIT, station_until)
                end = t + SORTIE - TRANSIT
                covered[start:min(end, M)] = True
                station_until = end
                ready[i] = t + SORTIE + RECHARGE
    return covered.mean()

res["dib_naive"] = [{"docks": n, "coverage": round(float(dib_coverage(n)), 4)} for n in range(1, 8)]
res["dib"] = [{"docks": n, "coverage": round(float(dib_relief(n)), 4)} for n in range(1, 8)]
res["dib_ideal"] = [{"docks": n, "coverage": round(min(1, n * (SORTIE - 2 * TRANSIT) / (SORTIE + RECHARGE)), 4)} for n in range(1, 8)]
res["perch"] = {"units": 1, "coverage": round(float(perch_coverage()), 4)}
res["n_min_continuous"] = math.ceil((SORTIE + RECHARGE) / (SORTIE - 2 * TRANSIT))
json.dump(res, open("out/dib_fleet.json", "w"), indent=1)
print(json.dumps(res, indent=1))