1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1111
  6rng = np.random.default_rng(SEED)
  7
  8class HystereticRouter:
  9    def __init__(self, eta_on=0.6, eta_off=0.4, active=0):
 10        self.on, self.off, self.active = eta_on, eta_off, active
 11        self.switches = []
 12    def step(self, h, t):
 13        # A switch is permitted only after the current attractor has exited.
 14        if h[self.active] <= self.off:
 15            candidates = [i for i in range(len(h)) if i != self.active and h[i] >= self.on]
 16            if candidates:
 17                old = self.active
 18                self.active = int(candidates[np.argmax(h[candidates])])
 19                self.switches.append((t, old, self.active))
 20        return self.active
 21
 22def dwell_sweep():
 23    # Scores are deliberately generated with a known bounded derivative.
 24    # Candidate rises from off to on, then the old score falls; traversal takes Delta/L.
 25    L = 0.017
 26    rows = []
 27    for delta in [0.04, 0.08, 0.14, 0.22]:
 28        off, on = 0.40, 0.40 + delta
 29        dt = 0.01
 30        n = 30000
 31        t = np.arange(n) * dt
 32        # candidate 1 rises while active score 0 falls, each slope <= L
 33        h0 = np.clip(0.75 - L*t, 0, 1)
 34        h1 = np.clip(0.20 + L*t, 0, 1)
 35        router = HystereticRouter(on, off)
 36        for k in range(n): router.step(np.array([h0[k], h1[k]]), t[k])
 37        assert router.switches
 38        first = router.switches[0][0]
 39        # Now explicitly measure the candidate's time from off to on.
 40        expected = delta / L
 41        observed = expected # continuous ramp identity, independently verified below
 42        # finite-difference derivative and actual threshold crossing times
 43        rise_on = np.where(h1 >= on)[0][0] * dt
 44        rise_off = np.where(h1 >= off)[0][0] * dt
 45        observed = rise_on - rise_off
 46        max_deriv = np.max(np.abs(np.diff(h1) / dt))
 47        rows.append({'delta': delta, 'predicted_dwell': expected,
 48                     'observed_dwell': observed, 'bound_ratio': observed / expected,
 49                     'max_fd_derivative': float(max_deriv), 'switch_time': float(first)})
 50    return rows
 51
 52def chatter_sweep():
 53    # A sinusoidal score of amplitude A around the threshold. Without hysteresis,
 54    # crossings occur repeatedly. If Delta > 2A, the score cannot traverse both thresholds.
 55    dt, duration, A, period = 0.01, 40.0, 0.08, 2.0
 56    t = np.arange(0, duration, dt)
 57    s = 0.5 + A*np.sin(2*np.pi*t/period)
 58    rows = []
 59    for delta in [0.00, 0.04, 0.10, 0.16, 0.20]:
 60        off, on = 0.5-delta/2, 0.5+delta/2
 61        # use complementary scores, which makes each threshold crossing meaningful
 62        r = HystereticRouter(on, off)
 63        for k, v in enumerate(s): r.step(np.array([1-v, v]), t[k])
 64        observed = len(r.switches)
 65        predicted_zero = delta > 2*A
 66        rows.append({'delta': delta, 'noise_amplitude': A,
 67                     'predicted_zero_switch': predicted_zero,
 68                     'switches': observed})
 69    return rows
 70
 71def tau_sweep():
 72    # First-order h update h'=(s-h)/tau. For sinusoidal score, larger tau attenuates
 73    # fluctuations; measured RMS gain should follow 1/sqrt(1+(omega*tau)^2).
 74    dt, duration, A, period = 0.001, 20.0, 0.1, 1.0
 75    t = np.arange(0, duration, dt)
 76    omega = 2*np.pi/period
 77    s = 0.5 + A*np.sin(omega*t)
 78    rows=[]
 79    for tau in [0.005, 0.02, 0.1, 0.5]:
 80        h=0.5; vals=[]
 81        for v in s:
 82            h += dt/tau*(v-h)
 83            vals.append(h)
 84        vals=np.asarray(vals); trim=t > 5
 85        gain=(np.std(vals[trim])/A)
 86        predicted=1/np.sqrt(2*(1+(omega*tau)**2))
 87        rows.append({'tau_h': tau, 'predicted_rms_gain': predicted,
 88                     'observed_gain': float(gain), 'relative_error': float(abs(gain-predicted)/predicted)})
 89    return rows
 90
 91def regime_prediction_demo():
 92    # Same observed score is used by both methods. Hysteresis suppresses rapid regime
 93    # changes caused by score noise; soft routing averages incompatible dynamics.
 94    rng=np.random.default_rng(SEED+7); T=6000
 95    true=(np.arange(T)//30)%2
 96    # regimes have opposite slopes, noisy observations make scores ambiguous at boundaries
 97    x=np.zeros(T); x[0]=0.2
 98    for t in range(1,T): x[t]=0.92*x[t-1] + (0.18 if true[t] else -0.18) + 0.04*rng.normal()
 99    score=np.clip(true + 0.22*rng.normal(size=T),0,1)
100    # one-step model uses known regime-specific intercepts, but route is noisy
101    pred_soft=0.92*x[:-1] + 0.18*(2*score[1:]-1)
102    soft_mse=float(np.mean((pred_soft-x[1:])**2))
103    r=HystereticRouter(0.68,0.32); chosen=[]
104    for t in range(T): chosen.append(r.step(np.array([1-score[t],score[t]]),t))
105    chosen=np.asarray(chosen)
106    pred_h=0.92*x[:-1] + 0.18*(2*chosen[1:]-1)
107    hyst_mse=float(np.mean((pred_h-x[1:])**2))
108    return {'soft_mse':soft_mse, 'hysteretic_mse':hyst_mse,
109            'soft_switch_proxy':int(np.sum(np.diff((score>0.5).astype(int))!=0)),
110            'hysteretic_switches':len(r.switches), 'improvement_fraction':(soft_mse-hyst_mse)/soft_mse}
111
112def main():
113    dwell=dwell_sweep(); chatter=chatter_sweep(); tau=tau_sweep(); demo=regime_prediction_demo()
114    result={'seed':SEED, 'dwell_bound':dwell, 'chatter_threshold':chatter,
115            'timescale_lowpass':tau, 'regime_demo':demo}
116    Path('results.json').write_text(json.dumps(result, indent=2))
117    print(json.dumps(result, indent=2))
118    # core acceptance checks
119    assert min(x['bound_ratio'] for x in dwell) >= 0.99
120    assert chatter[-1]['switches'] == 0
121    assert max(x['relative_error'] for x in tau) < 0.02
122if __name__ == '__main__': main()