Recipe 12.2: Python Implementation Example

pip install boto3

import json
import logging
import math
import random
import uuid
from collections import Counter, defaultdict
from datetime import date, datetime, timedelta, timezone
from decimal import Decimal
from statistics import mean, stdev

import boto3
from botocore.config import Config

# Structured logging. In production, ship JSON-formatted records
# to CloudWatch Logs Insights for cross-call investigation.
# Consumption data is PHI by association: a daily SKU count for
# orthopedic implants is a derived signal about which surgical
# cases happened on which days. Log structural metadata only
# (run_id, facility_id, segment, sku_count, mean_error, runtime_ms),
# never raw consumption rows tied to identifiable cases, never
# the SKU master with patient or case identifiers attached.
# Full consumption history lives in the encrypted S3 prefix with
# appropriate access controls; only the modeling-ready
# de-identified aggregate flows past the Glue boundary.
logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

# Adaptive retry handles throttling from S3, DynamoDB, EventBridge,
# CloudWatch, and SageMaker. The supply-forecast pipeline is a
# scheduled batch job with no human waiting on the result, so
# longer retries are acceptable than in an interactive bot. Cap
# them anyway so a stuck dependency does not balloon a 30-minute
# weekly run into a multi-hour incident.
BOTO3_RETRY_CONFIG = Config(
    retries={"max_attempts": 5, "mode": "adaptive"})

# Module-level clients. Reused across function calls within the
# pipeline so each call does not pay the connection cost. The
# demo wires up MockTable / MockEventBus / MockCloudWatch via
# run_demo() and never touches these real handles; they are
# staged here so production wiring is a one-line swap. boto3
# client and resource construction is lazy (no network call
# until first use), so the unused handles are free at import.
REGION = "us-east-1"
s3_client          = boto3.client("s3",
                                  region_name=REGION,
                                  config=BOTO3_RETRY_CONFIG)
dynamodb           = boto3.resource("dynamodb",
                                    region_name=REGION,
                                    config=BOTO3_RETRY_CONFIG)
eventbridge_client = boto3.client("events",
                                  region_name=REGION,
                                  config=BOTO3_RETRY_CONFIG)
cloudwatch_client  = boto3.client("cloudwatch",
                                  region_name=REGION,
                                  config=BOTO3_RETRY_CONFIG)
sagemaker_client   = boto3.client("sagemaker",
                                  region_name=REGION,
                                  config=BOTO3_RETRY_CONFIG)

# --- Resource Names ---
# Fill these in with your actual resource names. The demo prints
# what it would write rather than failing if the resources do
# not exist; see run_demo() at the bottom.
CONSUMPTION_HISTORY_BUCKET = "supply-forecast-consumption-history"
SKU_MASTER_BUCKET          = "supply-forecast-sku-master"
MODEL_ARTIFACT_BUCKET      = "supply-forecast-models"
FORECAST_OUTPUT_BUCKET     = "supply-forecast-outputs"
SKU_FORECASTS_TABLE        = "supply-forecast-sku-forecasts"
FORECAST_EVENT_BUS_NAME    = "supply-forecast-events-bus"
CLOUDWATCH_NAMESPACE       = "SupplyForecast"

# --- Versioning ---
# Every forecast record carries the model version active at the
# time of generation. This is how a future audit reconstructs
# which calibration produced which reorder point.
PIPELINE_VERSION       = "supply-forecast-v1.3"
SEGMENTATION_VERSION   = "syntetos-boylan-2025-04"
SAFETY_STOCK_VERSION   = "classical-z-score-v1"

# --- Segmentation Thresholds ---
# Standard four-corner classification from Syntetos, Boylan, and
# Croston. ADI = average demand interval (days between non-zero
# demand observations). CV² = coefficient of variation squared
# of non-zero demand size. The cutoffs (1.32 and 0.49) are the
# commonly cited defaults from the Syntetos-Boylan-Croston
# research. Tune them to your portfolio if your demand shapes
# are unusual.
ADI_THRESHOLD = Decimal("1.32")
CV2_THRESHOLD = Decimal("0.49")

# --- Forecast Horizon ---
# 90 days is a typical operational horizon: long enough to span
# most lead times plus a buffer, short enough that the forecast
# does not lose accuracy from extrapolating too far.
FORECAST_HORIZON_DAYS = 90

# --- Service Level Defaults ---
# Per-SKU target service levels in production come from the SKU
# master and reflect clinical criticality. The defaults below
# stand in for that lookup. A 99% service level is appropriate
# for emergency drugs and critical implants; 95% is fine for
# most consumables; 90% is acceptable for slow-moving non-clinical
# items. Below 90%, you are essentially accepting routine
# stockouts as a cost center.
DEFAULT_SERVICE_LEVEL_BY_SEGMENT = {
    "smooth":             Decimal("0.95"),
    "intermittent":       Decimal("0.97"),
    "erratic":            Decimal("0.95"),
    "lumpy":              Decimal("0.97"),
    "procedure_driven":   Decimal("0.99"),
}

# Z-score lookup for common service levels. Production reads
# the per-SKU service-level target from the SKU master and
# computes the z-score from the inverse normal CDF.
Z_SCORE_BY_SERVICE_LEVEL = {
    Decimal("0.90"): Decimal("1.282"),
    Decimal("0.95"): Decimal("1.645"),
    Decimal("0.97"): Decimal("1.881"),
    Decimal("0.99"): Decimal("2.326"),
}

# --- Lead Time Defaults ---
# Per-SKU lead times in production come from the vendor contract
# and the SKU master. The default below is a placeholder used
# when the SKU master lookup is missing.
DEFAULT_LEAD_TIME_DAYS = 7

# --- Order Quantity Defaults ---
# A simple cycle-stock heuristic stand-in for full Economic Order
# Quantity (EOQ) calculations. Production uses real holding-cost
# and order-cost data per SKU.
DEFAULT_ORDER_CYCLE_DAYS = 14

