Open the walk-count sub-problem (route-b combinatorial core) + W2

telperion/docs/bg_c_convexity.py

@@ -0,0 +1,46 @@
+"""Probe (c): can degree-sequence convexity prove m_2 >= phi(m_1) near the extremum band?
+Derived crude bound (Cauchy-Schwarz + Jensen):  m_2 >= 2 m_1^2 - m_1.  Test tightness vs caterpillar.
+Also fit best m_2 >= a*m_1^2+b*m_1+c and check the gap to the caterpillar boundary near m_1~0.52."""
+import sys; sys.path.insert(0,'telperion/src')
+from fractions import Fraction as F
+import numpy as np, networkx as nx
+def dm(n,e):
+    d=[0]*n;adj=[[] for _ in range(n)]
+    for a,b in e:d[a]+=1;d[b]+=1;adj[a].append(b);adj[b].append(a)
+    return d,adj
+def m12(n,e):
+    d,adj=dm(n,e);m1=F(0);m2=F(0)
+    for v in range(n):
+        dv=d[v];S=sum(F(1,d[a]) for a in adj[v]);Q=sum(F(1,d[a]**2) for a in adj[v])
+        m1+=S/dv;m2+=2*S*S/(dv*dv)-Q/(dv*dv)
+    return float(m1/n),float(m2/n)
+def cat(sp,a,L):
+    e=[];nid=sp
+    for i in range(sp-1):e.append((i,i+1))
+    for i in range(sp):
+        for _ in range(a):
+            p=i
+            for _ in range(L):e.append((p,nid));p=nid;nid+=1
+    return nid,e
+# crude bound test over all trees
+viol=0;tot=0;minslack=9
+for n in range(2,14):
+    for T in nx.nonisomorphic_trees(n):
+        idx={v:i for i,v in enumerate(T.nodes())};e=[(idx[a],idx[b]) for a,b in T.edges()]
+        m1,m2=m12(n,e);slack=m2-(2*m1*m1-m1)
+        tot+=1
+        if slack<-1e-12:viol+=1
+        minslack=min(minslack,slack)
+print(f"crude bound m_2 >= 2 m_1^2 - m_1: violations={viol}/{tot}, min slack={minslack:.5f}")
+# caterpillar boundary values near m1~0.52 and the crude-bound value there
+print("\ncaterpillar boundary vs crude bound vs true m_2:")
+for a in [5,7,9]:
+    n,e=cat(40,a,2);m1,m2=m12(n,e)
+    print(f"  a={a}: m1={m1:.4f} true m2={m2:.4f}  crude 2m1^2-m1={2*m1*m1-m1:.4f}  (crude is {'USELESS' if 2*m1*m1-m1 < m2-0.1 else 'ok'})")
+# how much of m_2 is 'variance' (2 avg x^2 - 2 m1^2) that the crude bound throws away?
+print("\nvariance content: m_2 = 2*avg(x^2) - avg(Q/d^2); caterpillar avg(x^2) vs m1^2:")
+for a in [7]:
+    n,e=cat(40,a,2);d,adj=dm(n,e)
+    xs=[float(sum(F(1,d[k]) for k in adj[v])/d[v]) for v in range(n)]
+    ax2=np.mean(np.array(xs)**2);m1=np.mean(xs)
+    print(f"  a={a}: avg(x^2)={ax2:.4f}  m1^2={m1*m1:.4f}  variance avg(x^2)-m1^2={ax2-m1*m1:.4f} (LOST by Jensen)")

telperion/docs/bg_cavity.py

