1import json
  2import math
  3import random
  4from pathlib import Path
  5
  6import numpy as np
  7from scipy.spatial import cKDTree
  8
  9
 10# Reachable set: unit square. The residual is exactly L-Lipschitz:
 11# r_L(x,y) = margin + L*distance((x,y), center), with an interior minimum.
 12# Thus the proposed implication can be checked without simulator/model noise.
 13
 14def make_cloud(n_side=101):
 15    z = np.linspace(0.0, 1.0, n_side)
 16    xx, yy = np.meshgrid(z, z, indexing="xy")
 17    return np.column_stack([xx.ravel(), yy.ravel()])
 18
 19
 20def residual(states, L, margin):
 21    center = np.array([0.5, 0.5])
 22    return margin + L * np.linalg.norm(states - center, axis=1)
 23
 24
 25def nearest_dist(query, samples):
 26    # Exact Euclidean nearest-neighbor distances, with a scalable spatial index.
 27    return cKDTree(samples).query(query, k=1)[0]
 28
 29
 30def estimate_L(cloud, margin, L):
 31    # A conservative finite-difference estimate on neighboring cloud points.
 32    vals = residual(cloud, L, margin)
 33    side = int(round(math.sqrt(len(cloud))))
 34    a = vals.reshape(side, side)
 35    dx = 1.0 / (side - 1)
 36    return float(np.max(np.abs(np.diff(a, axis=1))) / dx)
 37
 38
 39def run_policy(policy, cloud, L, margin, target_n, seed, gamma=2.0):
 40    rng = np.random.default_rng(seed)
 41    # Four corners provide a reproducible seed set; subsequent points are selected
 42    # from the validation/reachable cloud and therefore are valid reachable states.
 43    side = int(round(math.sqrt(len(cloud))))
 44    idx = [0, side - 1, len(cloud) - side, len(cloud) - 1]
 45    samples = cloud[idx].copy()
 46    chosen = np.zeros(len(cloud), dtype=bool); chosen[idx] = True
 47    available = np.flatnonzero(~chosen)
 48    history = []
 49    eps = 0.04  # below the true minimum margin; threshold is meaningful
 50    while len(samples) < target_n:
 51        d = nearest_dist(cloud, samples)
 52        r = residual(cloud, L, margin)
 53        if policy == "random":
 54            pick = int(rng.choice(available))
 55        elif policy == "adaptive":
 56            # q is the prescribed coverage plus low-margin score. Normalize distance
 57            # only so the two terms have comparable numerical scale.
 58            q = d + gamma * np.maximum(0.0, eps - r)
 59            q[d == 0] = -np.inf
 60            pick = int(np.argmax(q))
 61        elif policy == "farthest":
 62            d[np.isin(np.arange(len(cloud)), idx)] = -np.inf
 63            pick = int(np.argmax(d))
 64        else:
 65            raise ValueError(policy)
 66        samples = np.vstack([samples, cloud[pick]])
 67        idx.append(pick); chosen[pick] = True
 68        available = available[available != pick]
 69        delta = float(nearest_dist(cloud, samples).max())
 70        min_sample = float(residual(samples, L, margin).min())
 71        certified = min_sample >= eps and delta < eps / L
 72        history.append((len(samples), delta, min_sample, certified))
 73        if certified:
 74            break
 75    return samples, history
 76
 77
 78def first_cert(history):
 79    for n, d, r, ok in history:
 80        if ok:
 81            return n, d, r
 82    return None
 83
 84
 85def main():
 86    np.random.seed(7)
 87    random.seed(7)
 88    cloud = make_cloud(31)
 89    margin = 0.05
 90    eps = margin + 0.10
 91
 92    # Prediction 1: the finite-sample inequality is respected for every tested L,
 93    # sample set: true min >= sampled min - L*delta.
 94    bound_rows = []
 95    for L in [0.25, 0.5, 1.0, 2.0, 4.0]:
 96        for seed in range(4):
 97            samples, hist = run_policy("random", cloud, L, margin, 25, seed)
 98            delta = float(nearest_dist(cloud, samples).max())
 99            sample_min = float(residual(samples, L, margin).min())
100            actual_min = margin
101            lower = sample_min - L * delta
102            bound_rows.append({"L": L, "delta": delta, "sample_min": sample_min,
103                               "bound": lower, "actual_min": actual_min,
104                               "holds": lower <= actual_min + 1e-10})
105
106    # Prediction 2: certification turns on around delta = eps/L, and the
107    # observed crossing is reported in units of the predicted threshold.
108    threshold_rows = []
109    for L in [0.25, 0.5, 1.0]:
110        for policy in ["random", "adaptive"]:
111            crossings = []
112            for seed in range(5):
113                _, hist = run_policy(policy, cloud, L, margin, 300, seed)
114                c = first_cert(hist)
115                crossings.append(None if c is None else {"n": c[0], "delta": c[1], "ratio": c[1] / (eps / L)})
116            valid = [x for x in crossings if x is not None]
117            threshold_rows.append({"L": L, "policy": policy,
118                                   "predicted_delta": eps / L,
119                                   "n_success": len(valid),
120                                   "median_n": None if not valid else float(np.median([x["n"] for x in valid])),
121                                   "median_ratio": None if not valid else float(np.median([x["ratio"] for x in valid]))})
122
123    # Prediction 3: at fixed geometry, the conservatism gap scales linearly with L.
124    scaling_rows = []
125    for L in [0.25, 0.5, 1.0, 2.0, 4.0]:
126        gaps, products = [], []
127        for seed in range(5):
128            samples, _ = run_policy("random", cloud, L, margin, 20, seed + 100)
129            d = float(nearest_dist(cloud, samples).max())
130            gap = float(residual(samples, L, margin).min() - margin)
131            gaps.append(gap)
132            products.append(L * d)
133        scaling_rows.append({"L": L, "median_gap": float(np.median(gaps)),
134                             "median_L_delta": float(np.median(products)),
135                             "gap_over_L": float(np.median(gaps) / L),
136                             "max_gap_over_L": float(np.max(gaps) / L)})
137
138    # Secondary practical comparison at equal sample budget.
139    compare = []
140    for L in [1.0]:
141        for n in [12, 25, 50, 80]:
142            row = {"L": L, "n": n}
143            for policy in ["random", "adaptive"]:
144                ds, mins = [], []
145                for seed in range(6):
146                    s, _ = run_policy(policy, cloud, L, margin, n, seed)
147                    ds.append(float(nearest_dist(cloud, s).max()))
148                    mins.append(float(residual(s, L, margin).min()))
149                row[policy + "_median_delta"] = float(np.median(ds))
150                row[policy + "_median_min_sample_residual"] = float(np.median(mins))
151            compare.append(row)
152
153    out = {"config": {"margin": margin, "epsilon": eps, "cloud_points": len(cloud), "seed": 7},
154           "predictions": {
155               "bound_inequality": "actual_min >= sample_min - L*delta (reported as lower bound <= actual)",
156               "threshold": "certification when delta < epsilon/L and sample_min >= epsilon",
157               "scaling": "sample_min - actual_min is bounded by, and scales with, L*delta"},
158           "bound_sweep": bound_rows, "threshold_sweep": threshold_rows,
159           "scaling_sweep": scaling_rows, "equal_budget_comparison": compare}
160    Path("results.json").write_text(json.dumps(out, indent=2))
161    print(json.dumps({"bound_holds_fraction": float(np.mean([x["holds"] for x in bound_rows])),
162                      "threshold": threshold_rows, "scaling": scaling_rows,
163                      "comparison": compare}, indent=2))
164
165
166if __name__ == "__main__":
167    main()