1import json
 2from pathlib import Path
 3import numpy as np
 4
 5# Convex information-coordinate toy for I_q(Y;Z|U)>=H.
 6# A binary symmetric channel supplies the attainable range, while r is a
 7# differentiable information-capacity coordinate in [I0,Imax].
 8def entropy_binary(q):
 9    q = np.clip(q, 1e-12, 1-1e-12)
10    return -(q*np.log(q)+(1-q)*np.log(1-q))
11
12def channel_info(p):
13    return float(np.log(2.)-entropy_binary(1-p))
14
15def run_dual(target, r0, rmax, k=2., eta_l=.05, lam_max=8., steps=3000):
16    # F(r)=k/2(r-r0)^2; exact inner minimizer is r=min(rmax,r0+lambda/k).
17    lam = 0.; ema = r0; trace=[]
18    for t in range(steps):
19        r = float(np.clip(r0 + lam/k, r0, rmax))
20        ema = .9*ema + .1*r
21        lam = float(np.clip(lam + eta_l*(target-ema), 0., lam_max))
22        if t % 10 == 0: trace.append((t,r,lam))
23    return r,lam,trace
24
25def fixed_information_drift(target, fixed, eta=.07, steps=100, clip=8.):
26    lam=0.
27    for _ in range(steps): lam=np.clip(lam+eta*(target-fixed),0.,clip)
28    return float(lam), float(steps*eta*(target-fixed))
29
30def main():
31    # p=.55 is the unconstrained channel; pmax defines attainable information.
32    r0, rmax = channel_info(.55), channel_info(.995)
33    targets=[.5*r0, r0+.35*(rmax-r0), rmax+.10]
34    rows=[]
35    for name,h in zip(('inactive','binding','saturated'),targets):
36        r,lam,_=run_dual(h,r0,rmax)
37        rows.append({'regime':name,'H':h,
38          'predicted_I':r0 if name=='inactive' else (h if name=='binding' else rmax),
39          'observed_I':r,'abs_I_error':abs(r-(r0 if name=='inactive' else (h if name=='binding' else rmax))),
40          'predicted_lambda':0. if name=='inactive' else (None if name=='binding' else 8.),
41          'observed_lambda':lam})
42    sweep=[]
43    for frac in (.15,.35,.55,.75,.90):
44        h=r0+frac*(rmax-r0); r,lam,_=run_dual(h,r0,rmax)
45        sweep.append({'H':h,'observed_I':r,'abs_error':abs(r-h),'lambda':lam})
46    drifts=[]
47    for delta in (-.08,.05,.30):
48        obs,pred=fixed_information_drift(r0+delta,r0)
49        drifts.append({'delta_H_minus_I':delta,'predicted_unclipped_lambda':pred,'observed_lambda':obs})
50    # Standard fixed coefficient: same primal with lambda=1, cannot track changing H.
51    fixed_r=min(rmax,r0+1/2.)
52    out={'channel':{'I0':r0,'Imax':rmax},
53      'predictions':{'P1':'H<I0 => lambda=0 and I=I0','P2':'I0<H<Imax => I tracks H','P3':'H>Imax => lambda clips and I=Imax','P4':'fixed-I lambda slope=eta*(H-I)'},
54      'regimes':rows,'feasible_target_sweep':sweep,'fixed_I_drift':drifts,
55      'baseline_fixed_lambda':{'lambda':1.,'I':fixed_r}}
56    Path('results.json').write_text(json.dumps(out,indent=2)); print(json.dumps(out,indent=2))
57if __name__=='__main__': main()