Central vs Local Differential Privacy for GIS: Calibration, Validation, and Incident Response

import numpy as np
from dataclasses import dataclass
from typing import Dict


class BudgetExhaustedError(Exception):
    """Raised when a query would exceed the allocated epsilon budget."""


@dataclass
class PrivacyBudgetTracker:
    """Tracks cumulative epsilon/delta consumption across spatial queries.

    Sequential composition is assumed: every query that touches the same
    individuals adds its epsilon to the running total. Crossing the cap is a
    hard failure, never a silent clamp -- silent clamping is how pipelines
    leak their formal guarantee without anyone noticing.
    """
    epsilon_total: float
    delta: float = 1e-5
    epsilon_spent: float = 0.0

    def allocate(self, eps: float) -> float:
        if eps <= 0:
            raise ValueError("epsilon per query must be positive")
        if self.epsilon_spent + eps > self.epsilon_total + 1e-12:
            raise BudgetExhaustedError(
                f"Budget exhausted: spent {self.epsilon_spent:.4f}, "
                f"requested {eps:.4f}, cap {self.epsilon_total:.4f}"
            )
        self.epsilon_spent += eps
        return self.epsilon_total - self.epsilon_spent


class SpatialDPEngine:
    """Central- and local-model differential privacy for geospatial workloads.

    The same epsilon buys very different protection in the two models. In the
    central model noise is added once to a trusted aggregate, so the Laplace
    scale is sensitivity / epsilon. In the local model every record is perturbed
    independently before it leaves the device, so the scale is bounded by the
    whole coordinate-domain diameter and is far larger for the same epsilon.
    """

    def __init__(self, budget: PrivacyBudgetTracker, grid_resolution: int = 8,
                 domain_diameter_deg: float = 0.5) -> None:
        # grid_resolution: H3 resolution used to bin coordinates (6-9 typical).
        # domain_diameter_deg: L1 extent of the bounded coordinate domain; this
        # IS the local-DP coordinate sensitivity, so it must be clipped to, not
        # guessed at -- an unbounded domain makes the guarantee meaningless.
        self.budget = budget
        self.grid_res = grid_resolution
        self.domain_diameter = domain_diameter_deg

    @staticmethod
    def _count_sensitivity() -> float:
        # One person changes a disjoint spatial-bin count by at most 1 (L1 = 1).
        return 1.0

    def noise_scale(self, model: str, eps: float) -> float:
        """Expose the Laplace scale so callers can reason about utility first."""
        if model == "central":
            return self._count_sensitivity() / eps
        if model == "local":
            return self.domain_diameter / eps
        raise ValueError("model must be 'central' or 'local'")

    def central_dp_counts(self, bin_counts: np.ndarray, eps: float) -> np.ndarray:
        """Laplace mechanism on PRE-AGGREGATED bins (trusted-curator model)."""
        self.budget.allocate(eps)
        scale = self.noise_scale("central", eps)
        noisy = bin_counts + np.random.laplace(0.0, scale, size=bin_counts.shape)
        # Counts are non-negative; clamp AFTER noise. Post-processing a DP output
        # is free -- it consumes no additional budget.
        return np.maximum(0.0, noisy)

    def local_dp_coords(self, coords: np.ndarray, eps: float) -> np.ndarray:
        """Per-record Laplace perturbation at the edge (untrusted-curator model)."""
        self.budget.allocate(eps)
        # Clip into the declared domain BEFORE perturbing so sensitivity is truly
        # bounded; clipping is part of the mechanism, not optional hardening.
        clipped = np.clip(coords, -self.domain_diameter, self.domain_diameter)
        scale = self.noise_scale("local", eps)  # pure epsilon-LDP, no delta term
        return clipped + np.random.laplace(0.0, scale, size=coords.shape)

    def local_dp_count(self, membership: np.ndarray, eps: float) -> float:
        """Estimate a bin count from per-user randomized response (untrusted model).

        Each user reports a noisy membership bit and the aggregate is de-biased.
        Because noise is injected n times instead of once, this estimator's error
        grows like sqrt(n) / epsilon, versus the central model's constant
        sqrt(2) / epsilon -- that gap IS the formal price of an untrusted curator.
        """
        self.budget.allocate(eps)
        p = 1.0 / (1.0 + np.exp(-eps))  # truthful-response probability under RR
        truthful = np.random.random(membership.shape) < p
        reported = np.where(truthful, membership, 1 - membership)
        # Unbiased de-bias: E[reported] = p*x + (1-p)*(1-x).
        return float((reported.sum() - (1 - p) * membership.size) / (2 * p - 1))

    def displacement_report(self, original: np.ndarray,
                            perturbed: np.ndarray) -> Dict[str, float]:
        """Deterministic utility report: spatial drift in coordinate units."""
        original = np.atleast_2d(np.asarray(original, dtype=float))
        perturbed = np.atleast_2d(np.asarray(perturbed, dtype=float))
        drift = np.linalg.norm(original - perturbed, axis=1)
        return {
            "mean_drift": float(np.mean(drift)),
            "p95_drift": float(np.percentile(drift, 95)),
            "max_drift": float(np.max(drift)),
            "budget_remaining": float(
                self.budget.epsilon_total - self.budget.epsilon_spent),
        }

