1import json
  2from pathlib import Path
  3import numpy as np
  4from scipy.optimize import least_squares
  5from scipy.signal import periodogram
  6
  7SEED = 1120
  8rng = np.random.default_rng(SEED)
  9
 10
 11def simulate(r, omega, sigma, B=128, T=3000, burn=500, seed=0):
 12    """Independent noisy trajectories of h[t+1]=r R(omega)h[t]+sigma eps."""
 13    g = np.random.default_rng(seed)
 14    R = np.array([[np.cos(omega), -np.sin(omega)],
 15                  [np.sin(omega),  np.cos(omega)]])
 16    h = g.normal(size=(B, 2))
 17    out = np.empty((B, T))
 18    for t in range(T + burn):
 19        h = h @ (r * R).T + sigma * g.normal(size=(B, 2))
 20        if t >= burn:
 21            out[:, t-burn] = h[:, 0]
 22    return out
 23
 24
 25def ensemble_corr(y, maxlag=80):
 26    # Average stationary autocovariance over trajectories, as in the proposed estimator.
 27    y = y - y.mean(axis=1, keepdims=True)
 28    c = np.array([(y[:, :y.shape[1]-k] * y[:, k:]).mean() for k in range(maxlag+1)])
 29    return c / max(c[0], 1e-12)
 30
 31
 32def fit_resonance(c, max_r=0.9999):
 33    k = np.arange(len(c), dtype=float)
 34    # A damped sinusoid has two real coefficients and the resonance radius/frequency.
 35    def residual(p):
 36        r, w, a, b = p
 37        return r**k * (a*np.cos(w*k) + b*np.sin(w*k)) - c
 38    best = None
 39    for w0 in np.linspace(.12, 1.25, 8):
 40        x = least_squares(residual, [0.9, w0, 1.0, 0.0],
 41                          bounds=([.01, .01, -3, -3], [max_r, 1.5, 3, 3]),
 42                          max_nfev=1000)
 43        if best is None or np.sum(x.fun*x.fun) < best[0]:
 44            best = (np.sum(x.fun*x.fun), x.x)
 45    return best[1]
 46
 47
 48def spectrum_peak(y, omega_expected):
 49    # Average periodograms; ignore DC and select the nearest local spectral maximum.
 50    f, p = periodogram(y, axis=1, detrend='constant')
 51    p = p.mean(axis=0)
 52    target = omega_expected / (2*np.pi)
 53    band = (f > max(1/len(y[0]), target-.12)) & (f < target+.12)
 54    ix = np.where(band)[0][np.argmax(p[band])]
 55    return 2*np.pi*f[ix], p[ix], f, p
 56
 57
 58def jacobian_lyapunov(r, omega, T=10000):
 59    # Exact tangent update: J=rR, hence every normalized tangent grows by r.
 60    # Written as a numerical normalized iteration to verify the estimator itself.
 61    R = np.array([[np.cos(omega), -np.sin(omega)],
 62                  [np.sin(omega),  np.cos(omega)]])
 63    J = r * R
 64    v = np.array([1., .37]); v /= np.linalg.norm(v)
 65    logs=[]
 66    for _ in range(T):
 67        q = J @ v
 68        n = np.linalg.norm(q); logs.append(np.log(n)); v=q/n
 69    return float(np.mean(logs))
 70
 71
 72def run():
 73    outdir = Path('results'); outdir.mkdir(exist_ok=True)
 74    omega = .42
 75    # Prediction 1: tau=-1/log(r), estimated correlation envelope follows it.
 76    radius_rows=[]
 77    for i, r in enumerate([.70, .80, .90, .95, .98]):
 78        y = simulate(r, omega, .12, B=128, T=3500, seed=100+i)
 79        c = ensemble_corr(y, 100)
 80        fitr, fitw, a, b = fit_resonance(c)
 81        tau_pred = -1/np.log(r)
 82        tau_fit = -1/np.log(fitr)
 83        lam = jacobian_lyapunov(r, omega)
 84        radius_rows.append(dict(r=r, fitted_r=float(fitr), tau_pred=float(tau_pred),
 85                                tau_fit=float(tau_fit), lyapunov=float(lam),
 86                                rel_tau_error=float(abs(tau_fit-tau_pred)/tau_pred)))
 87
 88    # Prediction 2: peak frequency is omega/(2pi), independently of radius/noise.
 89    freq_rows=[]
 90    for i, w in enumerate([.25, .42, .70, 1.00]):
 91        y = simulate(.90, w, .12, B=128, T=3500, seed=200+i)
 92        c = ensemble_corr(y, 80)
 93        fitr, fitw, a, b = fit_resonance(c)
 94        peak, _, _, _ = spectrum_peak(y, w)
 95        freq_rows.append(dict(omega=w, fitted_omega=float(fitw), psd_peak=float(peak),
 96                              predicted=float(w), peak_abs_error=float(abs(peak-w)),
 97                              relative_error=float(abs(peak-w)/w)))
 98
 99    # Prediction 3: additive noise changes variance but not resonance or Lyapunov exponent.
