1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1222
  6OUT = Path('results.json')
  7
  8
  9def trajectory(c0, lambdas, eta, gamma=0.0, n=300):
 10    c0 = np.asarray(c0, dtype=float)
 11    factors = 1.0 - eta * np.asarray(lambdas, dtype=float)
 12    factors[0] -= gamma
 13    t = np.arange(n + 1)[:, None]
 14    return c0[None, :] * factors[None, :] ** t, factors
 15
 16
 17def A(c):
 18    return 0.5 * float(np.sum(np.asarray(c) ** 2))
 19
 20
 21def crossing_time_discrete(ca, cb, factors):
 22    # Exact continuous-time interpolation of the discrete recurrence.
 23    r1, r2 = abs(factors[0]), abs(factors[1])
 24    numerator = cb[0] ** 2 - ca[0] ** 2
 25    denominator = ca[1] ** 2 - cb[1] ** 2
 26    if numerator <= 0 or denominator <= 0 or r1 <= 0 or r2 <= 0 or r1 == r2:
 27        return np.nan
 28    return np.log(denominator / numerator) / (2 * np.log(r1 / r2))
 29
 30
 31def observed_cross(ca, cb, factors, n=100000):
 32    t = np.arange(n + 1)
 33    aa = .5 * (ca[0] ** 2 * factors[0] ** (2*t) + ca[1] ** 2 * factors[1] ** (2*t))
 34    ab = .5 * (cb[0] ** 2 * factors[0] ** (2*t) + cb[1] ** 2 * factors[1] ** (2*t))
 35    d = aa - ab
 36    ix = np.where((d[:-1] * d[1:] <= 0) & (t[:-1] > 0))[0]
 37    return float(ix[0]) if len(ix) else np.nan
 38
 39
 40def log_slope(values, start):
 41    t = np.arange(len(values))
 42    mask = (t >= start) & (values > 1e-250)
 43    return float(np.polyfit(t[mask], np.log(values[mask]), 1)[0])
 44
 45
 46def main():
 47    # Prediction 1: t* scales as 1/(lambda2-lambda1), and the exact discrete formula
 48    # agrees with the first integer step at which the ordering reverses.
 49    ca, cb, eta, lam1 = np.array([.20, 1.00]), np.array([.60, .20]), 1e-3, 1.
 50    crossing = []
 51    for gap in [.5, 1., 2., 4.]:
 52        lambdas = np.array([lam1, lam1 + gap])
 53        factors = 1 - eta * lambdas
 54        continuous = np.log((ca[1]**2-cb[1]**2)/(cb[0]**2-ca[0]**2)) / (2*gap)
 55        pred = crossing_time_discrete(ca, cb, factors)
 56        obs = observed_cross(ca, cb, factors)
 57        crossing.append({'gap': gap, 'continuous_time': continuous,
 58                         'discrete_formula_step': pred, 'observed_step': obs,
 59                         'relative_error': abs(pred-obs)/obs})
 60
 61    # Prediction 2: for c_{n+1}=(I-eta L-gamma e1e1^T)c_n, the slow-mode
 62    # stability edge is |1-eta lambda1-gamma|=1, i.e. gamma=2-eta lambda1.
 63    lambdas, eta2 = np.array([1., 3.]), .2
 64    gamma_star = 2 - eta2*lambdas[0]
 65    probes = np.linspace(gamma_star-.3, gamma_star+.3, 13)
 66    stability = []
 67    for gamma in probes:
 68        factors = np.array([1-eta2*lambdas[0]-gamma, 1-eta2*lambdas[1]])
 69        stable_exact = bool(np.max(np.abs(factors)) < 1)
 70        c, _ = trajectory([1., .2], lambdas, eta2, gamma=gamma, n=100)
 71        growth = float(np.max(np.linalg.norm(c, axis=1))/np.linalg.norm(c[0]))
 72        stability.append({'gamma': float(gamma), 'max_factor_abs': float(np.max(np.abs(factors))),
 73                          'predicted_stable': stable_exact, 'observed_bounded_100_steps': growth < 1.000001,
 74                          'norm_growth': growth})
 75
 76    # Prediction 3: late A2 slope is 2 log|r_j| for the slowest surviving mode.
 77    eta3, lambdas3 = .05, np.array([1., 2.5])
 78    slope = []
 79    for gamma in [0., .2, .6, .95]:
 80        c, factors = trajectory([.6, 1.], lambdas3, eta3, gamma=gamma, n=300)
 81        a2 = .5*np.sum(c*c, axis=1)
 82        dominant = int(np.argmax(np.abs(factors)))
 83        slope.append({'gamma': gamma, 'factors': factors.tolist(),
 84                      'predicted_slope': 2*np.log(abs(factors[dominant])),
 85                      'observed_slope': log_slope(a2, 100), 'dominant_mode': dominant+1})
 86    c, factors = trajectory([0., 1.], lambdas3, eta3, n=300)
 87    a2 = .5*np.sum(c*c, axis=1)
 88    ideal_filter = {'expected_fast_slope': 2*np.log(abs(factors[1])),
 89                    'observed_fast_slope': log_slope(a2, 50)}
 90
 91    # Training analogue: diagonal quadratic objective. The filtered run deliberately
 92    # begins with larger A2, but suppresses its slow coefficient and overtakes baseline.
 93    train = []
 94    for name, x0, gamma in [('baseline', [.60, .20], 0.), ('mode_filtered', [.20, 1.00], .75)]:
 95        x, factors = trajectory(x0, np.array([1., 3.]), .1, gamma=gamma, n=80)
 96        a2 = .5*np.sum(x*x, axis=1)
 97        train.append({'name': name, 'initial_A2': float(a2[0]), 'final_A2': float(a2[-1]),
 98                      'factors': factors.tolist(), 'a2_first_10': a2[:10].tolist()})
 99    # Find the crossing over full saved trajectories.
100    trajectories = []
101    for x0, gamma in [([.60,.20],0.),([.20,1.],.75)]:
102        x,_ = trajectory(x0, np.array([1.,3.]), .1, gamma=gamma, n=80)
103        trajectories.append(.5*np.sum(x*x,axis=1))
104    cross = next((i for i in range(1,81) if trajectories[1][i] < trajectories[0][i]), None)
105    result = {'seed': SEED, 'predictions': {'crossing': crossing,
106        'stability': {'eta': eta2, 'lambda1': 1., 'predicted_upper_gamma': gamma_star, 'sweep': stability},
107        'late_slope': slope, 'ideal_zero_slow_mode': ideal_filter},
108        'training_analogue': {'runs': train, 'first_filtered_below_baseline_step': cross}}
109    OUT.write_text(json.dumps(result, indent=2))
110    summary = {'crossing': [{'gap': x['gap'], 'pred': round(x['discrete_formula_step'],3), 'obs': x['observed_step'], 'relerr': round(x['relative_error'],6)} for x in crossing],
111               'stability_edge_predicted': gamma_star,
112               'stability_edge_observed_bracket': [x['gamma'] for x in stability if x['predicted_stable']][-1:],
113               'late_slope': slope, 'ideal_filter': ideal_filter,
114               'training': {'initial_A2': [x['initial_A2'] for x in train], 'final_A2': [x['final_A2'] for x in train], 'cross_step': cross}}
115    print(json.dumps(summary, indent=2))
116
117if __name__ == '__main__':
118    main()