1import json
  2import math
  3import numpy as np
  4
  5
  6def detailed_balance_check(seed=0):
  7    rng = np.random.default_rng(seed)
  8    beta, nu = 1.7, 0.8
  9    vals = []
 10    for _ in range(1000):
 11        work, du = rng.normal(size=2)
 12        kf = nu * np.exp(beta * (work - du) / 2)
 13        # Reverse work is -work and reverse energy change is -du.
 14        kr = nu * np.exp(beta * (-work + du) / 2)
 15        vals.append(math.log(kf / kr) - beta * (work - du))
 16    return float(np.max(np.abs(vals)))
 17
 18
 19def exact_response(eta, lam, T, delta, P):
 20    """Expected excess response for x'=a*x+sqrt(2 eta T) z, q=dx^2/eta."""
 21    a = 1.0 - eta * lam
 22    v0 = 2.0 * eta * T / (1.0 - a*a)
 23    v = v0
 24    q0 = 4.0 * T / (1.0 + a)
 25    excess = 0.0
 26    for _ in range(P):
 27        v = a*a * v + 2.0 * eta * (T + delta)
 28        q = (1.0 - a)**2 * v / eta + 2.0 * (T + delta)
 29        excess += q - q0
 30    return excess / delta
 31
 32
 33def monte_carlo_response(eta, lam, T, delta, P, ntraj=30000, seed=1):
 34    rng = np.random.default_rng(seed)
 35    a = 1.0 - eta * lam
 36    v0 = 2.0 * eta * T / (1.0 - a*a)
 37    x = rng.normal(0.0, math.sqrt(v0), size=ntraj)
 38    # Estimate baseline q from independent stationary transitions.
 39    xb = rng.normal(0.0, math.sqrt(v0), size=ntraj)
 40    qss = np.mean(((a * xb + math.sqrt(2*eta*T)*rng.normal(size=ntraj) - xb)**2) / eta)
 41    total = 0.0
 42    for _ in range(P):
 43        xn = a*x + math.sqrt(2*eta*(T+delta))*rng.normal(size=ntraj)
 44        total += np.mean((xn-x)**2 / eta) - qss
 45        x = xn
 46    return total / delta
 47
 48
 49def run():
 50    lam, T, P = 1.0, 1.0, 20
 51    # Prediction 1: explicit Euler SGD is stable iff 0 < eta*lambda < 2.
 52    etas = np.array([*np.linspace(0.2, 1.8, 9), 1.9, 1.95, 1.99, 2.01, 2.1])
 53    rows = []
 54    for eta in etas:
 55        a = 1 - eta*lam
 56        C = exact_response(eta, lam, T, 0.05*T, P)
 57        tau = -1.0 / math.log(abs(a)) if abs(a) not in (0, 1) else (0.0 if a == 0 else float('inf'))
 58        rows.append({'eta': float(eta), 'eta_lambda': float(eta*lam), 'stable': bool(abs(a)<1), 'C': float(C), 'tau': float(tau)})
 59    # Predicted pole C ~ const/(2-eta lambda): test scaling at the stable edge.
 60    edge = [r for r in rows if 1.8 <= r['eta_lambda'] < 2.0]
 61    scaled = [r['C'] * (2-r['eta_lambda']) for r in edge]
 62    edge_cv = float(np.std(scaled) / np.mean(scaled))
 63    response_monotone_stable = bool(np.all(np.diff([r['C'] for r in rows if r['stable']]) > 0))
 64    # Prediction 2: relaxation time increases toward the stability boundary.
 65    stable_near = [r for r in rows if r['stable'] and r['eta_lambda'] > 1.2]
 66    taus = np.array([r['tau'] for r in stable_near])
 67    monotone_tau = bool(np.all(np.diff(taus) > 0))
 68    # Prediction 3: differential response is linear in pulse size for small pulses.
 69    deltas = np.array([0.01, 0.02, 0.05, 0.1]) * T
 70    Cs = np.array([exact_response(1.5, lam, T, d, P) for d in deltas])
 71    # Compare to infinitesimal numerical derivative and coefficient of variation.
 72    C0 = exact_response(1.5, lam, T, 1e-5, P)
 73    linear_rel_error = float(np.max(np.abs(Cs/C0 - 1)))
 74    mc = monte_carlo_response(1.5, lam, T, 0.05, P, ntraj=100000, seed=4)
 75    exact = exact_response(1.5, lam, T, 0.05, P)
 76
 77    # Small optimizer comparison: noisy quadratic, fixed eta vs calorimetric trigger.
 78    def train(adaptive, seed):
 79        rng = np.random.default_rng(seed); x = 4.0; eta = 1.5; losses=[]; triggers=0
 80        for step in range(300):
 81            # Probe every 30 updates, using a short temperature pulse around current state.
 82            if adaptive and step > 0 and step % 30 == 0:
 83                c = exact_response(eta, lam, T, 0.05*T, 20)
 84                if c > 35.0:  # response threshold selected before trajectory, near edge
 85                    eta *= 0.7; triggers += 1
 86            grad = lam*x + rng.normal(0, math.sqrt(2*T/eta))
 87            x -= eta*grad
 88            losses.append(0.5*lam*x*x)
 89            if abs(x) > 1e6: return float('inf'), triggers
 90        return float(np.mean(losses[-30:])), triggers
 91    base = [train(False, s)[0] for s in range(20)]
 92    idea = [train(True, s)[0] for s in range(20)]
 93    triggers = [train(True, s)[1] for s in range(20)]
 94    out = {
 95      'detailed_balance_max_log_error': detailed_balance_check(),
 96      'prediction_stability_boundary': {'predicted_eta_lambda': 2.0, 'sweep_last_stable': max(r['eta_lambda'] for r in rows if r['stable']), 'sweep_first_unstable': min(r['eta_lambda'] for r in rows if not r['stable']), 'all_stable_below_2': all(r['stable'] for r in rows if r['eta_lambda'] < 2.0), 'first_unstable_is_above_2': all(not r['stable'] for r in rows if r['eta_lambda'] > 2.0)},
 97      'prediction_response_pole': {'predicted': 'C*(2-eta*lambda) approximately constant near edge', 'edge_cv': edge_cv, 'response_monotone_on_stable_sweep': response_monotone_stable, 'edge_rows': edge},
 98      'prediction_critical_slowing': {'predicted': 'tau increases as eta*lambda approaches 2 from below', 'monotone_on_1.2_to_edge': monotone_tau, 'taus': [float(x) for x in taus]},
 99      'prediction_pulse_linearity': {'deltas_over_T': [float(x) for x in deltas], 'C_values': [float(x) for x in Cs], 'max_relative_error_vs_infinitesimal': linear_rel_error},
100      'monte_carlo_check': {'exact_C': exact, 'mc_C': mc, 'relative_error': abs(mc-exact)/abs(exact)},
101      'mini_experiment': {'baseline_final_loss_mean': float(np.mean(base)), 'idea_final_loss_mean': float(np.mean(idea)), 'baseline_divergences': int(sum(np.isinf(base))), 'idea_divergences': int(sum(np.isinf(idea))), 'idea_trigger_count_mean': float(np.mean(triggers))},
102      'rows': rows
103    }
104    with open('results.json','w') as f: json.dump(out, f, indent=2)
105    print(json.dumps(out, indent=2))
106
107if __name__ == '__main__': run()