"""Exact two-neutral-atom obstruction table; stdlib, rational intervals.
The accompanying proof gives local two-atom stability, not global optimality.
"""
from pathlib import Path
from itertools import combinations
from fractions import Fraction as F
import runpy,sys,json
ROOT=Path(__file__).resolve().parent;OLD=ROOT
sys.argv=['verify_algebraic.py','--algebra-only'];z=runpy.run_path(str(OLD/'verify_algebraic.py'));D=json.loads((OLD/'stationarity-input.json').read_text());L=D['labels'];P2=D['P2'];C4=D['C4'];ai=z['ai'];li=z['li'];ee=list(combinations(range(6),2))+[(i,i)for i in range(6)]
rows=[r for r in json.loads((ROOT/'atom-certificate.json').read_text())['profiles']if 'role'in r]
# The cumulative birth-profile derivative D[k,i,label], multiplied by 2.
for row in rows:
 DD=[]
 for k in range(3):
  partner=row['neutral_partners'][k];dp=[[0]*5 for _ in range(7)]
  if partner<6:
   for e,flowrow in zip(ee,P2[k]):
    for i in set(e):
     for t in range(L[e[0]][e[1]],5):dp[i][t]-=flowrow[partner]
  birth=row['diagonal']if partner==6 else row['profile'][partner]
  for t in range(birth,5):
   dp[6][t]+=2
   if partner<6:dp[partner][t]+=2
  DD.append(dp)
 row['DD']=DD
count=negative=zero=positive=0;mingap=None;output=[]
for ii,r in enumerate(rows):
 for jj in range(ii,len(rows)):
  s=rows[jj]
  for label in range(1,max(r['diagonal'],s['diagonal'])+1):
   count+=1;LL=[L[i]+[r['profile'][i],s['profile'][i]]for i in range(6)]+[r['profile']+[r['diagonal'],label],s['profile']+[label,s['diagonal']]]
   coeff=[]
   for k in range(3):
    dr=r['DD'][k]+[[0]*5];ds=s['DD'][k][:6]+[[0]*5,s['DD'][k][6]]
    cross4=sum(dr[i][LL[i][j]]*ds[j][LL[i][j]]for i in range(8)for j in range(8));coeff.append(cross4-C4[k][r['role']][s['role']])
   kval=sum(li[k]*F(coeff[k],4)for k in range(3))
   if all(v==0 for v in coeff):zero+=1;continue
   if kval.l>=0:positive+=1;continue
   assert kval.h<0
   # Direct matching between new atoms, replacing their old-partner
   # extension on equal min(mass1,mass2) portions, has positive
   # first derivative in at least one old matching strategy.
   candidates=[]
   for k in range(3):
    hg=[sum(label<=LL[i][j]and L[e[0]][e[1]]<=LL[i][j]for i in(6,7)for j in set(e))for e in ee]
    cc=[sum(hg[e]*P2[k][e][j]for e in range(len(ee)))-F(C4[k][r['role']][j]+C4[k][s['role']][j],2)for j in range(6)]
    gain=sum(cc[j]*ai[j]for j in range(6))
    if gain.l>0:candidates.append((gain.l*li[k].l,k))
   assert candidates,(ii,jj,label,coeff)
   margin,k=max(candidates);assert margin>F(19,500);negative+=1;mingap=margin if mingap is None else min(mingap,margin)
   output.append(dict(i=ii,j=jj,label=label,negative_coefficient=coeff,improving_strategy=k))
assert count==979
print('PASS: neutral two-atom interaction table:',count,'cases;',negative,'negative cross terms all have a positive first-order direct-matching response;',zero,'zero;',positive,'positive.')
print('Every required weighted first-order response is >19/500; actual lower endpoint minimum >',float(mingap))
(ROOT/'neutral-pair-exact-results.json').write_text(json.dumps(dict(cases=count,negative=negative,zero=zero,positive=positive,strict_weighted_margin='19/500',responses=output),indent=2))