"""Focused reproducible checks of the Math page's precise statements.

This is not a formal proof checker. Exact enumeration, symbolic identities,
floating-point examples and counterexamples are labeled separately. Original
manuscripts and author scripts remain unchanged in ../sources.
Run: python research/math/validation/recheck.py
"""
from pathlib import Path
from fractions import Fraction as F
from itertools import product, combinations
import importlib.util
import hashlib
import json
import math
import platform
import sys
import numpy as np
import scipy
from scipy.linalg import expm
import sympy as sp

HERE = Path(__file__).resolve().parent
records = []

def record(id, title, kind, observations, counterexample=False):
    records.append(dict(id=id, title=title, kind=kind,
        outcome='counterexample_confirmed' if counterexample else 'check_satisfied',
        observations=observations))

# A telescoping product bound, with exact arithmetic on a small grid.
vectors=list(product([F(0),F(1,2),F(1)], repeat=4))
def budget(e): return 1-math.prod(1-x for x in e)
for a,b in product(vectors,repeat=2):
    assert abs(budget(a)-budget(b)) <= sum(abs(x-y) for x,y in zip(a,b))
record('R1','Component error-budget transfer','exact rational enumeration',
       {'pairs':len(vectors)**2,'bound':'|B(e)-B(f)| <= sum_i |e_i-f_i| for B(e)=1-product_i(1-e_i)'})

# Small uniform score error is not sufficient for decision-probability stability.
threshold_examples=[]
for eps in [F(1,10),F(1,100),F(1,10000)]:
    threshold_examples.append({'epsilon':str(eps),'score_difference':str(2*eps),
        'first_alert_probability':0,'second_alert_probability':1})
record('R2','Threshold mass defeats unconditional first-order decision transfer',
       'exact counterexample',{'construction':'All mass at scores -epsilon and +epsilon; threshold 0.',
       'examples':threshold_examples},True)

# Levi-Civita metric compatibility is not a bound on extrinsic bending.
arc=math.pi;chord=2.
assert arc>chord
record('R3','A constant induced metric does not make arc and chord equal',
       'analytic example with numerical values',{'embedding':'f(t)=(cos t,sin t)',
       'induced_metric':1,'covariant_derivative_of_metric':0,'jacobian_norm':1,
       'semicircle_arc':arc,'endpoint_chord':chord,'relative_distortion':1-chord/arc},True)

# Population Cantelli bound with priors: exact discrete distributions, not fitted estimates.
cantelli_cases=0
weight_sets=[(1,2,1),(2,1,3),(3,3,1)]
for weights0,weights1,shift,prior in product(weight_sets,weight_sets,[2,3,5],[F(1,100),F(1,2),F(9,10)]):
    values0=[F(-1),F(0),F(1)];values1=[x+shift for x in values0]
    probs0=[F(w,sum(weights0)) for w in weights0];probs1=[F(w,sum(weights1)) for w in weights1]
    mu0=sum(x*p for x,p in zip(values0,probs0));mu1=sum(x*p for x,p in zip(values1,probs1))
    var0=sum(p*(x-mu0)**2 for x,p in zip(values0,probs0));var1=sum(p*(x-mu1)**2 for x,p in zip(values1,probs1))
    threshold=(mu0+mu1)/2;a=(mu1-mu0)/2
    c0=var0/(var0+a*a);c1=var1/(var1+a*a)
    fp=sum(p for x,p in zip(values0,probs0) if x>threshold)
    tp=sum(p for x,p in zip(values1,probs1) if x>threshold)
    assert fp<=c0 and 1-tp<=c1
    lower=prior*(1-c1)/((1-prior)*c0+prior*(1-c1))
    actual=prior*tp/((1-prior)*fp+prior*tp)
    assert lower<=actual
    cantelli_cases+=1
record('R4','Population Cantelli posterior bound including class priors',
       'exact finite-distribution check',{'cases':cantelli_cases,'scope':'Population moments only; no plug-in confidence claim.'})

# Symbolic algebra and independently constructed spin matrices.
t=sp.symbols('t',real=True)
assert sp.trigsimp(sp.cos(t/2)**2+sp.sin(t/2)**2-1)==0
assert sp.trigsimp((1-sp.cos(t))**2+sp.sin(t)**2-4*sp.sin(t/2)**2)==0
spin_error=0.
for twice_j in range(1,6):
    j=twice_j/2;ms=np.arange(-j,j+1);raising=np.zeros((len(ms),len(ms)),complex)
    for i,m in enumerate(ms[:-1]):raising[i+1,i]=np.sqrt((j-m)*(j+m+1))
    jy=(raising-raising.T.conj())/(2j)
    for theta in np.linspace(0,2*np.pi,41):
        overlap=abs(expm(1j*theta*jy)[-1,-1])**2
        spin_error=max(spin_error,abs(overlap-np.cos(theta/2)**(2*twice_j)))