@@ -0,0 +1,94 @@
+"""Exact cavity (Bethe) free energy for the matching / monomer-dimer model on trees.
+
+Edge weight w_{uv} = 1/(d_u d_v).  Z = sum_matchings prod_{e in M} w_e = per(L)/prod(deg) = prod_{lam>0}(1+lam^2).
+Cavity messages on directed edges:
+    x_{u->v} = sum_{c ~ u, c != v} w_{uc} / (1 + x_{c->u})            (q_{c->u} = 1/(1+x_{c->u}))
+Exact Bethe free energy (exact on trees):
+    log Z = sum_v log(1 + sum_{a~v} w_{va} q_{a->v})  -  sum_{(u,v) in E} log(1 + w_{uv} q_{u->v} q_{v->u}).
+Verify log Z_Bethe == log(per/prod deg) on all trees n<=8, then compute the caterpillar's fixed point and
+its density, confirming F -> log rho*.
+"""
+import sys, math
+from fractions import Fraction as F
+sys.path.insert(0, 'telperion/src')
+import networkx as nx
+from telperion.girardeau import hard_core_boson_partition
+
+RHO = 1.2276458
+LOG_RHO = math.log(RHO)
+
+
+def cavity_messages(n, edges, iters=2000, tol=1e-14):
+    d = [0]*n; adj = [[] for _ in range(n)]
+    for a, b in edges:
+        d[a] += 1; d[b] += 1; adj[a].append(b); adj[b].append(a)
+    w = {}
+    for a, b in edges:
+        w[(a, b)] = w[(b, a)] = 1.0/(d[a]*d[b])
+    x = {}
+    for a, b in edges:
+        x[(a, b)] = x[(b, a)] = 0.0    # directed messages x[(u,v)] = message u->v
+    for _ in range(iters):
+        mx = 0.0
+        newx = {}
+        for (u, v) in x:
+            s = 0.0
+            for c in adj[u]:
+                if c == v: continue
+                s += w[(u, c)] / (1.0 + x[(c, u)])
+            newx[(u, v)] = s
+        for k in x:
+            mx = max(mx, abs(newx[k]-x[k]))
+        x = newx
+        if mx < tol: break
+    return d, adj, w, x
+
+
+def bethe_logZ(n, edges):
+    d, adj, w, x = cavity_messages(n, edges)
+    q = {k: 1.0/(1.0+x[k]) for k in x}
+    vsum = 0.0
+    for v in range(n):
+        a_v = 1.0 + sum(w[(v, a)]*q[(a, v)] for a in adj[v])
+        vsum += math.log(a_v)
+    esum = 0.0
+    seen = set()
+    for a, b in edges:
+        e = (min(a, b), max(a, b))
+        if e in seen: continue
+        seen.add(e)
+        esum += math.log(1.0 + w[(a, b)]*q[(a, b)]*q[(b, a)])
+    return vsum - esum
+
+
+def edges_of(T):
+    idx = {v: i for i, v in enumerate(T.nodes())}
+    return T.number_of_nodes(), [(idx[a], idx[b]) for a, b in T.edges()]
+
+
+print("=== verify Bethe logZ == log(per/prod deg) over all trees n<=8 ===")
+worst = 0.0; tot = 0
+for nn in range(2, 9):
+    for T in nx.nonisomorphic_trees(nn):
+        m, e = edges_of(T); tot += 1
+        exact = math.log(float(hard_core_boson_partition(m, e)))
+        beth = bethe_logZ(m, e)
+        worst = max(worst, abs(exact-beth))
+print(f"  {tot} trees, max |logZ_Bethe - log(per/prod)| = {worst:.2e}")
+
+
+def caterpillar_legs(sp, a, L):
+    e = []; nid = sp
+    for i in range(sp-1): e.append((i, i+1))
+    for i in range(sp):
+        for _ in range(a):
+            p = i
+            for _ in range(L): e.append((p, nid)); p = nid; nid += 1
+    return nid, e
+
+
+print("\n=== caterpillar cavity density (bulk) vs log rho* ===")
+for a in (5, 7, 9):
+    n, e = caterpillar_legs(60, a, 2)
+    F_density = bethe_logZ(n, e)/n
+    print(f"  a={a}: F = {F_density:.6f}   (log rho* = {LOG_RHO:.6f})")

telperion/docs/bg_cavity_caterpillar.py