# --- Synthetic Data ---
# Knobs for the demo's synthetic-consumption generator. The
# ranges produce a small but realistic-shaped portfolio across
# the four demand-pattern segments plus the procedure-driven
# branch.
SYNTHETIC_HISTORY_DAYS = 730     # two years of daily history
SYNTHETIC_FACILITY_ID  = "main-hospital-001"

# Deploy-time guardrail. Any blank resource name is a deploy-time
# bug, not a runtime surprise.
for _name, _value in [
    ("CONSUMPTION_HISTORY_BUCKET", CONSUMPTION_HISTORY_BUCKET),
    ("SKU_MASTER_BUCKET",          SKU_MASTER_BUCKET),
    ("MODEL_ARTIFACT_BUCKET",      MODEL_ARTIFACT_BUCKET),
    ("FORECAST_OUTPUT_BUCKET",     FORECAST_OUTPUT_BUCKET),
    ("SKU_FORECASTS_TABLE",        SKU_FORECASTS_TABLE),
    ("FORECAST_EVENT_BUS_NAME",    FORECAST_EVENT_BUS_NAME),
    ("CLOUDWATCH_NAMESPACE",       CLOUDWATCH_NAMESPACE),
]:
    assert _value, f"{_name} must be set before deploying."

def _to_decimal(value):
    """Convert numeric values to Decimal for DynamoDB-safe writes.

    DynamoDB rejects Python float at the SDK boundary because
    floating-point cannot represent every decimal value exactly.
    Pass everything numeric through this helper before any
    PutItem, BatchWriteItem, or UpdateItem call. Pandas and
    numpy types are not Decimal-friendly out of the box, so the
    helper covers the common cases (int, float, str-of-number,
    None) and lets exotic types fail loudly rather than silently.
    """
    if value is None:
        return None
    if isinstance(value, Decimal):
        return value
    if isinstance(value, bool):
        # bool is a subclass of int in Python; route it explicitly
        # so a True does not become Decimal('1') unexpectedly.
        return value
    if isinstance(value, (int, float)):
        # Round float through string to avoid the float->Decimal
        # repr surprise (Decimal(0.1) is not Decimal('0.1')).
        return Decimal(str(value))
    if isinstance(value, str):
        return Decimal(value)
    raise TypeError(f"Cannot convert {type(value).__name__} to Decimal")

class MockTable:
    """In-memory stand-in for a DynamoDB table.

    Production uses boto3.resource('dynamodb').Table(name). The
    mock supports the operations the demo calls: batch_writer,
    put_item, query, get_item. It is not a complete DynamoDB
    emulation; it covers what this pipeline needs.
    """

    def __init__(self, name):
        self.name           = name
        self.items          = {}    # (pk, sk) -> item dict
        self.write_count    = 0

    class _BatchWriter:
        def __init__(self, table):
            self.table = table

        def __enter__(self):
            return self

        def __exit__(self, *exc):
            return False

        def put_item(self, Item):
            pk = Item["facility_sku"]
            sk = Item["generated_at"]
            self.table.items[(pk, sk)] = dict(Item)
            self.table.write_count += 1

    def batch_writer(self):
        return self._BatchWriter(self)

    def put_item(self, Item):
        pk = Item["facility_sku"]
        sk = Item["generated_at"]
        self.items[(pk, sk)] = dict(Item)
        self.write_count += 1

class MockEventBus:
    """In-memory stand-in for EventBridge.

    Production uses boto3.client('events').put_events(...). The
    mock accumulates events so the demo can show what would have
    been emitted to the chat-events bus at each pipeline stage.
    """

    def __init__(self, name):
        self.name   = name
        self.events = []

    def put_events(self, Entries):
        self.events.extend(Entries)
        return {"FailedEntryCount": 0}

class MockCloudWatch:
    """In-memory stand-in for CloudWatch metrics.

    Production uses boto3.client('cloudwatch').put_metric_data.
    The mock collects metrics by name so the demo can summarize
    what would have been emitted.
    """

    def __init__(self):
        self.metrics = defaultdict(list)

    def put_metric_data(self, Namespace, MetricData):
        for m in MetricData:
            key = f"{Namespace}/{m['MetricName']}"
            self.metrics[key].append(m["Value"])

