1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1595
  6np.random.seed(SEED)
  7random.seed(SEED)
  8
  9
 10def fit_log_slope(values, start=5, stop=None):
 11    v = np.maximum(np.asarray(values, float), 1e-30)
 12    stop = len(v) if stop is None else stop
 13    x = np.arange(start, stop)
 14    return float(np.polyfit(x, np.log(v[start:stop]), 1)[0])
 15
 16
 17def trajectory(aq, af, c=0.7, n=80, z0=1.0, y0=0.0):
 18    z = np.zeros(n); y = np.zeros(n)
 19    z[0], y[0] = z0, y0
 20    for t in range(n - 1):
 21        z[t + 1] = aq * z[t]
 22        y[t + 1] = af * y[t] + c * z[t]
 23    return z, y
 24
 25
 26def max_rate_sweep():
 27    rows = []
 28    pairs = [(0.98,.30),(.70,.90),(.55,.55),(.90,.40),(.40,.90),(.80,.60)]
 29    for aq, af in pairs:
 30        z, y = trajectory(aq, af)
 31        norm = np.sqrt(z*z + y*y)
 32        observed = math.exp(fit_log_slope(norm, 15, 55))
 33        predicted = max(aq, af)
 34        rows.append(dict(aq=aq, af=af, predicted=predicted, observed=observed,
 35                         relative_error=abs(observed-predicted)/predicted))
 36    return rows
 37
 38
 39def slow_branch_invariance():
 40    rows = []
 41    aq = .92
 42    for af in [.10,.30,.60,.80]:
 43        z, y = trajectory(aq, af)
 44        observed = math.exp(fit_log_slope(np.sqrt(z*z+y*y), 15, 55))
 45        rows.append(dict(aq=aq, af=af, predicted=aq, observed=observed,
 46                         relative_error=abs(observed-aq)/aq))
 47    return rows
 48
 49
 50def stability_sweep():
 51    # For this triangular linear transition, eigenvalues are aq and af.
 52    # Thus the asymptotic norm slope is log(max(|aq|,|af|)); its sign changes at 1.
 53    rows = []
 54    for aq, af in [(0.95,.5),(.99,.99),(1.0,.5),(1.01,.5),(.5,1.01),(.5,.98)]:
 55        z, y = trajectory(aq, af, n=300)
 56        norm = np.sqrt(z*z+y*y)
 57        slope = fit_log_slope(norm, 100, 250)
 58        predicted_rate = max(abs(aq), abs(af))
 59        predicted_stable = predicted_rate < 1
 60        observed_rate = math.exp(slope)
 61        observed_stable = slope < -1e-5
 62        rows.append(dict(aq=aq, af=af, predicted_rate=predicted_rate,
 63                         observed_rate=observed_rate, predicted_stable=predicted_stable,
 64                         observed_stable=observed_stable,
 65                         slope=slope))
 66    return rows
 67
 68
 69def learning_comparison():
 70    # Fixed linear system identification: split model has independently constrained rates;
 71    # monolithic matched-size model has unconstrained 2x2 transition.
 72    rng = np.random.default_rng(SEED)
 73    T, N = 25, 1600
 74    X = rng.normal(size=(N,T,1))
 75    Z = np.zeros((N,T)); Y = np.zeros((N,T))
 76    for i in range(N):
 77        for t in range(T-1):
 78            Z[i,t+1] = .96*Z[i,t] + .08*X[i,t,0] + .03*rng.normal()
 79            Y[i,t+1] = .42*Y[i,t] + .65*Z[i,t] + .08*X[i,t,0] + .03*rng.normal()
 80    H = np.stack([Z,Y], axis=-1)
 81    train = np.arange(1200); test = np.arange(1200,N)
 82    # Least-squares one-step fit with states supplied, then multi-step free rollout.
 83    def mse(a,b): return float(np.mean((a-b)**2))
 84    inp = np.concatenate([H[train,:-1], X[train,:-1]], axis=-1).reshape(-1,3)
 85    target = H[train,1:].reshape(-1,2)
 86    W = np.linalg.lstsq(inp, target, rcond=None)[0]
 87    # split: z uses z,x; y uses y,z,x, preserving triangular structure
 88    wz = np.linalg.lstsq(inp[:,[0,2]], target[:,0], rcond=None)[0]
 89    wy = np.linalg.lstsq(inp[:,[1,0,2]], target[:,1], rcond=None)[0]
 90    one_m = mse(inp @ W, target)
 91    one_s = mse(np.column_stack([inp[:,[0,2]] @ wz, inp[:,[1,0,2]] @ wy]), target)
 92    def rollout(kind, idx):
 93        out=[]
 94        for i in idx:
 95            h = H[i,0].copy(); pred=[]
 96            for t in range(T-1):
 97                q = X[i,t,0]
 98                if kind=='mono': h = np.array([h[0],h[1],q]) @ W
 99                else: h = np.array([np.array([h[0],q]) @ wz, np.array([h[1],h[0],q]) @ wy])
100                pred.append(h.copy())
101            out.append(pred)
102        return np.asarray(out)
103    pm, ps = rollout('mono', test), rollout('split', test)
104    truth = H[test,1:]
105    return dict(monolithic_one_step=one_m, split_one_step=one_s,
106                monolithic_rollout=mse(pm,truth), split_rollout=mse(ps,truth),
107                n_test=len(test))
108
109
110def main():
111    result = dict(seed=SEED,
112                  prediction_1_max_rate=max_rate_sweep(),
113                  prediction_2_slow_branch_invariance=slow_branch_invariance(),
114                  prediction_3_stability_boundary=stability_sweep(),
115                  learning_comparison=learning_comparison())
116    Path('results.json').write_text(json.dumps(result, indent=2))
117    print(json.dumps(result, indent=2))
118
119if __name__ == '__main__':
120    main()