"""Reproduce the September 28, 2026 LSM geometric foundation checks.

Run in the site checkout: python scripts/check-foundations.py
Standalone: save this file as scripts/check-foundations.py, and save
/lsm-foundations-inputs.json as dist/lsm-foundations-inputs.json under the
same parent directory. Install numpy and scipy, then run the same command.
FCC incidences use exact doubled integer coordinates, not rounded meshes.
The independent HCP/Voronoi cross-check uses NumPy/SciPy and tolerance 1e-8.
The universal electron minimum follows from the proof on the site, not a
finite enumeration. These checks do not simulate particle dynamics.
"""
from itertools import product, combinations, permutations
from collections import Counter
from pathlib import Path
from fractions import Fraction
import json, math, hashlib

ROOT = Path(__file__).resolve().parents[1]

def cell(kind, p):
    p = tuple(p)  # coordinates in doubled units; hosts have even coordinate values
    if kind == 'O':
        hosts = [tuple(x + (s*2 if j == a else 0) for j, x in enumerate(p))
                 for a in range(3) for s in (-1, 1)]
    else:
        hosts = [q for q in product(*[(x-1, x+1) for x in p]) if (sum(q)//2) % 2 == 1]
    return {'type': kind, 'center': p, 'hosts': frozenset(hosts)}

def edges(cells):
    return [(i,j) for i,j in combinations(range(len(cells)),2)
            if len(cells[i]['hosts'] & cells[j]['hosts']) == 3]

def audit(cells):
    ts=sum(c['type']=='T' for c in cells); os=len(cells)-ts; es=edges(cells)
    reached={0}
    while True:
        nxt=reached | {b for a,b in es if a in reached} | {a for a,b in es if b in reached}
        if nxt==reached:break
        reached=nxt
    return {'T':ts,'O':os,'chargeNumerator':2*os-ts,'chargeDenominator':6,
            'faceContacts':len(es),'exposedApertures':4*ts+8*os-2*len(es),
            'connected':len(reached)==len(cells)}

def hcp_cross_check():
    import numpy as np
    from scipy.spatial import Delaunay, Voronoi, ConvexHull
    pts=np.array([(i+j/2+(k%2)/2,np.sqrt(3)*j/2+(k%2)*np.sqrt(3)/6,k*np.sqrt(2/3))
                  for i,j,k in product(range(-3,4),range(-3,4),range(-4,5))])
    groups={}
    for ids in Delaunay(pts).simplices:
        p=pts[ids]; a=2*(p[1:]-p[0]); b=(p[1:]**2).sum(1)-(p[0]**2).sum()
        if abs(np.linalg.det(a))<1e-10:continue
        c=np.linalg.solve(a,b)
        if sum((p[0]-c)**2)>.51:continue
        groups.setdefault(tuple(c.round(8)),set()).update(map(int,ids))
    result={}
    for typ in [4,6]:
        _,hosts=min([(k,v) for k,v in groups.items() if len(v)==typ],key=lambda kv:np.linalg.norm(kv[0]))
        ns=Counter(len(v) for v in groups.values() if len(v & hosts)==3)
        result['T' if typ==4 else 'O']={'T':ns[4],'O':ns[6]}
    assert result=={'T':{'T':1,'O':3},'O':{'T':6,'O':2}}
    vor=Voronoi(pts); origin=int(np.argmin((pts**2).sum(1)))
    verts=vor.vertices[vor.regions[vor.point_region[origin]]]
    faces=[ids for pair,ids in zip(vor.ridge_points,vor.ridge_vertices) if origin in pair]
    volume=float(ConvexHull(verts).volume)
    assert len(verts)==14 and len(faces)==12 and all(len(f)==4 for f in faces)
    assert abs(volume-1/math.sqrt(2))<1e-10
    return {'neighborCounts':result,'voronoi':{'vertices':14,'edges':24,'faces':12,
            'volumeOverD3':volume},'method':'Finite patch Delaunay cells grouped by circumcenter; numerical cross-check.'}

def run():
    electron=[cell('O',(0,0,0))]+[cell('T',p) for p in product((-1,1),repeat=3)]
    ea=audit(electron)
    assert ea==dict(T=8,O=1,chargeNumerator=-6,chargeDenominator=6,faceContacts=8,exposedApertures=24,connected=True)
    # Every nearest T shares exactly three hosts with the central O. Check all 48 signed permutations.
    sig={(c['type'],c['center'],c['hosts']) for c in electron}
    for perm in permutations(range(3)):
        for signs in product((-1,1),repeat=3):
            f=lambda p:tuple(signs[i]*p[perm[i]] for i in range(3))
            assert {(c['type'],f(c['center']),frozenset(f(h) for h in c['hosts'])) for c in electron}==sig
    # A full first shell around host (2,0,0), followed by one outside O.
    host=(2,0,0)
    cage=[cell('T',(host[0]+x,y,z)) for x,y,z in product((-1,1),repeat=3)]
    cage += [cell('O',tuple(x+(2*s if j==axis else 0) for j,x in enumerate(host)))
             for axis in range(3) for s in (-1,1)]
    assert all(host in c['hosts'] for c in cage)
    assert audit(cage)['faceContacts']==24
    smaller=cage+[cell('O',(4,2,2))]; sa=audit(smaller)
    assert sa==dict(T=8,O=7,chargeNumerator=6,chargeDenominator=6,faceContacts=25,exposedApertures=38,connected=True)
    omega_t=3*math.acos(1/3)-math.pi; omega_o=4*math.asin(1/3)
    assert abs(8*omega_t+6*omega_o-4*math.pi)<1e-12
    aperture=math.sqrt(3)-math.pi/2; st=4*omega_t; so=6*omega_o
    vt=4/(3*math.sqrt(2))-4*omega_t/3; vo=8*math.sqrt(2)/3-2*omega_o
    assert min(aperture,st,so,vt,vo)>0
    # Composite Monad: exact incidence inventory plus the analytic volume identity.
    assert 8*Fraction(1,4)==2 and 6*Fraction(1,6)==1
    assert 2*Fraction(-1,6)+Fraction(1,3)==0
    sphere_volume=4*math.pi/3
    assert abs(sphere_volume+2*vt+vo-4*math.sqrt(2))<1e-12
    for o in range(1,5):
        assert 2*o+o >= 3  # nonempty neutral mixed inventory T=2O

    area=lambda a:a['T']*st+a['O']*so+a['exposedApertures']*aperture
    source=ROOT/'scripts/public-particle-layouts.json'
    if source.exists():
        raw=json.loads(source.read_text())
        layouts={k:{'sites':raw[k]['sites']} for k in ['electron','positron','proton']}
        (ROOT/'dist/lsm-foundations-inputs.json').write_text(json.dumps(layouts,indent=2)+'\n')
    else:
        layouts=json.loads((ROOT/'dist/lsm-foundations-inputs.json').read_text())
    published={}
    for key in ['electron','positron','proton']:
        cs=[cell(s['type'],tuple(round(v*2) for v in s['pos'])) for s in layouts[key]['sites']]
        a=audit(cs);assert a['connected'];a['cavityAreaOverR2']=area(a);published[key]=a
    assert published['electron']['faceContacts']==8
    assert published['proton']['faceContacts']==32
    assert area(sa)<published['proton']['cavityAreaOverR2']
    # Six mutually separated T cells have the same total charge and smaller area;
    # this is outside the connected theorem and exposes the need for confinement.
    separated={'T':6,'O':0,'exposedApertures':24}
    assert area(separated)<area(ea)
    report={'edition':'2026-09-28','status':'Geometric checks; conditional surface theorem; no physical stability result',
            'monadAccounting':{'sphere':1,'T':2,'O':1,'geometricTypes':3,'pieces':4,'charge':0,'hostChargeAssumed':0,'volumeOverR3':4*math.sqrt(2),'exactVolumeOverR3':'4*sqrt(2)','hcpPrimitiveSphereCount':2},
            'electron':ea,'electronSymmetryOrder':48,'cagePlusOctah':sa,
            'cagePlusOctahCoordinates':[{k:list(v) if k=='center' else v for k,v in c.items() if k!='hosts'} for c in smaller],
            'published':published,'unitHostRadius':{'apertureArea':aperture,'tetrahWallArea':st,
            'octahWallArea':so,'tetrahVoidVolume':vt,'octahVoidVolume':vo,'electronCavityArea':area(ea),
            'cagePlusOctahCavityArea':area(sa),'sixSeparatedTetrahCavityArea':area(separated)},
            'hcp':hcp_cross_check(),
            'proofScope':'All finite nonempty face-connected unions of whole FCC interstitial cavities on a fixed undeformed host lattice, q/e=(2O-T)/6=-1. No dynamical or fragmentation conclusion.',
            'sourceSha256':{p:hashlib.sha256((ROOT/p).read_bytes()).hexdigest() for p in ['scripts/check-foundations.py','dist/lsm-foundations-inputs.json']}}
    (ROOT/'dist/lsm-foundations-checks.json').write_text(json.dumps(report,indent=2)+'\n')
    (ROOT/'dist/lsm-foundations-checker.txt').write_text(Path(__file__).read_text())
    print(json.dumps({'electron':ea,'cagePlusOctah':sa,'hcp':report['hcp'],'passed':True}))

if __name__=='__main__':run()