1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6SEED = 2136
  7rng = np.random.default_rng(SEED)
  8
  9
 10def sample_deficits(n_trials, n_candidates, dim, kappa=2.0, curvature=1.0, batch=2000):
 11    """Q(u*)-max Q(U_i) for U uniform in [-1,1]^d and Q=-curvature*||u||^k."""
 12    out = np.empty(n_trials)
 13    for lo in range(0, n_trials, batch):
 14        m = min(batch, n_trials - lo)
 15        u = rng.uniform(-1.0, 1.0, size=(m, n_candidates, dim))
 16        r = np.linalg.norm(u, axis=-1)
 17        out[lo:lo+m] = curvature * np.min(r ** kappa, axis=1)
 18    return out
 19
 20
 21def fit_log_slope(ns, means):
 22    return float(np.polyfit(np.log(ns), np.log(means), 1)[0])
 23
 24
 25def main():
 26    # Prediction 1: E[Delta_N] scales as N^(-kappa/d).
 27    # Prediction 2: doubling N multiplies error by 2^(-kappa/d).
 28    # We check both d=1 and d=2, with kappa=2.
 29    n_trials = 12000
 30    kappa = 2.0
 31    scaling = {}
 32    for dim in (1, 2):
 33        ns = np.array([4, 8, 16, 32, 64, 128, 256])
 34        means = np.array([sample_deficits(n_trials, int(n), dim, kappa).mean() for n in ns])
 35        slope = fit_log_slope(ns, means)
 36        ratios = means[1:] / means[:-1]
 37        predicted_slope = -kappa / dim
 38        predicted_ratio = 2 ** predicted_slope
 39        scaling[str(dim)] = {
 40            "N": ns.tolist(), "mean_deficit": means.tolist(),
 41            "observed_loglog_slope": slope, "predicted_loglog_slope": predicted_slope,
 42            "observed_doubling_ratio_mean": float(ratios.mean()),
 43            "predicted_doubling_ratio": predicted_ratio,
 44            "doubling_ratios": ratios.tolist(),
 45        }
 46
 47    # Prediction 3: scheduler N = ceil((C / ((1-beta) eps))^(d/kappa))
 48    # has N proportional to eps^(-d/kappa), and achieves the requested error scale.
 49    dim = 2
 50    beta = 0.9
 51    epsilons = np.array([0.04, 0.02, 0.01, 0.005])
 52    n_cal = 64
 53    cal = sample_deficits(30000, n_cal, dim, kappa)
 54    C_hat = float(cal.mean() * n_cal ** (kappa / dim))
 55    scheduled_ns = np.ceil((C_hat / ((1-beta) * epsilons)) ** (dim/kappa)).astype(int)
 56    scheduled_ns = np.maximum(scheduled_ns, 1)
 57    measured = np.array([sample_deficits(16000, int(n), dim, kappa).mean() for n in scheduled_ns])
 58    value_bound_scale = measured / (1-beta)
 59    eps_slope = fit_log_slope(epsilons, scheduled_ns)
 60    # For this geometry the asymptotic constant is known: area(pi r^2)/(area(box)) = pi r^2/4,
 61    # hence E[min ||U||^2] ~ 4/(pi N). This is also a check that C_hat is calibrated.
 62    C_exact = 4.0 / math.pi
 63    expected_at_schedule = C_exact * scheduled_ns ** (-dim / kappa) / (1-beta)
 64    scheduler = {
 65        "beta": beta, "C_hat": C_hat, "C_exact_asymptotic": C_exact,
 66        "C_hat_over_C_exact": C_hat / C_exact, "epsilons": epsilons.tolist(),
 67        "scheduled_N": scheduled_ns.tolist(), "measured_mean_delta": measured.tolist(),
 68        "measured_delta_over_1_minus_beta": value_bound_scale.tolist(),
 69        "expected_value_proxy_at_schedule": expected_at_schedule.tolist(),
 70        "mean_expected_proxy_over_epsilon": float(np.mean(expected_at_schedule / epsilons)),
 71        "observed_loglog_N_vs_epsilon_slope": eps_slope,
 72        "predicted_loglog_N_vs_epsilon_slope": -dim/kappa,
 73        "fraction_meeting_epsilon_bound": float(np.mean(value_bound_scale <= epsilons)),
 74    }
 75
 76    # Mini comparison: fixed pools versus the adaptive pool over heterogeneous states.
 77    # Curvature changes the local C, so fixed N is deliberately wasteful on easy states.
 78    n_states = 3000
 79    curvatures = rng.uniform(0.35, 2.0, size=n_states)
 80    eps = 0.015
 81    fixed_ns = [8, 32, 118, 128]
 82    rows = []
 83    # Use common state curvature and independent candidate draws; report action error
 84    # and Bellman-value error proxy Delta/(1-beta), with mean Q evaluations.
 85    for n in fixed_ns:
 86        ds = np.array([sample_deficits(1, n, dim, kappa, float(a))[0] for a in curvatures])
 87        rows.append({"method": f"fixed_N={n}", "mean_evaluations": n,
 88                     "mean_delta": float(ds.mean()), "mean_value_error_proxy": float((ds/(1-beta)).mean()),
 89                     "fraction_value_proxy_le_eps": float(np.mean(ds/(1-beta) <= eps))})
 90    # Estimate C per state from a cheap N=8 pilot, then use the formula and cap at 128.
 91    pilot_n = 8
 92    pilot = np.array([sample_deficits(1, pilot_n, dim, kappa, float(a))[0] for a in curvatures])
 93    local_C = pilot * pilot_n ** (kappa/dim)
 94    adaptive_n = np.ceil((local_C / ((1-beta)*eps)) ** (dim/kappa)).astype(int)
 95    adaptive_n = np.clip(adaptive_n, pilot_n, 128)
 96    adaptive_ds = np.array([sample_deficits(1, int(n), dim, kappa, float(a))[0]
 97                            for a, n in zip(curvatures, adaptive_n)])
 98    rows.append({"method": "adaptive_scheduler", "mean_evaluations": float(adaptive_n.mean()),
 99                 "mean_delta": float(adaptive_ds.mean()), "mean_value_error_proxy": float((adaptive_ds/(1-beta)).mean()),
100                 "fraction_value_proxy_le_eps": float(np.mean(adaptive_ds/(1-beta) <= eps)),
101                 "pool_min": int(adaptive_n.min()), "pool_median": float(np.median(adaptive_n)),
102                 "pool_max": int(adaptive_n.max())})
103
104    result = {"seed": SEED, "theory": {"action": "uniform box, Q(u)=-||u||^2", "kappa": kappa},
105              "scaling_checks": scaling, "scheduler_scaling_check": scheduler,
106              "mini_comparison": rows}
107    Path("results.json").write_text(json.dumps(result, indent=2))
108    print(json.dumps(result, indent=2))
109
110
111if __name__ == "__main__":
112    main()