1import json
  2import numpy as np
  3from scipy.optimize import root
  4
  5SEED = 1128
  6A, W0, H = 0.04, 1.0, 0.01
  7
  8# Modal oscillator: z' = (-A + i W0) z - q z(t-tau).
  9def char_root(q, tau, guess=None):
 10    if guess is None:
 11        guess = complex(-A, W0)
 12    def fun(v):
 13        s = v[0] + 1j*v[1]
 14        f = s + A - 1j*W0 + q*np.exp(-s*tau)
 15        return [f.real, f.imag]
 16    sol = root(fun, [guess.real, guess.imag], method='hybr')
 17    if not sol.success:
 18        raise RuntimeError(sol.message)
 19    return complex(sol.x[0], sol.x[1])
 20
 21
 22def compensated_root(q, tau, guess=None):
 23    """Root for the implemented predictor: x[t-d]+d*(x[t-d]-x[t-d-1])."""
 24    if guess is None:
 25        guess = complex(-A, W0)
 26    d = int(round(tau / H))
 27    def fun(v):
 28        s = v[0] + 1j*v[1]
 29        # predictor transfer factor for exp(s t), with d samples of delay
 30        pred = np.exp(-s*d*H) * (1.0 + d*(1.0-np.exp(-s*H)))
 31        f = s + A - 1j*W0 + q*pred
 32        return [f.real, f.imag]
 33    sol = root(fun, [guess.real, guess.imag], method='hybr')
 34    if not sol.success:
 35        raise RuntimeError(sol.message)
 36    return complex(sol.x[0], sol.x[1])
 37
 38
 39def simulate_mode(q, tau, compensated=False, T=240.0):
 40    n, d = int(T/H), int(round(tau/H))
 41    t = np.arange(n+1)*H
 42    z = np.exp((-A+1j*W0)*(t-tau-2*H)).astype(complex)
 43    for k in range(1, n+1):
 44        idx = max(0, k-d)
 45        src = z[idx]
 46        if compensated and d > 0:
 47            src = src + d*(src-z[max(0, idx-1)])
 48        z[k] = z[k-1] + H*((-A+1j*W0)*z[k-1]-q*src)
 49    return z
 50
 51
 52def dominant_freq(z, discard=.45):
 53    """Phase-slope estimator avoids FFT-bin quantization for a single analytic mode."""
 54    x = z[int(len(z)*discard):]
 55    phase = np.unwrap(np.angle(x))
 56    slope = np.polyfit(np.arange(len(x))*H, phase, 1)[0]
 57    return float(slope)
 58
 59
 60def main():
 61    # Prediction 1: delay phase has slope -omega0.
 62    taus = np.array([0., .1, .2, .3, .4, .5])
 63    phases = np.unwrap(np.angle(np.exp(-1j*W0*taus)))
 64    phase_slope = float(np.polyfit(taus, phases, 1)[0])
 65
 66    # Prediction 2: lambda=0 is delay invariant; a nonzero mode changes.
 67    z0, z0d = simulate_mode(0., .4), simulate_mode(0., .4, True)
 68    zero_err = float(np.max(abs(z0-z0d)))
 69    zbase, zdel = simulate_mode(.06, 0.), simulate_mode(.06, .4)
 70    nonzero_change = float(np.sqrt(np.mean(abs(zdel-zbase)**2)) /
 71                           np.sqrt(np.mean(abs(zbase)**2)))
 72
 73    # Prediction 3: first-order frequency shift is linear in q.
 74    tau = .30
 75    qs = np.array([.005, .01, .02, .04, .06])
 76    roots = np.array([char_root(float(q), tau) for q in qs])
 77    shifts = roots.imag-W0
 78    fit_slope = float(np.polyfit(qs, shifts, 1)[0])
 79    pred_slope = float(np.exp(A*tau)*np.sin(W0*tau))
 80
 81    # Four branches on a ring: Laplacian eigenvalues 0,2,2,4.
 82    lambdas, gamma, delay = [0.,2.,2.,4.], .03, .40
 83    rows = []
 84    for lam in lambdas:
 85        q = gamma*lam
 86        zd, zc = simulate_mode(q, delay), simulate_mode(q, delay, True)
 87        rb, rc = char_root(q, delay), compensated_root(q, delay)
 88        rows.append({'lambda':lam, 'baseline_freq':dominant_freq(zd),
 89            'compensated_freq':dominant_freq(zc),
 90            'baseline_root_freq':float(rb.imag), 'compensated_root_freq':float(rc.imag),
 91            'baseline_freq_error':abs(dominant_freq(zd)-W0),
 92            'compensated_freq_error':abs(dominant_freq(zc)-W0)})
 93    nz = [r for r in rows if r['lambda'] > 0]
 94    base_mean = float(np.mean([r['baseline_freq_error'] for r in nz]))
 95    comp_mean = float(np.mean([r['compensated_freq_error'] for r in nz]))
 96    result = {
 97      'seed':SEED, 'h':H, 'omega0':W0, 'damping':A,
 98      'prediction_1_phase': {'taus':taus.tolist(), 'observed_phase_slope':phase_slope,
 99                             'predicted_phase_slope':-W0},
100      'prediction_2_zero_mode': {'max_abs_error_lambda0':zero_err,
101                                 'relative_change_nonzero_q006_tau04':nonzero_change},
102      'prediction_3_frequency_scaling': {'tau':tau, 'q_values':qs.tolist(),
103        'measured_root_shifts':shifts.tolist(), 'fitted_shift_per_q':fit_slope,
104        'predicted_first_order_shift_per_q':pred_slope,
105        'relative_slope_error':abs(fit_slope-pred_slope)/abs(pred_slope)},
106      'mini_experiment': {'gamma':gamma, 'delay':delay, 'rows':rows,
107        'mean_nonzero_frequency_error_baseline':base_mean,
108        'mean_nonzero_frequency_error_compensated':comp_mean,
109        'relative_improvement':(base_mean-comp_mean)/base_mean},
110      'checks': {'phase_ok':bool(abs(phase_slope+W0)<1e-10),
111        'zero_mode_ok':bool(zero_err<1e-10 and nonzero_change>1e-3),
112        'scaling_ok':bool(abs(fit_slope-pred_slope)/abs(pred_slope)<.08),
113        'compensation_better':bool(comp_mean<base_mean)}
114    }
115    with open('results.json','w') as f: json.dump(result,f,indent=2)
116    print(json.dumps(result,indent=2))
117
118if __name__ == '__main__': main()