1import json
  2import math
  3import random
  4from pathlib import Path
  5import numpy as np
  6
  7EPS = 1e-12
  8
  9def sigmoid(z):
 10    return 1.0 / (1.0 + np.exp(-np.clip(z, -40.0, 40.0)))
 11
 12def binary_logloss(p, y):
 13    p = float(np.clip(p, 1e-12, 1.0 - 1e-12))
 14    return float(-y * np.log(p) - (1.0 - y) * np.log1p(-p))
 15
 16class LagNewton:
 17    def __init__(self, r, horizon, eta=1.0, eps=1e-8):
 18        self.r = np.asarray(r, dtype=float)
 19        self.T = len(r)
 20        n = np.maximum(np.arange(horizon, horizon - self.T, -1), 1.0)
 21        self.prior = 1.0 / (n * self.r * self.r + eps)
 22        self.H = self.prior.copy()
 23        self.theta = np.zeros(self.T)
 24        self.eta = eta
 25
 26    def step(self, x, y):
 27        p = float(sigmoid(np.dot(self.theta, x)))
 28        g = (p - y) * x
 29        h = p * (1.0 - p) * x * x + self.prior
 30        self.H += h
 31        self.theta = np.clip(self.theta - self.eta * g / self.H, -self.r, self.r)
 32        return p
 33
 34class AdaGrad:
 35    def __init__(self, r, lr=0.25, eps=1e-8):
 36        self.r = np.asarray(r, dtype=float)
 37        self.theta = np.zeros(len(r))
 38        self.G = np.zeros(len(r))
 39        self.lr, self.eps = lr, eps
 40
 41    def step(self, x, y):
 42        p = float(sigmoid(np.dot(self.theta, x)))
 43        g = (p - y) * x
 44        self.G += g * g
 45        self.theta = np.clip(self.theta - self.lr * g / (np.sqrt(self.G) + self.eps), -self.r, self.r)
 46        return p
 47
 48class AdamW:
 49    def __init__(self, r, lr=0.03, wd=1e-4, eps=1e-8):
 50        self.r = np.asarray(r, dtype=float)
 51        self.theta = np.zeros(len(r))
 52        self.m = np.zeros(len(r)); self.v = np.zeros(len(r)); self.t = 0
 53        self.lr, self.wd, self.eps = lr, wd, eps
 54
 55    def step(self, x, y):
 56        p = float(sigmoid(np.dot(self.theta, x)))
 57        g = (p - y) * x
 58        self.t += 1
 59        self.m = 0.9 * self.m + 0.1 * g
 60        self.v = 0.999 * self.v + 0.001 * g * g
 61        mh = self.m / (1.0 - 0.9 ** self.t)
 62        vh = self.v / (1.0 - 0.999 ** self.t)
 63        self.theta *= (1.0 - self.lr * self.wd)
 64        self.theta = np.clip(self.theta - self.lr * mh / (np.sqrt(vh) + self.eps), -self.r, self.r)
 65        return p
 66
 67def make_stream(seed, T, kind):
 68    rng = np.random.default_rng(seed)
 69    u = rng.choice([-1.0, 1.0], size=T + 40)
 70    j = np.arange(1, 21, dtype=float)
 71    if kind == 'exponential':
 72        true = 0.75 * np.exp(-j / 5.0)
 73    else:
 74        true = 0.95 / (j ** 0.8)
 75    true *= min(1.0, 2.0 / true.sum())
 76    y = np.zeros(T, dtype=int)
 77    for t in range(T):
 78        x = u[t + 20 - j.astype(int)]
 79        y[t] = int(rng.random() < sigmoid(np.dot(true, x)))
 80    return u, y, true
 81
 82def run_one(seed, kind, T=3000, d=20):
 83    u, y, true = make_stream(seed, T, kind)
 84    r = np.maximum(0.08, 1.35 * np.exp(-np.arange(1, d + 1) / 7.0))
 85    # Ensure the envelope contains the generating filter, while retaining lag weighting.
 86    r = np.maximum(r, np.abs(true) + 0.03)
 87    methods = {
 88        'lag_newton': LagNewton(r, T, eta=1.0),
 89        'adagrad': AdaGrad(r, lr=0.25),
 90        'adamw': AdamW(r, lr=0.03),
 91    }
 92    losses = {k: [] for k in methods}
 93    for t in range(T):
 94        x = u[t + 20 - np.arange(1, d + 1)]
 95        for name, opt in methods.items():
 96            p = opt.step(x, y[t])
 97            losses[name].append(binary_logloss(p, y[t]))
 98    out = {}
 99    for name, opt in methods.items():
100        a = np.asarray(losses[name])
101        out[name] = {
102            'cum_loss': float(a.sum()),
103            'mean_last_1000': float(a[-1000:].mean()),
104            'mean_first_500': float(a[:500].mean()),
105            'theta_mse': float(np.mean((opt.theta - true) ** 2)),
106            'max_abs_theta': float(np.max(np.abs(opt.theta))),
107        }
108    out['true_l1'] = float(np.sum(np.abs(true)))
109    out['newton_old_lag_abs_mean'] = float(np.mean(np.abs(methods['lag_newton'].theta[10:])))
110    return out
111
112def math_checks():
113    T, d = 100, 10
114    r = np.linspace(0.1, 0.5, d)
115    n = np.arange(T, T - d, -1.0)
116    gamma = np.sum(np.log1p(n * r * r))
117    # The claimed prior scale decreases as available rounds increase; verify exact formula.
118    prior = 1.0 / (n * r * r + 1e-8)
119    # For x^2 <= 1 and logistic curvature <= 1/4, each increment is bounded as stated.
120    p = np.linspace(0.01, 0.99, 50)
121    hmax = np.max(p * (1 - p) + prior[0])
122    H0 = prior.copy()
123    H = H0 + np.ones(d) * 0.25
124    return {
125        'gamma': float(gamma),
126        'gamma_nonnegative': bool(gamma >= 0),
127        'prior_positive': bool(np.all(prior > 0)),
128        'prior_matches_formula': bool(np.allclose(prior, 1.0 / (n * r * r + 1e-8))),
129        'curvature_upper_bound_check': bool(hmax <= 0.25 + prior[0] + 1e-12),
130        'H_monotone_check': bool(np.all(H >= H0)),
131        'recent_prior': float(prior[0]),
132        'old_prior': float(prior[-1]),
133        'recent_has_more_rounds': bool(n[0] > n[-1]),
134    }
135
136def main():
137    random.seed(7); np.random.seed(7)
138    result = {'math_checks': math_checks(), 'runs': {}}
139    for kind in ['exponential', 'polynomial']:
140        result['runs'][kind] = [run_one(100 + i, kind) for i in range(3)]
141    # Aggregate means for direct comparison.
142    result['means'] = {}
143    for kind, rows in result['runs'].items():
144        result['means'][kind] = {}
145        for method in ['lag_newton', 'adagrad', 'adamw']:
146            result['means'][kind][method] = {k: float(np.mean([x[method][k] for x in rows])) for k in rows[0][method]}
147    Path('results.json').write_text(json.dumps(result, indent=2))
148    print(json.dumps(result['math_checks'], indent=2))
149    print(json.dumps(result['means'], indent=2))
150
151if __name__ == '__main__':
152    main()