def generate_synthetic_consumption(seed=42):
    """Generate two years of daily SKU consumption with mixed shapes.

    Real consumption history comes from the materials management
    ledger, the OR case-cart system, and the inpatient pharmacy
    system. The demo replaces that extract with a deterministic
    generator that produces five SKUs covering each demand
    segment, plus an OR case-volume series for the procedure-
    driven branch. Seeding makes the run reproducible.
    """
    rng       = random.Random(seed)
    end_date  = date(2026, 4, 14)
    start_date = end_date - timedelta(days=SYNTHETIC_HISTORY_DAYS - 1)
    rows      = []

    # --- Smooth, high-volume: examination gloves ---
    # Steady weekday volume, lighter weekends, mild trend up.
    for i in range(SYNTHETIC_HISTORY_DAYS):
        d         = start_date + timedelta(days=i)
        dow       = d.weekday()
        base      = 95 + (i / SYNTHETIC_HISTORY_DAYS) * 10  # mild trend
        weekend   = 0.55 if dow >= 5 else 1.0
        noise     = rng.gauss(0, 5)
        qty       = max(0, int(base * weekend + noise))
        rows.append({
            "facility_id": SYNTHETIC_FACILITY_ID,
            "sku_id":      "GLOVE-NITRILE-MED-100CT",
            "date":        d.isoformat(),
            "quantity":    qty,
        })

    # --- Erratic, smooth interval but variable size: IV solution ---
    # Daily demand but with spiky size due to inpatient surges.
    for i in range(SYNTHETIC_HISTORY_DAYS):
        d         = start_date + timedelta(days=i)
        base      = 30
        spike     = 25 if rng.random() < 0.10 else 0
        noise     = rng.gauss(0, 8)
        qty       = max(0, int(base + spike + noise))
        rows.append({
            "facility_id": SYNTHETIC_FACILITY_ID,
            "sku_id":      "IV-SOLUTION-NS-1L",
            "date":        d.isoformat(),
            "quantity":    qty,
        })

    # --- Intermittent: niche surgical kit ---
    # Many zero days with small bursts, 2-3 cases per week.
    for i in range(SYNTHETIC_HISTORY_DAYS):
        d         = start_date + timedelta(days=i)
        if rng.random() < 0.30:
            qty = rng.randint(1, 3)
        else:
            qty = 0
        rows.append({
            "facility_id": SYNTHETIC_FACILITY_ID,
            "sku_id":      "KIT-VASCULAR-ACCESS",
            "date":        d.isoformat(),
            "quantity":    qty,
        })

    # --- Lumpy: specialty respiratory item ---
    # Long zero stretches, occasional batch use.
    for i in range(SYNTHETIC_HISTORY_DAYS):
        d         = start_date + timedelta(days=i)
        if rng.random() < 0.08:
            qty = rng.randint(5, 25)
        else:
            qty = 0
        rows.append({
            "facility_id": SYNTHETIC_FACILITY_ID,
            "sku_id":      "RESPIRATORY-FILTER-SPEC",
            "date":        d.isoformat(),
            "quantity":    qty,
        })

    # --- Procedure-driven: orthopedic staple cartridge ---
    # Consumption tracks scheduled OR cases (generated below).
    case_volume_by_date = {}
    for i in range(SYNTHETIC_HISTORY_DAYS):
        d         = start_date + timedelta(days=i)
        dow       = d.weekday()
        # ORs typically run Mon-Fri; weekend cases are urgent only.
        weekday_cases = max(0, int(rng.gauss(28, 4))) if dow < 5 else \
                        max(0, int(rng.gauss(4, 2)))
        case_volume_by_date[d.isoformat()] = weekday_cases

        # Roughly 0.4 staple cartridges per ortho case;
        # ~25% of total cases are ortho.
        ortho_cases = int(weekday_cases * 0.25)
        usage_per_case = 0.4
        noise          = rng.gauss(0, 0.5)
        qty            = max(0, int(round(ortho_cases * usage_per_case + noise)))
        rows.append({
            "facility_id": SYNTHETIC_FACILITY_ID,
            "sku_id":      "STAPLE-ORTHO-LARGE",
            "date":        d.isoformat(),
            "quantity":    qty,
        })

    return rows, case_volume_by_date

def generate_synthetic_sku_master():
    """Generate a small SKU master with successor maps and lead times.

    Real SKU master data comes from the institutional item catalog
    (typically GHX or the ERP item master) and includes vendor,
    contract, lead time, holding cost, order cost, clinical-
    criticality flag, procedure-driven flag, successor map, and
    per-SKU service-level target. The demo includes the fields
    the pipeline actually reads.
    """
    return {
        "GLOVE-NITRILE-MED-100CT": {
            "vendor_id":            "MEDLINE",
            "lead_time_days":       Decimal("5"),
            "is_procedure_driven":  False,
            "service_level_target": Decimal("0.95"),
            "holding_cost_per_unit": Decimal("0.02"),
            "order_cost":           Decimal("25.00"),
            "successor_sku_id":     None,
        },
        "IV-SOLUTION-NS-1L": {
            "vendor_id":            "BAXTER",
            "lead_time_days":       Decimal("7"),
            "is_procedure_driven":  False,
            "service_level_target": Decimal("0.97"),
            "holding_cost_per_unit": Decimal("0.10"),
            "order_cost":           Decimal("40.00"),
            "successor_sku_id":     None,
        },
        "KIT-VASCULAR-ACCESS": {
            "vendor_id":            "BD",
            "lead_time_days":       Decimal("10"),
            "is_procedure_driven":  False,
            "service_level_target": Decimal("0.97"),
            "holding_cost_per_unit": Decimal("4.50"),
            "order_cost":           Decimal("60.00"),
            "successor_sku_id":     None,
        },
        "RESPIRATORY-FILTER-SPEC": {
            "vendor_id":            "PHILIPS",
            "lead_time_days":       Decimal("14"),
            "is_procedure_driven":  False,
            "service_level_target": Decimal("0.95"),
            "holding_cost_per_unit": Decimal("3.00"),
            "order_cost":           Decimal("75.00"),
            "successor_sku_id":     None,
        },
        "STAPLE-ORTHO-LARGE": {
            "vendor_id":            "STRYKER",
            "lead_time_days":       Decimal("3"),
            "is_procedure_driven":  True,
            "service_level_target": Decimal("0.99"),
            "holding_cost_per_unit": Decimal("18.50"),
            "order_cost":           Decimal("100.00"),
            "successor_sku_id":     None,
        },
    }

