Detecting harsh braking and acceleration events

import numpy as np
import pandas as pd
from scipy.signal import savgol_filter


def detect_harsh_events(
    df: pd.DataFrame,
    time_col: str = "t",
    speed_col: str = "speed_ms",
    smooth_window: int = 5,
    min_duration_s: float = 0.8,
    a_ref: float = 3.6,
    v_ref: float = 8.0,
) -> pd.DataFrame:
    """
    Detect sustained harsh braking and acceleration events.

    The threshold is a_ref at v_ref and falls with speed as
    a_thresh(v) = a_ref * v_ref / max(v, v_ref), which approximates the
    achievable envelope of a passenger vehicle and, just as importantly,
    rises above the low-speed noise floor.

    Parameters
    ----------
    df : pd.DataFrame
        One vehicle's fixes, tz-aware time_col, speed in m/s from PROJECTED
        coordinates.
    min_duration_s : float
        An event must be sustained for this long. Single-sample spikes are
        measurement artefacts, not driving.

    Returns
    -------
    pd.DataFrame
        One row per event: start, end, duration, peak acceleration, mean
        speed, and the threshold it crossed.

    Raises
    ------
    ValueError
        On missing columns, an empty frame, or an even smoothing window.
    """
    required = {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 smooth_window % 2 == 0:
        raise ValueError("smooth_window must be odd.")

    d = df.sort_values(time_col).reset_index(drop=True).copy()
    t = d[time_col].to_numpy()
    dt = np.r_[np.nan, np.diff(t).astype("timedelta64[ns]").astype(float) / 1e9]
    if np.nanmedian(dt) <= 0:
        raise ValueError("Non-positive sampling interval — fix the clock first.")

    v = d[speed_col].to_numpy(dtype=float)
    if len(v) >= smooth_window:
        # Smooth the SPEED, not the acceleration: differentiating a smoothed
        # signal preserves the event edge, smoothing a differentiated one
        # spreads a spike across the whole window.
        v = savgol_filter(v, smooth_window, 2, mode="interp")

    a = np.gradient(v, np.nanmedian(dt))
    thresh = a_ref * v_ref / np.maximum(v, v_ref)

    over = np.abs(a) > thresh
    d["accel_ms2"] = a
    d["threshold_ms2"] = thresh

    # ── group consecutive over-threshold samples into events ──────────
    grp = (over != np.r_[False, over[:-1]]).cumsum()
    events = []
    for _, run in d[over].groupby(grp[over]):
        dur = (run[time_col].iloc[-1] - run[time_col].iloc[0]).total_seconds()
        if dur < min_duration_s:
            continue                     # impulsive: an artefact, not a manoeuvre
        peak = run["accel_ms2"].abs().max()
        events.append({
            "start": run[time_col].iloc[0],
            "end": run[time_col].iloc[-1],
            "duration_s": dur,
            "kind": "braking" if run["accel_ms2"].mean() < 0 else "acceleration",
            "peak_ms2": float(peak),
            "mean_speed_ms": float(run[speed_col].mean()),
            "threshold_ms2": float(run["threshold_ms2"].mean()),
        })
    return pd.DataFrame(events)


def events_per_100km(events: pd.DataFrame, distance_m: float) -> dict:
    """
    Normalise the count by exposure. A raw count compares nothing: a van that
    drove 400 km will out-score one that drove 40 no matter how it was driven.
    """
    if distance_m <= 0:
        raise ValueError("distance_m must be positive to normalise.")
    km = distance_m / 1000.0
    out = {"distance_km": km}
    for kind in ("braking", "acceleration"):
        n = int((events.get("kind", pd.Series(dtype=str)) == kind).sum())
        out[f"{kind}_per_100km"] = round(n * 100.0 / km, 2)
    return out

def validate_events(events: pd.DataFrame, df: pd.DataFrame,
                    distance_m: float, speed_col: str = "speed_ms") -> None:
    """Assertions that catch a detector measuring noise rather than driving."""
    rates = events_per_100km(events, distance_m)

    # 1. Plausible rate. Above ~40 per 100 km the detector is firing on noise;
    #    below ~0.2 it is almost certainly never firing at all.
    total = rates["braking_per_100km"] + rates["acceleration_per_100km"]
    assert 0.2 < total < 40, f"{total:.1f} events per 100 km is not plausible driving"

    # 2. No zero-speed events. Braking at a standstill is a receiver artefact.
    if len(events):
        assert (events["mean_speed_ms"] > 1.0).all(), (
            "events detected below walking pace — the low-speed noise floor "
            "is above the threshold; raise v_ref"
        )
        # 3. Durations must exceed the minimum by construction.
        assert (events["duration_s"] >= 0.79).all(), "sub-threshold duration leaked through"

    # 4. Speed must actually change across a braking event.
    for _, e in events[events["kind"] == "braking"].head(20).iterrows():
        win = df[(df["t"] >= e["start"]) & (df["t"] <= e["end"])][speed_col]
        assert win.iloc[0] - win.iloc[-1] > 1.0, (
            "braking event with no net speed loss — the derivative is noise"
        )
    print(f"OK — {len(events)} events, {total:.1f} per 100 km")