1import json, math
  2import numpy as np
  3
  4SEED = 937
  5rng = np.random.default_rng(SEED)
  6
  7
  8def burnin_prediction(r, eps, E0=1.0):
  9    return int(math.ceil(math.log(eps / E0) / math.log(r)))
 10
 11
 12def scalar_verification():
 13    rs = [0.50, 0.70, 0.85, 0.95, 0.99]
 14    epses = [1e-2, 1e-4, 1e-6]
 15    rows = []
 16    for r in rs:
 17        for eps in epses:
 18            pred = burnin_prediction(r, eps)
 19            e, k = 1.0, 0
 20            while e > eps and k < 100000:
 21                e *= r
 22                k += 1
 23            rows.append({'r': r, 'epsilon': eps, 'predicted_k': pred,
 24                         'observed_k': k, 'abs_error': abs(pred-k)})
 25    scaling = []
 26    for r in rs:
 27        x = np.log(1 / np.asarray(epses))
 28        y = np.asarray([burnin_prediction(r, e) for e in epses])
 29        slope = float(np.polyfit(x, y, 1)[0])
 30        target = 1 / abs(math.log(r))
 31        scaling.append({'r': r, 'observed_slope': slope,
 32                        'predicted_slope': target,
 33                        'relative_error': abs(slope-target)/target})
 34    boundary = []
 35    for r in [0.98, 0.999, 1.0, 1.001, 1.02]:
 36        e = 1.0
 37        for _ in range(1000):
 38            e *= r
 39        boundary.append({'r': r, 'E_1000': float(e),
 40                         'stable_observed': bool(e < 1),
 41                         'stable_predicted': bool(r < 1)})
 42    # Cost crossover directly evaluates the two asserted asymptotic forms.
 43    crossover = []
 44    for r in [0.5, 0.8, 0.95]:
 45        for eps in [1e-2, 1e-4, 1e-6]:
 46            obs = math.log(1/eps) / (1-r*r)
 47            rec = eps**(-1.0)  # d=1 toy attractor
 48            crossover.append({'r':r, 'epsilon':eps,
 49                              'observer_cost':obs, 'retrieval_cost':rec,
 50                              'observer_cheaper':obs < rec})
 51    return {'burnin_rows': rows, 'log_scaling': scaling,
 52            'stability_boundary': boundary, 'cost_scaling': crossover}
 53
 54
 55def lorenz_step(s, dt=0.01, sigma=10., rho=28., beta=8/3):
 56    def f(v):
 57        x, y, z = v
 58        return np.array([sigma*(y-x), x*(rho-z)-y, x*y-beta*z])
 59    k1 = f(s); k2 = f(s + dt*k1/2); k3 = f(s + dt*k2/2); k4 = f(s + dt*k3)
 60    return s + dt*(k1 + 2*k2 + 2*k3 + k4)/6
 61
 62
 63def generate(n=2600, dt=0.01):
 64    x = np.array([1., 1., 1.])
 65    truth = []
 66    for _ in range(n):
 67        x = lorenz_step(x, dt)
 68        truth.append(x.copy())
 69    truth = np.asarray(truth)
 70    # Scalar observation with moderate noise and occasional mismatch bursts.
 71    noise = rng.normal(0, 0.7, n)
 72    y = truth[:, 0] + noise
 73    bad = np.zeros(n, dtype=bool)
 74    bad[800:900] = True; bad[1700:1780] = True
 75    y[bad] += rng.normal(0, 7.0, bad.sum())
 76    return truth, y, bad
 77
 78
 79def run_policy(truth, y, mode, threshold=2.2, rmax=.94, eps=1e-3):
 80    # Deliberately small linear observer: x prediction plus scalar x correction.
 81    # Retrieval is a nearest scalar-observation context from a bank of latent states.
 82    z = truth[0].copy() + np.array([5., -3., 2.])
 83    bank_y, bank_z = [], []
 84    sqerr, retrievals, residuals, costs = [], 0, [], 0.0
 85    bad_retrievals = 0
 86    for t in range(len(y)):
 87        # One cheap model rollout. Its local contraction proxy is fixed but
 88        # varies by phase, making the switch's decision rule testable.
 89        r = .78 if (t % 700) < 500 else .985
 90        pred = lorenz_step(z)
 91        q = abs(y[t] - pred[0])
 92        use_retrieval = mode == 'retrieve' or (mode == 'adaptive' and (r > rmax or q > threshold))
 93        if mode == 'observer' or (mode == 'adaptive' and not use_retrieval):
 94            z = pred.copy()
 95            # Stable scalar innovation update, a simple observer correction.
 96            z[0] += .18 * (y[t] - z[0])
 97            costs += 1.0
 98        else:
 99            if bank_y:
100                j = int(np.argmin(np.abs(np.asarray(bank_y[-700:]) - y[t])))
101                z = bank_z[-700:][j].copy()
102            else:
103                z = pred.copy()
104            z[0] += .10 * (y[t] - z[0])
105            retrievals += 1; bad_retrievals += int(abs(y[t] - truth[t, 0]) > 4.0); costs += 7.0
106        bank_y.append(float(y[t])); bank_z.append(z.copy())
107        residuals.append(q); sqerr.append(float(np.mean((z-truth[t])**2)))
108    return {'mse':float(np.mean(sqerr)), 'rmse':float(np.sqrt(np.mean(sqerr))),
109            'retrievals':retrievals, 'retrieval_rate':retrievals/len(y),
110            'relative_cost':costs/len(y), 'mean_residual':float(np.mean(residuals)),
111            'bad_period_retrievals': int(bad_retrievals)}
112
113
114def main():
115    scalar = scalar_verification()
116    truth, y, bad = generate()
117    results = {}
118    for mode in ['observer', 'retrieve', 'adaptive']:
119        results[mode] = run_policy(truth, y, mode)
120    # Useful diagnostic specifically measures event sensitivity.
121    z = truth[0] + np.array([5., -3., 2.])
122    hits_good = hits_bad = 0
123    for t in range(len(y)):
124        r = .78 if (t % 700) < 500 else .985
125        q = abs(y[t] - lorenz_step(z)[0])
126        decision = r > .94 or q > 2.2
127        if bad[t]: hits_bad += int(decision)
128        else: hits_good += int(decision)
129        z = lorenz_step(z); z[0] += .18*(y[t]-z[0])
130    results['adaptive_event_diagnostics'] = {
131        'bad_steps': int(bad.sum()), 'bad_switches': hits_bad,
132        'bad_recall': hits_bad/max(1,int(bad.sum())),
133        'good_switch_rate': hits_good/max(1,int((~bad).sum()))}
134    out = {'seed': SEED, 'scalar_verification': scalar,
135           'lorenz_results': results,
136           'note': 'Toy observer/retrieval costs are normalized units; retrieval cost is 7x observer.'}
137    with open('results.json','w') as f: json.dump(out,f,indent=2)
138    print(json.dumps(out, indent=2))
139
140if __name__ == '__main__': main()