def prepare_consumption_data(raw_rows, sku_master, case_volume_by_date=None):
    """Step 1: Clean, fill gaps, reconcile successors, add features.

    See pseudocode Step 1 in the main recipe. The output is a list
    of rows, one per (facility_id, sku_id, date) triple, with
    quantity (zero-filled), calendar features (day_of_week,
    month, is_holiday placeholder, scheduled_cases,
    flu_season_index), and a flag for procedure-driven SKUs that
    the segmentation step uses to override the quantitative
    classification. The scheduled_cases and flu_season_index
    columns teach the row-as-feature-vector pattern the
    pseudocode advocates: the modeling-ready table is visibly
    self-contained with all exogenous signals attached per row.
    """
    if case_volume_by_date is None:
        case_volume_by_date = {}
    # 1a. Group raw transactions by (facility, sku, date) and sum.
    grouped = defaultdict(int)
    for row in raw_rows:
        key = (row["facility_id"], row["sku_id"], row["date"])
        grouped[key] += row["quantity"]

    # 1b. Determine the date range so we can fill gaps. A missing
    # day in the source extract is not the same as zero
    # consumption, but for forecasting purposes the safer
    # default is to assume the SKU was available and not used
    # rather than to leave a gap that the model interprets as
    # continuity. Production has a separate pipeline branch
    # that flags facility-wide outages (system downtime, holiday
    # closures) and excludes those days rather than zero-filling.
    all_dates = sorted({k[2] for k in grouped})
    if not all_dates:
        return []
    start_d = date.fromisoformat(all_dates[0])
    end_d   = date.fromisoformat(all_dates[-1])
    full_dates = [
        (start_d + timedelta(days=i)).isoformat()
        for i in range((end_d - start_d).days + 1)
    ]

    # 1c. Identify the (facility, sku) pairs we need to emit rows for.
    pairs = sorted({(k[0], k[1]) for k in grouped})

    # 1d. Build the cleaned, filled, feature-rich output.
    output = []
    for facility_id, sku_id in pairs:
        # Apply the SKU master successor map: if this SKU has been
        # superseded by a successor, attribute its history to the
        # successor so the new SKU has continuous history. The
        # demo SKUs have no successors; production routinely sees
        # 5-15% of the catalog change identifiers per year.
        master_entry        = sku_master.get(sku_id, {})
        successor_sku_id    = master_entry.get("successor_sku_id")
        target_sku_id       = successor_sku_id if successor_sku_id else sku_id
        is_procedure_driven = bool(master_entry.get("is_procedure_driven", False))

        for d_iso in full_dates:
            qty       = grouped.get((facility_id, sku_id, d_iso), 0)
            d         = date.fromisoformat(d_iso)
            # Scheduled cases: for procedure-driven SKUs, the
            # forecasted OR case volume is the primary signal.
            # For other SKUs it is an informational feature that
            # can improve a multivariate model. Production reads
            # this from the surgical scheduling system's forward-
            # looking extract.
            scheduled_cases = case_volume_by_date.get(d_iso, 0)
            # Flu/respiratory season index: a simple seasonal
            # indicator (0.0-1.0) for respiratory SKUs.
            # Peaks Nov-Feb in the Northern Hemisphere. Production
            # reads from a configurable seasonal-indicator table.
            month_num = d.month
            flu_season_index = (
                1.0 if month_num in (12, 1, 2) else
                0.5 if month_num in (11, 3) else
                0.0
            )
            output.append({
                "facility_id":         facility_id,
                "sku_id":              target_sku_id,
                "date":                d_iso,
                "quantity":            qty,
                "day_of_week":         d.weekday(),
                "month":               d.month,
                "is_holiday":          False,   # production reads from a holiday calendar table
                "scheduled_cases":     scheduled_cases,
                "flu_season_index":    flu_season_index,
                "is_procedure_driven": is_procedure_driven,
            })

    logger.info(
        "Prepared consumption data: %d rows across %d SKUs and %d days",
        len(output), len(pairs), len(full_dates))
    return output

def segment_skus(prepared_rows):
    """Step 2: Classify each SKU by demand pattern.

    Returns a mapping sku_id -> segment_label and a small
    diagnostics dict per SKU (adi, cv2, non_zero_count,
    total_days, mean_size). Production writes the segment
    assignments to a small DynamoDB table so downstream stages
    can read them without recomputing.
    """
    by_sku = defaultdict(list)
    procedure_driven_skus = set()
    for row in prepared_rows:
        by_sku[row["sku_id"]].append(row)
        if row["is_procedure_driven"]:
            procedure_driven_skus.add(row["sku_id"])

    segments    = {}
    diagnostics = {}

    for sku_id, rows in by_sku.items():
        # Procedure-driven override: route to the two-stage model
        # regardless of the quantitative classification. The
        # historical demand pattern of an implant SKU often looks
        # erratic or lumpy, but its underlying signal is OR case
        # volume, which is far more stable.
        if sku_id in procedure_driven_skus:
            segments[sku_id] = "procedure_driven"
            diagnostics[sku_id] = {
                "adi":            None,
                "cv2":            None,
                "total_days":     len(rows),
                "non_zero_count": sum(1 for r in rows if r["quantity"] > 0),
                "mean_size":      None,
                "override_reason": "procedure_driven_flag_in_master",
            }
            continue

        non_zero = [r["quantity"] for r in rows if r["quantity"] > 0]
        total_days     = len(rows)
        non_zero_count = len(non_zero)

        # If a SKU has fewer than ~12 non-zero observations we
        # cannot reliably classify it; route to "lumpy" as the
        # safest default (which uses category-level aggregation)
        # and let the master-data team flag it for review.
        if non_zero_count < 12:
            segments[sku_id] = "lumpy"
            diagnostics[sku_id] = {
                "adi":            None,
                "cv2":            None,
                "total_days":     total_days,
                "non_zero_count": non_zero_count,
                "mean_size":      None,
                "override_reason": "insufficient_non_zero_history",
            }
            continue

        # Average Demand Interval: total days / count of non-zero days.
        adi = Decimal(total_days) / Decimal(non_zero_count)

        # Coefficient of Variation Squared on non-zero demand sizes.
        avg_size = mean(non_zero)
        if avg_size <= 0:
            cv2 = Decimal("0")
        else:
            sd_size = stdev(non_zero) if len(non_zero) > 1 else 0
            cv2     = Decimal(str((sd_size / avg_size) ** 2))

        # Standard four-corner classification.
        if   adi <  ADI_THRESHOLD and cv2 <  CV2_THRESHOLD: label = "smooth"
        elif adi >= ADI_THRESHOLD and cv2 <  CV2_THRESHOLD: label = "intermittent"
        elif adi <  ADI_THRESHOLD and cv2 >= CV2_THRESHOLD: label = "erratic"
        else:                                                label = "lumpy"

        segments[sku_id]    = label
        diagnostics[sku_id] = {
            "adi":            adi,
            "cv2":            cv2,
            "total_days":     total_days,
            "non_zero_count": non_zero_count,
            "mean_size":      _to_decimal(avg_size),
            "override_reason": None,
        }

    seg_counts = Counter(segments.values())
    logger.info(
        "Segmented %d SKUs: %s",
        len(segments),
        ", ".join(f"{seg}={n}" for seg, n in sorted(seg_counts.items())))
    return segments, diagnostics

