1import json
  2import math
  3import random
  4from pathlib import Path
  5import numpy as np
  6from scipy.optimize import brentq, minimize
  7from scipy.special import logsumexp
  8
  9SEED = 17
 10rng = np.random.default_rng(SEED)
 11
 12
 13def logcomb(n, k):
 14    return math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1)
 15
 16
 17def epsilon_exact(N, k, beta):
 18    """Theorem 1 epsilon: solve beta/N * sum(...) = 1."""
 19    if k >= N:
 20        return 1.0
 21    logs = np.array([logcomb(m, k) - logcomb(N, k) for m in range(k, N)])
 22    powers = np.arange(1, N - k + 1)  # N-m
 23    def f(e):
 24        z = logs - powers * math.log1p(-e)
 25        log_value = math.log(beta / N) + logsumexp(z)
 26        # Monotone residual in log-space; root is log_value == 0.
 27        return log_value
 28    lo = k / N
 29    # For sufficiently large k, the monotone equation reaches one at the
 30    # admissible endpoint; the interval bound is then the active certificate.
 31    a = lo + 1e-12
 32    if f(a) >= 0:
 33        return float(lo)
 34    return float(brentq(f, a, 1 - 1e-12, xtol=1e-12))
 35
 36
 37def calibrate(A, y, ids):
 38    """Minimum L2-norm theta satisfying selected linear constraints A theta >= y."""
 39    if len(ids) == 0:
 40        return np.zeros(A.shape[1])
 41    As, ys = A[ids], y[ids]
 42    fun = lambda t: 0.5 * float(t @ t)
 43    jac = lambda t: t
 44    cons = {'type': 'ineq', 'fun': lambda t: As @ t - ys,
 45            'jac': lambda t: As}
 46    out = minimize(fun, np.zeros(A.shape[1]), jac=jac, constraints=cons,
 47                   method='SLSQP', options={'maxiter': 300, 'ftol': 1e-11})
 48    if not out.success:
 49        raise RuntimeError(out.message)
 50    return out.x
 51
 52
 53def compress(A, y, k, mode, seed):
 54    local = np.random.default_rng(seed)
 55    chosen = []
 56    for _ in range(k):
 57        t = calibrate(A, y, chosen)
 58        violation = y - A @ t
 59        candidates = np.array([i for i in range(len(y)) if i not in chosen])
 60        if mode == 'greedy':
 61            i = int(candidates[np.argmax(violation[candidates])])
 62        else:
 63            i = int(local.choice(candidates))
 64        chosen.append(i)
 65    return np.array(chosen), calibrate(A, y, chosen)
 66
 67
 68def run_toy():
 69    # Scenario uncertainty is iid: each constraint is a random affine requirement.
 70    n, nt, k = 400, 50000, 2
 71    train = rng.uniform(-1, 1, size=n)
 72    test = rng.uniform(-1, 1, size=nt)
 73    # Convex envelope-like constraints; two extreme slopes are informative.
 74    A = np.column_stack([train, np.ones(n)])
 75    y = 0.34 + 0.42 * np.abs(train) + rng.normal(0, 0.018, n)
 76    At = np.column_stack([test, np.ones(nt)])
 77    yt = 0.34 + 0.42 * np.abs(test) + rng.normal(0, 0.018, nt)
 78    rows = []
 79    for kk in [2, 4, 8]:
 80        for mode in ['greedy', 'random']:
 81            chosen, theta = compress(A, y, kk, mode, SEED + kk)
 82            train_rate = float(np.mean(y - A @ theta > 1e-8))
 83            fresh_rate = float(np.mean(yt - At @ theta > 0))
 84            rows.append({'k': kk, 'method': mode, 'train_violation': train_rate,
 85                         'fresh_violation': fresh_rate, 'theta': theta.tolist(),
 86                         'selected': chosen.tolist()})
 87    # Non-iid stress test: shifted fresh distribution should not be called certified.
 88    shifted = rng.uniform(0.55, 1.0, size=nt)
 89    Ash = np.column_stack([shifted, np.ones(nt)])
 90    ysh = 0.34 + 0.42 * np.abs(shifted) + rng.normal(0, 0.018, nt)
 91    chosen, theta = compress(A, y, 2, 'greedy', SEED)
 92    shifted_rate = float(np.mean(ysh - Ash @ theta > 0))
 93    return rows, shifted_rate
 94
 95
 96def main():
 97    # Cheap math verification and explicit quantitative predictions.
 98    checks = []
 99    for N in [100, 200, 400, 800]:
100        checks.append({'N': N, 'k': 2, 'beta': 1e-5, 'epsilon': epsilon_exact(N, 2, 1e-5)})
101    beta_rows = [{'beta': b, 'epsilon': epsilon_exact(400, 2, b)}
102                 for b in [1e-3, 1e-5, 1e-7]]
103    k_rows = [{'k': k, 'epsilon': epsilon_exact(400, k, 1e-5)}
104              for k in [0, 1, 2, 4, 8, 16]]
105    toy, shifted = run_toy()
106    result = {'seed': SEED, 'theorem_checks': checks, 'beta_sweep': beta_rows,
107              'k_sweep': k_rows, 'toy': toy, 'shifted_fresh_violation_k2_greedy': shifted}
108    Path('results.json').write_text(json.dumps(result, indent=2))
109    print(json.dumps(result, indent=2))
110
111if __name__ == '__main__':
112    main()