Tuning DBSCAN eps and min_samples for stay points

import numpy as np
import pandas as pd
from sklearn.neighbors import NearestNeighbors


def suggest_eps(
    xy: np.ndarray,
    min_samples: int = 5,
    position_error_m: float = 8.0,
    merge_ceiling_m: float = 70.0,
) -> dict:
    """
    Derive an eps candidate from the k-distance knee and bound-check it.

    Parameters
    ----------
    xy : np.ndarray
        (n, 2) PROJECTED metric coordinates of the candidate stationary fixes.
        Run this on slow fixes only — including motorway travel drowns the
        knee in a uniform distance distribution.
    position_error_m : float
        1-sigma horizontal error. eps below about 2x this fragments a single
        stationary vehicle, so the floor is a property of the receiver rather
        than of the tuning.

    Returns
    -------
    dict
        Suggested eps, the raw knee, and whether it violates either bound.

    Raises
    ------
    ValueError
        On too few points or wrongly shaped input.
    """
    xy = np.asarray(xy, dtype=float)
    if xy.ndim != 2 or xy.shape[1] != 2:
        raise ValueError("xy must be an (n, 2) array of projected coordinates.")
    if len(xy) <= min_samples:
        raise ValueError(f"Need more than {min_samples} points to compute a k-distance.")

    nn = NearestNeighbors(n_neighbors=min_samples).fit(xy)
    d, _ = nn.kneighbors(xy)
    kdist = np.sort(d[:, -1])                    # distance to the k-th neighbour

    # Knee by maximum distance from the chord joining the curve's endpoints —
    # the standard Kneedle construction, and robust enough here.
    i = np.arange(len(kdist), dtype=float)
    i_n = (i - i[0]) / max(i[-1] - i[0], 1e-9)
    k_n = (kdist - kdist[0]) / max(kdist[-1] - kdist[0], 1e-9)
    knee_idx = int(np.argmax(k_n - i_n))
    knee = float(kdist[knee_idx])

    floor = 2.0 * position_error_m
    suggested = float(np.clip(knee, floor, merge_ceiling_m))
    return {
        "knee_m": knee,
        "suggested_eps_m": suggested,
        "below_floor": knee < floor,
        "above_ceiling": knee > merge_ceiling_m,
        "floor_m": floor,
        "ceiling_m": merge_ceiling_m,
    }


def suggest_min_samples(median_interval_s: float, min_dwell_s: float = 180.0,
                        safety: float = 0.6) -> int:
    """
    min_samples from the sampling rate and the dwell you intend to detect.

    A stop of min_dwell_s at this cadence contains min_dwell_s / interval
    fixes. Requiring all of them is brittle — a couple will be rejected by
    quality filters — so take a fraction of that count.
    """
    if median_interval_s <= 0:
        raise ValueError("median_interval_s must be positive.")
    expected = min_dwell_s / median_interval_s
    return max(3, int(round(expected * safety)))

def validate_parameters(sweep: pd.DataFrame, truth_stops: int | None = None) -> dict:
    """
    Prefer a plateau over a best cell.

    sweep : DataFrame with columns eps, min_samples, n_clusters.
    A parameter set that only works at one value has been fitted to the
    sweep rather than to the behaviour, and will drift with the next month.
    """
    piv = sweep.pivot(index="eps", columns="min_samples", values="n_clusters")

    # 1. Find the cell whose 3x3 neighbourhood varies least — that is the plateau.
    var = piv.rolling(3, center=True, min_periods=3).std().T.rolling(
        3, center=True, min_periods=3).std().T
    if var.notna().sum().sum() == 0:
        raise ValueError("Sweep grid too small to identify a plateau; widen it.")
    eps_star, ms_star = var.stack().idxmin()

    out = {"eps": float(eps_star), "min_samples": int(ms_star),
           "n_clusters": int(piv.loc[eps_star, ms_star])}

    # 2. The plateau should be within reach of the labelled count, if known.
    if truth_stops:
        err = abs(out["n_clusters"] - truth_stops) / truth_stops
        assert err < 0.25, (
            f"plateau gives {out['n_clusters']} stops against {truth_stops} "
            "labelled — the stop definition, not the parameters, is the problem"
        )
        out["count_error"] = err

    # 3. A sharp optimum is a warning sign, not a success.
    neighbourhood = piv.loc[
        piv.index[max(piv.index.get_loc(eps_star) - 1, 0):piv.index.get_loc(eps_star) + 2]
    ]
    spread = neighbourhood.std().std()
    assert spread < 0.25 * out["n_clusters"], (
        "cluster count varies sharply around the chosen cell — no plateau exists"
    )
    return out