Wiki:Packs/Probabilistic Staffing/Modules
From WFM Labs
Modules for the Probabilistic Staffing pack (CP-WFM-018). Save each block under the filename in its heading and upload it as project knowledge with the rest of the pack. Dependencies: numpy, scipy, pandas, matplotlib, pyyaml.
The pack page carries the method, the parameters, the intake schema and the daily protocol; this page carries only the code those describe. Read Wiki:Packs/Probabilistic Staffing first.
Block 4 — 03-simulator.md
# 03 — Simulator
The model itself: distributions, the Erlang helpers, the two-pass demand calculation, the
supply pipeline, variance attribution and charts. The posterior-update and scoring helpers it
imports live in `03b-updating.md`; the intake and reconstruction layer is `05b-intake.md`, and the cycle that drives
them is `05c-cycle.md`.
Save the fenced block below as `mc_staffing.py`. Dependencies: `numpy`, `scipy`, `pandas`,
`matplotlib`, `pyyaml`.
Run `python mc_staffing.py` with no arguments for a self-test on demonstration parameters —
useful for checking the environment before real data exists. With `params.yaml` present it
uses that instead.
---
## `mc_staffing.py`
```python
"""Probabilistic staffing model: two-sided Monte Carlo over a weekly horizon.
Demand -> required productive hours (workload / occupancy)
Supply -> available productive hours (heads x ramp x hours x (1 - shrink))
Answer -> P(supply >= demand), shortfall distribution, variance attribution.
Shrinkage is applied ONCE, on the supply side. Demand is never grossed up for it.
"""
from __future__ import annotations
import json
import os
import warnings
from dataclasses import dataclass, field
import numpy as np
import pandas as pd
from scipy import optimize, stats
SEED = 20260919
Z90 = 1.6448536269514722
# --------------------------------------------------------------------------
# distribution helpers
# --------------------------------------------------------------------------
def lognormal_from_ci(lo: float, hi: float, conf: float = 0.90):
"""Lognormal (mu, sigma) on the log scale matching a central interval."""
if lo <= 0 or hi <= lo:
raise ValueError(f"bad lognormal interval: ({lo}, {hi})")
z = stats.norm.ppf(0.5 + conf / 2)
mu = (np.log(lo) + np.log(hi)) / 2
sigma = (np.log(hi) - np.log(lo)) / (2 * z)
return mu, sigma
def beta_from_mean_ci(mean: float, lo: float, hi: float, conf: float = 0.90):
"""Beta(a, b) with the given mean whose central interval width matches (lo, hi).
Parameterized as Beta(mean*k, (1-mean)*k) and solved for the concentration k.
The search is DELIBERATELY restricted to k > 1/min(mean, 1-mean), the region
where both shape parameters exceed 1 and the density is unimodal. Interval
width is monotone decreasing in k only inside that region; below it the width
rises, peaks and collapses again as the density piles up on the boundaries.
An unrestricted search therefore has two roots, and for a small mean it finds
the wrong one: a degenerate U-shaped Beta that reproduces the requested mean
while drawing almost every sample at 0 with rare samples at 1. That failure is
silent -- the mean is right and every quantile is wrong.
"""
if not 0 < mean < 1:
raise ValueError(f"beta mean must be in (0,1), got {mean}")
if not lo < mean < hi:
raise ValueError(f"mean {mean} must lie inside its interval ({lo}, {hi})")
target = hi - lo
ql, qh = (1 - conf) / 2, 0.5 + conf / 2
def width(k):
a, b = mean * k, (1 - mean) * k
return stats.beta.ppf(qh, a, b) - stats.beta.ppf(ql, a, b)
k_min = (1.0 / min(mean, 1 - mean)) * (1 + 1e-9) # a > 1 and b > 1
k_max = 1e9
if width(k_min) <= target:
# No unimodal Beta with this mean is that wide. Use the widest one and say so.
warnings.warn(
f"interval ({lo}, {hi}) is wider than any unimodal Beta with mean {mean} "
f"can produce; using the widest available (k={k_min:.3g}). Reconsider "
f"whether a Beta is the right family here.", RuntimeWarning)
k = k_min
else:
k = np.exp(optimize.brentq(lambda lk: width(np.exp(lk)) - target,
np.log(k_min), np.log(k_max)))
return mean * k, (1 - mean) * k
def check_beta(mean: float, lo: float, hi: float, conf: float = 0.90, tol: float = 0.15):
"""Round-trip a Beta spec and confirm the realized mean and interval match.
Run this over every rate in the config before trusting a single result. A
distribution can reproduce a requested mean while being the wrong shape
entirely; only the quantiles catch it.
"""
a, b = beta_from_mean_ci(mean, lo, hi, conf)
ql, qh = (1 - conf) / 2, 0.5 + conf / 2
got_lo, got_hi = stats.beta.ppf(ql, a, b), stats.beta.ppf(qh, a, b)
got_mean = a / (a + b)
ok = (abs(got_mean - mean) <= tol * mean
and abs((got_hi - got_lo) - (hi - lo)) <= tol * (hi - lo)
and a > 1 and b > 1)
return {"a": a, "b": b, "mean": got_mean, "lo": got_lo, "hi": got_hi, "ok": bool(ok)}
def pert_ab(lo: float, mode: float, hi: float, lam: float = 4.0):
"""Beta-PERT shape parameters for a (min, most likely, max) elicitation."""
if hi <= lo:
raise ValueError(f"bad PERT range: ({lo}, {mode}, {hi})")
a = 1 + lam * (mode - lo) / (hi - lo)
b = 1 + lam * (hi - mode) / (hi - lo)
return a, b
def draw_beta(rng, spec, size):
"""spec: {'mean':, 'ci':[lo,hi]} or {'a':, 'b':} (a posterior)."""
if "a" in spec:
a, b = spec["a"], spec["b"]
else:
a, b = beta_from_mean_ci(spec["mean"], *spec["ci"])
return rng.beta(a, b, size=size)
def draw_lognormal(rng, spec, size):
"""spec: {'ci':[lo,hi]} or {'mu':, 'sigma':} (a posterior)."""
if "mu" in spec:
mu, sigma = spec["mu"], spec["sigma"]
else:
mu, sigma = lognormal_from_ci(*spec["ci"])
return rng.lognormal(mu, sigma, size=size)
def _draw_curve(rng, spec_list, N, monotone: bool = False) -> np.ndarray:
"""A tenure curve as (N, L). Entries may be fixed floats or Beta specs.
Ramp and early attrition were fixed arrays in an earlier version while the
parameter register described both as distributions. That is the wrong way
round: `01-method.md` calls the ramp multiplier the place most plans lie to
themselves, and a curve carrying no uncertainty cannot express the failure
mode `06-interpretation.md` warns about -- a ramp reaching 1.0 sooner than
the tenure data supports.
With `monotone`, each draw is made non-decreasing across tenure. Sampling
each week independently can put week 3 below week 2, which would mean agents
getting less proficient with practice.
"""
cols = []
for spec in spec_list:
if isinstance(spec, dict):
cols.append(draw_beta(rng, spec, N))
else:
cols.append(np.full(N, float(spec)))
out = np.column_stack(cols)
return np.maximum.accumulate(out, axis=1) if monotone else out
# --------------------------------------------------------------------------
# occupancy as a function of load — the alternative to a fixed assumption
# --------------------------------------------------------------------------
def erlang_b(n: int, a: float) -> float:
"""Erlang B by the stable recurrence. No factorials, so no overflow at
realistic group sizes -- the direct formula fails above roughly n=170."""
inv = 1.0
for i in range(1, n + 1):
inv = 1.0 + inv * i / a
return 1.0 / inv
def erlang_c(n: int, a: float) -> float:
if a >= n:
return 1.0
b = erlang_b(n, a)
rho = a / n
return b / (1.0 - rho * (1.0 - b))
def service_level(n: int, a: float, asa: float, aht: float) -> float:
if a >= n:
return 0.0
return 1.0 - erlang_c(n, a) * np.exp(-(n - a) * asa / aht)
def agents_for_sl(a: float, target: float, asa: float, aht: float,
cap_mult: float = 3.0) -> int:
n = int(np.floor(a)) + 1
limit = max(int(a * cap_mult) + 10, n + 10)
while n < limit:
if service_level(n, a, asa, aht) >= target:
return n
n += 1
return limit
def erlang_a(n: int, a: float, aht: float, patience: float, jmax: int = 4000) -> dict:
"""M/M/n+M — Erlang C plus impatience. Returns abandonment and delay rates.
Erlang C assumes infinite patience: every caller waits forever, so once
staffing falls below offered load the queue grows without bound and the
model reports total collapse. Real callers leave, and that abandonment is
what stabilizes the system. Using Erlang C to describe an understaffed
operation therefore overstates what happens -- the failure is severe but it
is finite, and it shows up as lost contacts rather than as infinite waits.
Computed from the birth-death chain in log space so that loads of several
hundred erlangs do not overflow. `patience` is mean time to abandon, in the
same units as `aht`.
"""
mu, theta = 1.0 / aht, 1.0 / patience
lam = a * mu
logr = [0.0]
for k in range(1, n + 1):
logr.append(logr[-1] + np.log(a) - np.log(k))
for j in range(1, jmax + 1):
logr.append(logr[-1] + np.log(lam) - np.log(n * mu + j * theta))
r = np.array(logr)
p = np.exp(r - r.max())
p /= p.sum()
pq = p[n:]
eq = float((np.arange(len(pq)) * pq).sum())
return {"p_abandon": float(theta * eq / lam), "p_delay": float(pq.sum()),
"mean_queue": eq}
def occupancy_curve(target_sl: float, asa: float, aht: float,
a_min: float = 0.5, a_max: float = 5000.0, points: int = 60):
"""Tabulate achievable occupancy against offered load, at a fixed service level.
This is the honest answer to the model's weakest joint. Occupancy is not a
property of an operation, it is a consequence of how much load is being
pooled: at 450s AHT and 80/20, roughly 83% at 30 erlangs and 97% at 480.
Holding it fixed while volume swings across the simulated range understates
the requirement in low draws and overstates it in high ones, compressing the
very spread the model exists to show.
Returns (a_grid, occ_grid) for interpolation. Built once per channel at the
channel's median AHT -- the curve's shape is driven by load, and re-deriving
it per draw would cost far more than the second-order accuracy it buys.
"""
a_grid = np.geomspace(a_min, a_max, points)
occ = np.array([a / agents_for_sl(a, target_sl, asa, aht) for a in a_grid])
return a_grid, occ
def _occupancy_for(spec, rng, N, W, workload_hours, aht_median):
"""Resolve an occupancy spec to an (N, W) array.
Accepts either a Beta spec (fixed occupancy, sampled) or
{kind: erlang_curve, target_sl, asa_seconds, operating_hours_per_week,
spread}, in which case occupancy tracks the load in every draw.
"""
if not (isinstance(spec, dict) and spec.get("kind") == "erlang_curve"):
return draw_beta(rng, spec, (N, 1)) * np.ones((1, W))
oh = float(spec["operating_hours_per_week"])
a_grid, occ_grid = occupancy_curve(
float(spec.get("target_sl", 0.80)), float(spec.get("asa_seconds", 20.0)),
aht_median)
load = np.clip(workload_hours / oh, a_grid[0], a_grid[-1]) # erlangs
occ = np.interp(load, a_grid, occ_grid)
# Residual uncertainty about the curve itself: the operation does not sit
# exactly on its theoretical curve, and pretending otherwise would replace
# one false certainty with another.
spread = float(spec.get("spread", 0.03))
if spread > 0:
occ = occ * rng.lognormal(0.0, spread, size=(N, 1))
return np.clip(occ, 0.30, 0.98)
# --------------------------------------------------------------------------
# results container
# --------------------------------------------------------------------------
@dataclass
class Result:
demand_interactive: np.ndarray # (N, W) productive hours, queued work
demand_deferrable: np.ndarray # (N, W) productive hours, email and similar
supply: np.ndarray # (N, W) available productive hours
shrink: np.ndarray # (N, W) sampled shrink, kept for the S-curve
hours_per_head: float
inputs: dict = field(default_factory=dict) # name -> (N,) scalar draws
endogenous: set = field(default_factory=set) # inputs that are model OUTPUTS too
@property
def demand(self) -> np.ndarray:
"""Total required productive hours.
Read with `coverage(lane="interactive")` as a priority model: queued work
has first call on capacity, deferrable work absorbs what is left. The
pooled figure is therefore the STRICTER test, not the flattering one --
the gap between the two is what the email backlog absorbs before the
queue is affected at all.
Summing does assume the hours are interchangeable. Where interactive and
deferrable work are staffed by different, non-cross-skilled people, the
residual is not actually available and both figures need a supply split
the model does not currently carry.
"""
return self.demand_interactive + self.demand_deferrable
@property
def shortfall(self) -> np.ndarray:
return self.demand - self.supply
def coverage(self, week=None, lane: str = "all") -> float:
"""lane: 'interactive' asks whether queued work is covered when it has
first call on capacity; 'all' asks whether everything is, deferrable
work included. The gap between them is the deferrable exposure."""
if lane == "all":
sf = self.shortfall
elif lane == "interactive":
sf = self.demand_interactive - self.supply
else:
raise ValueError(f"lane must be 'all' or 'interactive', got {lane!r}")
sf = sf if week is None else sf[:, week]
if week is None:
return float((sf.max(axis=1) <= 0).mean()) # every week covered
return float((sf <= 0).mean())
def productive_hours_per_head(self, week: int) -> np.ndarray:
"""What one head actually delivers, per draw: scheduled hours net of shrink.
Converting a shortfall to FTE by dividing by SCHEDULED hours is the
shrink-once rule broken in the other direction -- the shortfall is in
productive hours, so its denominator must be productive hours too.
Dividing by 37.5 rather than 37.5 x (1 - shrink) understates the gap by
the whole shrink percentage.
"""
return self.hours_per_head * (1.0 - self.shrink[:, week])
def shortfall_fte(self, week: int) -> np.ndarray:
"""Per-draw shortfall in FTE. Converted inside the draw, then summarized
-- not a percentile of hours divided by an average head."""
return self.shortfall[:, week] / self.productive_hours_per_head(week)
def summary(self) -> pd.DataFrame:
rows = []
for w in range(self.demand.shape[1]):
d, s = self.demand[:, w], self.supply[:, w]
sf = d - s
fte = self.shortfall_fte(w)
rows.append({
"week": w,
"demand_p10": np.percentile(d, 10),
"demand_p50": np.percentile(d, 50),
"demand_p90": np.percentile(d, 90),
"supply_p50": np.percentile(s, 50),
"coverage": float((sf <= 0).mean()),
"coverage_interactive": self.coverage(w, lane="interactive"),
"deferrable_share": float(np.mean(self.demand_deferrable[:, w]
/ np.maximum(d, 1e-9))),
"shortfall_p50_hrs": np.percentile(sf, 50),
"shortfall_p90_hrs": np.percentile(sf, 90),
"shortfall_p50_fte": np.percentile(fte, 50),
"shortfall_p90_fte": np.percentile(fte, 90),
})
return pd.DataFrame(rows)
# --------------------------------------------------------------------------
# demand
# --------------------------------------------------------------------------
def simulate_demand(cfg, rng, N, W, inputs, endogenous):
"""Two passes: accumulate workload per channel, then convert to productive hours.
The passes are separate because occupancy may be a function of LOAD, and the
load on a channel is a portfolio quantity -- it is not known until every
segment has contributed. Converting inside the segment loop would apply an
occupancy derived from one segment's volume to the whole channel.
Returns (interactive, deferrable), each (N, W) productive hours.
"""
chans = cfg["channels"]
# Channel parameters are properties of the PORTFOLIO, not of a segment, and
# are drawn once per simulation before the segment loop.
#
# Drawing them inside the loop -- as an earlier version did -- gives voice AHT
# an independent value for each segment, so a high draw in one is offset by a
# low draw in another and the portfolio-level risk is diversified away. It is
# the same error the forecast-error comment below warns about, one level up.
# Where a channel genuinely behaves differently by segment, give that segment
# its own channel entry rather than relying on resampling to represent it.
cpar = {}
for c, ch in chans.items():
if ch["kind"] == "deferrable":
tph = draw_lognormal(rng, ch["items_per_productive_hour"], (N, 1))
cpar[c] = {"tph": tph}
inputs[f"tph[{c}]"] = tph[:, 0]
continue
aht = draw_lognormal(rng, ch["aht_seconds"], (N, 1))
conc = (draw_lognormal(rng, ch["concurrency"], (N, 1))
if isinstance(ch.get("concurrency"), dict) else
np.full((N, 1), float(ch.get("concurrency", 1.0))))
cpar[c] = {"aht": aht, "conc": conc}
inputs[f"AHT[{c}]"] = aht[:, 0]
if isinstance(ch.get("concurrency"), dict):
inputs[f"conc[{c}]"] = conc[:, 0]
# ---- pass 1: workload (interactive) and item counts (deferrable) --------
workload = {c: np.zeros((N, W)) for c in chans}
for sname, seg in cfg["segments"].items():
lo, mode, hi = seg["transactions_per_week"]
a, b = pert_ab(lo, mode, hi)
base = lo + (hi - lo) * rng.beta(a, b, size=(N, 1))
index = np.asarray(seg.get("weekly_index", [1.0] * W), float)
if index.size != W:
raise ValueError(f"{sname}: weekly_index length {index.size} != horizon {W}")
# Forecast error has TWO components, and the model needs both.
#
# The systematic part is drawn once per simulation: forecast bias
# persists, and redrawing it weekly would average it away. But drawing
# ONLY that leaves every week perfectly correlated -- measured at 0.9998
# before `forecast_error_weekly` existed -- which makes "the probability
# we cover every week" identical to "the probability we cover the worst
# week", and quietly removes week-to-week forecast error, a large and
# entirely real source of staffing risk.
err_sys = draw_lognormal(rng, seg["forecast_error"], (N, 1))
wspec = seg.get("forecast_error_weekly")
err_wk = draw_lognormal(rng, wspec, (N, W)) if wspec else 1.0
txn = base * index[None, :] * err_sys * err_wk
inputs[f"txn[{sname}]"] = txn.mean(axis=1)
inputs[f"fcst_err[{sname}]"] = err_sys[:, 0]
cr = draw_beta(rng, seg["contact_rate"], (N, 1))
inputs[f"CR[{sname}]"] = cr[:, 0]
contacts = txn * cr
alpha = np.asarray([seg["channel_split"][c] for c in chans], float)
split = rng.dirichlet(alpha, size=N) # (N, n_chan)
for i, c in enumerate(chans):
inputs[f"split[{sname},{c}]"] = split[:, i]
for i, (c, ch) in enumerate(chans.items()):
n_c = contacts * split[:, i : i + 1]
if ch["kind"] == "deferrable":
workload[c] += n_c # items, not hours
else:
workload[c] += n_c * cpar[c]["aht"] / 3600.0 / cpar[c]["conc"]
# ---- pass 2: workload -> productive hours -------------------------------
interactive = np.zeros((N, W))
deferrable = np.zeros((N, W))
for c, ch in chans.items():
if ch["kind"] == "deferrable":
deferrable += workload[c] / cpar[c]["tph"]
continue
aht_med = float(np.median(cpar[c]["aht"]))
occ = _occupancy_for(ch["occupancy"], rng, N, W, workload[c], aht_med)
inputs[f"occ[{c}]"] = occ.mean(axis=1)
if isinstance(ch["occupancy"], dict) and ch["occupancy"].get("kind") == "erlang_curve":
endogenous.add(f"occ[{c}]")
interactive += workload[c] / occ
return interactive, deferrable
# --------------------------------------------------------------------------
# supply
# --------------------------------------------------------------------------
def simulate_supply(cfg, rng, N, W, inputs):
sup = cfg["supply"]
hours = float(sup["scheduled_hours_per_head"])
ramp = _draw_curve(rng, sup["ramp_curve"], N, monotone=True) # (N, L)
early = _draw_curve(rng, sup["early_attrition_weekly"], N) # (N, Le)
inputs["ramp_wk1"] = ramp[:, 0]
inputs["attrition_early_wk1"] = early[:, 0]
tenured_rate = draw_beta(rng, sup["tenured_attrition_weekly"], (N, 1))[:, 0]
inputs["attrition_tenured"] = tenured_rate
# --- starting tenured pool -------------------------------------------
tenured = np.full(N, float(sup["starting_productive_heads"]))
# --- pipeline rates: drawn ONCE, shared across every cohort -----------
# These are properties of the labor market and the training operation, not
# of an individual class. Drawing them per class let a bad class be canceled
# by a good one, which averaged away the pipeline risk this whole side of the
# model exists to represent -- and left the attribution table reporting only
# the first cohort's draw.
#
# What this deliberately omits is per-class EXECUTION noise on top of the
# shared rate (a binomial draw around it). That is a second-order effect next
# to not knowing the rate, and adding it without evidence would manufacture
# precision. Revisit once several cohorts have been observed.
fill = draw_beta(rng, sup["requisition_fill_prob"], N)
cf = draw_beta(rng, sup["class_fill_rate"], N)
grad = draw_beta(rng, sup["graduation_rate"], N)
ttf = draw_lognormal(rng, sup["time_to_fill_weeks"], N)
inputs["req_fill_prob"] = fill
inputs["class_fill_rate"] = cf
inputs["graduation_rate"] = grad
inputs["time_to_fill_wks"] = ttf
# --- build cohorts from the class plan --------------------------------
grad_week, coh_heads = [], []
for cls in sup.get("class_plan", []):
offers = float(cls["requisitions_opened"]) * fill
seats = np.minimum(float(cls["seats_planned"]) * cf, offers)
heads = seats * grad
slip = np.ceil(np.maximum(ttf - float(cls.get("lead_weeks", 0)), 0)).astype(int)
gw = int(cls["planned_start_week"]) + int(sup["training_weeks"]) + slip
grad_week.append(gw) # (N,) — slip is per draw
coh_heads.append(heads)
grad_week = np.array(grad_week).T if grad_week else np.zeros((N, 0), int) # (N, C)
coh_heads = np.array(coh_heads).T if coh_heads else np.zeros((N, 0)) # (N, C)
# --- shrinkage --------------------------------------------------------
shrink_p = draw_beta(rng, sup["shrinkage_planned"], (N, 1))
shrink_u = draw_beta(rng, sup["shrinkage_unplanned"], (N, 1))
shrink = np.clip(shrink_p + shrink_u, 0.0, 0.95) * np.ones((1, W))
inputs["shrink_planned"] = shrink_p[:, 0]
inputs["shrink_unplanned"] = shrink_u[:, 0]
# --- walk the horizon -------------------------------------------------
supply = np.zeros((N, W))
for w in range(W):
tenured *= (1.0 - tenured_rate)
productive = tenured.copy()
if coh_heads.shape[1]:
k = w - grad_week # (N, C) weeks since grad
live = k >= 0
if live.any():
rows = np.arange(coh_heads.shape[0])[:, None]
L, Le = ramp.shape[1], early.shape[1]
eidx = np.clip(k, 0, Le - 1)
haz = np.where(k < Le, early[rows, eidx], tenured_rate[:, None])
coh_heads = np.where(live, coh_heads * (1.0 - haz), coh_heads)
ridx = np.clip(k, 0, L - 1)
mult = np.where(k < L, ramp[rows, ridx], 1.0)
productive += np.where(live, coh_heads * mult, 0.0).sum(axis=1)
supply[:, w] = productive * hours * (1.0 - shrink[:, w])
return supply, shrink, hours
# --------------------------------------------------------------------------
# run
# --------------------------------------------------------------------------
def run(cfg, n_draws: int = 20000, seed: int = SEED) -> Result:
rng = np.random.default_rng(seed)
W = int(cfg["horizon_weeks"])
inputs: dict = {}
endogenous: set = set()
inter, defer = simulate_demand(cfg, rng, n_draws, W, inputs, endogenous)
supply, shrink, hours = simulate_supply(cfg, rng, n_draws, W, inputs)
return Result(inter, defer, supply, shrink, hours, inputs, endogenous)
# --------------------------------------------------------------------------
# variance attribution
# --------------------------------------------------------------------------
def attribution(res: Result, week: int | None = None,
include_endogenous: bool = False) -> pd.DataFrame:
"""Rank inputs by share of shortfall variance, via squared rank correlation.
SCREENING TOOL ONLY. It assumes each input's effect is monotone and roughly
additive. It will mislead where two inputs interact (occupancy and volume do).
Confirm the top two or three by pinning them at their mean and re-running.
ENDOGENOUS inputs are flagged and, by default, excluded. Under `erlang_curve`
occupancy is a FUNCTION of load, so it inherits volume's correlation with
shortfall and ranks near the top -- while actually reducing the requirement.
Left in, it reads as the second-largest driver and sends the gap register
chasing a quantity that is an output of the model, not an input to it.
"""
y = res.shortfall.max(axis=1) if week is None else res.shortfall[:, week]
rows = []
for name, x in res.inputs.items():
if np.std(x) == 0:
continue
endo = name in res.endogenous
if endo and not include_endogenous:
continue
r = stats.spearmanr(x, y).statistic
rows.append({"input": name, "rho": r, "rho2": r * r, "endogenous": endo})
df = pd.DataFrame(rows)
if df.empty:
return df
df["share"] = df["rho2"] / df["rho2"].sum()
return df.sort_values("share", ascending=False).reset_index(drop=True)
def pin_check(cfg, name_to_fixed: dict, n_draws=20000, seed=SEED, week=None):
"""Exact contribution: re-run with a parameter's spread removed.
Pass e.g. {'channels.chat.aht_seconds': {'mu': m, 'sigma': 0.0}} and compare
the shortfall variance against the full run.
"""
import copy
c = copy.deepcopy(cfg)
for dotted, val in name_to_fixed.items():
node = c
*path, leaf = dotted.split(".")
for p in path:
node = node[p]
node[leaf] = val
res = run(c, n_draws, seed)
y = res.shortfall.max(axis=1) if week is None else res.shortfall[:, week]
return float(np.var(y))
# --------------------------------------------------------------------------
# the decision curve
# --------------------------------------------------------------------------
def coverage_curve(res: Result, extra_heads=range(0, 121, 5)) -> pd.DataFrame:
"""P(covered every week) against additional fully-ramped heads from week 0.
Fully-ramped is deliberately optimistic: it is the ceiling on what extra
hiring can buy. Real heads arrive late and ramped, so read this as an upper
bound and use the pipeline in the config to model the achievable version.
"""
rows = []
for h in extra_heads:
added = h * res.hours_per_head * (1.0 - res.shrink)
sf = res.demand - (res.supply + added)
by_week = (sf <= 0).mean(axis=0) # coverage in each week
rows.append({"extra_heads": h,
"coverage_all_weeks": float((sf.max(axis=1) <= 0).mean()),
"coverage_worst_week": float(by_week.min()),
"coverage_mean_week": float(by_week.mean())})
return pd.DataFrame(rows)
# --------------------------------------------------------------------------
# charts
# --------------------------------------------------------------------------
def charts(res: Result, week: int | None = None, outdir: str = "."):
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
w = week if week is not None else res.demand.shape[1] - 1
paths = []
# 1. demand vs supply for the chosen week
fig, ax = plt.subplots(figsize=(9, 4.5))
ax.hist(res.demand[:, w], bins=80, alpha=0.6, label="required productive hours")
ax.hist(res.supply[:, w], bins=80, alpha=0.6, label="available productive hours")
ax.set_title(f"Week {w}: demand vs supply | coverage = {res.coverage(w):.1%}")
ax.set_xlabel("productive hours"); ax.legend()
p = f"{outdir}/01-demand-vs-supply.png"; fig.savefig(p, dpi=140, bbox_inches="tight")
plt.close(fig); paths.append(p)
# 2. coverage by week
s = res.summary()
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(s["week"], s["coverage"], marker="o")
ax.axhline(0.8, ls="--", lw=1); ax.set_ylim(0, 1)
ax.set_title("Coverage probability by week"); ax.set_xlabel("week")
p = f"{outdir}/02-coverage-by-week.png"; fig.savefig(p, dpi=140, bbox_inches="tight")
plt.close(fig); paths.append(p)
# 3. requirement fan
fig, ax = plt.subplots(figsize=(9, 4))
ax.fill_between(s["week"], s["demand_p10"], s["demand_p90"], alpha=0.25,
label="demand P10-P90")
ax.plot(s["week"], s["demand_p50"], label="demand P50")
ax.plot(s["week"], s["supply_p50"], ls="--", label="supply P50")
ax.set_title("Required vs available productive hours"); ax.legend()
p = f"{outdir}/03-fan.png"; fig.savefig(p, dpi=140, bbox_inches="tight")
plt.close(fig); paths.append(p)
# 4. attribution tornado
a = attribution(res).head(12)
if not a.empty:
fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(a["input"][::-1], a["share"][::-1])
ax.set_title("Share of shortfall variance (screening)")
ax.set_xlabel("share")
p = f"{outdir}/04-attribution.png"; fig.savefig(p, dpi=140, bbox_inches="tight")
plt.close(fig); paths.append(p)
# 5. the decision curve
cc = coverage_curve(res)
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(cc["extra_heads"], cc["coverage_all_weeks"], marker="o")
ax.axhline(0.8, ls="--", lw=1); ax.set_ylim(0, 1)
ax.set_title("Coverage probability vs additional heads (upper bound)")
ax.set_xlabel("additional fully-ramped heads")
p = f"{outdir}/05-decision-curve.png"; fig.savefig(p, dpi=140, bbox_inches="tight")
plt.close(fig); paths.append(p)
return paths
# --------------------------------------------------------------------------
# demonstration parameters — REPLACE with params.yaml
# --------------------------------------------------------------------------
DEMO = {
"horizon_weeks": 13,
"channels": {
"voice": {"kind": "interactive", "aht_seconds": {"ci": [380, 520]},
"concurrency": 1.0,
"occupancy": {"kind": "erlang_curve", "target_sl": 0.80,
"asa_seconds": 20, "operating_hours_per_week": 168,
"spread": 0.03}},
"chat": {"kind": "interactive", "aht_seconds": {"ci": [600, 950]},
"concurrency": {"ci": [1.4, 2.2]},
"occupancy": {"kind": "erlang_curve", "target_sl": 0.80,
"asa_seconds": 30, "operating_hours_per_week": 168,
"spread": 0.03}},
"email": {"kind": "deferrable", "items_per_productive_hour": {"ci": [4.0, 7.0]}},
},
"segments": {
"SEG-A": {"transactions_per_week": [950_000, 1_100_000, 1_320_000],
"weekly_index": [1.00, 1.02, 1.05, 1.09, 1.12, 1.10, 1.06,
1.03, 1.00, 0.97, 0.95, 0.94, 0.96],
"forecast_error": {"ci": [0.90, 1.12]},
"forecast_error_weekly": {"ci": [0.93, 1.08]},
"contact_rate": {"mean": 0.031, "ci": [0.024, 0.040]},
"channel_split": {"voice": 55, "chat": 30, "email": 15}},
"SEG-B": {"transactions_per_week": [300_000, 375_000, 525_000],
"weekly_index": [1.00, 1.01, 1.03, 1.04, 1.06, 1.05, 1.03,
1.01, 1.00, 0.99, 0.98, 0.98, 0.99],
"forecast_error": {"ci": [0.85, 1.20]},
"forecast_error_weekly": {"ci": [0.90, 1.11]},
"contact_rate": {"mean": 0.052, "ci": [0.038, 0.070]},
"channel_split": {"voice": 40, "chat": 42, "email": 18}},
},
"supply": {
"starting_productive_heads": 300,
"scheduled_hours_per_head": 37.5,
"ramp_curve": [{"mean": 0.45, "ci": [0.30, 0.60]},
{"mean": 0.70, "ci": [0.55, 0.82]},
{"mean": 0.85, "ci": [0.73, 0.93]}, 1.00],
"early_attrition_weekly": [{"mean": 0.030, "ci": [0.015, 0.050]},
{"mean": 0.025, "ci": [0.012, 0.043]},
{"mean": 0.018, "ci": [0.008, 0.033]},
{"mean": 0.012, "ci": [0.005, 0.024]}],
"tenured_attrition_weekly": {"mean": 0.006, "ci": [0.003, 0.010]},
"shrinkage_planned": {"mean": 0.18, "ci": [0.15, 0.21]},
"shrinkage_unplanned": {"mean": 0.09, "ci": [0.05, 0.15]},
"requisition_fill_prob": {"mean": 0.70, "ci": [0.50, 0.86]},
"time_to_fill_weeks": {"ci": [3.0, 9.0]},
"class_fill_rate": {"mean": 0.88, "ci": [0.72, 0.97]},
"graduation_rate": {"mean": 0.84, "ci": [0.70, 0.93]},
"training_weeks": 5,
"class_plan": [
{"planned_start_week": 0, "seats_planned": 24,
"requisitions_opened": 40, "lead_weeks": 6},
{"planned_start_week": 4, "seats_planned": 24,
"requisitions_opened": 40, "lead_weeks": 6},
],
},
}
def load_config(path: str = "params.yaml"):
if os.path.exists(path):
import yaml
with open(path) as fh:
return yaml.safe_load(fh)
print(f"[warn] {path} not found — running on DEMONSTRATION parameters. "
f"Results are illustrative only.")
return DEMO
if __name__ == "__main__":
cfg = load_config()
res = run(cfg)
pd.set_option("display.width", 160)
print(res.summary().to_string(index=False, float_format=lambda v: f"{v:,.2f}"))
print(f"\nCoverage, every week: {res.coverage():.1%}")
print("\nVariance attribution (screening):")
print(attribution(res).head(10).to_string(index=False))
print("\nDecision curve:")
print(coverage_curve(res).to_string(index=False))
print("\nCharts:", charts(res))
```
Block 5 — 03b-updating.md
# 03b — Updating and scoring
The posterior-update, forgetting, quantile-archive and calibration helpers. `cycle.py`
imports from this module; `mc_staffing.py` does not. Save all four modules side by side.
Split from `03-simulator.md` for a mechanical reason worth knowing: a single block carrying
both modules ran past 1,000 lines, which is the point at which the wiki's syntax highlighter
silently falls back to unhighlighted plain text and the block loses its visible boundary.
---
## `update.py`
```python
"""Posterior updating, calibration scoring and regime detection."""
from __future__ import annotations
import numpy as np
from scipy import stats
# --------------------------------------------------------------------------
# forgetting
# --------------------------------------------------------------------------
def discount_beta(a, b, lam):
"""Power-discount a Beta toward its uniform base so old evidence decays."""
return 1 + lam * (a - 1), 1 + lam * (b - 1)
def discount_dirichlet(alpha, lam):
return [1 + lam * (ai - 1) for ai in alpha]
# --------------------------------------------------------------------------
# conjugate updates
# --------------------------------------------------------------------------
def update_beta(a, b, successes, trials, lam=1.0, effective_n=None):
"""Beta-Binomial update.
For genuinely independent events (graduations out of starts, offers out of
requisitions) pass the raw counts and leave `effective_n` alone.
For quantities measured in hours or on correlated units -- shrinkage above
all -- pass `effective_n` as the number of INDEPENDENT units behind the
observation (agent-weeks, not hours). Successes are rescaled to match, which
is the whole point: treating 150,000 shrinkage hours as 150,000 independent
trials produces a posterior so tight it stops responding to reality.
Rescaling here rather than at the call site is deliberate. Passing a raw
success count against a reduced trial count is an easy mistake to make and
silently inflates the rate.
"""
if trials <= 0:
raise ValueError(f"trials must be positive, got {trials}")
if successes < 0 or successes > trials:
raise ValueError(f"successes {successes} outside [0, {trials}]")
if effective_n is not None:
if effective_n <= 0:
raise ValueError(f"effective_n must be positive, got {effective_n}")
successes = successes / trials * effective_n
trials = effective_n
a, b = discount_beta(a, b, lam)
return a + successes, b + (trials - successes)
def update_dirichlet(alpha, counts, lam=1.0):
alpha = discount_dirichlet(alpha, lam)
return [ai + ci for ai, ci in zip(alpha, counts)]
def update_normal_invgamma(mu0, kappa0, alpha0, beta0, x, lam=1.0):
"""Normal-Inverse-Gamma update on log-scale observations (AHT, time-to-fill,
forecast error). Updates BOTH location and spread — never freeze the variance.
Returns (mu, kappa, alpha, beta); the predictive is Student-t.
"""
x = np.log(np.asarray(x, float))
x = x[np.isfinite(x)]
n, xbar = len(x), x.mean() if len(x) else 0.0
# Power-discount ALL THREE accumulating parameters. Discounting kappa and
# alpha while leaving beta untouched makes alpha reach a steady state while
# beta grows without bound, so sqrt(beta/(alpha-1)) inflates forever: against
# a true sigma of 0.20 at lam=0.98, one observation a day, it drifts past 0.5
# inside a year. That failure is insidious because it errs toward WIDE
# intervals, which read as humility rather than as a bug.
kappa0, alpha0, beta0 = lam * kappa0, lam * alpha0, lam * beta0
kappa = kappa0 + n
mu = (kappa0 * mu0 + n * xbar) / kappa
alpha = alpha0 + n / 2
ss = ((x - xbar) ** 2).sum() if n else 0.0
beta = beta0 + 0.5 * ss + (kappa0 * n * (xbar - mu0) ** 2) / (2 * kappa)
return mu, kappa, alpha, beta
def lognormal_params(mu, kappa, alpha, beta):
"""Point (mu, sigma) for sampling, from a Normal-Inverse-Gamma posterior."""
return mu, float(np.sqrt(beta / max(alpha - 1, 1e-6)))
def ci_from_beta(a, b, conf=0.90):
ql, qh = (1 - conf) / 2, 0.5 + conf / 2
return float(stats.beta.ppf(ql, a, b)), float(stats.beta.ppf(qh, a, b))
# --------------------------------------------------------------------------
# calibration
# --------------------------------------------------------------------------
def pit(ensemble, actual):
"""Probability integral transform: where the actual fell in the forecast.
A well-calibrated model produces PIT values uniform on [0,1]. Clustering
near 1 means the model forecasts LOW; near 0, HIGH; clustering in the
middle means the intervals are too wide.
"""
e = np.asarray(ensemble, float)
return float((e <= actual).mean())
def crps(ensemble, actual):
"""Continuous ranked probability score. Lower is better; units are the
units of the forecast. Sorted-sample form, O(n log n)."""
x = np.sort(np.asarray(ensemble, float))
n = len(x)
term1 = np.abs(x - actual).mean()
i = np.arange(1, n + 1)
term2 = (2.0 / (n * n)) * np.sum((2 * i - n - 1) * x)
return float(term1 - 0.5 * term2)
def interval_hit(ensemble, actual, conf=0.80):
lo = np.percentile(ensemble, (1 - conf) / 2 * 100)
hi = np.percentile(ensemble, (0.5 + conf / 2) * 100)
return bool(lo <= actual <= hi)
# --------------------------------------------------------------------------
# forecast archiving — the mechanism that makes scoring possible at all
# --------------------------------------------------------------------------
# Midpoint quadrature over [0,1]: tau_i = (i - 0.5)/M. The midpoint rule is used
# rather than an evenly-spaced grid from 0.025 to 0.975 because CRPS is an
# integral over the WHOLE unit interval -- a grid that stops short of the tails
# omits that mass and biases every score in the same direction. Measured against
# the ensemble CRPS: an 0.025-0.975 grid runs 2.2-2.6% high, this one within
# 0.2%, at identical storage cost.
QM = 40
QGRID = (np.arange(1, QM + 1) - 0.5) / QM
def archive_quantiles(ensemble, grid=QGRID):
"""Compress a predictive ensemble to a quantile grid for storage.
A session cannot keep 40,000 draws between cycles, and without the previous
cycle's forecast there is nothing to score -- which would leave the model
making claims it cannot be held to. A 40-point grid is small enough to sit
in MODEL-STATE.md and dense enough for PIT, interval hits and CRPS.
"""
return np.percentile(np.asarray(ensemble, float), grid * 100).tolist()
def pit_from_quantiles(q, actual, grid=QGRID):
"""PIT by interpolating the stored grid. Returns 0.0 or 1.0 beyond the ends.
Values pinned at exactly 0 or 1 are real information -- the actual fell
outside everything the model considered -- and must not be discarded.
"""
q = np.asarray(q, float)
if actual <= q[0]:
return 0.0
if actual >= q[-1]:
return 1.0
return float(np.interp(actual, q, grid))
def crps_from_quantiles(q, actual, grid=QGRID):
"""CRPS via the pinball-loss identity: CRPS = 2 * mean over tau of pinball.
Agrees with the ensemble `crps()` to within 0.2% on the 40-point midpoint
grid, measured on both lognormal and normal ensembles.
"""
q = np.asarray(q, float)
pin = np.where(actual >= q, (actual - q) * grid, (q - actual) * (1 - grid))
return float(2.0 * pin.mean())
def interval_hit_from_quantiles(q, actual, conf=0.80, grid=QGRID):
q = np.asarray(q, float)
lo = float(np.interp((1 - conf) / 2, grid, q))
hi = float(np.interp(0.5 + conf / 2, grid, q))
return bool(lo <= actual <= hi)
# --------------------------------------------------------------------------
# aging — forgetting is a function of ELAPSED TIME, not of call count
# --------------------------------------------------------------------------
def lam_for(lam_daily: float, days_elapsed: float) -> float:
"""Convert a daily forgetting factor to the factor for a gap of N days.
Every update function discounts once per call. If a parameter is updated
weekly and lambda is quoted per day, calling update once a week ages it at
0.98 per WEEK instead of 0.98^7 = 0.868 -- roughly seven times too little
forgetting, silently. Pass `lam=lam_for(0.98, days_since_last_update)`.
This is also how a missing day is honoured: age the posterior without
feeding it evidence, via `discount_only` below. A gap is information about
how stale the model is; skipping it pretends no time passed.
"""
if not 0 < lam_daily <= 1:
raise ValueError(f"lam_daily must be in (0,1], got {lam_daily}")
if days_elapsed < 0:
raise ValueError(f"days_elapsed must be >= 0, got {days_elapsed}")
return float(lam_daily ** days_elapsed)
def discount_only_beta(a, b, lam_daily, days):
"""Age a Beta across days with no observation. Widens; never shifts the mean."""
return discount_beta(a, b, lam_for(lam_daily, days))
# --------------------------------------------------------------------------
# effective sample size
# --------------------------------------------------------------------------
def effective_n(raw_denominator: float, kind: str) -> float:
"""Independent units behind an observation, which is rarely the raw count.
Contact rate is the trap. Its denominator is transactions -- of the order of
a million a week -- but those transactions are not a million independent
trials of a stable propensity to make contact. They share a day, a
disruption, a marketing campaign, a system outage. Updated on the raw count,
the posterior concentration reaches the millions within days and the 90%
interval collapses to a few thousandths of a percentage point, on a
parameter that is the single largest contributor to the model's output
variance -- 44% across both segments on the demo, over exogenous inputs.
The model stops learning because it has decided it already knows.
The divisors below are deliberately crude and deliberately conservative.
They are `[estimated]` and should be replaced with a measured intra-day
correlation as soon as enough daily history exists to compute one.
"""
DIVISOR = {
"events": 1.0, # genuinely independent: graduations out of starts
"contacts": 20.0, # contacts within a day share conditions [estimated]
"transactions": 50.0, # bulk volume, heavily correlated [estimated]
"hours": 40.0, # ~ one agent-week per 40 hours [estimated]
}
if kind not in DIVISOR:
raise ValueError(f"unknown kind {kind!r}; expected one of {sorted(DIVISOR)}")
return max(raw_denominator / DIVISOR[kind], 1.0)
def regime_flag(pit_history, window=10, threshold=0.8, trigger=7):
"""Distinguish bad luck from a biased model.
Returns 'HIGH', 'LOW' or None. If 7 of the last 10 actuals landed above the
80th percentile of their forecast, the model is biased low — a structural
change, not a run of bad luck. Do not absorb this silently into a posterior:
surface it, and consider lowering the forgetting factor.
"""
h = np.asarray(pit_history[-window:], float)
if len(h) < window:
return None
if (h > threshold).sum() >= trigger:
return "LOW" # model forecasting low
if (h < 1 - threshold).sum() >= trigger:
return "HIGH" # model forecasting high
return None
```
Block 8 — 05b-intake.md
# 05b — Intake and reconstruction
`intake.py`: loading the files, asserting they are usable, and reconstructing what a closed
week actually required. The cycle that drives it is `05c-cycle.md`.
**This block is where the pack stopped auditing columns and started asserting properties.**
Four evaluation rounds fixed reconstruction defects one column at a time, and each round a
sibling column failed the same way. `assert_reconstructable` states the preconditions instead
— one row per (date, segment, channel), every configured channel and segment actually
reporting, every value inside a declared physical range, no blanks in any column the
arithmetic consumes — and raises `Unscoreable` for anything outside them.
---
## `intake.py`
```python
"""Intake and reconstruction: getting data in, asserting it is usable, and
rebuilding what a closed week actually required.
Split from `cycle.py` because this layer grew to carry the precondition
assertions -- row-set identity, config-data parity, value ranges -- that four
earlier rounds kept rediscovering one column at a time. It is now the larger
half and it stands on its own: nothing here touches MODEL-STATE.md.
"""
from __future__ import annotations
import os
from datetime import date, datetime
import numpy as np
import pandas as pd
class Unscoreable(Exception):
"""A week cannot be reconstructed honestly from the data supplied.
Raised rather than returned as a number, because every previous version of
this module answered a missing or partial column with a plausible-looking
figure instead of a complaint. Measured on the demonstration intake, the
silent answers were: concurrency absent +20.5%, AHT absent -81.3%, AHT
present for only part of the week -0.1%, a supply column all-NaN -100%, and
two missing supply days -28.3%. Every one of those passed validation and
landed in the calibration log as though it were evidence.
The rule this class enforces: a reconstruction either uses the data it
claims to use, or it says why it cannot.
"""
# `aht_seconds` is required: without it `realized_required_hours` computes 0.0
# for every interactive channel, `score()` drops the row, and the calibration
# log stays empty forever -- the pack's central claim quietly never engaging.
# A loud failure at load is strictly better than a silent one at scoring.
REQUIRED_DEMAND = ["date", "segment", "channel", "transactions",
"contacts_offered", "contacts_handled", "productive_hours",
"aht_seconds"]
REQUIRED_SUPPLY = ["date", "scheduled_hours", "productive_hours",
"shrink_planned_hours", "shrink_unplanned_hours"]
def required_demand_columns(cfg) -> list[str]:
"""`REQUIRED_DEMAND`, plus `concurrency_effective` when any configured
channel is concurrent.
This has to fail at LOAD. Making it a per-week refusal produced the round-3
failure in a new costume: a deployment built to the documented minimum with
a chat channel scored nothing, ever, while `MODEL-STATE.md` reported a
comfortable "Scored periods 0" and the same truncated flag was re-appended
every cycle. Refusing once, loudly, at the door is the only version an operator
can act on.
"""
cols = list(REQUIRED_DEMAND)
if any(ch.get("kind") != "deferrable" and _planned_concurrency(ch) > 1 + 1e-9
for ch in cfg["channels"].values()):
cols.append("concurrency_effective")
# The same argument applied to the column beside it. A file built to the
# documented minimum with a deferrable channel and no planned throughput
# used to load clean and then fail at scoring, every week, forever --
# refusing once at the door is the only version an operator can act on.
if any(ch.get("kind") == "deferrable" and _planned_throughput(ch) is None
for ch in cfg["channels"].values()):
cols.append("items_resolved")
return cols
def load_intake(directory: str, cfg=None):
def rd(name, required):
p = os.path.join(directory, name)
if not os.path.exists(p):
return None
df = pd.read_csv(p, parse_dates=["date"])
missing = [c for c in required if c not in df.columns]
if missing:
raise ValueError(f"{name} missing required columns: {missing}")
return df
dem_req = required_demand_columns(cfg) if cfg is not None else REQUIRED_DEMAND
return (rd("daily_demand.csv", dem_req),
rd("daily_supply.csv", REQUIRED_SUPPLY),
rd("pipeline_events.csv", ["date", "event", "count"]))
def validate(demand, supply, cfg, tol_reconcile: float = 0.02) -> list[str]:
"""Checks 1-4 and 6 from `04`. Check 5, the `check_beta` round-trip over
`params.yaml`, is a config check rather than a data check and is run
separately — see `check_beta` in `03-simulator.md`.
Returns failures; an empty list means clean.
Failures are returned rather than raised so a cycle can report ALL of them
at once. A planner fixing one thing at a time per run will stop bothering.
"""
fails = []
if demand is None:
return ["daily_demand.csv not found"]
known_seg, known_chan = set(cfg["segments"]), set(cfg["channels"])
for col, known, label in (("segment", known_seg, "segment"),
("channel", known_chan, "channel")):
unknown = set(demand[col].unique()) - known
if unknown:
fails.append(f"unknown {label}(s) not in params.yaml: {sorted(unknown)}")
bad = demand["contacts_handled"] > demand["contacts_offered"]
if bad.any():
fails.append(f"contacts_handled > contacts_offered on {int(bad.sum())} row(s)")
dup = demand.duplicated(subset=["date", "segment", "channel"]).sum()
if dup:
fails.append(f"{dup} duplicate (date, segment, channel) row(s) in daily_demand.csv "
f"— the file may have been loaded twice")
if supply is not None:
key = [c for c in ("date", "site", "cohort_id") if c in supply]
sdup = supply.duplicated(subset=key).sum()
if sdup:
fails.append(f"{sdup} duplicate row(s) in daily_supply.csv on {key}")
if "occupancy_achieved" in demand:
oob = demand[(demand["occupancy_achieved"] > 1.0)
| (demand["occupancy_achieved"] <= 0)]
if len(oob):
# Occupancy above 1 is the most common occupancy data defect -- after-
# call-work definitions routinely push it over 100% -- and it used to
# pass here and then kill the cycle from inside update_beta with an
# unhandled ValueError, after scoring, with no state written.
fails.append(f"{len(oob)} occupancy_achieved value(s) outside (0, 1] "
f"(max {demand['occupancy_achieved'].max():.3f}) — check the "
f"definition, particularly whether after-call work is counted")
if supply is not None:
sh = supply["shrink_planned_hours"] + supply["shrink_unplanned_hours"]
over = sh > supply["scheduled_hours"]
if over.any():
fails.append(f"shrink hours exceed scheduled hours on {int(over.sum())} row(s)")
d_hrs = demand.groupby("date")["productive_hours"].sum()
s_hrs = supply.groupby("date")["productive_hours"].sum()
common = d_hrs.index.intersection(s_hrs.index)
if len(common):
gap = ((d_hrs[common] - s_hrs[common]).abs()
/ s_hrs[common].replace(0, np.nan))
worst = float(gap.max())
if worst > tol_reconcile:
fails.append(
f"productive hours reconcile to {worst:.1%} at worst "
f"(tolerance {tol_reconcile:.0%}) — demand and supply files "
f"disagree on what a productive hour is")
ops, closed = operating_days(cfg), closed_dates(cfg)
for frame, label in ((demand, "daily_demand.csv"), (supply, "daily_supply.csv")):
if frame is None:
continue
days = pd.to_datetime(sorted(frame["date"].unique()))
if len(days) > 1:
span = pd.date_range(days.min(), days.max(), freq="D")
# Only days the calendar says the operation runs. Checking all seven
# made every weekend a validation failure for a five-day operation,
# which halted the cycle permanently -- and the remedy the message
# implied (interpolate) is the one `04` forbids.
expect = [d for d in span if d.weekday() in ops
and d.normalize() not in closed]
gaps = pd.DatetimeIndex(expect).difference(days)
if len(gaps):
fails.append(f"{label}: {len(gaps)} missing operating date(s), first "
f"{gaps[0].date()} — record as missing, do not interpolate; "
f"if the operation was shut, declare it in "
f"calendar.closed_dates")
if supply is not None:
d_days = set(pd.to_datetime(demand["date"].unique()))
s_days = set(pd.to_datetime(supply["date"].unique()))
only_d, only_s = sorted(d_days - s_days), sorted(s_days - d_days)
# Dates present in one file and not the other used to be dropped by the
# reconciliation intersection rather than reported, which is how a
# supply file missing two days passed validation clean.
if only_d:
fails.append(f"{len(only_d)} date(s) in demand but not supply, first "
f"{only_d[0].date()}")
if only_s:
fails.append(f"{len(only_s)} date(s) in supply but not demand, first "
f"{only_s[0].date()}")
fails += check_calendar(cfg)
fails += check_closed_dates(cfg, demand, supply)
fails += _prior_band_check(demand, cfg)
fails += [m for m in _consistency_residual(demand, cfg) if not m.startswith(ADVISORY)]
return fails
def advisories(demand, cfg) -> list[str]:
"""Notices that reduce assurance without making a week unscoreable.
Kept out of `validate()` because a non-empty failure list halts the cycle,
and halting on "this check is providing less assurance than you think" is
the over-refusal trap again: it stops the model running over a message about
the model's own limits. These are reported in the cycle result and carried
in the state file instead, where an operator can see them every cycle
without being blocked by them.
"""
return [m.removeprefix(ADVISORY) for m in _consistency_residual(demand, cfg)
if m.startswith(ADVISORY)]
def _prior_band_check(demand, cfg, factor: float = 3.0) -> list[str]:
"""Compare each observed quantity against the prior it was elicited from.
This exists because the consistency residual below cannot be relied on
alone. That residual compares `handled x AHT / conc / occupancy` against
`productive_hours` -- and `04` defines occupancy as productive-working over
productive-available, which makes the identity true BY CONSTRUCTION whenever
occupancy is derived from those same columns, as most BI layers derive it.
Measured: with derived occupancy, AHT supplied in minutes instead of seconds
passes the residual silently. The residual has detection power only when
occupancy is independently sourced.
This check has no such dependency. A units error shows up as a 60x
discrepancy against the elicited prior, and the prior comes from a different
place entirely -- a person, before the data existed.
"""
if demand is None:
return []
out = []
for c, ch in cfg["channels"].items():
g = demand[demand["channel"] == c]
if g.empty:
continue
for col, spec_key, label in (("aht_seconds", "aht_seconds", "AHT"),
("concurrency_effective", "concurrency", "concurrency"),
("items_resolved", None, None)):
if spec_key is None or col not in g or g[col].isna().all():
continue
spec = ch.get(spec_key)
if not isinstance(spec, dict):
continue
if "ci" in spec:
lo, hi = float(spec["ci"][0]), float(spec["ci"][1])
elif "mu" in spec:
mid = float(np.exp(spec["mu"]))
lo = hi = mid
else:
continue
obs = float(g[col].median())
if obs < lo / factor or obs > hi * factor:
out.append(
f"{c} {label} median {obs:,.4g} is more than {factor:g}x outside its "
f"elicited prior [{lo:,.4g}, {hi:,.4g}] — check units (seconds vs "
f"minutes) and the column definition before trusting any score")
return out
def _consistency_residual(demand, cfg, tol: float = 0.05) -> list[str]:
"""Cross-check the columns against each other, not just against a range.
`realized_required_hours` derives, and correctly refuses to SCORE on, the
identity
handled x AHT / concurrency / occupancy_achieved == productive_hours
because it is circular. Circularity is exactly what makes it a good
VALIDATION residual: the two sides come from different columns, so they only
agree if the columns agree with each other.
Range checks cannot do this. A value can sit inside its physical range, be
non-blank, belong to a complete row set, and still be wrong by an order of
magnitude -- AHT supplied in minutes rather than seconds passes every range
check and understates the requirement by about 80%. This residual catches
that class: unit errors, upstream double-counts, and a segment zeroed in
place while its hours stay put.
It is a CHECK, never an input to the answer.
"""
need = ["contacts_handled", "aht_seconds", "occupancy_achieved", "productive_hours"]
if demand is None:
return []
missing = [c for c in need if c not in demand]
if missing:
# Saying nothing here was the defect: on a file built to the documented
# minimum the check simply did not run, and reported no opinion, so an
# operator had no way to know the value-axis check was absent.
return [ADVISORY + f"consistency residual INACTIVE — {missing} not supplied. The value-axis "
f"check is not running; a units error or a double-count will not be "
f"caught by it"]
d = demand[demand["channel"].isin(
[c for c, ch in cfg["channels"].items() if ch.get("kind") != "deferrable"])]
d = d.dropna(subset=need)
if d.empty:
return []
# PER WEEK, not pooled over all history. Pooled, a bad week is diluted by
# every good one, so the check got weaker the longer the model ran -- the
# opposite of what a calibration instrument should do. A single bad week in
# thirteen fell under tolerance; at twenty-six it vanished.
out, gaps = [], []
conc_all = d.get("concurrency_effective")
for wk, g in d.assign(_w=d["date"].map(week_key)).groupby("_w"):
conc = g["concurrency_effective"].fillna(1.0) if conc_all is not None else 1.0
implied = (g["contacts_handled"] * g["aht_seconds"] / 3600.0
/ conc / g["occupancy_achieved"]).sum()
actual = g["productive_hours"].sum()
if actual <= 0:
continue
gap = implied / actual - 1.0
gaps.append(gap)
if abs(gap) > tol:
out.append(f"consistency residual {gap:+.1%} in {wk} (tolerance {tol:.0%}): "
f"handled x AHT / concurrency / occupancy implies {implied:,.0f} "
f"productive hours against {actual:,.0f} reported. Check units on "
f"aht_seconds, the concurrency definition, and for double-counted "
f"contacts")
# A residual that is identically zero every week means occupancy was derived
# from the very columns being compared, which makes the identity true by
# construction and the check worthless. Say so rather than reporting health.
if len(gaps) >= 3 and max(abs(g) for g in gaps) < 1e-6:
out.append(ADVISORY + "consistency residual is identically zero every week — "
"occupancy_achieved appears to be DERIVED from handled x AHT / "
"concurrency / productive_hours, which makes this check a tautology. "
"It is providing no assurance. Source occupancy independently from the "
"ACD, or rely on the prior-band check instead")
return out
ADVISORY = "[advisory] "
DAY_NAMES = ["mon", "tue", "wed", "thu", "fri", "sat", "sun"]
def operating_days(cfg) -> set[int]:
"""Weekday indices the operation actually runs, from `calendar.operating_days`.
Defaults to all seven. Without this the pack refused every real operating
shape it was handed: a Monday-to-Friday operation halted permanently on the
date-contiguity check, and nothing in the pack mentioned five-day weeks at
all. A refusal with no documented remedy blocks a deployment exactly as
effectively as a silent wrong number.
"""
cal = (cfg or {}).get("calendar") or {}
days = cal.get("operating_days")
if not days:
return set(range(7))
out = set()
for d in days:
key = str(d).strip().lower()[:3]
if key not in DAY_NAMES:
raise ValueError(f"calendar.operating_days: unknown day {d!r}")
out.add(DAY_NAMES.index(key))
return out
def check_calendar(cfg) -> list[str]:
"""The calendar must agree with the rest of the config, and be believable.
Declaring five operating days changed what the INTAKE accepted and nothing
else: `operating_hours_per_week` stayed at 168, the volume priors were
untouched, and `target_occupancy` kept dividing realized load by 168. A
five-day deployment therefore ran to completion with every actual in the
bottom 5% of its own interval and no flag raised — a refusal replaced by a
silent wrong number, which is the trade this pack exists to refuse.
Completion is not the bar. The calibration log coming back honest is.
"""
fails = []
ops = operating_days(cfg)
if len(ops) < 7:
implied = 24.0 * len(ops)
for c, ch in (cfg.get("channels") or {}).items():
spec = ch.get("occupancy")
if not (isinstance(spec, dict) and spec.get("kind") == "erlang_curve"):
continue
oh = float(spec.get("operating_hours_per_week", 168))
if oh > implied * 1.01:
fails.append(
f"calendar declares {len(ops)} operating day(s) but channel {c} "
f"spreads load over operating_hours_per_week={oh:g}. At most "
f"{implied:g} is consistent. Occupancy, and therefore the whole "
f"requirement, is computed against the wrong denominator")
return fails
def check_closed_dates(cfg, demand, supply) -> list[str]:
"""`closed_dates` suppresses gap detection, so it needs its own guard.
An outage declared as a closure scored clean and about 13% low, straight
into the calibration log as evidence. The whole value axis this pack spent
two rounds hardening was moved into an unvalidated operator-supplied config
field. A closure is a claim about the world and it should cost something to
make: it has to fall inside the data span, and it has to actually have no
data behind it.
"""
fails = []
closed = closed_dates(cfg)
if not closed:
return fails
frames = [f for f in (demand, supply) if f is not None and len(f)]
if not frames:
return fails
lo = min(pd.Timestamp(f["date"].min()).normalize() for f in frames)
hi = max(pd.Timestamp(f["date"].max()).normalize() for f in frames)
span_days = max((hi - lo).days + 1, 1)
inside = {d for d in closed if lo <= d <= hi}
for f, label in ((demand, "daily_demand.csv"), (supply, "daily_supply.csv")):
if f is None:
continue
have = {pd.Timestamp(x).normalize() for x in f["date"].unique()}
contradicted = sorted(inside & have)
if contradicted:
fails.append(
f"{len(contradicted)} closed_date(s) have rows in {label}, first "
f"{contradicted[0].date()} — a declared closure with data behind it is "
f"either a mislabelled operating day or an undeclared change")
share = len(inside) / span_days
if share > 0.2:
fails.append(
f"{len(inside)} of {span_days} days in the data span are declared closed "
f"({share:.0%}) — above 20% this stops being a holiday list and starts "
f"suppressing gap detection wholesale")
return fails
def closed_dates(cfg) -> set:
"""Dates the operation was deliberately shut — holidays, shutdowns.
A closed day is DECLARED, not inferred. Declaring it is what distinguishes
a holiday from a data outage, and the two need opposite responses: one is
expected and the other is a defect. Previously neither encoding worked —
omitting the day tripped contiguity, zeroing it tripped the zero check, and
`04` forbade the only remedy the code would accept.
"""
cal = (cfg or {}).get("calendar") or {}
return {pd.Timestamp(d).normalize() for d in (cal.get("closed_dates") or [])}
def expected_dates(cfg, week: str):
"""The dates a given week should carry, given the calendar."""
# date.fromisocalendar, not a parsed string: pandas does not read "2026-W37-1",
# and week_key emits ISO year/week, which is not always the calendar year.
iso_year, iso_week = int(week[:4]), int(week.split("W")[1])
monday = pd.Timestamp(date.fromisocalendar(iso_year, iso_week, 1))
ops, closed = operating_days(cfg), closed_dates(cfg)
return [d for i in range(7)
for d in [monday + pd.Timedelta(days=i)]
if d.weekday() in ops and d.normalize() not in closed]
def week_key(d) -> str:
iso = pd.Timestamp(d).isocalendar()
return f"{iso[0]}-W{iso[1]:02d}"
def complete_weeks(df, cfg=None) -> list[str]:
"""Weeks carrying every date the calendar expects. Partial weeks never score.
With no calendar this is all seven days, as before. With one it is the
operating days minus declared closures, so a five-day operation and a
holiday week both complete normally instead of never completing at all.
"""
w = df.assign(_w=df["date"].map(week_key))
have = w.groupby("_w")["date"].apply(lambda s: {pd.Timestamp(x).normalize()
for x in s.unique()})
out = []
for wk, dates in have.items():
want = {d.normalize() for d in expected_dates(cfg, wk)} if cfg is not None \
else None
if want is None:
if len(dates) == 7:
out.append(wk)
elif want and dates == want:
# Equality, not issubset. A lower bound let a week reporting mon-sun
# and a week reporting mon-fri both count as complete under a 5-day
# calendar, while the reconstruction summed every row present -- the
# two scored 30% apart against the same forecast.
out.append(wk)
return sorted(out)
# Physical ranges for every value the reconstruction reads. A quantity outside
# its range is a data defect, not a small number -- negative concurrency once
# produced 4.1 billion hours through the divisor floor, and a negative AHT
# produced a negative requirement. Both passed validation.
RANGES = {
"contacts_offered": (0.0, 1e9),
"contacts_handled": (0.0, 1e9),
"aht_seconds": (1.0, 36000.0),
"concurrency_effective": (1.0, 20.0),
"productive_hours": (0.0, 1e7),
"occupancy_achieved": (0.01, 1.0),
"items_resolved": (0.0, 1e9),
}
def assert_reconstructable(demand, cfg, week: str) -> None:
"""Assert the PRECONDITIONS of a reconstruction, rather than probing for
known failures one column at a time.
Four rounds of this pack fixed columns that had been found to fail, and each
round a sibling column failed the same way. Column-by-column auditing closes
the cases you thought of. These assertions close the properties instead:
1. ROW-SET IDENTITY — exactly one row per (date, segment, channel).
Double-loading an intake file inflated the requirement 88% and the
supply 100%, and validation was clean.
2. CONFIG-DATA PARITY — every configured channel and segment actually
reports. The reconstruction iterates the channels in the DATA, so a
channel that stops reporting silently removes its whole requirement:
-33% with nothing raised. A model cannot notice demand that never
arrives in its own input.
3. VALUE RANGES — every quantity inside a declared physical range.
4. COMPLETENESS — no blanks in any column the arithmetic consumes.
Anything outside these raises `Unscoreable`, which is caught by `score()`
and reported. This is falsifiable in a way "we probed thirty-one cases" is
not: to break it you must violate a stated property.
"""
d = demand[demand["date"].map(week_key) == week]
if d.empty:
raise Unscoreable(f"{week}: no demand rows")
# 1. row-set identity
key = ["date", "segment", "channel"]
dup = d.duplicated(subset=key).sum()
if dup:
raise Unscoreable(f"{week}: {dup} duplicate (date, segment, channel) row(s) — "
f"the file may have been loaded twice")
want = {x.normalize() for x in expected_dates(cfg, week)}
have = {pd.Timestamp(x).normalize() for x in d["date"].unique()}
if want - have:
missing = sorted(want - have)
raise Unscoreable(f"{week}: demand missing {len(missing)} expected day(s), "
f"first {missing[0].date()} — declare it in "
f"calendar.closed_dates if the operation was shut")
if have - want:
extra = sorted(have - want)
raise Unscoreable(f"{week}: demand carries {len(extra)} day(s) the calendar does "
f"not expect, first {extra[0].date()} — either the calendar is "
f"wrong or the operation ran when it says it was shut")
# 2. config-data parity, both directions
want_ch, have_ch = set(cfg["channels"]), set(d["channel"].unique())
if want_ch - have_ch:
raise Unscoreable(f"{week}: configured channel(s) absent from the data: "
f"{sorted(want_ch - have_ch)}")
if have_ch - want_ch:
raise Unscoreable(f"{week}: channel(s) in the data but not in params.yaml: "
f"{sorted(have_ch - want_ch)}")
want_sg, have_sg = set(cfg["segments"]), set(d["segment"].unique())
if want_sg - have_sg:
raise Unscoreable(f"{week}: configured segment(s) absent from the data: "
f"{sorted(want_sg - have_sg)}")
# 3 and 4. values present and in range, for the columns each lane consumes
for c, g in d.groupby("channel"):
kind = cfg["channels"][c]["kind"]
needed = ["contacts_offered"]
if kind == "deferrable":
# Only demanded when the reconstruction will actually read them.
# Requiring columns the arithmetic never touches contradicts this
# function's own stated property, and `04` tells an operator that
# leaving an unsourced column empty is legitimate.
if _planned_throughput(cfg["channels"][c]) is None:
needed += ["items_resolved", "productive_hours"]
else:
needed.append("aht_seconds")
if kind != "deferrable" and _planned_concurrency(cfg["channels"][c]) > 1 + 1e-9:
needed.append("concurrency_effective")
for col in needed:
if col not in g:
raise Unscoreable(f"{week}: {c} has no {col} column")
n_na = int(g[col].isna().sum())
if n_na:
raise Unscoreable(f"{week}: {c} missing {col} on {n_na} of {len(g)} rows")
lo, hi = RANGES[col]
bad = g[(g[col] < lo) | (g[col] > hi)]
if len(bad):
raise Unscoreable(f"{week}: {c} has {len(bad)} {col} value(s) outside "
f"[{lo}, {hi}] (min {g[col].min():.4g}, "
f"max {g[col].max():.4g})")
def assert_supply_reconstructable(supply, week: str, cfg=None) -> None:
"""Row-set and value preconditions for the supply side."""
s = supply[supply["date"].map(week_key) == week]
if s.empty:
raise Unscoreable(f"{week}: no supply rows")
want = {x.normalize() for x in expected_dates(cfg, week)}
have = {pd.Timestamp(x).normalize() for x in s["date"].unique()}
if want - have:
missing = sorted(want - have)
raise Unscoreable(f"{week}: supply missing {len(missing)} expected day(s), "
f"first {missing[0].date()} — declare it in "
f"calendar.closed_dates if the operation was shut")
key = [c for c in ("date", "site", "cohort_id") if c in s]
if s.duplicated(subset=key).sum():
raise Unscoreable(f"{week}: duplicate supply rows on {key} — "
f"the file may have been loaded twice")
if "productive_hours" not in s:
raise Unscoreable(f"{week}: supply has no productive_hours column")
n_na = int(s["productive_hours"].isna().sum())
if n_na:
raise Unscoreable(f"{week}: supply productive_hours blank on {n_na} row(s)")
lo, hi = RANGES["productive_hours"]
bad = s[(s["productive_hours"] < lo) | (s["productive_hours"] > hi)]
if len(bad):
raise Unscoreable(f"{week}: {len(bad)} supply productive_hours outside [{lo}, {hi}]")
total = float(s["productive_hours"].sum())
if total <= 0:
raise Unscoreable(f"{week}: supply productive_hours totals zero")
closed = closed_dates(cfg)
working = s[~s["date"].map(lambda x: pd.Timestamp(x).normalize()).isin(closed)]
zeroed = int((working["productive_hours"] == 0).sum())
if zeroed:
# Zero is inside the declared range, so the range check alone lets a
# non-reporting day through and scores the week low -- which is the
# missing-days defect with the rows present and zeroed instead. A day
# that genuinely delivered nothing has to be stated deliberately, not
# encoded as a zero among working days.
raise Unscoreable(
f"{week}: {zeroed} supply row(s) report zero productive hours on days the "
f"calendar says were open — if the operation was shut, declare the date in "
f"calendar.closed_dates; otherwise this is a reporting gap, not a zero")
def _planned_throughput(ch: dict):
"""Planned items per productive hour from params.yaml, or None."""
spec = ch.get("items_per_productive_hour")
if isinstance(spec, dict):
if "mu" in spec:
return float(np.exp(spec["mu"]))
if "ci" in spec:
return float(np.exp(np.mean(np.log(spec["ci"]))))
return None
def _planned_concurrency(ch: dict) -> float:
"""The concurrency the plan assumed, from params.yaml."""
spec = ch.get("concurrency", 1.0)
if isinstance(spec, dict):
if "mu" in spec:
return float(np.exp(spec["mu"]))
if "ci" in spec:
return float(np.exp(np.mean(np.log(spec["ci"]))))
return max(float(spec), 1e-6)
def target_occupancy(cfg, channel: str, workload_hours: float, aht: float) -> float:
"""The occupancy the plan was built on, at the realized load.
Deliberately NOT the achieved occupancy. See `realized_required_hours`.
"""
spec = cfg["channels"][channel]["occupancy"]
if isinstance(spec, dict) and spec.get("kind") == "erlang_curve":
from mc_staffing import occupancy_curve
a_grid, occ_grid = occupancy_curve(
float(spec.get("target_sl", 0.80)), float(spec.get("asa_seconds", 20.0)), aht)
load = workload_hours / float(spec["operating_hours_per_week"])
return float(np.interp(np.clip(load, a_grid[0], a_grid[-1]), a_grid, occ_grid))
# A fixed occupancy arrives in either form: elicited {mean, ci}, or {a, b}
# once `apply_posteriors` has written a learned posterior back. Reading only
# `mean` raised KeyError on the second cycle of every deployment using the
# documented fixed-occupancy fallback -- the two headline fixes of the
# previous round, each correct alone, meeting each other.
if "a" in spec and "b" in spec:
return float(spec["a"] / (spec["a"] + spec["b"]))
return float(spec["mean"])
def realized_required_hours(demand, cfg, week: str) -> float:
"""What the week WOULD have required to serve everything offered at target.
This is a COUNTERFACTUAL RECONSTRUCTION, and it has to be, for a reason that
is easy to get wrong -- an earlier version of this function got it wrong.
The obvious formulation uses handled contacts at ACHIEVED occupancy. That is
circular. `04-intake-schema.md` defines achieved occupancy as productive-
working over productive-available, and productive-working is exactly
handled x AHT / concurrency. So the expression reduces algebraically to the
productive hours the floor actually delivered. Scored that way, the
"requirement" can never exceed what was staffed: in an understaffed week
contacts go unhandled, occupancy rises, and the computed requirement stays
pinned to delivered hours. The calibration log would then be incapable of
detecting a demand miss -- the one failure the model exists to price -- while
appearing perfectly well calibrated.
So: OFFERED contacts, at the occupancy the plan assumed (or the curve's
occupancy at the realized load). The assumption being made is explicit --
that the work offered should have been served at target service -- and it is
a different quantity from `realized_supply_hours`, not an algebraic
restatement of it.
The two scored quantities are still not independent: both are driven by the
same week. Read their PIT histories together, not as separate evidence.
"""
assert_reconstructable(demand, cfg, week)
d = demand[demand["date"].map(week_key) == week]
total = 0.0
for c, g in d.groupby("channel"):
ch = cfg["channels"][c]
offered = float(g["contacts_offered"].sum())
if offered <= 0:
continue
if ch["kind"] == "deferrable":
# Planned throughput, for the same reason the interactive lane uses
# target occupancy: dividing offered work by the throughput the week
# actually achieved reintroduces the circularity in miniature -- a
# week that cleared its backlog by working faster would score as
# having required less. Observed throughput is used only if the
# config carries no planned figure.
spec = ch.get("items_per_productive_hour")
tph = None
if isinstance(spec, dict):
if "mu" in spec:
tph = float(np.exp(spec["mu"]))
elif "ci" in spec:
tph = float(np.exp(np.mean(np.log(spec["ci"]))))
if tph is None:
hrs = float(g["productive_hours"].sum())
items = float(g.get("items_resolved", pd.Series(dtype=float)).sum())
tph = (items / hrs) if hrs > 0 and items > 0 else None
if not tph or tph <= 0:
# The interactive lane refuses when it cannot reconstruct; this
# lane used to return a smaller number instead, which is the
# same asymmetry between the lanes that an earlier round fixed
# in the opposite direction.
raise Unscoreable(
f"{week}: {c} has no planned items_per_productive_hour and no "
f"usable observed throughput")
total += offered / tph
continue
if "aht_seconds" not in g:
raise Unscoreable(f"{week}: {c} has no aht_seconds column")
missing = int(g["aht_seconds"].isna().sum())
if missing:
# Partial coverage is the dangerous case: dropping the blank rows
# quietly scores the week on a subset of its own contacts and looks
# almost right (-0.1% on the demo), which is far harder to notice
# than a column that is absent altogether.
raise Unscoreable(
f"{week}: {c} missing aht_seconds on {missing} of {len(g)} rows")
wt = g["contacts_offered"].clip(lower=0)
if wt.sum() <= 0:
raise Unscoreable(f"{week}: {c} has no offered contacts")
aht = float((g["aht_seconds"] * wt).sum() / wt.sum())
# Concurrency falls back to the PLANNED value from params.yaml, never to
# 1.0. Defaulting to 1.0 silently inflated a chat channel's requirement
# by the whole concurrency factor -- and it is the one value guaranteed
# to be wrong for any channel configured with concurrency at all.
planned_conc = _planned_concurrency(ch)
conc_col = g.get("concurrency_effective", pd.Series(dtype=float))
has_obs = len(conc_col) > 0 and conc_col.notna().any()
if has_obs:
cc = g.loc[conc_col.notna()]
cw = cc["contacts_offered"].clip(lower=0)
conc = (max(float((cc["concurrency_effective"] * cw).sum() / cw.sum()), 1e-6)
if cw.sum() > 0 else planned_conc)
elif planned_conc > 1.0 + 1e-9:
# A channel the plan treats as concurrent cannot be reconstructed
# from a file that never measured its concurrency. Falling back to
# the plan here would score the channel against its own assumption
# and report a difference of zero, which is worse than not scoring.
raise Unscoreable(
f"{week}: {c} is configured with concurrency {planned_conc:.2f} "
f"but the intake carries no concurrency_effective")
else:
conc = planned_conc # single-session channel; plan and floor agree
workload = offered * aht / 3600.0 / conc
occ = target_occupancy(cfg, c, workload, aht)
if occ > 0:
total += workload / occ
return float(total)
def realized_supply_hours(supply, week: str, cfg=None) -> float:
"""Available productive hours for a COMPLETE week of supply data.
Week completeness was previously checked on the demand frame only, so a
supply file missing two days scored 28% low with nothing flagged -- and
because `run_cycle` only ever ran the regime check on required hours, the
detector was pointed at the other series.
"""
assert_supply_reconstructable(supply, week, cfg)
s = supply[supply["date"].map(week_key) == week]
return float(s["productive_hours"].sum())
```
Block 9 — 05c-cycle.md
# 05c — The cycle module
`cycle.py`: state, scoring, posterior updates and the run loop. It implements the protocol in
`05-daily-protocol.md` and sits on the intake layer in `05b-intake.md`.
Reads and rewrites `MODEL-STATE.md`. Run it against a directory holding `daily_demand.csv`,
`daily_supply.csv` and optionally `pipeline_events.csv`.
---
## `cycle.py`
```python
"""The cycle: state, scoring, posterior updates, and the run loop.
Reads and rewrites MODEL-STATE.md. The intake and reconstruction layer it sits
on is `intake.py`.
"""
from __future__ import annotations
import json
import re
from datetime import date, datetime
import numpy as np
import pandas as pd
from intake import (Unscoreable, advisories, complete_weeks, load_intake,
realized_required_hours, realized_supply_hours, validate,
week_key)
from update import (archive_quantiles, crps_from_quantiles, discount_only_beta,
effective_n, interval_hit_from_quantiles, lam_for,
pit_from_quantiles, regime_flag, update_beta,
update_dirichlet, update_normal_invgamma)
# NB: the fence marker is built rather than written literally. A literal triple
# backtick inside this file would close the very code block that carries it, and
# anything extracting the block would silently truncate here.
FENCE = "`" * 3
STATE_RE = re.compile(FENCE + r"json\s*\n(.*?)\n" + FENCE, re.S)
def read_state(path: str = "MODEL-STATE.md") -> dict:
with open(path) as fh:
m = STATE_RE.search(fh.read())
if not m:
raise ValueError(f"{path} contains no json state block")
return json.loads(m.group(1))
def blank_state(revision: int = 0) -> dict:
return {"revision": revision, "covers_through": None, "lambda_daily": {},
"posteriors": {}, "forecast_archive": [], "calibration": [],
"flags": [], "definitions": {}}
def bootstrap_state(cfg, default_lambda: float = 0.98) -> dict:
"""Build revision 0 from the elicited priors in params.yaml.
Every posterior starts marked `prior_only` with n_obs = 0. The mark is
removed by the first real observation and by nothing else -- not by age, and
not by a parameter having been carried a long time.
"""
from mc_staffing import beta_from_mean_ci, lognormal_from_ci
st = blank_state(0)
def beta_node(spec):
a, b = beta_from_mean_ci(spec["mean"], *spec["ci"])
return {"family": "beta", "a": a, "b": b, "n_obs": 0, "prior_only": True}
def nig_node(spec, kappa=5.0, alpha=3.0):
mu, sigma = lognormal_from_ci(*spec["ci"])
# beta chosen so E[sigma^2] = beta/(alpha-1) matches the elicited sigma.
return {"family": "nig", "mu": mu, "kappa": kappa, "alpha": alpha,
"beta": sigma ** 2 * (alpha - 1), "n_obs": 0, "prior_only": True}
for sname, seg in cfg["segments"].items():
st["posteriors"][f"CR[{sname}]"] = beta_node(seg["contact_rate"])
st["posteriors"][f"fcst_err[{sname}]"] = nig_node(seg["forecast_error"])
st["posteriors"][f"split[{sname}]"] = {
"family": "dirichlet",
"alpha": [float(seg["channel_split"][c]) for c in cfg["channels"]],
"channels": list(cfg["channels"]), "n_obs": 0, "prior_only": True}
for c, ch in cfg["channels"].items():
if ch["kind"] == "deferrable":
continue
st["posteriors"][f"AHT[{c}]"] = nig_node(ch["aht_seconds"])
occ_spec = ch["occupancy"]
if isinstance(occ_spec, dict) and occ_spec.get("kind") == "erlang_curve":
# The curve supplies occupancy to the simulation, so there is no
# occupancy parameter to learn. Achieved occupancy is still worth
# tracking: it is the empirical test of whether this operation sits
# on its theoretical curve, and it is what a fitted replacement
# curve would eventually be built from. Kept as a DIAGNOSTIC -- it
# is recorded and reported, and it does not feed the model.
node = beta_node({"mean": 0.80, "ci": [0.60, 0.93]})
node["role"] = "diagnostic"
node["note"] = ("occupancy comes from erlang_curve; this tracks "
"achieved occupancy to test the curve, and is not a "
"model input")
st["posteriors"][f"occ[{c}]"] = node
else:
st["posteriors"][f"occ[{c}]"] = beta_node(occ_spec)
sup = cfg["supply"]
for name, key in (("shrink_planned", "shrinkage_planned"),
("shrink_unplanned", "shrinkage_unplanned"),
("req_fill_prob", "requisition_fill_prob"),
("class_fill_rate", "class_fill_rate"),
("graduation_rate", "graduation_rate"),
("attrition_tenured", "tenured_attrition_weekly")):
st["posteriors"][name] = beta_node(sup[key])
SLOW = {"attrition_tenured": 0.995, "graduation_rate": 0.99,
"class_fill_rate": 0.99, "req_fill_prob": 0.99}
for p in st["posteriors"]:
st["lambda_daily"][p] = SLOW.get(p, default_lambda)
st["definitions"] = {"chat_aht": "UNCONFIRMED", "productive_hours": "UNCONFIRMED",
"occupancy": "UNCONFIRMED"}
return st
def frozen_parameters(cfg, state) -> list[str]:
"""Sampled parameters that carry no node in MODEL-STATE.md.
These neither learn nor age: they sit in params.yaml at whatever was
elicited, permanently. The pack's prose stated this count as four, then as
five, and both were wrong -- so it is computed here and reported in the
state file rather than asserted anywhere. A number that has been wrong
twice should not be written down a third time.
"""
have = set(state.get("posteriors", {}))
configured = []
for s in cfg.get("segments", {}):
configured += [f"CR[{s}]", f"fcst_err[{s}]", f"split[{s}]"]
if cfg["segments"][s].get("forecast_error_weekly"):
configured.append(f"fcst_err_weekly[{s}]")
for c, ch in cfg.get("channels", {}).items():
if ch.get("kind") == "deferrable":
configured.append(f"tph[{c}]")
continue
configured += [f"AHT[{c}]", f"occ[{c}]"]
if isinstance(ch.get("concurrency"), dict):
configured.append(f"conc[{c}]")
sup = cfg.get("supply", {})
for name, key in (("shrink_planned", "shrinkage_planned"),
("shrink_unplanned", "shrinkage_unplanned"),
("req_fill_prob", "requisition_fill_prob"),
("class_fill_rate", "class_fill_rate"),
("graduation_rate", "graduation_rate"),
("attrition_tenured", "tenured_attrition_weekly"),
("time_to_fill", "time_to_fill_weeks")):
if key in sup:
configured.append(name)
if any(isinstance(x, dict) for x in sup.get("ramp_curve", [])):
configured.append("ramp_curve")
if any(isinstance(x, dict) for x in sup.get("early_attrition_weekly", [])):
configured.append("early_attrition")
return sorted(p for p in configured if p not in have)
def write_state(state: dict, path: str = "MODEL-STATE.md", readout: str = "",
cfg=None) -> None:
"""Emit MODEL-STATE.md: rendered tables for humans, JSON block for the code."""
cal = calibration_summary(state)
prior_only = [p for p, v in state["posteriors"].items() if v.get("prior_only")]
diagnostic = [p for p, v in state["posteriors"].items()
if v.get("role") == "diagnostic"]
lines = [f"# MODEL-STATE.md", "",
f"**Revision:** {state['revision']} · "
f"**Covers data through:** {state['covers_through'] or '—'} · "
f"**Emitted:** {date.today()}", ""]
if prior_only:
lines += [f"**{len(prior_only)} of {len(state['posteriors'])} parameters are "
f"`[prior only]`:** {', '.join(sorted(prior_only))}", ""]
if cfg is not None:
frozen = frozen_parameters(cfg, state)
if frozen:
lines += [f"**{len(frozen)} sampled parameter(s) carry no state node** and so "
f"neither learn nor age: {', '.join(frozen)}. They sit in "
f"`params.yaml` at whatever was elicited.", ""]
if diagnostic:
lines += [f"**Tracked but not model inputs:** {', '.join(sorted(diagnostic))} — "
f"compare against what the occupancy curve predicts at the observed "
f"load; a persistent gap means the curve needs refitting.", ""]
sup = calibration_summary(state, quantity="available_productive_hours")
lines += ["## Calibration", "",
"| Metric | Required hours | Available hours |", "|---|---|---|",
f"| Scored periods | {cal.get('n', 0)} | {sup.get('n', 0)} |",
f"| 80% interval hit rate | {cal.get('hit80', '—')} | {sup.get('hit80', '—')} |",
f"| Mean CRPS | {cal.get('crps_mean', '—')} | {sup.get('crps_mean', '—')} |",
f"| PIT mean (0.5 = unbiased) | {cal.get('pit_mean', '—')} | {sup.get('pit_mean', '—')} |",
f"| regime_flag | {cal.get('regime') or 'none'} | {sup.get('regime') or 'none'} |", ""]
if state.get("covers_through") and cal.get("n", 0) == 0:
lines += ["> **The required-hours log is empty.** Forecasts are being archived and not "
"scored. Check Open flags below — this is the state in which the model "
"cannot be shown to be wrong, which is the one state it exists to avoid.", ""]
if state["flags"]:
lines += ["## Open flags", ""] + [f"- {f}" for f in state["flags"][-10:]] + [""]
if readout:
lines += ["## This cycle", "", readout, ""]
lines += ["## Machine state", "",
"Do not hand-edit. Parameters move through `cycle.py` or not at all.", "",
FENCE + "json", json.dumps(state, indent=1, sort_keys=True), FENCE, ""]
with open(path, "w") as fh:
fh.write("\n".join(lines))
def score(state: dict, demand, supply, cfg):
"""Score every archived forecast whose target week has now closed.
Returns (scored, skipped). A tuple rather than a list with a sentinel dict
in it: `05-daily-protocol.md` tells an operator to call this directly, and a
sentinel would poison `state["calibration"]` on the next cycle.
"""
if demand is None:
return [], ["no demand file"]
done = complete_weeks(demand, cfg)
already = {(c["target_week"], c["quantity"]) for c in state["calibration"]}
new, skipped = [], []
for fc in state["forecast_archive"]:
key = (fc["target_week"], fc["quantity"])
if fc["target_week"] not in done or key in already:
continue
try:
if fc["quantity"] == "required_productive_hours":
actual = realized_required_hours(demand, cfg, fc["target_week"])
elif fc["quantity"] == "available_productive_hours":
if supply is None:
raise Unscoreable(f"{fc['target_week']}: no supply file")
actual = realized_supply_hours(supply, fc["target_week"], cfg)
else:
continue
except Unscoreable as exc:
skipped.append(f"{fc['quantity']}: {exc}")
continue
if not np.isfinite(actual) or actual <= 0:
skipped.append(f"{fc['quantity']}: {fc['target_week']} computed as {actual}")
continue
q = fc["q"]
new.append({"target_week": fc["target_week"], "quantity": fc["quantity"],
"issued": fc["issued"], "actual": round(actual, 2),
"pit": round(pit_from_quantiles(q, actual), 4),
"crps": round(crps_from_quantiles(q, actual), 2),
"hit80": interval_hit_from_quantiles(q, actual, 0.80)})
return new, skipped
def calibration_summary(state: dict, quantity="required_productive_hours") -> dict:
rows = [c for c in state["calibration"] if c["quantity"] == quantity]
if not rows:
return {"n": 0}
pits = [r["pit"] for r in rows]
return {"n": len(rows),
"hit80": round(float(np.mean([r["hit80"] for r in rows])), 3),
"crps_mean": round(float(np.mean([r["crps"] for r in rows])), 2),
"pit_mean": round(float(np.mean(pits)), 3),
"regime": regime_flag(pits)}
def daily_statistics(demand, supply, events, cfg) -> list[dict]:
"""Turn raw rows into (parameter, successes, trials, kind) records.
Every rate names the KIND of its denominator so effective_n can be applied.
Contact rate's denominator is transactions -- a bulk correlated count, not
independent trials -- and updating it raw collapses the interval on the
single largest contributor to the model's output variance (44% on the demo,
over exogenous inputs).
"""
stats = []
if demand is not None:
per_seg = demand.groupby(["date", "segment"]).agg(
transactions=("transactions", "max"),
contacts=("contacts_offered", "sum")).reset_index()
for _, r in per_seg.iterrows():
if r["transactions"] > 0:
stats.append({"param": f"CR[{r['segment']}]", "family": "beta",
"successes": float(r["contacts"]),
"trials": float(r["transactions"]),
"kind": "transactions"})
for (dt, seg), g in demand.groupby(["date", "segment"]):
# Carry the channel NAMES with the counts. Passing a bare positional
# list aligned them to whatever order the alpha happened to be
# stored in at bootstrap, so reordering channels in params.yaml sent
# each count to the wrong component silently, and adding one
# truncated the last count through zip.
names = list(cfg["channels"])
counts = [float(g.loc[g["channel"] == c, "contacts_offered"].sum())
for c in names]
if sum(counts) > 0:
stats.append({"param": f"split[{seg}]", "family": "dirichlet",
"counts": counts, "channels": names,
"trials": sum(counts), "kind": "contacts"})
for c in cfg["channels"]:
sub = demand[(demand["channel"] == c) & demand.get(
"aht_seconds", pd.Series(dtype=float)).notna()]
vals = [v for v in sub.get("aht_seconds", []) if v > 0]
if vals:
stats.append({"param": f"AHT[{c}]", "family": "nig", "obs": vals})
# Per DATE, not pooled across the batch. Pooling several days into one
# update makes effective_n treat a fortnight as a single observation
# and leaves n_obs reporting 1 when twenty days were seen.
if cfg["channels"][c].get("kind") == "deferrable":
# Deferrable channels have no occupancy in the model, so an
# occ[<channel>] statistic had nowhere to land: it created an
# untyped node in state that never aged and never cleared, and
# then killed apply_posteriors with KeyError: 'occupancy'.
continue
occ = demand[demand["channel"] == c]
if "occupancy_achieved" in occ and occ["occupancy_achieved"].notna().any():
for _, o in occ.dropna(subset=["occupancy_achieved"]).groupby("date"):
hrs = float(o["productive_hours"].sum())
if hrs > 0:
worked = float((o["occupancy_achieved"] * o["productive_hours"]).sum())
stats.append({"param": f"occ[{c}]", "family": "beta",
"successes": worked, "trials": hrs, "kind": "hours"})
if "transactions_forecast" in demand:
f = demand.groupby(["date", "segment"]).agg(
a=("transactions", "max"), f=("transactions_forecast", "max")).reset_index()
for seg, g in f.groupby("segment"):
ratios = [float(x.a / x.f) for x in g.itertuples()
if pd.notna(x.f) and x.f > 0 and x.a > 0]
if ratios:
stats.append({"param": f"fcst_err[{seg}]", "family": "nig",
"obs": ratios})
if supply is not None:
for _, s in supply.groupby("date"): # per date, not pooled
sched = float(s["scheduled_hours"].sum())
if sched <= 0:
continue
for name, col in (("shrink_planned", "shrink_planned_hours"),
("shrink_unplanned", "shrink_unplanned_hours")):
stats.append({"param": name, "family": "beta",
"successes": float(s[col].sum()),
"trials": sched, "kind": "hours"})
if events is not None and "planned_count" in events.columns:
# `planned_count` is optional in the schema, so its absence must not be
# a crash -- a sparse file written to the documented minimum used to
# raise KeyError and kill the whole cycle.
ev = events.groupby("event")[["count", "planned_count"]].sum(min_count=1)
# The DENOMINATOR for each rate, pinned. graduation_rate is graduates
# over class STARTS, never over seats planned: the simulator composes
# class fill and graduation multiplicatively (`seats x grad`), so using
# seats planned here would charge the same class-fill shortfall twice.
DENOM = {"offer_accepted": ("req_fill_prob", "requisitions opened"),
"class_started": ("class_fill_rate", "seats planned"),
"graduated": ("graduation_rate", "class starts")}
for e, (param, _) in DENOM.items():
if e not in ev.index:
continue
planned = ev.loc[e, "planned_count"]
if pd.notna(planned) and planned > 0:
stats.append({"param": param, "family": "beta",
"successes": float(ev.loc[e, "count"]),
"trials": float(planned), "kind": "events"})
if supply is not None and "separations" in supply.columns \
and "heads_on_roll" in supply.columns:
# `tenured_attrition_weekly` is a WEEKLY hazard and must be learned
# weekly. Grouping by date fed it a DAILY rate, which settled roughly
# 3x low with an interval confidently excluding the truth -- and it errs
# toward inflating supply, the direction that hides a shortfall.
#
# The denominator is agent-weeks at risk: mean heads on roll across the
# week, not the sum of daily headcounts, which counts the same person
# seven times and is this pack's own contact-rate trap reproduced inside
# the extractor that was added to fix the last round's finding.
sup_w = supply.assign(_w=supply["date"].map(week_key))
for wk in complete_weeks(supply, cfg):
s = sup_w[sup_w["_w"] == wk]
days = s["date"].nunique()
heads = float(s["heads_on_roll"].sum()) / max(days, 1) # agent-weeks
sep = float(s["separations"].sum())
if heads > 0 and 0 <= sep <= heads:
stats.append({"param": "attrition_tenured", "family": "beta",
"successes": sep, "trials": heads, "kind": "events"})
return stats
def apply_updates(state: dict, stats: list[dict], as_of: str,
default_lambda: float = 0.98, run_date: str | None = None) -> list[str]:
"""Update every posterior; age the ones with no observation today.
Aging is by ELAPSED DAYS, not per call. A weekly-updated parameter must
discount by lambda^7, not lambda. Parameters with no data in this batch are
aged anyway -- time passed for them too, and a posterior that never widens
while the world moves is the overconfidence this model exists to avoid.
"""
run_date = run_date or as_of
notes = []
seen = set()
# Aging keys off CALENDAR time since the last run, not off the data date.
#
# Keying it off `covers_through` made idle-day aging unreachable: with no
# new data the data date cannot advance, so the gap was always zero and the
# cycle reported "posteriors aged only" while aging nothing. That silently
# canceled the whole point of aging -- a model going stale because no data
# arrived is exactly when its intervals should widen.
#
# Zero still means zero, so running twice in one day changes nothing.
prev_run = state.get("last_run") or state.get("covers_through")
gap = 1.0 if not prev_run else float(max(
(datetime.fromisoformat(run_date).date()
- datetime.fromisoformat(prev_run).date()).days, 0))
for s in stats:
p = s["param"]
seen.add(p)
if p not in state["posteriors"]:
# setdefault wrote an empty {} that had no family, so the ageing
# loop skipped it forever while it inflated the [prior only] count.
notes.append(f"{p}: no node in MODEL-STATE.md, skipped "
f"(bootstrap does not carry this parameter)")
continue
post = state["posteriors"][p]
lam_d = state["lambda_daily"].get(p, default_lambda)
lam = lam_for(lam_d, gap)
if gap == 0 and not s.get("_force"):
# Same day, already absorbed: apply evidence without aging.
lam = 1.0
if s["family"] == "beta":
if "a" not in post:
notes.append(f"{p}: no prior in state, skipped")
continue
n_eff = effective_n(s["trials"], s["kind"])
post["a"], post["b"] = update_beta(post["a"], post["b"],
s["successes"], s["trials"],
lam=lam, effective_n=n_eff)
post["n_obs"] = post.get("n_obs", 0) + 1
elif s["family"] == "dirichlet":
if "alpha" not in post:
notes.append(f"{p}: no prior in state, skipped")
continue
order = post.get("channels")
if order and s.get("channels"):
if set(order) != set(s["channels"]):
notes.append(f"{p}: channel set changed "
f"({sorted(order)} -> {sorted(s['channels'])}), skipped")
continue
idx = {c: i for i, c in enumerate(s["channels"])}
counts = [s["counts"][idx[c]] for c in order] # realign by name
else:
counts = s["counts"]
scale = effective_n(s["trials"], s["kind"]) / max(s["trials"], 1e-9)
post["alpha"] = update_dirichlet(post["alpha"],
[c * scale for c in counts], lam=lam)
post["n_obs"] = post.get("n_obs", 0) + 1
elif s["family"] == "nig":
need = ("mu", "kappa", "alpha", "beta")
if not all(k in post for k in need):
notes.append(f"{p}: no prior in state, skipped")
continue
post["mu"], post["kappa"], post["alpha"], post["beta"] = \
update_normal_invgamma(post["mu"], post["kappa"], post["alpha"],
post["beta"], s["obs"], lam=lam)
post["n_obs"] = post.get("n_obs", 0) + len(s["obs"])
post["last_update"] = as_of
post.pop("prior_only", None)
# Age every unobserved posterior, whatever its family. An earlier version
# guarded on `"a" not in post`, which silently skipped NIG and Dirichlet
# nodes -- so AHT and forecast-error intervals never widened with age, which
# is precisely the overconfidence aging exists to prevent.
if gap > 0:
for p, post in state["posteriors"].items():
if p in seen:
continue
lam_d = state["lambda_daily"].get(p, default_lambda)
lam = lam_for(lam_d, gap)
fam = post.get("family")
if fam == "beta" and "a" in post:
post["a"], post["b"] = discount_only_beta(
post["a"], post["b"], lam_d, gap)
elif fam == "dirichlet" and "alpha" in post:
post["alpha"] = [1 + lam * (a - 1) for a in post["alpha"]]
elif fam == "nig" and "kappa" in post:
post["kappa"] = lam * post["kappa"]
post["alpha"] = lam * post["alpha"]
post["beta"] = lam * post["beta"]
else:
continue
notes.append(f"{p}: no observation, aged {gap:.0f} day(s) only")
return notes
def archive_forecast(state: dict, res, issued: str, first_week_start: str,
weeks: int = 4) -> None:
"""Store a quantile grid for the next `weeks` weeks, both quantities.
Only the near horizon is archived. Weeks further out will be re-forecast
several times before they close, and storing them all would bloat the state
file without adding anything scoreable sooner.
"""
start = pd.Timestamp(first_week_start)
for w in range(min(weeks, res.demand.shape[1])):
wk = week_key(start + pd.Timedelta(weeks=w))
for name, arr in (("required_productive_hours", res.demand[:, w]),
("available_productive_hours", res.supply[:, w])):
state["forecast_archive"] = [
f for f in state["forecast_archive"]
if not (f["target_week"] == wk and f["quantity"] == name)]
state["forecast_archive"].append(
{"issued": issued, "target_week": wk, "quantity": name,
"q": [round(float(x), 2) for x in archive_quantiles(arr)]})
def prune_archive(state: dict, keep_weeks: int = 26) -> None:
"""Drop scored forecasts older than the calibration window."""
scored = {(c["target_week"], c["quantity"]) for c in state["calibration"]}
keep = sorted({f["target_week"] for f in state["forecast_archive"]})[-keep_weeks:]
state["forecast_archive"] = [
f for f in state["forecast_archive"]
if f["target_week"] in keep or (f["target_week"], f["quantity"]) not in scored]
# two scored quantities per week, so keep_weeks x 2 rows
state["calibration"] = state["calibration"][-(keep_weeks * 2):]
def apply_posteriors(cfg: dict, state: dict) -> dict:
"""Push learned posteriors back into a copy of the config, for step 6.
Without this the pack's central claim -- parameters move on evidence, not on
argument -- is a manual splice the operator has to invent every cycle.
Nodes marked `role: diagnostic` are SKIPPED. Achieved occupancy under an
`erlang_curve` is tracked to test the curve, not to feed the model; writing
it back would silently replace the curve with a fixed occupancy and undo the
mechanism it exists to validate.
"""
import copy
from update import lognormal_params
cfg = copy.deepcopy(cfg)
applied, skipped = [], []
for name, post in state.get("posteriors", {}).items():
if post.get("role") == "diagnostic":
skipped.append(name)
continue
fam = post.get("family")
target = None
if name.startswith("CR[") and name[3:-1] in cfg["segments"]:
target = cfg["segments"][name[3:-1]]["contact_rate"]
elif name.startswith("fcst_err[") and name[9:-1] in cfg["segments"]:
target = cfg["segments"][name[9:-1]]["forecast_error"]
elif name.startswith("AHT[") and name[4:-1] in cfg["channels"]:
target = cfg["channels"][name[4:-1]]["aht_seconds"]
elif name.startswith("occ[") and name[4:-1] in cfg["channels"]:
target = cfg["channels"][name[4:-1]]["occupancy"]
elif name.startswith("split[") and name[6:-1] in cfg["segments"]:
seg = cfg["segments"][name[6:-1]]
order = post.get("channels") or list(cfg["channels"])
seg["channel_split"] = {c: float(a) for c, a
in zip(order, post["alpha"])}
applied.append(name)
continue
else:
target = {"shrink_planned": ("supply", "shrinkage_planned"),
"shrink_unplanned": ("supply", "shrinkage_unplanned"),
"req_fill_prob": ("supply", "requisition_fill_prob"),
"class_fill_rate": ("supply", "class_fill_rate"),
"graduation_rate": ("supply", "graduation_rate"),
"attrition_tenured": ("supply", "tenured_attrition_weekly"),
}.get(name)
target = cfg["supply"][target[1]] if target else None
if target is None or not isinstance(target, dict):
skipped.append(name)
continue
if fam == "beta" and "a" in post:
target.clear(); target.update({"a": post["a"], "b": post["b"]})
applied.append(name)
elif fam == "nig" and "mu" in post:
mu, sigma = lognormal_params(post["mu"], post["kappa"],
post["alpha"], post["beta"])
target.clear(); target.update({"mu": mu, "sigma": sigma})
applied.append(name)
else:
skipped.append(name)
cfg["_posteriors_applied"] = applied
cfg["_posteriors_skipped"] = skipped
return cfg
def run_cycle(directory: str, cfg, state_path: str = "MODEL-STATE.md",
as_of: str | None = None, run_date: str | None = None) -> dict:
"""Steps 1-5 and 8. Returns a report; does not itself re-run the model.
`as_of` is the last DATA date; `run_date` is the calendar day the cycle is
run, defaulting to today. They are separate on purpose: evidence is bounded
by the data, aging is bounded by the calendar, and conflating them means a
model that never widens while no data arrives.
"""
state = read_state(state_path)
demand, supply, events = load_intake(directory, cfg)
as_of = as_of or str(demand["date"].max().date())
run_date = run_date or str(date.today())
fails = validate(demand, supply, cfg)
if fails:
return {"as_of": as_of, "halted": True, "validation_failures": fails,
"advisories": advisories(demand, cfg)}
advice = advisories(demand, cfg)
for a in advice:
line = f"{as_of}: {a}"
if line not in state["flags"]:
state["flags"].append(line)
# Only dates not already absorbed may update a posterior. The intake files
# grow -- a daily cycle points at the same directory every day -- so without
# this the same evidence is applied again on every run and the posterior
# concentrates on nothing but repetition. Scoring, by contrast, reads the
# FULL history: it needs completed weeks, and `score` already skips any
# (week, quantity) pair it has recorded before.
ingested = state.get("ingested_through")
if ingested:
cut = pd.Timestamp(ingested)
fresh = [None if d is None else d[d["date"] > cut]
for d in (demand, supply, events)]
if all(d is None or d.empty for d in fresh):
new_demand, new_supply, new_events = None, None, None
else:
new_demand, new_supply, new_events = fresh
else:
new_demand, new_supply, new_events = demand, supply, events
scored, unscored = score(state, demand, supply, cfg)
state["calibration"].extend(scored)
cal = calibration_summary(state)
# Check the regime on BOTH scored series. Running it on required hours only
# left the available-hours PIT unmonitored, which is where a supply-data
# problem shows up.
cal_sup = calibration_summary(state, quantity="available_productive_hours")
for label, summ in (("required", cal), ("available", cal_sup)):
if summ.get("regime"):
state["flags"].append(
f"{as_of}: regime_flag={summ['regime']} on {label} hours — seven of the "
f"last ten actuals in one tail. Surface before updating further; "
f"consider lowering lambda.")
if unscored:
msg = f"{len(unscored)} forecast(s) could not be scored: " + "; ".join(unscored[:4])
if len(unscored) > 4:
msg += f"; and {len(unscored) - 4} more"
line = f"{as_of}: {msg}"
# Do not re-append an identical flag every cycle; ten repeats of one
# message used to evict every other flag from the rendered list.
if line not in state["flags"]:
state["flags"].append(line)
if new_demand is None:
notes = [f"no dates after {ingested} — nothing new to absorb"]
notes += apply_updates(state, [], as_of, run_date=run_date)
else:
notes = apply_updates(
state, daily_statistics(new_demand, new_supply, new_events, cfg), as_of,
run_date=run_date)
state["revision"] += 1
state["covers_through"] = as_of
state["ingested_through"] = as_of
state["last_run"] = run_date
prune_archive(state)
return {"as_of": as_of, "run_date": run_date, "halted": False,
"revision": state["revision"], "scored": scored, "unscored": unscored,
"advisories": advice, "calibration": cal, "calibration_supply": cal_sup,
"update_notes": notes, "state": state}
```
See also
- Wiki:Packs/Probabilistic Staffing — the pack these modules belong to
- Wiki:Packs — the pack index