class SmoothModel:
    """Weekly-seasonal baseline plus level + linear trend.

    Production replaces this with Prophet, ETS, or SARIMA. The
    interface (fit, predict_horizon, sigma) is what matters: any
    model that can produce point forecasts and a residual standard
    deviation drops in.
    """

    def __init__(self):
        self.dow_factor      = {}     # day-of-week -> seasonal factor
        self.level           = 0.0    # final level
        self.trend_per_day   = 0.0    # estimated daily trend
        self.sigma           = 0.0    # residual standard deviation

    def fit(self, daily_quantities, dates):
        # 3a. Estimate level + trend by ordinary least squares on the
        # de-seasonalized series. Compute weekly seasonality first
        # by averaging quantity within each day-of-week.
        n = len(daily_quantities)
        if n == 0:
            return self

        dow_buckets = defaultdict(list)
        for q, d_iso in zip(daily_quantities, dates):
            dow = date.fromisoformat(d_iso).weekday()
            dow_buckets[dow].append(q)

        overall_mean = mean(daily_quantities) if daily_quantities else 0
        self.dow_factor = {
            dow: (mean(qs) / overall_mean) if overall_mean > 0 else 1.0
            for dow, qs in dow_buckets.items()
        }

        # 3b. Linear trend on the de-seasonalized series.
        deseasonalized = [
            q / self.dow_factor.get(date.fromisoformat(d).weekday(), 1.0)
            for q, d in zip(daily_quantities, dates)
        ]
        if n >= 2:
            xs   = list(range(n))
            x_m  = sum(xs) / n
            y_m  = sum(deseasonalized) / n
            num  = sum((xs[i] - x_m) * (deseasonalized[i] - y_m) for i in range(n))
            den  = sum((xs[i] - x_m) ** 2 for i in range(n))
            self.trend_per_day = num / den if den > 0 else 0.0
            self.level         = y_m + self.trend_per_day * (n - 1 - x_m)
        else:
            self.trend_per_day = 0.0
            self.level         = deseasonalized[0]

        # 3c. Residuals against the in-sample fit produce the
        # forecast-error standard deviation we feed into the
        # safety-stock calculation. A real production model uses
        # rolling-origin cross-validation on the validation
        # window, not in-sample residuals.
        residuals = []
        for i, (q, d) in enumerate(zip(daily_quantities, dates)):
            f = self._fitted_value(i, d, n)
            residuals.append(q - f)
        if len(residuals) > 1:
            self.sigma = stdev(residuals)
        else:
            self.sigma = 0.0
        return self

    def _fitted_value(self, t, d_iso, n):
        # Reconstruct the in-sample fitted value for residuals.
        anchor = self.level - self.trend_per_day * (n - 1)
        ds_val = anchor + self.trend_per_day * t
        dow    = date.fromisoformat(d_iso).weekday()
        return ds_val * self.dow_factor.get(dow, 1.0)

    def predict_horizon(self, last_date_iso, horizon_days):
        last_d   = date.fromisoformat(last_date_iso)
        forecasts = []
        for h in range(1, horizon_days + 1):
            future_d = last_d + timedelta(days=h)
            ds_val   = self.level + self.trend_per_day * h
            dow_f    = self.dow_factor.get(future_d.weekday(), 1.0)
            point    = max(0.0, ds_val * dow_f)
            forecasts.append({"date": future_d.isoformat(), "point": point})
        return forecasts

class SBAModel:
    """Syntetos-Boylan Approximation for intermittent demand.

    Croston decomposes intermittent demand into a non-zero
    demand size series and an inter-arrival time series, then
    forecasts each with simple exponential smoothing. Croston's
    original estimator is biased upward; SBA applies a
    correction factor of (1 - alpha/2) to remove it. The
    forecast for any future day is constant: size / interval *
    (1 - alpha/2).
    """

    def __init__(self, alpha=0.1):
        self.alpha           = alpha
        self.size_smoothed   = 0.0  # exp-smoothed non-zero size
        self.interval_smoothed = 0.0  # exp-smoothed inter-arrival time
        self.daily_forecast  = 0.0  # point forecast per day
        self.sigma           = 0.0  # residual std dev

    def fit(self, daily_quantities, dates):
        if not daily_quantities:
            return self

        # Walk through the series tracking days since last non-zero
        # demand. Apply exponential smoothing to the size series
        # and the inter-arrival series, but only update at the
        # non-zero points (this is what makes it Croston-style).
        size_s     = None
        interval_s = None
        gap_count  = 0
        for q in daily_quantities:
            gap_count += 1
            if q > 0:
                if size_s is None:
                    size_s     = float(q)
                    interval_s = float(gap_count)
                else:
                    size_s     = self.alpha * float(q)        + (1 - self.alpha) * size_s
                    interval_s = self.alpha * float(gap_count) + (1 - self.alpha) * interval_s
                gap_count = 0

        if size_s is None or interval_s is None or interval_s == 0:
            self.daily_forecast = 0.0
            self.sigma          = 0.0
            return self

        self.size_smoothed     = size_s
        self.interval_smoothed = interval_s

        # SBA bias correction: multiply by (1 - alpha/2) to
        # remove Croston's upward bias. Without this, intermittent
        # SKUs systematically over-stock.
        self.daily_forecast = (size_s / interval_s) * (1 - self.alpha / 2.0)

        # Residual standard deviation against the constant forecast.
        # This is the per-day demand variability, which the
        # safety-stock formula consumes via z * sqrt(lead_time)
        # * sigma. For accuracy reporting, production uses MASE
        # rather than MAPE because MAPE is undefined on zero-
        # demand days; that is a separate metric from the std-dev
        # used for reorder-point calculation.
        residuals = [q - self.daily_forecast for q in daily_quantities]
        if len(residuals) > 1:
            self.sigma = stdev(residuals)
        else:
            self.sigma = 0.0
        return self

    def predict_horizon(self, last_date_iso, horizon_days):
        last_d = date.fromisoformat(last_date_iso)
        return [
            {
                "date":  (last_d + timedelta(days=h)).isoformat(),
                "point": max(0.0, self.daily_forecast),
            }
            for h in range(1, horizon_days + 1)
        ]

