"""Exact first-order/neutral-second-order one-atom stability certificate.
Only standard-library arithmetic. The old six blocks are held fixed;
their masses and a seventh mass may vary in a neighborhood of (a*,0).
This does not prove global graphon optimality or multi-atom stability.
"""
from pathlib import Path
from fractions import Fraction as F
from itertools import combinations,product
import json,runpy,sys
ROOT=Path(__file__).resolve().parent;OLD=ROOT
oldargv=sys.argv[:];sys.argv=['verify_algebraic.py','--algebra-only'];z=runpy.run_path(str(OLD/'verify_algebraic.py'));sys.argv=oldargv
D=json.loads((OLD/'stationarity-input.json').read_text());data=json.loads((ROOT/'atom-certificate.json').read_text());L=D['labels'];P2=D['P2'];C4=D['C4'];ai=z['ai'];li=z['li'];Ti=1-z['gI'];weighted=[[li[k]*ai[j]/2 for j in range(6)]for k in range(3)]
edges=list(combinations(range(6),2))+[(i,i)for i in range(6)];seen=set();positive=neutral=0;worst=None
# Strict positive definiteness of the original mixed quadratic.
Hm=[[sum(li[k]*F(C4[k][i][j],4)for k in range(3))for j in range(6)]for i in range(6)]
LD=[[z['I'](0)for _ in range(6)]for _ in range(6)];piv=[]
for i in range(6):
 v=Hm[i][i]-sum(LD[i][k]*LD[i][k]*piv[k]for k in range(i));assert v.l>0;piv.append(v);LD[i][i]=z['I'](1)
 for j in range(i+1,6):LD[j][i]=(Hm[j][i]-sum(LD[j][k]*LD[i][k]*piv[k]for k in range(i)))/v

# All active flow maps use strictly positive coordinates at a*. They have
# full endpoint rank and hence remain feasible after small perturbations.
for k in range(3):
 assert all(sum(P2[k][e][j]for e,pair in enumerate(edges)if i in pair)==2*int(i==j)for i in range(6)for j in range(6))
 assert all(sum(P2[k][e][j]*ai[j]/2 for j in range(6)).l>0 for e in range(len(edges))if any(P2[k][e]))

def coeff(LL,k,partner):
 e=(partner,6);hg=[sum(LL[e[0]][e[1]]<=LL[i][j]and LL[f[0]][f[1]]<=LL[i][j]for i in set(e)for j in set(f))for f in edges]
 return [2*sum(hg[e]*P2[k][e][j]for e in range(len(edges)))-(C4[k][partner][j]if partner<6 else 0)for j in range(6)]
for row in data['profiles']:
 diag=row['diagonal'];profile=row['profile'];key=(diag,tuple(profile));assert key not in seen;seen.add(key)
 assert 1<=diag<=4 and len(profile)==6 and all(1<=profile[i]<=max(diag,L[i][i])for i in range(6))
 LL=[r+[profile[i]]for i,r in enumerate(L)]+[profile+[diag]]
 assert len(row['partners'])==3 and all(0<=v<=6 for v in row['partners'])
 cc=[coeff(LL,k,row['partners'][k])for k in range(3)]
 if 'role'not in row:
  gap=sum(cc[k][j]*weighted[k][j]for k in range(3)for j in range(6))-2*Ti
  assert gap.l>F(19,1000),(key,float(gap.l))
  worst=gap.l if worst is None else min(worst,gap.l);positive+=1;continue
 role=row['role'];assert 0<=role<6 and cc==[C4[k][role]for k in range(3)]
 # The weighted first-order price is exactly 2T by the already checked
 # mass-stationarity identity; no decimal tolerance is used here.
 ee=list(combinations(range(7),2))+[(i,i)for i in range(7)];m=len(ee)
 H=[[sum(LL[e[0]][e[1]]<=LL[i][j]and LL[f[0]][f[1]]<=LL[i][j]for i in set(e)for j in set(f))for f in ee]for e in ee]
 for k in range(3):
  pp=[[0]*6 for _ in ee]
  for j,e in enumerate(edges):pp[ee.index(e)]=P2[k][j]
  partner=row['neutral_partners'][k];assert 0<=partner<=6
  dv=[-rr[partner]if partner<6 else 0 for rr in pp];dv[ee.index((partner,6))]+=2
  Hd=[sum(H[i][j]*dv[j]for j in range(m))for i in range(m)]
  cross=[sum(Hd[e]*pp[e][j]for e in range(m))for j in range(6)];square=sum(dv[e]*Hd[e]for e in range(m))
  assert cross==C4[k][role]
  assert square-C4[k][role][role]==row['neutral_excess4'][k]>=0
 neutral+=1
expected={(diag,profile)for diag in range(1,5)for profile in product(*(range(1,max(diag,L[i][i])+1)for i in range(6)))}
assert seen==expected and len(seen)==5824 and positive==5801 and neutral==23
print('PASS:',positive,'profiles have first-order weighted matching price gap >19/1000; actual minimum lower endpoint >',float(worst))
print('PASS:',neutral,'neutral profiles have exact old-role quadratic extensions plus nonnegative square terms.')
print('PASS: one-atom local stability of the exact six-block construction (with old label cells fixed).')