100    noise_rows=[]
101    for i, sig in enumerate([.03, .12, .40, .80]):
102        y = simulate(.93, omega, sig, B=128, T=3500, seed=300+i)
103        c=ensemble_corr(y, 80)
104        fitr, fitw, a, b=fit_resonance(c)
105        noise_rows.append(dict(sigma=sig, fitted_r=float(fitr), fitted_omega=float(fitw),
106                               lyapunov=float(jacobian_lyapunov(.93,omega)),
107                               variance=float(y.var())))
108
109    # Controller: resonance penalty clips an estimated radius above r_max.
110    # Compare the uncontrolled high-Q oscillator with the controlled one.
111    r_raw=.975; r_max=.90
112    yc=simulate(r_raw, omega, .12, B=128, T=4000, seed=500)
113    # In this linear MVP, applying the penalty is equivalent to the prescribed gain control.
114    yctrl=simulate(min(r_raw,r_max), omega, .12, B=128, T=4000, seed=501)
115    def summarize(y, r):
116        c=ensemble_corr(y,100); fr,fw,_,_=fit_resonance(c); peak,powr,_,_=spectrum_peak(y,omega)
117        return dict(fitted_r=float(fr), fitted_omega=float(fw), decay_time=float(-1/np.log(fr)),
118                    peak_frequency=float(peak), peak_power=float(powr), long_lag_corr=float(abs(c[80])),
119                    lyapunov=float(jacobian_lyapunov(r,omega)))
120    control = {'uncontrolled':summarize(yc,r_raw), 'resonance_control':summarize(yctrl,r_max),
121               'target_radius':r_max}
122
123    decay_pass = all(z['rel_tau_error'] < .05 for z in radius_rows)
124    freq_pass = all(z['relative_error'] < .10 for z in freq_rows)
125    noise_r_span = max(z['fitted_r'] for z in noise_rows) - min(z['fitted_r'] for z in noise_rows)
126    noise_w_span = max(z['fitted_omega'] for z in noise_rows) - min(z['fitted_omega'] for z in noise_rows)
127    noise_pass = noise_r_span < .01 and noise_w_span < .01
128    control_gain = control['uncontrolled']['peak_power'] / control['resonance_control']['peak_power']
129    report={'seed':SEED, 'model':'h[t+1]=r R(omega) h[t]+sigma Normal(0,I), y=h[0]',
130            'predictions':{
131              'decay_time':'tau=-1/log(r); fitted resonance radius predicts correlation decay',
132              'frequency':'PSD peak and fitted omega equal injected omega',
133              'independence':'additive noise changes variance, not J-based Lambda or resonance'},
134            'radius_sweep':radius_rows, 'frequency_sweep':freq_rows,
135            'noise_sweep':noise_rows, 'control':control,
136            'quantitative_checks': {
137              'decay_time_max_relative_error': max(z['rel_tau_error'] for z in radius_rows),
138              'decay_time_pass_under_5pct': decay_pass,
139              'frequency_max_relative_error': max(z['relative_error'] for z in freq_rows),
140              'frequency_pass_under_10pct': freq_pass,
141              'noise_fitted_radius_span': noise_r_span,
142              'noise_fitted_frequency_span': noise_w_span,
143              'noise_invariance_pass_under_0.01': noise_pass,
144              'baseline_peak_power': control['uncontrolled']['peak_power'],
145              'idea_peak_power': control['resonance_control']['peak_power'],
146              'peak_power_reduction_factor': control_gain,
147              'baseline_long_lag_abs_corr': control['uncontrolled']['long_lag_corr'],
148              'idea_long_lag_abs_corr': control['resonance_control']['long_lag_corr']
149            }}
150    Path('results/report.json').write_text(json.dumps(report, indent=2))
151    # Compact CSV-like human-readable summary for quick inspection.
152    print(json.dumps(report, indent=2))
153
154if __name__ == '__main__':
155    run()