assert spin_error<1e-12
record('D1','Chord identity and transverse highest-weight spin overlaps',
       'symbolic identities and matrix-exponential checks',{'spin_cases':205,'max_absolute_error':spin_error,
       'scope':'d=|sin(theta/2)|, K=2d; specified transverse rotation, not arbitrary dynamics.'})
assert sp.cos(0)**2==1 and sp.cos(sp.pi/2)**2==0
record('D2','Half-angle overlap does not descend to the arrow-axis quotient',
       'exact counterexample',{'q(0)':'[1:0]','q(pi)':'[-1:0]=[1:0]',
       'S(0)':1,'S(pi)':0,'scope':'Keep the quotient angle and the Bloch/rotation angle distinct.'},True)

# Finite aromatic-ring spectrum: exact formula precedes its large-N limit.
ring_error=0.;deviations={}
for n in [6,10,14,50]:
    deviations[n]=0.
    for theta in np.linspace(0,np.pi,21):
        H=np.zeros((n,n),complex)
        for i in range(n):H[i,(i+1)%n]=-np.exp(1j*theta/n);H[(i+1)%n,i]=H[i,(i+1)%n].conjugate()
        ev=np.linalg.eigvalsh(H)
        gap=ev[n//2]-ev[n//2-1]
        ratio=gap/(4*np.sin(np.pi/n))
        formula=np.sin((np.pi-theta)/n)/np.sin(np.pi/n)
        ring_error=max(ring_error,abs(ratio-formula))
        deviations[n]=max(deviations[n],abs(formula-(1-theta/np.pi)))
assert ring_error<1e-12
assert all(deviations[a]>deviations[b] for a,b in zip([6,10,14],[10,14,50]))
record('D3','Finite aromatic-ring gap and the linear large-ring limit',
       'Hermitian eigensolver versus analytic formula',{'cases':84,'max_absolute_error':ring_error,
       'formula':'sin((pi-theta)/N)/sin(pi/N), N=4k+2, 0<=theta<=pi',
       'max_deviations_from_linear_limit':deviations})
record('D4','Commutativity alone does not force linear gap closure',
       'exact counterexample',{'family':'diag(0, 1+theta^2)','gap':'1+theta^2',
       'commutator':'0 for every pair of parameter values'},True)

holonomy=np.exp(2j*np.pi*.5)
assert abs(holonomy+1)<1e-12
record('C1','Flatness need not imply trivial global holonomy',
       'analytic example with numerical value',{'base':'circle','connection':'A=(1/2) dtheta in a U(1) trivialization',
       'curvature':0,'loop_holonomy_real':float(holonomy.real),'loop_holonomy_imag':float(holonomy.imag),
       'scope':'The manuscript main theorem adds simple connectivity; its earlier unqualified lemma needs that restriction.'},True)
record('C2','No uniform derivative L2 bound from a connection L2 norm',
       'analytic counterexample family',{'base':'flat torus [0,2pi]^2','connection':'A_n=sin(nx)dy/n',
       'curvature':'cos(nx) dx wedge dy','norm_A_squared':'2*pi^2/n^2','YM_energy':'pi^2',
       'ratios':{n:2/(n*n) for n in [1,2,4,8,16,32]},
       'scope':'Refutes the unrestricted compact-base L2 lemma; the simply-connected main claim also needs a defined, justified norm comparison.'},True)

# Interpolation is exact at nodes even where the Hessian is nonzero.
record('A1','A positive universal pointwise lower bound fails at interpolation nodes',
       'exact counterexample',{'local_function':'f(x)=x^2','linear_interpolant_on_0_h':'h*x',
       'hessian':2,'error_at_0':0,'error_at_h':0,'maximum_error':'h^2/4',
       'compact_manifold_example':'f(t)=1-cos(t) on S1 also has error 0 and Hessian 1 at an interpolated node t=0.'},True)
interpolation=[]
for h in [1.,.5,.25,.125]:
    x=np.linspace(0,h,1001);error=float(np.max(abs(x*x-h*x)))
    assert abs(error-h*h/4)<1e-14
    interpolation.append({'h':h,'maximum_error':error,'sup_second_derivative':2.,'bound':2*h*h/8})
record('A2','A concrete uniform-norm interpolation bound with quadratic scaling',
       'numerical check of exact polynomial formula',{'examples':interpolation,
       'statement':'For C2 f on [a,b], ||f-I_h f||_infinity <= h^2 sup|f_second|/8.'})

# Lift identities over exact finite joint distributions, including zero cells.
tables=0
for a,b,c,d in product(range(6),repeat=4):
    total=a+b+c+d
    if not total:continue
    pa=F(a+b,total);pb=F(a+c,total);pab=F(a,total)
    if not pa or not pb:continue
    lift=pab/(pa*pb)
    assert pab/pb==lift*pa and pab/pa==lift*pb
    assert (lift==1)==(pab==pa*pb)
    assert lift==1+(pab-pa*pb)/(pa*pb)
    tables+=1
record('P1','Lift, conditioning, symmetry and independence',
       'exact rational enumeration',{'joint_tables':tables,'cell_counts':'0..5; exclude undefined marginal denominators',
       'scope':'Ordinary probability identities; no Riemannian curvature is inferred.'})

# An independent restricted-growth enumeration versus the source bitmask DP.
spec=importlib.util.spec_from_file_location('source_partition_counter',HERE.parent/'sources/test2b_fixed_n_topology.py')
source=importlib.util.module_from_spec(spec);spec.loader.exec_module(source)
import networkx as nx
def partitions(n):
    def extend(labels):
        if len(labels)==n:
            yield tuple(sum(1<<i for i,x in enumerate(labels) if x==k) for k in range(max(labels)+1));return
        for value in range(max(labels)+2):yield from extend(labels+[value])
    yield from extend([0])
def connected(mask,adj):
    seen=mask&-mask;front=seen
    while front:
        reach=0
        for i in range(len(adj)):
            if front&(1<<i):reach|=adj[i]
        front=reach&mask&~seen;seen|=front
    return seen==mask
graph_count=edge_additions=connected_graphs=0
family_table=[]
for n in range(1,6):
    edges=list(combinations(range(n),2));parts=list(partitions(n));counts={}
    for emask in range(1<<len(edges)):
        adj=[0]*n;G=nx.Graph();G.add_nodes_from(range(n))
        for k,(a,b) in enumerate(edges):
            if emask&(1<<k):adj[a]|=1<<b;adj[b]|=1<<a;G.add_edge(a,b)
        count=sum(all(connected(block,adj) for block in part) for part in parts)
        assert count==source.count_connected_partitions(G)
        counts[emask]=count;graph_count+=1
        if nx.is_connected(G):
            connected_graphs+=1;assert 2**(n-1)<=count<=len(parts)
    for emask,count in counts.items():
        for k in range(len(edges)):
            if not emask&(1<<k):
                assert count<=counts[emask|(1<<k)];edge_additions+=1
for n in range(1,11):
    tree=source.count_connected_partitions(nx.path_graph(n));complete=source.count_connected_partitions(nx.complete_graph(n))
    assert tree==2**(n-1) and complete==int(sp.bell(n))
    cycle=source.count_connected_partitions(nx.cycle_graph(n)) if n>=3 else None
    if n>=3:assert cycle==2**n-n
    family_table.append({'n':n,'tree':tree,'cycle':cycle,'complete':complete})
record('Z1','Connected partitions: independent enumeration and edge-addition monotonicity',
       'exact exhaustive enumeration',{'all_labeled_graphs_n_1_to_5':graph_count,
       'single_edge_additions':edge_additions,'connected_graphs_with_sandwich_bound':connected_graphs,
       'family_counts':family_table,'scope':'Finite checks supplement the injection and counting proofs; this is the established connected-set-partition object.'})

output={'date':'2026-10-07','scope':'Focused claim checks, not formal verification or blanket manuscript validation.',
        'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__,'sympy':sp.__version__,
        'script_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),'checks':records}
(HERE/'results.json').write_text(json.dumps(output,indent=2,ensure_ascii=False)+'\n',encoding='utf-8')
print(json.dumps({'checks_completed':len(records),'counterexamples_confirmed':sum(x['outcome']=='counterexample_confirmed' for x in records),
                  'graphs_checked':graph_count,'edge_additions_checked':edge_additions},indent=2))