import json, math from pathlib import Path import numpy as np SEED = 1049 class RankController: def __init__(self, dmax, alpha=0.08, con=2.0, coff=1.2, M=3, cooldown=0): self.dmax, self.alpha, self.con, self.coff, self.M = dmax, alpha, con, coff, M self.C = np.zeros((dmax, dmax)); self.active = np.zeros(dmax, dtype=bool) self.above = np.zeros(dmax, dtype=int); self.below = np.zeros(dmax, dtype=int) self.cooldown = cooldown; self.cool = 0; self.events = [] def update(self, y, noise_var): self.C = (1-self.alpha)*self.C + self.alpha*np.outer(y, y) S = (self.C - noise_var*np.eye(self.dmax) + (self.C - noise_var*np.eye(self.dmax)).T)/2 vals, vecs = np.linalg.eigh(S); vals = vals[::-1]; vecs = vecs[:, ::-1] # In this diagonal toy, sorting is equivalent to ordered subspace eigenvalues. on = vals > self.con*noise_var; off = vals < self.coff*noise_var if self.cool: self.cool -= 1 for i in range(self.dmax): if on[i]: self.above[i] += 1 else: self.above[i] = 0 if off[i]: self.below[i] += 1 else: self.below[i] = 0 if not self.cool: for i in range(self.dmax): if (not self.active[i]) and self.above[i] >= self.M: self.active[i] = True; self.events.append((len(self.events), 'on', i)); self.cool=self.cooldown elif self.active[i] and self.below[i] >= self.M: self.active[i] = False; self.events.append((len(self.events), 'off', i)); self.cool=self.cooldown return vals, self.active.copy() def crossing_sweep(): # Prediction: stationary corrected population eigenvalue is signal variance s; # activation should occur iff s > c_on*sigma^2, approximately at ratio c_on. rng=np.random.default_rng(SEED); alpha=.08; con=2.; coff=1.2; M=3 ratios=np.array([0.5,1.,1.5,1.9,2.0,2.1,2.5,4.]) outcomes=[] for q in ratios: hits=[] for rep in range(20): ctl=RankController(1,alpha,con,coff,M) for _ in range(500): ctl.update(np.array([rng.normal(0, math.sqrt(q)) + rng.normal()]),1.) hits.append(bool(ctl.active[0])) outcomes.append(float(np.mean(hits))) # Estimate transition by interpolating the 50% point across a dense sweep. dense=np.linspace(.5,4,36); probs=[] for q in dense: h=[] for rep in range(12): ctl=RankController(1,alpha,con,coff,M) for _ in range(400): ctl.update(np.array([rng.normal(0,math.sqrt(q))+rng.normal()]),1.) h.append(ctl.active[0]) probs.append(np.mean(h)) trans=float(dense[np.argmin(np.abs(np.array(probs)-.5))]) # Deterministic population check: corrected covariance is exactly q, so crossing is q > c_on. deterministic = [] for q in ratios: ctl=RankController(1,alpha,con,coff,M) # Feed a constant-magnitude sequence whose sample variance is q+noise. for t in range(200): y=np.array([math.sqrt(q+1.) if t%2==0 else -math.sqrt(q+1.)]) ctl.update(y,1.) deterministic.append(bool(ctl.active[0])) 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]]} def delay_sweep(): # After an abrupt jump, EWMA expectation is C_t = s_new+(C0-s_new)(1-a)^t. # Crossing delay is the first t with C_t-noise > con*noise. rng=np.random.default_rng(SEED+1); alpha=.1; noise=1.; con=2.; M=1 s0=.2; s1=5.; target=(1+con) predicted=math.ceil(math.log((target-s1)/(s0+1-s1))/math.log(1-alpha)) delays=[] for rep in range(30): ctl=RankController(1,alpha,con,1.2,M) for _ in range(100): ctl.update(np.array([rng.normal(0,math.sqrt(s0))+rng.normal()]),noise) delay=None for t in range(1,150): _,a=ctl.update(np.array([rng.normal(0,math.sqrt(s1))+rng.normal()]),noise) if a[0]: delay=t; break delays.append(delay if delay is not None else 150) return {'predicted_first_crossing_steps':predicted,'observed_median_steps':float(np.median(delays)),'observed_mean_steps':float(np.mean(delays)),'delays':delays} def chatter_test(): # Near-threshold noisy signal: hysteresis should reduce toggles relative to one threshold. rng=np.random.default_rng(SEED+2); n=3000 def run(coff): ctl=RankController(1,.12,2.,coff,3) toggles=0; prev=False for _ in range(n): _,a=ctl.update(np.array([rng.normal(0,math.sqrt(1.9))+rng.normal()]),1.) toggles += int(a[0] != prev); prev=bool(a[0]) return toggles, len(ctl.events), bool(prev) # same random distribution, reset stream for fair comparison rng=np.random.default_rng(SEED+2); h=run(1.2) rng=np.random.default_rng(SEED+2); no=run(2.0) return {'hysteresis_coff_1.2':h,'single_threshold_coff_2.0':no,'toggle_reduction_fraction':1-h[0]/max(no[0],1)} def switched_rank_demo(): # Minimal order-adaptive reconstruction: observations are independent latent channels + sensor noise. # Fixed rank-2, fixed rank-5, and controller rank; report MSE on clean latent reconstruction. rng=np.random.default_rng(SEED+3); d=5; T=2400; noise=.35 ys=[]; xs=[]; true=[] for t in range(T): r=2 if t<700 or t>=1700 else 5 x=np.zeros(d); x[:r]=rng.normal(size=r) y=x+rng.normal(0,noise,size=d); ys.append(y); xs.append(x); true.append(r) ys=np.array(ys); xs=np.array(xs) ctl=RankController(d,.1,2.,1.2,3); pred=[]; ranks=[] for y in ys: _,a=ctl.update(y,noise**2); ranks.append(int(a.sum())); pred.append(y*a) # oracle fixed masks are idealized standard fixed-width references mse_ad=float(np.mean((np.array(pred)-xs)**2)) mse2=float(np.mean((ys*np.array([1,1,0,0,0])-xs)**2)) mse5=float(np.mean((ys-xs)**2)) 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))]} def main(): out={'seed':SEED,'crossing_sweep':crossing_sweep(),'ewma_delay':delay_sweep(),'hysteresis':chatter_test(),'switched_demo':switched_rank_demo()} Path('results.json').write_text(json.dumps(out,indent=2)) print(json.dumps(out,indent=2)) if __name__=='__main__': main()