1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1049
  6
  7class RankController:
  8    def __init__(self, dmax, alpha=0.08, con=2.0, coff=1.2, M=3, cooldown=0):
  9        self.dmax, self.alpha, self.con, self.coff, self.M = dmax, alpha, con, coff, M
 10        self.C = np.zeros((dmax, dmax)); self.active = np.zeros(dmax, dtype=bool)
 11        self.above = np.zeros(dmax, dtype=int); self.below = np.zeros(dmax, dtype=int)
 12        self.cooldown = cooldown; self.cool = 0; self.events = []
 13    def update(self, y, noise_var):
 14        self.C = (1-self.alpha)*self.C + self.alpha*np.outer(y, y)
 15        S = (self.C - noise_var*np.eye(self.dmax) + (self.C - noise_var*np.eye(self.dmax)).T)/2
 16        vals, vecs = np.linalg.eigh(S); vals = vals[::-1]; vecs = vecs[:, ::-1]
 17        # In this diagonal toy, sorting is equivalent to ordered subspace eigenvalues.
 18        on = vals > self.con*noise_var; off = vals < self.coff*noise_var
 19        if self.cool: self.cool -= 1
 20        for i in range(self.dmax):
 21            if on[i]: self.above[i] += 1
 22            else: self.above[i] = 0
 23            if off[i]: self.below[i] += 1
 24            else: self.below[i] = 0
 25        if not self.cool:
 26            for i in range(self.dmax):
 27                if (not self.active[i]) and self.above[i] >= self.M:
 28                    self.active[i] = True; self.events.append((len(self.events), 'on', i)); self.cool=self.cooldown
 29                elif self.active[i] and self.below[i] >= self.M:
 30                    self.active[i] = False; self.events.append((len(self.events), 'off', i)); self.cool=self.cooldown
 31        return vals, self.active.copy()
 32
 33def crossing_sweep():
 34    # Prediction: stationary corrected population eigenvalue is signal variance s;
 35    # activation should occur iff s > c_on*sigma^2, approximately at ratio c_on.
 36    rng=np.random.default_rng(SEED); alpha=.08; con=2.; coff=1.2; M=3
 37    ratios=np.array([0.5,1.,1.5,1.9,2.0,2.1,2.5,4.])
 38    outcomes=[]
 39    for q in ratios:
 40        hits=[]
 41        for rep in range(20):
 42            ctl=RankController(1,alpha,con,coff,M)
 43            for _ in range(500): ctl.update(np.array([rng.normal(0, math.sqrt(q)) + rng.normal()]),1.)
 44            hits.append(bool(ctl.active[0]))
 45        outcomes.append(float(np.mean(hits)))
 46    # Estimate transition by interpolating the 50% point across a dense sweep.
 47    dense=np.linspace(.5,4,36); probs=[]
 48    for q in dense:
 49        h=[]
 50        for rep in range(12):
 51            ctl=RankController(1,alpha,con,coff,M)
 52            for _ in range(400): ctl.update(np.array([rng.normal(0,math.sqrt(q))+rng.normal()]),1.)
 53            h.append(ctl.active[0])
 54        probs.append(np.mean(h))
 55    trans=float(dense[np.argmin(np.abs(np.array(probs)-.5))])
 56    # Deterministic population check: corrected covariance is exactly q, so crossing is q > c_on.
 57    deterministic = []
 58    for q in ratios:
 59        ctl=RankController(1,alpha,con,coff,M)
 60        # Feed a constant-magnitude sequence whose sample variance is q+noise.
 61        for t in range(200):
 62            y=np.array([math.sqrt(q+1.) if t%2==0 else -math.sqrt(q+1.)])
 63            ctl.update(y,1.)
 64        deterministic.append(bool(ctl.active[0]))
 65    return {'predicted_ratio':con,'swept_ratios':ratios.tolist(),'activation_probability':outcomes,'observed_50pct_ratio':trans,'deterministic_activation':deterministic,'deterministic_boundary_between': [float(ratios[i]) for i in range(len(ratios)-1) if deterministic[i]!=deterministic[i+1]]}
 66
 67def delay_sweep():
 68    # After an abrupt jump, EWMA expectation is C_t = s_new+(C0-s_new)(1-a)^t.
 69    # Crossing delay is the first t with C_t-noise > con*noise.
 70    rng=np.random.default_rng(SEED+1); alpha=.1; noise=1.; con=2.; M=1
 71    s0=.2; s1=5.; target=(1+con)
 72    predicted=math.ceil(math.log((target-s1)/(s0+1-s1))/math.log(1-alpha))
 73    delays=[]
 74    for rep in range(30):
 75        ctl=RankController(1,alpha,con,1.2,M)
 76        for _ in range(100): ctl.update(np.array([rng.normal(0,math.sqrt(s0))+rng.normal()]),noise)
 77        delay=None
 78        for t in range(1,150):
 79            _,a=ctl.update(np.array([rng.normal(0,math.sqrt(s1))+rng.normal()]),noise)
 80            if a[0]: delay=t; break
 81        delays.append(delay if delay is not None else 150)
 82    return {'predicted_first_crossing_steps':predicted,'observed_median_steps':float(np.median(delays)),'observed_mean_steps':float(np.mean(delays)),'delays':delays}
 83
 84def chatter_test():
 85    # Near-threshold noisy signal: hysteresis should reduce toggles relative to one threshold.
 86    rng=np.random.default_rng(SEED+2); n=3000
 87    def run(coff):
 88        ctl=RankController(1,.12,2.,coff,3)
 89        toggles=0; prev=False
 90        for _ in range(n):
 91            _,a=ctl.update(np.array([rng.normal(0,math.sqrt(1.9))+rng.normal()]),1.)
 92            toggles += int(a[0] != prev); prev=bool(a[0])
 93        return toggles, len(ctl.events), bool(prev)
 94    # same random distribution, reset stream for fair comparison
 95    rng=np.random.default_rng(SEED+2); h=run(1.2)
 96    rng=np.random.default_rng(SEED+2); no=run(2.0)
 97    return {'hysteresis_coff_1.2':h,'single_threshold_coff_2.0':no,'toggle_reduction_fraction':1-h[0]/max(no[0],1)}
 98
 99def switched_rank_demo():