@@ -0,0 +1,48 @@
+"""Exact infinite-caterpillar cavity fixed point + free-energy density F(a), maximised over arm-count a.
+
+Cell = 1 hub (deg a+2) + a arm-mids (deg 2) + a leaves (deg 1); weights w=1/(d_u d_v).
+Messages (fixed point): x_leaf->AM=0; x_AM->H=1/2; X:=x_H->H and x_H->AM solve the hub cavity;
+x_AM->L closes the arm.  Density F(a) = per-cell Bethe free energy / (2a+1).  Confirm max_a F(a)=log rho*.
+"""
+import math
+from scipy.optimize import brentq, minimize_scalar
+
+LOG_RHO = math.log(1.2276458)
+
+
+def fixed_point(a):
+    wHH = 1.0/(a+2)**2
+    wHA = 1.0/(2*(a+2))
+    wAL = 0.5
+    xAMH = wAL/(1.0+0.0)          # arm-mid -> hub  = 1/2
+    # X = x_H->H  solves  X = wHH/(1+X) + (a * wHA)/(1+xAMH)
+    c = a*wHA/(1.0+xAMH)
+    X = brentq(lambda X: wHH/(1.0+X) + c - X, 0.0, 2.0)
+    xHAM = 2*wHH/(1.0+X) + (a-1)*wHA/(1.0+xAMH)   # hub -> arm-mid (excludes 1 arm-mid)
+    xAML = wHA/(1.0+xHAM)                          # arm-mid -> leaf
+    return dict(wHH=wHH, wHA=wHA, wAL=wAL, X=X, xAMH=xAMH, xHAM=xHAM, xAML=xAML)
+
+
+def density(a):
+    p = fixed_point(a)
+    q = lambda x: 1.0/(1.0+x)
+    qHH = q(p['X']); qAMH = q(p['xAMH']); qHAM = q(p['xHAM']); qAML = q(p['xAML']); qLAM = 1.0
+    A_H = 1.0 + 2*p['wHH']*qHH + a*p['wHA']*qAMH
+    A_AM = 1.0 + p['wHA']*qHAM + p['wAL']*qLAM
+    A_L = 1.0 + p['wAL']*qAML
+    B_HH = 1.0 + p['wHH']*qHH*qHH
+    B_HA = 1.0 + p['wHA']*qHAM*qAMH
+    B_AL = 1.0 + p['wAL']*qAML*qLAM
+    vsum = math.log(A_H) + a*math.log(A_AM) + a*math.log(A_L)
+    esum = math.log(B_HH) + a*math.log(B_HA) + a*math.log(B_AL)
+    return (vsum - esum)/(2*a+1)
+
+
+print(" a  |   F(a)     | vs log rho*")
+for a in (5, 6, 7, 8, 9, 10):
+    print(f" {a:2d} | {density(a):.6f} | {density(a)-LOG_RHO:+.6f}")
+
+res = minimize_scalar(lambda a: -density(a), bounds=(4, 12), method='bounded')
+astar = res.x
+print(f"\ncontinuous max: a* = {astar:.4f}   F(a*) = {density(astar):.6f}   log rho* = {LOG_RHO:.6f}")
+print(f"  => cavity infinite-caterpillar density maximised at a*~{astar:.2f}, equals log rho* to {abs(density(astar)-LOG_RHO):.2e}")

telperion/docs/bg_cavity_contraction.py

