import os, math, heapq
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
from pathlib import Path
warnings.filterwarnings("ignore")
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#a0aec0","#6b46c1"
here = Path.cwd(); DATA = next(p for p in [here/"data", here.parent/"data"] if (p/"ed_hourly.parquet").exists())
h = pd.read_parquet(DATA/"ed_hourly.parquet")["arrivals"]; h.index = pd.to_datetime(h.index)
lam_hour = h.groupby(h.index.hour).mean().values          # mean arrivals/hour by time-of-day (the rate lambda(t))
MEAN_SERVICE_H = 3.0                                       # ASSUMED: mean ED length-of-stay per treatment space (hours)
MU = 1.0/MEAN_SERVICE_H
print(f"{len(h):,} hourly obs | peak rate {lam_hour.max():.1f}/hr (hr {lam_hour.argmax()}), trough {lam_hour.min():.1f}/hr (hr {lam_hour.argmin()})")
print(f"service assumption: mean {MEAN_SERVICE_H} h per space -> mu = {MU:.3f}/hr per server")

33,984 hourly obs | peak rate 6.0/hr (hr 17), trough 1.6/hr (hr 4)
service assumption: mean 3.0 h per space -> mu = 0.333/hr per server

fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].plot(range(24), lam_hour, "o-", color=BLUE); ax[0].fill_between(range(24), lam_hour, alpha=.1, color=BLUE)
ax[0].set_xlabel("hour of day"); ax[0].set_ylabel("arrivals/hour  (λ)"); ax[0].set_title("Arrival rate λ(t) — real diurnal profile")
ax[1].plot(range(24), lam_hour*MEAN_SERVICE_H, "o-", color=PURP); ax[1].set_xlabel("hour of day"); ax[1].set_ylabel("offered load a = λ/μ (Erlangs)")
ax[1].set_title("Offered load through the day"); fig.tight_layout(); plt.show()

def erlang_c(c, a):
    c = int(c)
    if c <= a: return 1.0
    s = sum(a**k/math.factorial(k) for k in range(c))
    last = a**c/math.factorial(c) * (c/(c-a))
    return last/(s+last)
def servers_for_sla(lam, mu, target=0.2):
    a = lam/mu; c = int(np.ceil(a)) + 1
    while erlang_c(c, a) > target: c += 1
    return c
lam_peak = lam_hour.max(); a_peak = lam_peak/MU
cs = np.arange(int(a_peak)+1, int(a_peak)+13)
pw = [erlang_c(c, a_peak) for c in cs]
ew = [erlang_c(c, a_peak)/(c*MU-lam_peak)*60 for c in cs]
c_star = servers_for_sla(lam_peak, MU, 0.2)
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].plot(cs, pw, "o-", color=BLUE); ax[0].axhline(0.2, color=RED, ls="--", label="SLA: P(wait)=0.20")
ax[0].axvline(c_star, color=GREEN, ls=":", label=f"c* = {c_star}"); ax[0].set_xlabel("servers c"); ax[0].set_ylabel("P(wait)"); ax[0].set_title(f"Peak hour (a={a_peak:.0f} Erlangs): probability of waiting"); ax[0].legend(fontsize=8)
ax[1].plot(cs, ew, "o-", color=PURP); ax[1].axvline(c_star, color=GREEN, ls=":"); ax[1].set_xlabel("servers c"); ax[1].set_ylabel("E[wait] (min)"); ax[1].set_title("Expected wait vs servers")
fig.tight_layout(); plt.show()
print(f"peak hour needs c* = {c_star} servers for P(wait)<0.20  (offered load {a_peak:.0f} Erlangs)")

peak hour needs c* = 23 servers for P(wait)<0.20  (offered load 18 Erlangs)

targets = {0.5:"lenient (P<0.50)", 0.2:"standard (P<0.20)", 0.1:"tight (P<0.10)"}
plan = {t: np.array([servers_for_sla(lam_hour[hr], MU, t) for hr in range(24)]) for t in targets}
fig, ax = plt.subplots(figsize=(12, 4.5))
for (t, lab), c in zip(targets.items(), [GREY, BLUE, RED]):
    ax.step(range(24), plan[t], where="mid", color=c, lw=2, label=lab)
ax.plot(range(24), lam_hour*MEAN_SERVICE_H, "o--", color=PURP, alpha=.6, label="offered load a(t)")
ax.set_xlabel("hour of day"); ax.set_ylabel("servers required"); ax.set_title("Hour-by-hour staffing curve (SIPP Erlang-C)"); ax.legend(fontsize=8)
fig.tight_layout(); plt.show()
print("standard-SLA plan (servers by hour):", plan[0.2].tolist())
print(f"total server-hours/day: lenient {plan[0.5].sum()}, standard {plan[0.2].sum()}, tight {plan[0.1].sum()}")

standard-SLA plan (servers by hour): [12, 10, 9, 8, 8, 8, 9, 12, 17, 20, 22, 22, 23, 22, 22, 22, 23, 23, 23, 22, 21, 20, 18, 14]
total server-hours/day: lenient 358, standard 410, tight 444

def simulate(lam_hour, c_hour, mu, days=120, seed=0):
    rng = np.random.default_rng(seed)
    H = []
    for d in range(days):
        for hr in range(24):
            for u in rng.uniform(0, 1, rng.poisson(lam_hour[hr])):
                H.append((d*24+hr+u, 0, d*24+hr+u))     # (event_time, type: 0=arrival, arrival_time)
    heapq.heapify(H); busy = 0; queue = []; waits = []
    while H:
        t, typ, at = heapq.heappop(H)
        c = c_hour[int(t) % 24]
        if typ == 0:                                     # arrival
            if busy < c:
                busy += 1; heapq.heappush(H, (t + rng.exponential(1/mu), 1, t)); waits.append((int(t)%24, 0.0))
            else:
                queue.append(t)
        else:                                            # service completion
            busy -= 1
            if queue and busy < c_hour[int(t) % 24]:
                a0 = queue.pop(0); busy += 1; heapq.heappush(H, (t + rng.exponential(1/mu), 1, a0)); waits.append((int(a0)%24, t - a0))
    return pd.DataFrame(waits, columns=["hour", "wait"])