def _run_validation() -> None:
    rng = np.random.default_rng(7)

    # 1. Same query, equal epsilon: local aggregate error must exceed central.
    #    Central injects one Laplace(1/eps) draw; local injects n randomized
    #    responses, so its error grows ~ sqrt(n)/eps -- the untrusted-curator tax.
    members = (rng.random(5000) < 0.3).astype(float)
    true_count = members.sum()
    eps = 1.0
    central_err = float(np.mean([
        abs((true_count + np.random.laplace(0.0, 1.0 / eps)) - true_count)
        for _ in range(200)]))
    edge = SpatialDPEngine(PrivacyBudgetTracker(epsilon_total=1e6))
    local_err = float(np.mean([
        abs(edge.local_dp_count(members, eps) - true_count) for _ in range(200)]))
    assert local_err > central_err, (local_err, central_err)

    # 2. Central counts stay non-negative and spend exactly the requested budget.
    engine = SpatialDPEngine(PrivacyBudgetTracker(epsilon_total=10.0),
                             domain_diameter_deg=0.5)
    counts = rng.integers(0, 50, size=64).astype(float)
    noisy = engine.central_dp_counts(counts, eps=0.5)
    assert (noisy >= 0).all(), "negative count leaked past the clamp"
    assert abs(engine.budget.epsilon_spent - 0.5) < 1e-9

    # 3. Smaller epsilon must increase expected drift (more privacy, less utility).
    coords = rng.normal(0.0, 0.05, size=(2000, 2))
    loose = SpatialDPEngine(PrivacyBudgetTracker(20.0), domain_diameter_deg=0.5)
    tight = SpatialDPEngine(PrivacyBudgetTracker(20.0), domain_diameter_deg=0.5)
    drift_loose = loose.displacement_report(coords, loose.local_dp_coords(coords, 3.0))
    drift_tight = tight.displacement_report(coords, tight.local_dp_coords(coords, 1.0))
    assert drift_tight["mean_drift"] > drift_loose["mean_drift"]

    # 4. Budget exhaustion is a hard failure, never a silent pass.
    small = SpatialDPEngine(PrivacyBudgetTracker(epsilon_total=1.0))
    small.central_dp_counts(counts, eps=0.9)
    try:
        small.central_dp_counts(counts, eps=0.5)
        raise AssertionError("over-spend was not rejected")
    except BudgetExhaustedError:
        pass

    print("spatial DP engine: all assertions passed")


if __name__ == "__main__":
    _run_validation()