1import json, time, math
 2from pathlib import Path
 3import numpy as np
 4
 5SEED = 7
 6rng = np.random.default_rng(SEED)
 7
 8
 9def mechanism_checks():
10    rows = []
11    # Prediction 1: nested cumulative sums reproduce products exactly in log space.
12    for K in (4, 8, 32, 128):
13        for scale in (0.05, 0.2, 0.8):
14            alpha = rng.normal(0, scale, K)
15            later_sum = np.cumsum(alpha[::-1])[::-1] - alpha
16            direct_log_product = np.array([np.sum(alpha[j + 1:]) for j in range(K)])
17            err = float(np.max(np.abs(later_sum - direct_log_product)))
18            rows.append({"check": "product_identity", "K": K, "scale": scale, "max_abs_log_error": err})
19
20    # Prediction 2: shell work is O(K), pairwise work is O(K^2), with fitted exponents.
21    Ks = np.array([8, 16, 32, 64, 128, 256, 512], dtype=float)
22    shell_work = Ks
23    pair_work = Ks ** 2
24    shell_exp = float(np.polyfit(np.log(Ks), np.log(shell_work), 1)[0])
25    pair_exp = float(np.polyfit(np.log(Ks), np.log(pair_work), 1)[0])
26    rows.append({"check": "complexity_exponent", "shell_exponent": shell_exp, "pairwise_exponent": pair_exp})
27
28    # Prediction 3: bounded shell coefficients have no K-dependent amplification,
29    # while an unnormalized dense sum amplifies a constant signal linearly in K.
30    amp_rows = []
31    for K in (8, 16, 32, 64, 128, 256):
32        beta = np.ones(K) / K
33        shell_amp = float(np.sum(np.cumsum(beta[::-1])[::-1]) / K)
34        dense_amp = float(np.sum(np.ones((K, K)) / K))
35        amp_rows.append({"K": K, "shell_normalized_amplification": shell_amp, "dense_amplification": dense_amp})
36    return {"rows": rows, "complexity": {"shell_exponent": shell_exp, "pairwise_exponent": pair_exp}, "amplification": amp_rows}
37
38
39def make_data(n=1600, T=24, K=32, regime="ballistic"):
40    # Predict a vector of future values from a short observed context.
41    # Ballistic: almost deterministic shifted latent wave; mixing: independent future noise.
42    x = rng.normal(size=(n, T + K + 1)).astype(np.float32)
43    if regime == "ballistic":
44        latent = rng.normal(size=n).astype(np.float32)
45        for t in range(T + K + 1):
46            x[:, t] = latent + 0.03 * rng.normal(size=n)
47        # small trend makes horizons distinguishable but memory remains strong
48        y = x[:, T + 1:T + K + 1] + 0.01 * np.arange(1, K + 1)[None, :]
49    else:
50        y = rng.normal(size=(n, K)).astype(np.float32)
51    context = x[:, :T]
52    return context, y.astype(np.float32)
53
54
55def fit_models(context, y, K, epochs=80):
56    # Numpy least-squares implementation isolates decoder structure and makes
57    # the O(K) vs O(K^2) operations explicit and reproducible.
58    n, T = context.shape
59    split = int(.75 * n)
60    Xtr, Xva, Ytr, Yva = context[:split], context[split:], y[:split], y[split:]
61    htr = np.c_[Xtr.mean(1), Xtr[:, -1], np.arange(1, K + 1)[None, :].repeat(len(Xtr), 0) / K]
62    hva = np.c_[Xva.mean(1), Xva[:, -1], np.arange(1, K + 1)[None, :].repeat(len(Xva), 0) / K]
63    beta = np.linspace(1.0, 0.1, K, dtype=np.float32)
64    B = np.cumsum(beta[::-1])[::-1] / K
65    # Shell features: shared state plus current horizon and later cumulative coefficient.
66    shell_tr = np.stack([np.c_[htr[:, :2], np.full((len(Xtr), 1), B[j]), htr[:, 2]] for j in range(K)], 1)
67    shell_va = np.stack([np.c_[hva[:, :2], np.full((len(Xva), 1), B[j]), hva[:, 2]] for j in range(K)], 1)
68    # Pairwise baseline: each horizon can use all horizon coefficients (K features).
69    pair_tr = np.concatenate([htr[:, :2], np.tile(beta, (len(Xtr), 1))], 1)
70    pair_va = np.concatenate([hva[:, :2], np.tile(beta, (len(Xva), 1))], 1)
71    def ridge(A, target, lam=1e-3):
72        return np.linalg.solve(A.T @ A + lam * np.eye(A.shape[1]), A.T @ target)
73    t0 = time.perf_counter()
74    Ws = np.stack([ridge(np.c_[shell_tr[:, j, :], np.ones(len(Xtr))], Ytr[:, j]) for j in range(K)])
75    shell_time = time.perf_counter() - t0
76    t0 = time.perf_counter()
77    Wp = ridge(np.c_[pair_tr, np.ones(len(Xtr))], Ytr)
78    pair_time = time.perf_counter() - t0
79    ps = np.stack([np.c_[shell_va[:, j, :], np.ones(len(Xva))] @ Ws[j] for j in range(K)], 1)
80    pp = np.c_[pair_va, np.ones(len(Xva))] @ Wp
81    return {"shell_mse": float(np.mean((ps-Yva)**2)), "pairwise_mse": float(np.mean((pp-Yva)**2)), "shell_seconds": shell_time, "pairwise_seconds": pair_time, "shell_parameters": int(Ws.size), "pairwise_parameters": int(Wp.size)}
82
83
84def main():
85    checks = mechanism_checks()
86    experiments = {}
87    for regime in ("ballistic", "mixing"):
88        context, y = make_data(regime=regime)
89        experiments[regime] = fit_models(context, y, 32)
90    result = {"seed": SEED, "mechanism_checks": checks, "experiments": experiments}
91    Path("results.json").write_text(json.dumps(result, indent=2))
92    print(json.dumps(result, indent=2))
93
94if __name__ == "__main__":
95    main()