class ProcedureDrivenModel:
    """Two-stage forecast: case_volume * per-case usage rate.

    Forecast the case volume with the SmoothModel, then apply a
    per-case usage rate derived from the historical relationship
    between cases and SKU consumption. The usage rate is the
    SKU master's per-case usage if available; otherwise it is
    estimated from the training window.
    """

    def __init__(self):
        self.case_model       = SmoothModel()
        self.usage_per_case   = 0.0
        self.sigma            = 0.0

    def fit(self, daily_quantities, dates, case_volume_by_date):
        # Estimate the per-case usage rate over the training window.
        # Only count days where there were both cases and consumption
        # to avoid division by zero.
        ratios = []
        for q, d in zip(daily_quantities, dates):
            cases = case_volume_by_date.get(d, 0)
            if cases > 0 and q > 0:
                ratios.append(q / cases)
        self.usage_per_case = mean(ratios) if ratios else 0.0

        # Fit the case-volume model on the case series.
        case_series = [case_volume_by_date.get(d, 0) for d in dates]
        self.case_model.fit(case_series, dates)

        # Residual standard deviation comes from comparing actual
        # SKU consumption to (case_volume * usage_per_case).
        residuals = []
        for q, d in zip(daily_quantities, dates):
            cases    = case_volume_by_date.get(d, 0)
            expected = cases * self.usage_per_case
            residuals.append(q - expected)
        if len(residuals) > 1:
            self.sigma = stdev(residuals)
        else:
            self.sigma = 0.0
        return self

    def predict_horizon(self, last_date_iso, horizon_days):
        case_forecasts = self.case_model.predict_horizon(last_date_iso, horizon_days)
        return [
            {"date": cf["date"], "point": max(0.0, cf["point"] * self.usage_per_case)}
            for cf in case_forecasts
        ]

def train_segment_models(prepared_rows, segments, case_volume_by_date):
    """Step 3: Fit a model per (sku, segment).

    Returns a dict sku_id -> trained_model. The held-out window
    in production is the last 90 days; the demo uses the full
    training history for simplicity. Production also runs a
    quality-gate check against the previous production model
    and rejects regressions, which the demo skips.
    """
    by_sku = defaultdict(list)
    for row in prepared_rows:
        by_sku[row["sku_id"]].append(row)

    trained = {}
    for sku_id, rows in by_sku.items():
        # Sort by date so the model sees a regular series.
        rows.sort(key=lambda r: r["date"])
        quantities = [r["quantity"] for r in rows]
        dates_list = [r["date"]     for r in rows]

        segment = segments[sku_id]

        if segment == "smooth":
            model = SmoothModel().fit(quantities, dates_list)

        elif segment == "procedure_driven":
            model = ProcedureDrivenModel().fit(
                quantities, dates_list, case_volume_by_date)

        else:
            # intermittent / erratic / lumpy all route to SBA in
            # this demo. Production uses Croston for one, SBA for
            # another, and TSB (Teunter-Syntetos-Babai) for SKUs
            # at risk of obsolescence; lumpy SKUs often forecast
            # best at the category level with hierarchical
            # reconciliation. Pick the right tool per segment.
            model = SBAModel(alpha=0.1).fit(quantities, dates_list)

        trained[sku_id] = {
            "model":          model,
            "segment":        segment,
            "last_date":      dates_list[-1] if dates_list else None,
            "model_version":  f"{PIPELINE_VERSION}-{segment}-{datetime.now(timezone.utc).strftime('%Y%m%d')}",
        }
        logger.info(
            "Trained %s model for SKU=%s sigma=%.3f",
            segment, sku_id, getattr(model, "sigma", 0.0))

    return trained

def _z_score_for_service_level(service_level):
    """Look up the z-score for a target service level.

    Production reads per-SKU service-level targets from the SKU
    master and computes the z-score from the inverse normal CDF.
    The lookup table here covers the common service levels and
    falls back to a linear interpolation for in-between values.
    """
    if service_level in Z_SCORE_BY_SERVICE_LEVEL:
        return Z_SCORE_BY_SERVICE_LEVEL[service_level]
    # Simple fallback for values not in the table; production
    # uses scipy.stats.norm.ppf(service_level).
    sorted_levels = sorted(Z_SCORE_BY_SERVICE_LEVEL.keys())
    for i in range(len(sorted_levels) - 1):
        lo, hi = sorted_levels[i], sorted_levels[i + 1]
        if lo <= service_level <= hi:
            frac = (service_level - lo) / (hi - lo)
            return Z_SCORE_BY_SERVICE_LEVEL[lo] + frac * (
                Z_SCORE_BY_SERVICE_LEVEL[hi] - Z_SCORE_BY_SERVICE_LEVEL[lo])
    # Default to 95% if outside the table range.
    return Z_SCORE_BY_SERVICE_LEVEL[Decimal("0.95")]

