1import json
  2import math
  3import numpy as np
  4
  5
  6def dg_moments(q, limit=200):
  7    z = np.arange(-limit, limit + 1)
  8    p = (1-q)/(1+q) * q**np.abs(z)
  9    p /= p.sum()
 10    t = np.abs(z)
 11    return p, z, float((p*t).sum()), float((p*t*t).sum() - (p*t).sum()**2)
 12
 13
 14def math_check():
 15    q = 0.37
 16    p, z, mu_emp, var_emp = dg_moments(q)
 17    mu = 2*q/(1-q*q)
 18    # Exact expectation and derivative for a quadratic objective at one coordinate.
 19    target, x = 2, -1
 20    def expected(eta):
 21        qq = math.exp(eta)
 22        pp = (1-qq)/(1+qq) * qq**np.abs(z)
 23        return float((pp * (x + z - target)**2).sum())
 24    eta = math.log(q)
 25    eps = 2e-5
 26    fd = (expected(eta+eps)-expected(eta-eps))/(2*eps)
 27    F = (x + z - target)**2
 28    score_cov = float((p*(F-(p*F).sum())*(np.abs(z)-mu)).sum())
 29    fisher = var_emp
 30    return {
 31        "q": q, "pmf_sum": float(p.sum()), "analytic_mu": mu,
 32        "empirical_mu": mu_emp, "empirical_fisher": fisher,
 33        "score_identity_cov": score_cov, "finite_difference_dE_deta": fd,
 34        "identity_abs_error": abs(score_cov-fd),
 35        "natural_gradient": score_cov/(fisher+1e-12),
 36    }
 37
 38
 39def utilities(losses):
 40    # Higher utility means better (lower) loss; centered rank utilities.
 41    ranks = np.argsort(np.argsort(losses))
 42    u = (len(losses)-1 - 2*ranks).astype(float) / max(1, len(losses)-1)
 43    return u - u.mean()
 44
 45
 46def dg_mutation(rng, q, d):
 47    mag = rng.geometric(1-q, size=d) - 1
 48    sign = np.where(rng.random(d) < .5, -1, 1)
 49    return mag * sign
 50
 51
 52def optimize(method, seed, targets, weights, iterations=90, K=24):
 53    rng = np.random.default_rng(seed)
 54    d = len(targets)
 55    x = np.zeros(d, dtype=int)
 56    q = np.full(d, .35)
 57    eta = np.log(q)
 58    h = np.zeros(d)
 59    beta, rho = .85, .22
 60    best = float(np.sum(weights*(x-targets)**2))
 61    curve = []
 62    for _ in range(iterations):
 63        zs = []
 64        if method == "double_geometric":
 65            for _k in range(K): zs.append(dg_mutation(rng, q, d))
 66        elif method == "gaussian_es":
 67            # Standard Gaussian ES, rounded to the same integer search space.
 68            zs = [np.rint(rng.normal(0, 1.5, size=d)).astype(int) for _k in range(K)]
 69        elif method == "fixed_radius":
 70            # Simple integer random search with independent {-1,0,+1} mutations.
 71            zs = [rng.integers(-1, 2, size=d) for _k in range(K)]
 72        losses = np.array([np.sum(weights*(np.clip(x+z, -10, 10)-targets)**2) for z in zs])
 73        i = int(np.argmin(losses))
 74        if losses[i] < best: best = float(losses[i])
 75        curve.append(best)
 76        # All methods use the same rank-based selection signal; only DG adapts q.
 77        if method == "double_geometric":
 78            U = utilities(losses)
 79            T = np.array([np.abs(z) for z in zs], dtype=float)
 80            mu = 2*q/(1-q*q)
 81            var = T.var(axis=0)
 82            # Cov(utility,T)/Var(T) estimates negative objective natural gradient.
 83            g_utility = np.mean(U[:,None]*(T-mu), axis=0)/(var+1e-6)
 84            h = beta*h + (1-beta)*g_utility
 85            eta = np.clip(eta + rho*h, math.log(.08), math.log(.80))
 86            q = np.exp(eta)
 87        # ES update: select the best candidate, matching a minimal zero-order optimizer.
 88        x = np.clip(x + zs[i], -10, 10).astype(int)
 89    return {"final_best": best, "curve": curve, "final_x": x.tolist(), "final_q": q.tolist()}
 90
 91
 92def main():
 93    check = math_check()
 94    # Six layerwise integer controls with heterogeneous scales/targets.
 95    targets = np.array([4, -3, 2, 5, -2, 3])
 96    weights = np.array([1., 1.4, .8, 1.2, .9, 1.1])
 97    methods = ["double_geometric", "gaussian_es", "fixed_radius"]
 98    results = {m: [optimize(m, s, targets, weights) for s in range(8)] for m in methods}
 99    summary = {}
100    for m, rr in results.items():
101        finals = np.array([r["final_best"] for r in rr])
102        # Equal evaluation budget: report final objective and fraction reaching exact target.
103        summary[m] = {"final_mean": float(finals.mean()), "final_std": float(finals.std()),
104                       "exact_fraction": float(np.mean(finals < 1e-9)),
105                       "median": float(np.median(finals))}
106    out = {"math_check": check, "summary": summary, "raw": results,
107           "evaluations_per_run": 90*24, "targets": targets.tolist()}
108    with open("results.json", "w") as f: json.dump(out, f, indent=2)
109    print(json.dumps({"math_check": check, "summary": summary}, indent=2))
110
111if __name__ == "__main__": main()