Hausdorff vs LCSS for partial trajectory matching

import numpy as np
from scipy.spatial.distance import cdist


def discrete_hausdorff(a: np.ndarray, b: np.ndarray) -> float:
    """
    Symmetric discrete Hausdorff distance between two projected paths.

    Parameters
    ----------
    a, b : np.ndarray
        (n, 2) and (m, 2) arrays of METRIC coordinates. Degrees give a
        latitude-dependent answer that cannot be compared between cities.

    Returns
    -------
    float
        The greatest distance from any vertex of either path to the nearest
        vertex of the other. Dominated by the single worst point, by design.

    Raises
    ------
    ValueError
        On empty or wrongly shaped input.
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    if a.ndim != 2 or a.shape[1] != 2 or b.ndim != 2 or b.shape[1] != 2:
        raise ValueError("Both paths must be (n, 2) coordinate arrays.")
    if a.size == 0 or b.size == 0:
        raise ValueError("Cannot compare an empty path.")

    d = cdist(a, b)                       # O(n·m) memory — block long trips
    return float(max(d.min(axis=1).max(), d.min(axis=0).max()))


def lcss(
    a: np.ndarray,
    b: np.ndarray,
    eps_m: float = 40.0,
    delta: int = 20,
) -> float:
    """
    Longest Common SubSequence similarity in [0, 1].

    Parameters
    ----------
    eps_m : float
        Spatial tolerance in metres. Two points match if they are within it.
        Set from positional error — roughly 2-3x the p95 fix error — not by
        eye, or the score measures your tolerance rather than the routes.
    delta : int
        Maximum index offset allowed when matching, which stops the algorithm
        pairing the start of one trip with the end of the other. Also caps
        the cost at O(n·delta) rather than O(n·m).

    Returns
    -------
    float
        Matched length divided by the shorter path's length. 1.0 means the
        shorter path is entirely contained in the longer one within eps.
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    if a.size == 0 or b.size == 0:
        raise ValueError("Cannot compare an empty path.")
    if eps_m <= 0:
        raise ValueError("eps_m must be positive.")

    n, m = len(a), len(b)
    prev = np.zeros(m + 1, dtype=np.int32)
    cur = np.zeros(m + 1, dtype=np.int32)

    for i in range(1, n + 1):
        lo, hi = max(1, i - delta), min(m, i + delta)
        cur[:] = 0
        for j in range(lo, hi + 1):
            if np.hypot(a[i - 1, 0] - b[j - 1, 0], a[i - 1, 1] - b[j - 1, 1]) <= eps_m:
                cur[j] = prev[j - 1] + 1
            else:
                cur[j] = max(prev[j], cur[j - 1])
        prev, cur = cur, prev

    # Normalise by the SHORTER path: a 2 km trip fully contained in a 20 km
    # one should score 1.0, not 0.1. Normalising by the longer path instead
    # makes every short trip look dissimilar to everything.
    return float(prev[m] / max(min(n, m), 1))

def validate_measures(a: np.ndarray, b: np.ndarray, eps_m: float = 40.0) -> dict:
    """Property tests that catch the usual implementation errors."""
    # 1. Identity. A path against itself is Hausdorff 0 and LCSS 1.
    assert discrete_hausdorff(a, a) == 0.0, "Hausdorff of a path with itself is not 0"
    assert abs(lcss(a, a, eps_m) - 1.0) < 1e-9, "LCSS of a path with itself is not 1"

    # 2. Symmetry. Both measures are symmetric; an asymmetric result means
    #    the one-sided Hausdorff was returned, or delta is applied unevenly.
    assert abs(discrete_hausdorff(a, b) - discrete_hausdorff(b, a)) < 1e-6
    assert abs(lcss(a, b, eps_m) - lcss(b, a, eps_m)) < 1e-9

    # 3. Tolerance monotonicity. LCSS must not fall as eps rises.
    lo, hi = lcss(a, b, eps_m), lcss(a, b, eps_m * 2)
    assert hi >= lo - 1e-9, "LCSS decreased with a looser tolerance"

    # 4. Containment. A strict prefix of a path must score LCSS 1.0.
    assert abs(lcss(a[: len(a) // 2], a, eps_m) - 1.0) < 1e-9, (
        "a prefix does not score 1.0 — the normalisation uses the longer "
        "path instead of the shorter one"
    )
    return {"hausdorff_m": discrete_hausdorff(a, b), "lcss": lcss(a, b, eps_m)}