1import json
  2import time
  3import numpy as np
  4
  5
  6def sigmoid(x):
  7    x = np.clip(x, -30.0, 30.0)
  8    return 1.0 / (1.0 + np.exp(-x))
  9
 10
 11def spectrum(T, r):
 12    j = np.arange(1, len(r) + 1)
 13    n = np.maximum(T - j + 1, 0)
 14    return float(np.sum(np.log1p(n * np.asarray(r) ** 2)))
 15
 16
 17def selection_objective(theta, rhat, t, lam):
 18    """Return C(h), h=0..H, using n_{t,j}=max(t-j+1,0)."""
 19    H = len(theta)
 20    j = np.arange(1, H + 1)
 21    n = np.maximum(t - j + 1, 0.0)
 22    price = lam * np.log1p(n * np.asarray(rhat) ** 2)
 23    tail = n * np.asarray(theta) ** 2
 24    # C(h)=sum_{j<=h} price_j + sum_{j>h} tail_j
 25    return np.r_[np.sum(tail), np.cumsum(price - tail) + np.sum(tail)]
 26
 27
 28def math_check(seed=7):
 29    rng = np.random.default_rng(seed)
 30    H, T = 18, 80
 31    theta = rng.normal(0, .12, H)
 32    rhat = np.abs(theta) + .2 / np.sqrt(np.arange(1, H + 1))
 33    lam, t = .35, 63
 34    C = selection_objective(theta, rhat, t, lam)
 35    j = np.arange(1, H + 1)
 36    n = np.maximum(t - j + 1, 0.)
 37    marginal = lam * np.log1p(n * rhat ** 2) - n * theta ** 2
 38    algebra_error = float(np.max(np.abs(np.diff(C) - marginal)))
 39
 40    # Directly verify the paper's weighted tail scale against logistic KL.
 41    # For independent Rademacher inputs, compare source and truncated logits.
 42    B = .75
 43    th = rng.normal(0, .035, 7)
 44    th *= min(1., B / max(np.sum(np.abs(th)), 1e-12))
 45    hh, TT = 2, 120
 46    V = sum((TT-j+1) * th[j-1] ** 2 for j in range(hh + 1, len(th) + 1))
 47    kl_sum = 0.0
 48    reps = 4000
 49    for _ in range(reps):
 50        u = rng.choice([-1., 1.], size=len(th))
 51        a = float(np.dot(th, u))
 52        b = float(np.dot(th[:hh], u[:hh]))
 53        p, q = sigmoid(a), sigmoid(b)
 54        kl = p * np.log(p / q) + (1-p) * np.log((1-p) / (1-q))
 55        kl_sum += kl
 56    # Same expected KL is incurred on each affected round, hence V is its scale.
 57    empirical_ratio = (TT-hh) * kl_sum / reps / max(V, 1e-12)
 58    return {"objective_marginal_max_error": algebra_error,
 59            "tail_energy": float(V), "empirical_logistic_KL_over_tail_energy": float(empirical_ratio),
 60            "math_check_pass": bool(algebra_error < 1e-10 and .05 < empirical_ratio < .2)}
 61
 62
 63def make_data(T, H, seed):
 64    rng = np.random.default_rng(seed)
 65    u = rng.choice([-1., 1.], size=T + H + 1)
 66    theta = np.zeros((T, H))
 67    # Abruptly switch from short exponential memory to a longer polynomial-like tail.
 68    a = .40 * (.55 ** np.arange(H))
 69    b = .17 / (np.arange(1, H + 1) ** .72)
 70    theta[:T//2] = a
 71    theta[T//2:] = b
 72    y = np.zeros(T, dtype=np.int8)
 73    for t in range(T):
 74        x = u[H + t - np.arange(1, H + 1)]
 75        y[t] = rng.random() < sigmoid(np.dot(theta[t], x))
 76    return u, y, theta
 77
 78
 79def run_model(u, y, true_theta, mode, H, lam=.012, k=.08, lr=.035, warmup=100):
 80    T = len(y)
 81    est = np.zeros(H)
 82    count = np.zeros(H)
 83    active_h = H if mode == "fixed_long" else (4 if mode == "fixed_short" else 1)
 84    losses, hs, ops = [], [], []
 85    t0 = time.perf_counter()
 86    for t in range(T):
 87        x = u[H + t - np.arange(1, H + 1)]
 88        h_used = H if mode == "fixed_long" else (4 if mode == "fixed_short" else active_h)
 89        logit = float(np.dot(est[:h_used], x[:h_used]))
 90        p = float(sigmoid(logit))
 91        losses.append(-(y[t] * np.log(p + 1e-12) + (1-y[t]) * np.log(1-p + 1e-12)))
 92        err = p - y[t]
 93        # Full estimates are maintained as a transparent online pilot for unseen lags;
 94        # only the active prefix contributes to prediction and its measured inference ops.
 95        est -= lr * err * x
 96        count += 1
 97        rhat = np.abs(est) + k / np.sqrt(count)
 98        if mode == "adaptive" and t >= warmup and t % 5 == 0:
 99            C = selection_objective(est, rhat, t + 1, lam)
100            best = np.min(C)
101            candidates = np.flatnonzero(C <= best * 1.02 + 1e-9)
102            active_h = int(candidates[0]) if len(candidates) else int(np.argmin(C))
103            active_h = max(1, min(H, active_h))
104        hs.append(h_used)
105        ops.append(h_used)
106    elapsed = time.perf_counter() - t0
107    losses = np.asarray(losses)
108    return {"cum_logloss": float(np.sum(losses)),
109            "first_half_logloss": float(np.sum(losses[:T//2])),
110            "second_half_logloss": float(np.sum(losses[T//2:])),
111            "mean_active_h": float(np.mean(hs)),
112            "second_half_active_h": float(np.mean(hs[T//2:])),
113            "inference_ops_relative": float(np.sum(ops) / (T * H)),
114            "wall_seconds": float(elapsed)}
115
116
117def main():
118    checks = math_check()
119    T, H = 7000, 24
120    u, y, truth = make_data(T, H, seed=19)
121    results = {}
122    for mode in ["fixed_long", "fixed_short", "adaptive"]:
123        results[mode] = run_model(u, y, truth, mode, H)
124    # Report the idea against the conventional fixed maximum context and a short fixed window.
125    out = {"math": checks, "experiment": results,
126           "settings": {"T": T, "H_max": H, "lambda": .012, "confidence_k": .08,
127                        "seed": 19, "data": "short exponential -> long polynomial abrupt switch"}}
128    with open("results.json", "w") as f:
129        json.dump(out, f, indent=2)
130    print(json.dumps(out, indent=2))
131
132
133if __name__ == "__main__":
134    main()