import json import math import random from pathlib import Path import numpy as np from scipy.spatial import cKDTree # Reachable set: unit square. The residual is exactly L-Lipschitz: # r_L(x,y) = margin + L*distance((x,y), center), with an interior minimum. # Thus the proposed implication can be checked without simulator/model noise. def make_cloud(n_side=101): z = np.linspace(0.0, 1.0, n_side) xx, yy = np.meshgrid(z, z, indexing="xy") return np.column_stack([xx.ravel(), yy.ravel()]) def residual(states, L, margin): center = np.array([0.5, 0.5]) return margin + L * np.linalg.norm(states - center, axis=1) def nearest_dist(query, samples): # Exact Euclidean nearest-neighbor distances, with a scalable spatial index. return cKDTree(samples).query(query, k=1)[0] def estimate_L(cloud, margin, L): # A conservative finite-difference estimate on neighboring cloud points. vals = residual(cloud, L, margin) side = int(round(math.sqrt(len(cloud)))) a = vals.reshape(side, side) dx = 1.0 / (side - 1) return float(np.max(np.abs(np.diff(a, axis=1))) / dx) def run_policy(policy, cloud, L, margin, target_n, seed, gamma=2.0): rng = np.random.default_rng(seed) # Four corners provide a reproducible seed set; subsequent points are selected # from the validation/reachable cloud and therefore are valid reachable states. side = int(round(math.sqrt(len(cloud)))) idx = [0, side - 1, len(cloud) - side, len(cloud) - 1] samples = cloud[idx].copy() chosen = np.zeros(len(cloud), dtype=bool); chosen[idx] = True available = np.flatnonzero(~chosen) history = [] eps = 0.04 # below the true minimum margin; threshold is meaningful while len(samples) < target_n: d = nearest_dist(cloud, samples) r = residual(cloud, L, margin) if policy == "random": pick = int(rng.choice(available)) elif policy == "adaptive": # q is the prescribed coverage plus low-margin score. Normalize distance # only so the two terms have comparable numerical scale. q = d + gamma * np.maximum(0.0, eps - r) q[d == 0] = -np.inf pick = int(np.argmax(q)) elif policy == "farthest": d[np.isin(np.arange(len(cloud)), idx)] = -np.inf pick = int(np.argmax(d)) else: raise ValueError(policy) samples = np.vstack([samples, cloud[pick]]) idx.append(pick); chosen[pick] = True available = available[available != pick] delta = float(nearest_dist(cloud, samples).max()) min_sample = float(residual(samples, L, margin).min()) certified = min_sample >= eps and delta < eps / L history.append((len(samples), delta, min_sample, certified)) if certified: break return samples, history def first_cert(history): for n, d, r, ok in history: if ok: return n, d, r return None def main(): np.random.seed(7) random.seed(7) cloud = make_cloud(31) margin = 0.05 eps = margin + 0.10 # Prediction 1: the finite-sample inequality is respected for every tested L, # sample set: true min >= sampled min - L*delta. bound_rows = [] for L in [0.25, 0.5, 1.0, 2.0, 4.0]: for seed in range(4): samples, hist = run_policy("random", cloud, L, margin, 25, seed) delta = float(nearest_dist(cloud, samples).max()) sample_min = float(residual(samples, L, margin).min()) actual_min = margin lower = sample_min - L * delta bound_rows.append({"L": L, "delta": delta, "sample_min": sample_min, "bound": lower, "actual_min": actual_min, "holds": lower <= actual_min + 1e-10}) # Prediction 2: certification turns on around delta = eps/L, and the # observed crossing is reported in units of the predicted threshold. threshold_rows = [] for L in [0.25, 0.5, 1.0]: for policy in ["random", "adaptive"]: crossings = [] for seed in range(5): _, hist = run_policy(policy, cloud, L, margin, 300, seed) c = first_cert(hist) crossings.append(None if c is None else {"n": c[0], "delta": c[1], "ratio": c[1] / (eps / L)}) valid = [x for x in crossings if x is not None] threshold_rows.append({"L": L, "policy": policy, "predicted_delta": eps / L, "n_success": len(valid), "median_n": None if not valid else float(np.median([x["n"] for x in valid])), "median_ratio": None if not valid else float(np.median([x["ratio"] for x in valid]))}) # Prediction 3: at fixed geometry, the conservatism gap scales linearly with L. scaling_rows = [] for L in [0.25, 0.5, 1.0, 2.0, 4.0]: gaps, products = [], [] for seed in range(5): samples, _ = run_policy("random", cloud, L, margin, 20, seed + 100) d = float(nearest_dist(cloud, samples).max()) gap = float(residual(samples, L, margin).min() - margin) gaps.append(gap) products.append(L * d) scaling_rows.append({"L": L, "median_gap": float(np.median(gaps)), "median_L_delta": float(np.median(products)), "gap_over_L": float(np.median(gaps) / L), "max_gap_over_L": float(np.max(gaps) / L)}) # Secondary practical comparison at equal sample budget. compare = [] for L in [1.0]: for n in [12, 25, 50, 80]: row = {"L": L, "n": n} for policy in ["random", "adaptive"]: ds, mins = [], [] for seed in range(6): s, _ = run_policy(policy, cloud, L, margin, n, seed) ds.append(float(nearest_dist(cloud, s).max())) mins.append(float(residual(s, L, margin).min())) row[policy + "_median_delta"] = float(np.median(ds)) row[policy + "_median_min_sample_residual"] = float(np.median(mins)) compare.append(row) out = {"config": {"margin": margin, "epsilon": eps, "cloud_points": len(cloud), "seed": 7}, "predictions": { "bound_inequality": "actual_min >= sample_min - L*delta (reported as lower bound <= actual)", "threshold": "certification when delta < epsilon/L and sample_min >= epsilon", "scaling": "sample_min - actual_min is bounded by, and scales with, L*delta"}, "bound_sweep": bound_rows, "threshold_sweep": threshold_rows, "scaling_sweep": scaling_rows, "equal_budget_comparison": compare} Path("results.json").write_text(json.dumps(out, indent=2)) print(json.dumps({"bound_holds_fraction": float(np.mean([x["holds"] for x in bound_rows])), "threshold": threshold_rows, "scaling": scaling_rows, "comparison": compare}, indent=2)) if __name__ == "__main__": main()