100    # Minimal order-adaptive reconstruction: observations are independent latent channels + sensor noise.
101    # Fixed rank-2, fixed rank-5, and controller rank; report MSE on clean latent reconstruction.
102    rng=np.random.default_rng(SEED+3); d=5; T=2400; noise=.35
103    ys=[]; xs=[]; true=[]
104    for t in range(T):
105        r=2 if t<700 or t>=1700 else 5
106        x=np.zeros(d); x[:r]=rng.normal(size=r)
107        y=x+rng.normal(0,noise,size=d); ys.append(y); xs.append(x); true.append(r)
108    ys=np.array(ys); xs=np.array(xs)
109    ctl=RankController(d,.1,2.,1.2,3); pred=[]; ranks=[]
110    for y in ys:
111        _,a=ctl.update(y,noise**2); ranks.append(int(a.sum())); pred.append(y*a)
112    # oracle fixed masks are idealized standard fixed-width references
113    mse_ad=float(np.mean((np.array(pred)-xs)**2))
114    mse2=float(np.mean((ys*np.array([1,1,0,0,0])-xs)**2))
115    mse5=float(np.mean((ys-xs)**2))
116    return {'adaptive_mse':mse_ad,'fixed_rank2_mse':mse2,'fixed_rank5_mse':mse5,'mean_active_rank':float(np.mean(ranks)),'segment_mean_ranks':[float(np.mean(ranks[:700])),float(np.mean(ranks[700:1700])),float(np.mean(ranks[1700:]))], 'segment_mse_adaptive':[float(np.mean((np.array(pred[:700])-xs[:700])**2)),float(np.mean((np.array(pred[700:1700])-xs[700:1700])**2)),float(np.mean((np.array(pred[1700:])-xs[1700:])**2))]}
117
118def main():
119    out={'seed':SEED,'crossing_sweep':crossing_sweep(),'ewma_delay':delay_sweep(),'hysteresis':chatter_test(),'switched_demo':switched_rank_demo()}
120    Path('results.json').write_text(json.dumps(out,indent=2))
121    print(json.dumps(out,indent=2))
122if __name__=='__main__': main()