"""Two-level Monte Carlo for W3 entry 044 (Flux-specific):
P(Flux Steel Works makes a first commercial GOES shipment from Louisville, Ohio by 2029-12-31).
Outer loop: uncertainty about the inputs. Inner loop: schedule outcomes given those inputs.
Time unit: months from 2029-01-01. Deadline = 12. All inputs are stated judgments over public facts."""
import numpy as np
rng = np.random.default_rng(20261111)
OUTER, INNER = 4000, 4000

def inner(p_fin, interp_commercial, p_ontime, slip_med, n=INNER):
    target = rng.uniform(0, 12, n)                                   # "2029", month unstated
    slip = np.where(rng.random(n) < p_ontime, 0.0, rng.lognormal(np.log(slip_med), 0.8, n))
    if interp_commercial:   # "commercial production in 2029" already includes qualification
        qual = rng.triangular(0, 0.5, 2, n)
    else:                   # "first coils in 2029", then customer qualification
        qual = rng.triangular(2, 5, 12, n)
    late_fin = rng.random(n) >= p_fin                                # money not in place in time
    fin_delay = np.where(late_fin, rng.uniform(12, 36, n), 0.0)
    return target + slip + qual + fin_delay

ps, all_ship = [], []
for _ in range(OUTER):
    p_fin = rng.beta(6, 4)                     # mean 0.60
    interp = rng.random() < 0.5                # what Flux's "2029" means: 50/50, unresolved in public record
    p_on = rng.beta(2, 8)                      # mean 0.20 chance of no schedule slip
    slip_med = rng.uniform(4, 15)              # median slip, months
    s = inner(p_fin, interp, p_on, slip_med)
    ps.append((s <= 12).mean()); all_ship.append(s[:250])
ps = np.array(ps); ship = np.concatenate(all_ship)

def ym(m):  # months from 2029-01 -> "Mon YYYY"
    y, mo = divmod(int(np.floor(m)), 12)
    return ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"][mo] + f" {2029+y}"

print(f"REGISTER PROBABILITY (mean, 2 dp): {ps.mean():.2f}   exact {ps.mean():.4f}")
print(f"P(YES by 2029-12-31): mean {ps.mean():.3f}   P10 {np.percentile(ps,10):.3f}   P50 {np.percentile(ps,50):.3f}   P90 {np.percentile(ps,90):.3f}")
for lab, m in [("end-2029", 12), ("end-2030", 24), ("end-2031", 36), ("end-2032", 48)]:
    print(f"  cumulative P(first commercial shipment by {lab}): {(ship<=m).mean():.3f}")
print(f"  median first commercial shipment: {ym(np.median(ship))}")

# Split by what Flux's "2029" means (separate seeded run so the main result above is unchanged)
rng = np.random.default_rng(2029)
for interp, lab in [(False, 'first coils'), (True, 'commercial sales')]:
    pp = np.array([(inner(rng.beta(6, 4), interp, rng.beta(2, 8), rng.uniform(4, 15)) <= 12).mean() for _ in range(1500)])
    print(f"  if '2029' means {lab}: mean {pp.mean():.3f}  P10 {np.percentile(pp,10):.3f}  P90 {np.percentile(pp,90):.3f}")

# Read by scripts/check_forecast_models.py (the main run's mean; the split above does not touch ps).
print(f"correlated_mean={ps.mean():.4f}")