Choosing KDE bandwidth for mobility hotspots

import numpy as np
from sklearn.model_selection import GridSearchCV
from sklearn.neighbors import KernelDensity
from scipy.stats import spearmanr


def silverman_bandwidth(xy: np.ndarray) -> float:
    """Rule-of-thumb bandwidth. Useful ONLY as an upper bound on clustered data."""
    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.")
    n = len(xy)
    if n < 10:
        raise ValueError("Too few points for a bandwidth estimate.")
    # Median absolute deviation is more robust than the standard deviation
    # here, which is itself dominated by inter-cluster separation.
    sigma = np.median(np.abs(xy - np.median(xy, axis=0))) * 1.4826
    return float(1.06 * sigma * n ** (-1 / 6))


def cv_bandwidth(
    xy: np.ndarray,
    candidates: np.ndarray | None = None,
    sample: int = 20_000,
    cv: int = 4,
    seed: int = 0,
) -> dict:
    """
    Likelihood cross-validated bandwidth on a stratified subsample.

    Parameters
    ----------
    xy : np.ndarray
        (n, 2) PROJECTED metric coordinates.
    sample : int
        Subsample size. Full-data CV is O(n²) per fold; 20 000 points reaches
        within a few per cent of the full-data optimum on typical urban sets
        and runs in seconds instead of hours.

    Returns
    -------
    dict
        Best bandwidth, the Silverman upper bound, and the score curve.

    Raises
    ------
    ValueError
        On bad shape or too few points.
    """
    xy = np.asarray(xy, dtype=float)
    if xy.ndim != 2 or xy.shape[1] != 2:
        raise ValueError("xy must be an (n, 2) array.")
    if len(xy) < 200:
        raise ValueError("Cross-validation needs at least a few hundred points.")

    rng = np.random.default_rng(seed)
    idx = rng.choice(len(xy), size=min(sample, len(xy)), replace=False)
    sub = xy[idx]

    if candidates is None:
        upper = silverman_bandwidth(xy)
        # Search well below the rule of thumb: on clustered data the optimum
        # is routinely a quarter of it, and a grid that starts at it misses.
        candidates = np.geomspace(upper / 20.0, upper, 18)

    search = GridSearchCV(
        KernelDensity(kernel="gaussian"),
        {"bandwidth": candidates},
        cv=cv,
        n_jobs=-1,
    ).fit(sub)

    best = float(search.best_params_["bandwidth"])
    return {
        "bandwidth_m": best,
        "silverman_m": silverman_bandwidth(xy),
        "candidates_m": candidates,
        "scores": search.cv_results_["mean_test_score"],
        "at_grid_edge": bool(best <= candidates[0] * 1.01 or best >= candidates[-1] * 0.99),
    }


def hotspot_stability(xy_a: np.ndarray, xy_b: np.ndarray, grid: np.ndarray,
                      bandwidth_m: float, top_k: int = 50) -> float:
    """
    Rank correlation of the top-k cells between two independent periods.

    This is the test that matters: a bandwidth whose hotspot ranking does not
    reproduce on a held-out week is describing the sample, not the city.
    """
    def dens(pts):
        kde = KernelDensity(bandwidth=bandwidth_m, kernel="gaussian").fit(pts)
        return np.exp(kde.score_samples(grid))

    da, db = dens(xy_a), dens(xy_b)
    top = np.argsort(da)[-top_k:]
    rho, _ = spearmanr(da[top], db[top])
    return float(rho)

def validate_bandwidth(res: dict, stability_rho: float,
                       cell_size_m: float, min_rho: float = 0.7) -> None:
    """Four assertions that catch the ways a bandwidth choice goes wrong."""
    bw = res["bandwidth_m"]

    # 1. The optimum must be interior to the search grid, or it is a boundary
    #    artefact and the grid needs widening.
    assert not res["at_grid_edge"], (
        f"CV optimum {bw:.0f} m sits at the edge of the search grid — widen it"
    )

    # 2. Cell size must resolve the kernel. Below ~3 cells per bandwidth the
    #    surface is a blocky sampling of the kernel rather than the kernel.
    assert cell_size_m <= bw / 3, (
        f"cell {cell_size_m:.0f} m is too coarse for bandwidth {bw:.0f} m"
    )

    # 3. Reproducibility on held-out data. This is the real test.
    assert stability_rho >= min_rho, (
        f"hotspot ranking correlates only {stability_rho:.2f} across periods "
        "— widen the bandwidth until it reproduces"
    )

    # 4. Sanity against the rule of thumb: CV should land BELOW it on
    #    clustered data. Landing above suggests the points are not clustered
    #    and a rule of thumb would have been fine.
    if bw > res["silverman_m"]:
        print(f"note: CV bandwidth {bw:.0f} m exceeds Silverman "
              f"{res['silverman_m']:.0f} m — is this data actually clustered?")
    print(f"OK — bandwidth {bw:.0f} m, stability ρ = {stability_rho:.2f}")