Differentially Private Stay-Point and POI Extraction

from __future__ import annotations

import math
import random
from dataclasses import dataclass
from typing import Iterable, Sequence


@dataclass(frozen=True)
class Fix:
    """One timestamped location fix in metres (projected, not degrees)."""

    t: float
    x: float
    y: float


@dataclass(frozen=True)
class Stay:
    """A detected stay point: centroid, duration, and the fixes it consumed."""

    x: float
    y: float
    start: float
    end: float

    @property
    def dwell(self) -> float:
        return self.end - self.start


def extract_stays(fixes: Sequence[Fix], *, radius_m: float, dwell_s: float) -> list[Stay]:
    """Detect maximal runs that stay within `radius_m` of their own centroid.

    Deliberately greedy and deterministic: a randomised extractor would make the
    per-device stay count data-dependent in a way that is hard to bound, and the
    bound is what the accounting rests on.
    """
    stays: list[Stay] = []
    i = 0
    n = len(fixes)
    while i < n:
        j = i + 1
        cx, cy = fixes[i].x, fixes[i].y
        while j < n:
            k = j - i + 1
            nx = (cx * (k - 1) + fixes[j].x) / k
            ny = (cy * (k - 1) + fixes[j].y) / k
            if max(math.hypot(f.x - nx, f.y - ny) for f in fixes[i:j + 1]) > radius_m:
                break
            cx, cy = nx, ny
            j += 1
        if fixes[j - 1].t - fixes[i].t >= dwell_s:
            stays.append(Stay(cx, cy, fixes[i].t, fixes[j - 1].t))
            i = j
        else:
            i += 1
    return stays


def bound_contribution(
    stays: Sequence[Stay], *, cap: int, seed: int
) -> list[Stay]:
    """Sample at most `cap` stays per device, uniformly and without replacement.

    Uniform sampling — rather than taking the longest or the first `cap` stays —
    keeps the released distribution unbiased with respect to dwell time, which a
    "keep the longest" rule would skew toward workplaces and homes.
    """
    if len(stays) <= cap:
        return list(stays)
    rng = random.Random(seed)
    return rng.sample(list(stays), cap)


def poi_histogram(
    per_device: Iterable[Sequence[Stay]],
    *,
    cell_of: "callable",
    cap: int,
    epsilon: float,
    min_devices: int,
    seed: int = 0,
) -> dict[str, float]:
    """Release a noisy per-cell stay count with a distinct-device floor.

    Sensitivity is `cap`, not 1: one device may contribute up to `cap` stays across
    the grid, so the Laplace scale must be cap/epsilon. Getting this wrong is the
    most common defect in stay-point releases.
    """
    rng = random.Random(seed)
    counts: dict[str, float] = {}
    devices: dict[str, set[int]] = {}
    for device_id, stays in enumerate(per_device):
        for stay in bound_contribution(stays, cap=cap, seed=seed + device_id):
            cell = cell_of(stay)
            counts[cell] = counts.get(cell, 0.0) + 1.0
            devices.setdefault(cell, set()).add(device_id)

    scale = cap / epsilon
    released: dict[str, float] = {}
    for cell, count in counts.items():
        if len(devices[cell]) < min_devices:
            continue                      # suppressed: too few distinct visitors
        u = rng.random() - 0.5
        noise = -scale * math.copysign(1.0, u) * math.log1p(-2 * abs(u))
        released[cell] = max(0.0, count + noise)
    return released

def _validate() -> None:
    # A synthetic device: 20 minutes parked, a drive, then 40 minutes parked.
    fixes = (
        [Fix(t=60.0 * i, x=100.0 + (i % 2), y=200.0) for i in range(20)]
        + [Fix(t=1_200.0 + 60.0 * i, x=100.0 + 400.0 * i, y=200.0 + 300.0 * i) for i in range(8)]
        + [Fix(t=1_700.0 + 60.0 * i, x=3_300.0, y=2_600.0 + (i % 2)) for i in range(40)]
    )
    stays = extract_stays(fixes, radius_m=150.0, dwell_s=300.0)

    # 1. Both stays are found, and the drive between them is not one.
    assert len(stays) == 2, [(s.x, s.dwell) for s in stays]
    assert all(s.dwell >= 300.0 for s in stays)

    # 2. The contribution cap is enforced exactly.
    many = [Stay(float(i), 0.0, 0.0, 600.0) for i in range(30)]
    assert len(bound_contribution(many, cap=8, seed=1)) == 8

    # 3. Sampling must be unbiased with respect to dwell, not "longest first".
    mixed = [Stay(float(i), 0.0, 0.0, 300.0 + 60.0 * i) for i in range(40)]
    picked = bound_contribution(mixed, cap=8, seed=2)
    assert min(s.dwell for s in picked) < 900.0, "sampling must reach short stays"

    # 4. Sensitivity must scale with the cap, not be fixed at 1.
    devices = [[Stay(0.0, 0.0, 0.0, 600.0)] * 5 for _ in range(400)]
    tight = poi_histogram(devices, cell_of=lambda s: "c0", cap=1,
                          epsilon=1.0, min_devices=5, seed=3)
    loose = poi_histogram(devices, cell_of=lambda s: "c0", cap=5,
                          epsilon=1.0, min_devices=5, seed=3)
    assert tight["c0"] < loose["c0"], "a larger cap admits more stays per device"

    # 5. A cell visited by too few distinct devices must be suppressed.
    sparse = [[Stay(0.0, 0.0, 0.0, 600.0)] for _ in range(3)]
    assert poi_histogram(sparse, cell_of=lambda s: "c9", cap=1,
                         epsilon=1.0, min_devices=5) == {}

    print("stay-point extraction: all assertions passed")


_validate()