@@ -0,0 +1,97 @@
+"""Lead 1 foundation: is the matching cavity map a CONTRACTION (Heilmann-Lieb), and does the
+finite-DEPTH truncation of the free energy converge GEOMETRICALLY to log rho*?
+
+Cavity map on a directed edge:  x_{u->v} = sum_{c~u, c!=v} w_{uc}/(1+x_{c->u}),  w=1/(d_u d_c).
+Jacobian magnitude  sum_c w_{uc}/(1+x_{c->u})^2 < (d_u-1)/d_u < 1  => contraction.
+
+Test A: propagate a perturbation of the leaf message up an arm/spine of the caterpillar; measure decay.
+Test B: compute the caterpillar free energy with messages truncated at depth d (set to 0 below depth d),
+        vs the exact fixed point; error should be ~ rho_contract^d (geometric).  If so, Lead 1 gives
+        F(T) <= log rho* + C rho^d -- a GEOMETRICALLY convergent bound (better than W9's slow moment hierarchy).
+"""
+import math
+import numpy as np
+from scipy.optimize import brentq
+
+LOG_RHO = math.log(1.2276458)
+
+
+def caterpillar_edges(spine_len, a, leg=2):
+    e = []; nid = spine_len
+    for i in range(spine_len - 1):
+        e.append((i, i + 1))
+    for i in range(spine_len):
+        for _ in range(a):
+            p = i
+            for _ in range(leg):
+                e.append((p, nid)); p = nid; nid += 1
+    return nid, e
+
+
+def cavity_fixed(n, edges, iters, init=0.0):
+    """Run the cavity recursion for `iters` sweeps from init; return messages + degrees/adj/w."""
+    d = [0] * n; adj = [[] for _ in range(n)]
+    for a, b in edges:
+        d[a] += 1; d[b] += 1; adj[a].append(b); adj[b].append(a)
+    w = {}
+    for a, b in edges:
+        w[(a, b)] = w[(b, a)] = 1.0 / (d[a] * d[b])
+    x = {(a, b): init for a, b in edges}
+    x.update({(b, a): init for a, b in edges})
+    for _ in range(iters):
+        nx = {}
+        for (u, v) in x:
+            nx[(u, v)] = sum(w[(u, c)] / (1.0 + x[(c, u)]) for c in adj[u] if c != v)
+        x = nx
+    return d, adj, w, x
+
+
+def bethe_density(n, edges, iters, init=0.0):
+    d, adj, w, x = cavity_fixed(n, edges, iters, init)
+    q = {k: 1.0 / (1.0 + x[k]) for k in x}
+    vs = sum(math.log(1.0 + sum(w[(v, a)] * q[(a, v)] for a in adj[v])) for v in range(n))
+    seen = set(); es = 0.0
+    for a, b in edges:
+        e = (min(a, b), max(a, b))
+        if e in seen: continue
+        seen.add(e)
+        es += math.log(1.0 + w[(a, b)] * q[(a, b)] * q[(b, a)])
+    return (vs - es) / n
+
+
+# Test A: contraction rate along a long arm/spine
+print("=== Test A: perturbation decay (contraction factor) ===")
+n, e = caterpillar_edges(30, 7, 2)
+d, adj, w, x0 = cavity_fixed(n, e, 400)          # converged
+# perturb one leaf's message, re-propagate a few sweeps, measure how the change shrinks per hop
+import copy
+xp = dict(x0)
+# find a leaf and its directed edge to arm-mid
+deg = d
+leaf = next(v for v in range(n) if deg[v] == 1)
+am = adj[leaf][0]
+xp[(leaf, am)] += 0.1
+# one synchronous sweep: measure max change at distance-k edges (crude contraction proxy)
+prev = 0.1
+for hop in range(1, 8):
+    nx = {}
+    for (u, v) in xp:
+        nx[(u, v)] = sum(w[(u, c)] / (1.0 + xp[(c, u)]) for c in adj[u] if c != v)
+    diffs = [abs(nx[k] - x0[k]) for k in nx]
+    mx = max(diffs)
+    print(f"  sweep {hop}: max|deviation from fixed pt| = {mx:.3e}  ratio={mx/prev:.3f}")
+    prev = mx if mx > 0 else prev
+    xp = nx
+
+# Test B: depth-truncated free energy convergence
+print("\n=== Test B: free energy vs cavity ITERATIONS (depth) -> geometric convergence to log rho* ===")
+n, e = caterpillar_edges(40, 7, 2)
+Fstar = bethe_density(n, e, 500)
+print(f"  converged F (a=7, spine 40) = {Fstar:.6f}   log rho* = {LOG_RHO:.6f}")
+prev_err = None
+for it in range(1, 12):
+    Fit = bethe_density(n, e, it, init=0.0)
+    err = abs(Fit - Fstar)
+    ratio = (err / prev_err) if prev_err and prev_err > 0 else float('nan')
+    print(f"  iters={it:2d}: F={Fit:.6f}  |F-F*|={err:.3e}  ratio={ratio:.3f}")
+    prev_err = err

telperion/docs/bg_cavity_potential.py