sim = simulate(lam_hour, plan[0.2], MU, days=120)
by = sim.groupby("hour")["wait"]
realized_pw = by.apply(lambda w: (w > 1e-9).mean()); realized_ew = by.mean()*60
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].step(range(24), realized_pw.reindex(range(24)).values, where="mid", color=BLUE, lw=2, label="simulated P(wait)")
ax[0].axhline(0.2, color=RED, ls="--", label="target 0.20"); ax[0].set_xlabel("hour"); ax[0].set_ylabel("P(wait)"); ax[0].set_title("Realized wait probability under the SIPP plan"); ax[0].legend(fontsize=8)
ax[1].step(range(24), realized_ew.reindex(range(24)).values, where="mid", color=PURP, lw=2); ax[1].set_xlabel("hour"); ax[1].set_ylabel("E[wait] (min)"); ax[1].set_title("Realized mean wait by hour")
fig.tight_layout(); plt.show()
# the hour-by-hour extremes are what this panel is about, so report them rather than leaving
# them to be read off the chart
_wh = int(realized_pw.idxmax()); _wm = int(realized_ew.idxmax())
print(f"worst hour {_wh:02d}:00 -> P(wait) {realized_pw.max():.2f}; longest mean wait "
      f"{realized_ew.max():.0f} min at {_wm:02d}:00")
print(f"overall simulated P(wait) {(sim.wait>1e-9).mean():.2f} | mean wait {sim.wait.mean()*60:.1f} min | "
      f"hours breaching target: {(realized_pw>0.22).sum()}")

worst hour 01:00 -> P(wait) 0.73; longest mean wait 110 min at 01:00
overall simulated P(wait) 0.20 | mean wait 17.7 min | hours breaching target: 11

import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS
hh = pd.DataFrame({"y": h.values, "hour": h.index.hour, "dow": h.index.dayofweek})
def model(hour, dow, y=None):
    a = numpyro.sample("a", dist.Normal(np.log(h.mean()+1e-3), 1.0))
    b_h = numpyro.sample("b_h", dist.Normal(0,1).expand([24]))
    b_d = numpyro.sample("b_d", dist.Normal(0,0.5).expand([7]))
    numpyro.sample("obs", dist.Poisson(jnp.exp(a + b_h[hour] + b_d[dow])), obs=y)
mcmc = MCMC(NUTS(model), num_warmup=400, num_samples=500, num_chains=1, progress_bar=False)
mcmc.run(jax.random.PRNGKey(0), jnp.array(hh.hour.values), jnp.array(hh.dow.values), y=jnp.array(hh.y.values, float))
post = mcmc.get_samples()
# posterior lambda(hour) for an average weekday, and required servers per draw
lam_draws = np.exp(post["a"][:,None] + post["b_h"] + post["b_d"].mean(1)[:,None])   # (S, 24)
S = lam_draws.shape[0]
req = np.empty((S, 24), int)
for s in range(S):
    for hr in range(24):
        req[s,hr] = servers_for_sla(lam_draws[s,hr], MU, 0.2)
lo, med, hi = np.percentile(req, [5,50,95], axis=0)
fig, ax = plt.subplots(figsize=(12, 4.5))
ax.fill_between(range(24), lo, hi, color=BLUE, alpha=.2, label="90% credible band")
ax.step(range(24), med, where="mid", color=BLUE, lw=2, label="posterior median requirement")
ax.step(range(24), plan[0.2], where="mid", color=GREY, lw=1.5, ls="--", label="point-estimate plan (§3)")
ax.set_xlabel("hour of day"); ax.set_ylabel("servers required (P(wait)<0.20)"); ax.set_title("Required servers with Bayesian uncertainty"); ax.legend(fontsize=8)
fig.tight_layout(); plt.show()
print(f"server-hours/day: point plan {plan[0.2].sum():.0f} | posterior-median {med.sum():.0f} | robust (95th pct) {hi.sum():.1f}")
print(f"the 95th-percentile plan adds {hi.sum()-med.sum():.1f} server-hours as an uncertainty margin")

server-hours/day: point plan 410 | posterior-median 409 | robust (95th pct) 417.0
the 95th-percentile plan adds 8.0 server-hours as an uncertainty margin

quantity	meaning	at the ED peak
$\lambda$	arrival rate (patients/hour)	~6 at 5 pm
$\mu$	service rate per server (patients/hour)	$1/3\approx0.33$ (3 h mean stay)
$a=\lambda/\mu$	offered load, in Erlangs	18

Capacity planning for an emergency department — a queueing view¶

Machine Learning in Operations Research · how many servers to hit a service target?¶

1 · The arrival process — real, and strongly time-varying¶

2 · Classical queueing — the M/M/c model and Erlang-C¶

What a queue is¶

Kendall's notation: the M/M/c queue¶

The three numbers that drive everything¶

Erlang-C — the probability an arrival has to wait¶

From P(wait) to the expected wait — and Little's law¶

Why the calculation beats intuition¶

Assumptions — and why we simulate later¶

3 · Time-varying staffing — Erlang-C hour by hour (SIPP)¶

4 · Does the plan hold up? — a discrete-event simulation¶

5 · Bayesian arrival rates → staffing under uncertainty¶

6 · Takeaways¶