"""Simulation behind Fig. 12 of the safety stock article.

Runs one million replenishment cycles for the worked example and reports, for
each method, the share of cycles that ended without a stockout. Needs numpy.

    python safety-stock-simulation.py
"""
import numpy as np

rng = np.random.default_rng(20261018)
N = 1_000_000
Z = 1.6448536                               # 95% cycle service level, one-sided

# Daily use after 2% scrap: three independent parents (180, 140, 80 a day)
d = 400 * 1.02
sd = (58**2 + 46**2 + 30**2) ** 0.5 * 1.02
peak = 560 * 1.02                           # peak daily use before scrap, grossed up

def service(lead_times, planned_lt, buffer):
    """Share of cycles whose demand over the realised lead time stays within the reorder point."""
    demand = rng.normal(d * lead_times, sd * np.sqrt(lead_times))
    return float((demand <= d * planned_lt + buffer).mean())

# Case 1: each delivery follows the quote, a triangle from 35 to 60 days, usually 45
a, c, b = 35, 45, 60
lt1 = rng.triangular(a, c, b, N)
mean1 = (a + b + c) / 3
sd1 = ((a*a + b*b + c*c - a*b - a*c - b*c) / 18) ** 0.5
p95_1 = np.percentile(lt1, 95)

# Case 2: 85% of deliveries follow a 35 to 55 day triangle, 15% are delayed to 60 to 80 days.
# Buffers use the measured mean and spread of these deliveries, as they would from your receipts.
late = rng.random(N) < 0.15
lt2 = np.where(late, rng.uniform(60, 80, N), rng.triangular(35, 45, 55, N))
mean2, sd2, p95_2 = lt2.mean(), lt2.std(), np.percentile(lt2, 95)

def methods(mean, sd_lt, longest, p95):
    return {
        "Combined": Z * np.sqrt(mean * sd**2 + d**2 * sd_lt**2),
        "Lead time only": Z * d * sd_lt,
        "Demand only": Z * sd * np.sqrt(mean),
        "Safety days (5 days)": 5 * d,
        "50% rule": 0.5 * d * mean,
        "Max/average": peak * longest - d * mean,
        "P95 lead time": d * (p95 - mean),
    }

m1, m2 = methods(mean1, sd1, 60, p95_1), methods(mean2, sd2, 80, p95_2)
print(f"{'Method':24}{'Case 1 buffer':>15}{'no stockout':>13}{'Case 2 buffer':>15}{'no stockout':>13}")
for name in m1:
    s1 = service(lt1, mean1, m1[name])
    s2 = service(lt2, mean2, m2[name])
    print(f"{name:24}{m1[name]:>15,.0f}{s1:>12.1%} {m2[name]:>15,.0f}{s2:>12.1%}")