Speed percentile profiles for corridor benchmarking

import numpy as np
import pandas as pd


def corridor_speed_profile(
    df: pd.DataFrame,
    corridor_col: str = "link_id",
    dir_col: str = "direction",
    time_col: str = "t",
    speed_col: str = "speed_ms",
    bin_minutes: int = 15,
    min_obs: int = 30,
    percentiles: tuple = (15, 50, 85),
) -> pd.DataFrame:
    """
    Percentile speed profile per corridor, direction and time-of-day bin.

    Parameters
    ----------
    min_obs : int
        Bins with fewer observations are emitted with NaN percentiles and a
        flag rather than dropped. A missing bin and a quiet bin look identical
        downstream unless the distinction is carried explicitly.

    Returns
    -------
    pd.DataFrame
        One row per (corridor, direction, bin) with the requested percentiles,
        the observation count and a sufficient-sample flag.

    Raises
    ------
    ValueError
        On missing columns, an empty frame, or non-positive speeds.
    """
    required = {corridor_col, dir_col, time_col, speed_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if df.empty:
        raise ValueError("Input DataFrame is empty.")
    if (df[speed_col] < 0).any():
        raise ValueError("Negative speeds present — check the derivation.")

    d = df.copy()
    # Bin edges anchored to midnight LOCAL time: a corridor profile is a
    # statement about the working day, and 08:00 has to mean 08:00 there.
    minutes = d[time_col].dt.hour * 60 + d[time_col].dt.minute
    d["tod_bin"] = (minutes // bin_minutes) * bin_minutes
    d["day_type"] = np.where(d[time_col].dt.dayofweek < 5, "weekday", "weekend")

    keys = [corridor_col, dir_col, "day_type", "tod_bin"]
    g = d.groupby(keys, observed=True)[speed_col]

    out = g.agg(n="size").reset_index()
    for p in percentiles:
        out[f"p{p}_ms"] = g.quantile(p / 100.0).to_numpy()
    out["sufficient"] = out["n"] >= min_obs
    out.loc[~out["sufficient"], [f"p{p}_ms" for p in percentiles]] = np.nan
    return out


def add_congestion_ratio(
    profile: pd.DataFrame,
    corridor_col: str = "link_id",
    dir_col: str = "direction",
    offpeak_hours: tuple = (22, 5),
) -> pd.DataFrame:
    """
    Express each bin's median as a fraction of the corridor's own free-flow.

    Free-flow is the p85 of the OFF-PEAK distribution for that corridor and
    direction, which is the standard definition and the only one that makes
    two different roads comparable.
    """
    lo, hi = offpeak_hours
    h = profile["tod_bin"] // 60
    off = profile[(h >= lo) | (h < hi)]
    if off.empty:
        raise ValueError(
            "No off-peak bins — widen offpeak_hours or the free-flow baseline "
            "will be taken from congested data and every ratio will look fine."
        )
    ff = (off.groupby([corridor_col, dir_col], observed=True)["p85_ms"]
             .max().rename("freeflow_ms"))
    out = profile.merge(ff, on=[corridor_col, dir_col], how="left")
    out["congestion_ratio"] = out["p50_ms"] / out["freeflow_ms"]
    return out

def validate_profile(profile: pd.DataFrame, min_obs: int = 30) -> dict:
    """Four checks that catch a profile built on too little or the wrong data."""
    # 1. Percentiles must be ordered. Out-of-order values mean the groupby
    #    keys and the quantile results were misaligned.
    ok = profile.dropna(subset=["p15_ms", "p50_ms", "p85_ms"])
    assert (ok["p15_ms"] <= ok["p50_ms"]).all() and (ok["p50_ms"] <= ok["p85_ms"]).all(), \
        "percentiles out of order — the aggregation is misaligned"

    # 2. Coverage. If most bins are under the sample floor, the corridor is
    #    too finely divided or the fleet is too small for this resolution.
    cov = profile["sufficient"].mean()
    assert cov > 0.6, (
        f"only {cov:.0%} of bins have ≥{min_obs} observations — widen the bin "
        "or aggregate corridors before believing this profile"
    )

    # 3. Congestion ratio in a sane range. Above 1.1 means the free-flow
    #    baseline came from congested data.
    if "congestion_ratio" in profile:
        r = profile["congestion_ratio"].dropna()
        assert r.max() < 1.15, f"ratio up to {r.max():.2f} — free-flow baseline is too low"
        assert r.min() > 0.02, "ratio near zero — stationary vehicles included in the median"

    # 4. Directions must both be present, or the profile is one carriageway.
    assert profile["direction"].nunique() >= 2, "only one direction in the profile"
    return {"bins": len(profile), "coverage": float(cov)}