@@ -0,0 +1,80 @@
+"""Cavity-potential bound: does a discharge potential P(x) on cavity messages prove F(T) <= log rho*?
+
+Per-vertex free energy (edges split 1/2 to each end):  pv(v) = log A_v - (1/2) sum_a log B_{va},
+  A_v = 1 + sum_a q_a/(d*d_a),  q_a = 1/(1+x_a),  x_a = incoming message from neighbour a (degree d_a),
+  x_{v->a} = sum_{c != a} q_c/(d*d_c)  (outgoing),  B_{va} = 1 + q_a q_{v->a}/(d*d_a).
+Discharge telescopes on trees:  sum_v sum_a [P(x_a) - P(x_{v->a})] = 0.
+Certificate: exists P and bound B with   pv(v) - sum_a [P(x_a) - P(x_{v->a})] <= B   for every local config
+(degree d, neighbours (d_a, x_a)).  If min B == log rho*, the cavity potential proves the density bound
+(tight at the caterpillar) -- the exact route, in message space, NOT moments.
+
+Config space: degrees 1..DMAX, messages on a grid; P piecewise-linear on the grid (linear interpolation ->
+constraints linear in P's grid values).  Cutting-plane over configs.
+"""
+import sys, itertools, math
+import numpy as np
+from scipy.optimize import linprog
+sys.path.insert(0, 'telperion/src')
+LOG_RHO = math.log(1.2276458)
+
+DMAX = 4
+GRID = np.linspace(0.0, 1.0, 11)      # message grid for P (piecewise linear)
+NG = len(GRID)
+
+def interp_row(x):
+    """row over grid s.t. row.P = P(x) by linear interpolation; x clamped to [0,1]."""
+    x = min(max(x, 0.0), 1.0)
+    j = min(int(x*(NG-1)), NG-2)
+    t = (x - GRID[j])/(GRID[j+1]-GRID[j])
+    r = np.zeros(NG); r[j] = 1-t; r[j+1] = t
+    return r
+
+# per-vertex config = (d, tuple of (d_a, x_a)).  message-grid values for neighbours' x_a
+MSG = [0.0, 0.1, 0.2, 1.0/3, 0.4, 0.5]   # candidate incoming message values (coarse)
+
+def pv_and_disc(d, nbrs):
+    """nbrs = list of (d_a, x_a). Returns (pv, discharge_row over P-grid)."""
+    q = [1.0/(1.0+x) for (_, x) in nbrs]
+    A = 1.0 + sum(q[i]/(d*nbrs[i][0]) for i in range(d))
+    disc = np.zeros(NG)
+    pv = math.log(A)
+    for i in range(d):
+        d_a, x_a = nbrs[i]
+        x_out = sum(q[c]/(d*nbrs[c][0]) for c in range(d) if c != i)
+        q_out = 1.0/(1.0+x_out)
+        B = 1.0 + q[i]*q_out/(d*d_a)
+        pv -= 0.5*math.log(B)
+        disc += interp_row(x_a) - interp_row(x_out)     # +P(x_a) - P(x_out)
+    return pv, disc
+
+# enumerate configs: degree d, multiset of (d_a in 1..DMAX, x_a in MSG)
+def configs():
+    nb_types = [(da, xa) for da in range(1, DMAX+1) for xa in MSG]
+    for d in range(1, DMAX+1):
+        for combo in itertools.combinations_with_replacement(nb_types, d):
+            yield d, list(combo)
+
+CONF = list(configs())
+print(f"DMAX={DMAX}, grid={NG}, {len(CONF)} configs")
+
+# LP: variables [P_0..P_{NG-1}, B, beta]; per-config  pv - disc.P - beta*d <= B
+# density bound (handshake sum d = 2n-2):  F(T) <= B + beta*(2 - 2/n) -> B + 2 beta (bulk).  Minimize B+2beta.
+NV = NG + 2
+BI, BETA = NG, NG+1
+A_ub = []; b_ub = []
+for d, nbrs in CONF:
+    pv, disc = pv_and_disc(d, nbrs)
+    row = np.zeros(NV); row[:NG] = -disc; row[BI] = -1.0; row[BETA] = -d   # -disc.P - B - beta*d <= -pv
+    A_ub.append(row); b_ub.append(-pv)
+A_eq = np.zeros((1, NV)); A_eq[0, 0] = 1.0     # gauge P(0)=0
+c = np.zeros(NV); c[BI] = 1.0; c[BETA] = 2.0    # minimize B + 2 beta (bulk density bound)
+res = linprog(c, A_ub=np.array(A_ub), b_ub=np.array(b_ub), A_eq=A_eq, b_eq=[0.0],
+              bounds=[(-10, 10)]*NG + [(-2, 2), (-2, 2)], method='highs')
+if res.success:
+    dens = res.x[BI] + 2*res.x[BETA]
+    print(f"density bound B+2beta = {dens:.6f}   log rho* = {LOG_RHO:.6f}   gap = {dens-LOG_RHO:+.6f}")
+    print(f"  (B={res.x[BI]:.5f}, beta={res.x[BETA]:.5f})")
+    print("=> cavity potential (with handshake)", "CLOSES to log rho*" if abs(dens-LOG_RHO) < 2e-3 else "plateaus ABOVE log rho*")
+    print(f"   P grid = {np.round(res.x[:NG], 4)}")
+else:
+    print("LP failed:", res.message)