def _suggest_order_quantity(mean_demand_per_day, sku_master_entry):
    """Heuristic order quantity. Production uses real EOQ math.

    EOQ = sqrt(2 * annual_demand * order_cost / holding_cost_per_unit_per_year).
    The demo uses a simple two-week cycle stock heuristic so the
    pseudocode is readable. The numbers it produces are
    operationally reasonable for steady-demand consumables; they
    are not appropriate for short-shelf-life or high-cost items
    where shelf-life or working-capital constraints dominate.
    """
    return max(1, int(round(mean_demand_per_day * DEFAULT_ORDER_CYCLE_DAYS)))

def generate_sku_forecasts_and_reorder_points(trained, sku_master, run_id):
    """Step 4: Forecast, then translate variance into reorder points.

    Returns a list of forecast records, one per SKU, ready to
    write to DynamoDB.
    """
    records = []
    now_utc = datetime.now(timezone.utc).isoformat()

    for sku_id, info in trained.items():
        model         = info["model"]
        segment       = info["segment"]
        last_date     = info["last_date"]
        model_version = info["model_version"]

        master_entry        = sku_master.get(sku_id, {})
        lead_time_days      = int(master_entry.get(
            "lead_time_days", Decimal(DEFAULT_LEAD_TIME_DAYS)))
        service_level       = master_entry.get(
            "service_level_target",
            DEFAULT_SERVICE_LEVEL_BY_SEGMENT.get(segment, Decimal("0.95")))
        z_score             = _z_score_for_service_level(service_level)

        # Run the model forward over the operational horizon.
        forecast = model.predict_horizon(last_date, FORECAST_HORIZON_DAYS)
        if not forecast:
            logger.warning("Empty forecast for SKU=%s; skipping", sku_id)
            continue

        # Aggregate to lead-time horizon for the reorder calculation.
        # The reorder calculation cares about demand during the
        # lead-time window only, not the full forecast horizon.
        lead_window_points = forecast[:lead_time_days]
        mean_demand_lead   = sum(p["point"] for p in lead_window_points)
        mean_demand_per_day = (
            sum(p["point"] for p in forecast) / len(forecast))

        # Safety stock: classical formula, z * sqrt(lead_time) * sigma_daily.
        # The square root reflects that the standard deviation of
        # the sum of independent daily demands grows with the
        # square root of the number of days, not linearly.
        sigma_daily   = Decimal(str(getattr(model, "sigma", 0.0)))
        safety_stock  = (z_score
                         * Decimal(str(math.sqrt(lead_time_days)))
                         * sigma_daily)

        # Compute reorder point in float (Decimal only at the
        # DynamoDB write boundary). The safety-stock formula
        # operates at unit-count precision; sub-unit decimals
        # are meaningless for physical inventory.
        reorder_point = int(round(
            mean_demand_lead + float(safety_stock)))

        order_quantity = _suggest_order_quantity(
            mean_demand_per_day, master_entry)

        # Aggregate forecast totals over the horizon for the
        # operational dashboard. Lower / upper bounds come from a
        # simple z-score band around the point forecast for the
        # demo; production uses the model's native prediction
        # interval (Prophet's yhat_lower / yhat_upper, the
        # statsmodels confidence interval, or DeepAR's quantile
        # outputs).
        horizon_total = sum(p["point"] for p in forecast)
        z_for_band    = float(z_score)
        sigma_total   = float(sigma_daily) * math.sqrt(FORECAST_HORIZON_DAYS)
        lower_bound   = max(0, int(round(horizon_total - z_for_band * sigma_total)))
        upper_bound   =          int(round(horizon_total + z_for_band * sigma_total))

        record = {
            "facility_sku":          f"{SYNTHETIC_FACILITY_ID}#{sku_id}",
            "facility_id":           SYNTHETIC_FACILITY_ID,
            "sku_id":                sku_id,
            "segment":               segment,
            "forecast_date_from":    forecast[0]["date"],
            "forecast_horizon_days": _to_decimal(FORECAST_HORIZON_DAYS),
            "mean_demand_horizon":   _to_decimal(int(round(horizon_total))),
            "lower_bound":           _to_decimal(lower_bound),
            "upper_bound":           _to_decimal(upper_bound),
            "lead_time_days":        _to_decimal(lead_time_days),
            "service_level_target":  _to_decimal(service_level),
            "reorder_point":         _to_decimal(reorder_point),
            "order_quantity":        _to_decimal(order_quantity),
            "sigma_daily":           _to_decimal(sigma_daily),
            "generated_at":          now_utc,
            "run_id":                run_id,
            "model_version":         model_version,
            "pipeline_version":      PIPELINE_VERSION,
            "safety_stock_version":  SAFETY_STOCK_VERSION,
            "segmentation_version":  SEGMENTATION_VERSION,
        }
        records.append(record)
        logger.info(
            "Forecast SKU=%s segment=%s reorder_point=%d order_qty=%d",
            sku_id, segment, reorder_point, order_quantity)

    return records

