1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2710
  6rng = np.random.default_rng(SEED)
  7
  8# Three-state clockwise cycle: i -> i+1 at a and i -> i-1 at b.
  9# Its stationary distribution is uniform, sigma=(a-b) log(a/b),
 10# and the nonconstant output mode has tau=1/(1.5*(a+b)).
 11def exact_stats(a, b):
 12    sigma = (a - b) * math.log(a / b)
 13    tau = 1.0 / (1.5 * (a + b))
 14    return sigma, tau
 15
 16def simulate_ctmc(a, b, dt=0.01, steps=300000, burn=30000):
 17    """Uniformized discrete simulation; estimates p, q, sigma and output tau."""
 18    z = 0
 19    states = np.empty(steps-burn, dtype=np.int8)
 20    # y is a centered observable with nonzero projection on both complex modes.
 21    yy = np.array([1.0, -0.5, -0.5])
 22    counts = np.zeros((3,3), dtype=np.int64)
 23    occup = np.zeros(3, dtype=np.int64)
 24    # Exact one-step transition probabilities to avoid multiple jumps.
 25    # dt is small for all rates used below.
 26    for t in range(steps):
 27        u = rng.random()
 28        if u < a*dt:
 29            zn = (z + 1) % 3
 30        elif u < (a+b)*dt:
 31            zn = (z - 1) % 3
 32        else:
 33            zn = z
 34        if t >= burn:
 35            states[t-burn] = z
 36            occup[z] += 1
 37            if zn != z:
 38                counts[z, zn] += 1
 39        z = zn
 40    p = occup / occup.sum()
 41    qhat = counts / (occup[:, None] * dt + 1e-30)
 42    sig = 0.0
 43    for i in range(3):
 44        for j in range(i+1, 3):
 45            f = p[i]*qhat[i,j]
 46            g = p[j]*qhat[j,i]
 47            if f > 0 and g > 0:
 48                sig += (f-g)*math.log(f/g)
 49    # Correlation integral via FFT, truncated where correlation is negligible.
 50    y = yy[states].astype(float)
 51    y -= y.mean()
 52    n = len(y)
 53    f = np.fft.rfft(y, 2*n)
 54    ac = np.fft.irfft(f*np.conjugate(f), 2*n)[:n]
 55    ac /= np.arange(n, 0, -1)
 56    ac /= ac[0]
 57    cutoff = np.where(ac < 0.02)[0]
 58    m = int(cutoff[0]) if len(cutoff) else min(n-1, int(10/(1.5*(a+b)*dt)))
 59    tauhat = max(0.0, np.trapz(ac[:m], dx=dt))
 60    return sig, tauhat, p
 61
 62def tracking_mse(rate, freq, amp=1.0, dt=0.002, duration=80.0):
 63    """Mean-field response of the stochastic mode to a changing target.
 64    This is the exact mean response of a symmetric reset-like Markov decoder.
 65    """
 66    n = int(duration/dt)
 67    t = np.arange(n)*dt
 68    target = amp*np.sin(2*np.pi*freq*t)
 69    out = 0.0
 70    err = np.empty(n)
 71    for k in range(n):
 72        out += dt*rate*(target[k]-out)
 73        err[k] = out-target[k]
 74    burn = n//4
 75    return float(np.mean(err[burn:]**2)), float(np.mean((2*np.pi*freq*amp*np.cos(2*np.pi*freq*t[burn:]))**2))
 76
 77def main():
 78    # Cheap sanity check of the exact entropy-production expression on random rates.
 79    random_checks = []
 80    for _ in range(1000):
 81        a, b = 10**rng.uniform(-1, 1, 2)
 82        s, _ = exact_stats(a,b)
 83        random_checks.append(s >= -1e-12)
 84
 85    # Prediction 1: at fixed bias ratio, sigma is linear in global activity.
 86    ratio = 3.0
 87    rates = np.array([0.25, 0.5, 1.0, 2.0, 4.0])
 88    sigmas = np.array([exact_stats(r*ratio/(ratio+1), r/(ratio+1))[0] for r in rates])
 89    slope_sigma = np.polyfit(np.log(rates), np.log(sigmas), 1)[0]
 90
 91    # Prediction 2: the output memory time is inverse in activity.
 92    taus = np.array([exact_stats(r*ratio/(ratio+1), r/(ratio+1))[1] for r in rates])
 93    slope_tau = np.polyfit(np.log(rates), np.log(taus), 1)[0]
 94
 95    # Estimate both expressions once, checking numerical estimators against exact values.
 96    a, b = ratio/(ratio+1), 1/(ratio+1)
 97    exact_s, exact_tau = exact_stats(a,b)
 98    est_s, est_tau, p = simulate_ctmc(a,b)
 99
100    # Prediction 3: claimed budget says epsilon*sigma >= kappa*tau*v^2.
101    # Test the dimensionless ratio with kappa=1 over rate and frequency sweeps.
102    budget_rows = []
103    for r in [0.5, 1.0, 2.0, 4.0]:
104        aa, bb = r*ratio/(ratio+1), r/(ratio+1)
105        s, tau = exact_stats(aa,bb)
106        for f in [0.01, 0.03, 0.1, 0.3]:
107            eps, v2 = tracking_mse(r, f)
108            budget_rows.append({"rate":r, "freq":f, "sigma":s, "tau":tau,
109                                "epsilon":eps, "v2":v2,
110                                "budget_ratio":eps*s/(tau*v2)})
111
112    # Tiny baseline comparison: fixed slow deterministic low-pass vs budget-selected activity.
113    # The proposed controller selects the smallest tested rate whose estimated budget ratio
114    # exceeds one; this is deliberately reported as a test, not a fitted claim.
115    base = []
116    idea = []
117    for f in [0.03, 0.1, 0.3]:
118        eb, _ = tracking_mse(1.0, f)
119        choices = [(r, tracking_mse(r,f)[0]) for r in [0.5,1,2,4]]
120        chosen = next((x for x in choices if next(q["budget_ratio"] for q in budget_rows if q["rate"]==x[0] and q["freq"]==f) >= 1), choices[-1])
121        base.append(eb); idea.append(chosen[1])
122
123    result = {
124      "seed": SEED,
125      "math_sanity": {"all_random_entropy_rates_nonnegative": all(random_checks),
126                       "exact_sigma": exact_s, "estimated_sigma": est_s,
127                       "exact_tau": exact_tau, "estimated_tau": est_tau,
128                       "stationary_p": p.tolist()},
129      "predictions": {
130        "P1_sigma_activity_exponent_predicted": 1.0,
131        "P1_sigma_activity_exponent_observed": float(slope_sigma),
132        "P2_tau_activity_exponent_predicted": -1.0,
133        "P2_tau_activity_exponent_observed": float(slope_tau),
134        "P3_budget_ratio_should_be_ge_kappa_1": True,
135        "P3_budget_ratio_min": float(min(x["budget_ratio"] for x in budget_rows)),
136        "P3_budget_ratio_median": float(np.median([x["budget_ratio"] for x in budget_rows])),
137        "P3_budget_ratio_max": float(max(x["budget_ratio"] for x in budget_rows))
138      },
139      "budget_sweep": budget_rows,
140      "mini_comparison": {"frequencies":[0.03,0.1,0.3], "baseline_fixed_rate1_mse":base,
141                          "budget_controller_mse":idea}
142    }
143    Path("results.json").write_text(json.dumps(result, indent=2))
144    print(json.dumps(result, indent=2))
145
146if __name__ == "__main__":
147    main()