def load_forecasts_to_dynamodb(records, table, event_bus, cloudwatch):
    """Step 5: Batched writes to the served table plus a CURRENT pointer.

    The boto3 resource-level batch_writer() chunks into 25-item
    batches and retries UnprocessedItems with exponential backoff
    internally. The explicit chunking below is for clarity in a
    pedagogical mock; production typically just hands the full
    list of records to a single batch_writer() context.
    """
    if not records:
        return 0

    written = 0
    chunk   = 25
    for i in range(0, len(records), chunk):
        batch = records[i:i + chunk]
        with table.batch_writer() as bw:
            for record in batch:
                bw.put_item(Item=record)
                written += 1

    # Upsert a CURRENT pointer per SKU so consumers can do a
    # single GetItem instead of querying-and-sorting client-side.
    # The CURRENT row references the latest generated_at and
    # carries the same reorder_point / order_quantity / bounds.
    for record in records:
        current_pointer = dict(record)
        current_pointer["generated_at"] = "CURRENT"
        current_pointer["points_to"]    = record["generated_at"]
        table.put_item(Item=current_pointer)

    # Emit an EventBridge event so downstream consumers (the ERP
    # integration job, the materials management dashboard, the
    # forecast-monitoring job) know a new run completed. The
    # event payload deliberately carries no PHI: just the run
    # identifier, the facility, the SKU count, and the timestamp.
    event_bus.put_events(Entries=[{
        "Source":        "supply.forecast",
        "DetailType":    "ForecastRunCompleted",
        "EventBusName":  FORECAST_EVENT_BUS_NAME,
        "Time":          datetime.now(timezone.utc),
        "Detail":        json.dumps({
            "facility_id":   SYNTHETIC_FACILITY_ID,
            "run_id":        records[0]["run_id"],
            "sku_count":     len(records),
            "pipeline_version": PIPELINE_VERSION,
        }),
    }])

    cloudwatch.put_metric_data(
        Namespace=CLOUDWATCH_NAMESPACE,
        MetricData=[
            {"MetricName": "RecordsWritten",
             "Value":      float(written),
             "Unit":       "Count"},
            {"MetricName": "SkuCount",
             "Value":      float(len(records)),
             "Unit":       "Count"},
        ])

    logger.info("Loaded %d forecast records into %s", written, table.name)
    return written

def run_supply_forecast_pipeline(table, event_bus, cloudwatch):
    """End-to-end pipeline orchestration.

    The demo wires up synthetic data; production starts with an
    S3 read of the consumption-history Parquet partition for the
    target facility and date range.
    """
    run_id = str(uuid.uuid4())
    print(f"\n=== Supply Forecast Pipeline run_id={run_id} ===\n")

    # --- Generate synthetic input data (production reads from S3) ---
    raw_rows, case_volume_by_date = generate_synthetic_consumption()
    sku_master                    = generate_synthetic_sku_master()
    print(f"[input] {len(raw_rows)} raw rows, "
          f"{len(case_volume_by_date)} case-volume days, "
          f"{len(sku_master)} SKUs in master")

    # --- Step 1: Prepare and clean the consumption data ---
    print("\n[step 1] prepare_consumption_data")
    prepared = prepare_consumption_data(raw_rows, sku_master, case_volume_by_date)
    print(f"  -> {len(prepared)} prepared rows")

    # --- Step 2: Segment SKUs by demand pattern ---
    print("\n[step 2] segment_skus")
    segments, diagnostics = segment_skus(prepared)
    for sku_id, segment in segments.items():
        diag = diagnostics[sku_id]
        adi  = f"{diag['adi']:.2f}" if diag["adi"] is not None else "n/a"
        cv2  = f"{diag['cv2']:.2f}" if diag["cv2"] is not None else "n/a"
        print(f"  SKU={sku_id:30s} segment={segment:18s} adi={adi:>5s} cv2={cv2:>5s}")

    # --- Step 3: Train a model per segment ---
    print("\n[step 3] train_segment_models")
    trained = train_segment_models(prepared, segments, case_volume_by_date)
    for sku_id, info in trained.items():
        sigma = getattr(info["model"], "sigma", 0.0)
        print(f"  SKU={sku_id:30s} segment={info['segment']:18s} "
              f"sigma={sigma:.3f} model_version={info['model_version']}")

    # --- Step 4: Generate forecasts and reorder points ---
    print("\n[step 4] generate_sku_forecasts_and_reorder_points")
    records = generate_sku_forecasts_and_reorder_points(trained, sku_master, run_id)
    for r in records:
        print(f"  SKU={r['sku_id']:30s} segment={r['segment']:18s} "
              f"reorder_point={int(r['reorder_point']):>5d} "
              f"order_qty={int(r['order_quantity']):>5d} "
              f"horizon={int(r['mean_demand_horizon']):>6d} "
              f"({int(r['lower_bound']):>5d}..{int(r['upper_bound']):>5d})")

    # --- Step 5: Load to DynamoDB ---
    print("\n[step 5] load_forecasts_to_dynamodb")
    written = load_forecasts_to_dynamodb(records, table, event_bus, cloudwatch)
    print(f"  -> wrote {written} records "
          f"(plus {len(records)} CURRENT pointers)")
    print(f"  -> emitted {len(event_bus.events)} EventBridge events")
    print(f"  -> emitted CloudWatch metrics: "
          f"{sorted(cloudwatch.metrics.keys())}")

    return records

def run_demo():
    """Run the pipeline end-to-end against the in-memory mocks.

    No AWS resources are touched; every external dependency is
    a mock. Useful for sanity-checking the forecasting logic and
    the reorder-point math before wiring to real services.
    """
    table      = MockTable(SKU_FORECASTS_TABLE)
    event_bus  = MockEventBus(FORECAST_EVENT_BUS_NAME)
    cloudwatch = MockCloudWatch()

    records = run_supply_forecast_pipeline(table, event_bus, cloudwatch)

    print("\n=== Sample CURRENT record ===")
    sample_pk = records[0]["facility_sku"]
    sample    = table.items[(sample_pk, "CURRENT")]
    # Convert Decimal to str for readable JSON output.
    def _decimalify(o):
        if isinstance(o, Decimal):
            return str(o)
        if isinstance(o, datetime):
            return o.isoformat()
        return o
    print(json.dumps(sample, default=_decimalify, indent=2))

    return records

if __name__ == "__main__":
    run_demo()