"""Exact rational certificates for integer-polyhedral teaching examples.
Python 3 standard library only. All LP optimization is exhaustive vertex enumeration.
"""
from fractions import Fraction as F
from itertools import combinations,product,permutations
import pathlib,json,sys,math
OUT=pathlib.Path(sys.argv[1]) if len(sys.argv)>1 else pathlib.Path.cwd()/'foundation-integer-capstone-results.json'
def dot(a,b):return sum((x*y for x,y in zip(a,b)),F(0))
def solve(A,b):
 n=len(b);M=[[F(x) for x in row]+[F(v)] for row,v in zip(A,b)]
 for j in range(n):
  p=next((i for i in range(j,n) if M[i][j]),None)
  if p is None:return None
  M[j],M[p]=M[p],M[j];v=M[j][j];M[j]=[x/v for x in M[j]]
  for i in range(n):
   if i!=j and M[i][j]:
    v=M[i][j];M[i]=[x-v*y for x,y in zip(M[i],M[j])]
 return tuple(M[i][-1] for i in range(n))
def determinant(A):
 n=len(A)
 if n==0:return F(1)
 M=[[F(x) for x in r] for r in A];d=F(1)
 for j in range(n):
  p=next((i for i in range(j,n) if M[i][j]),None)
  if p is None:return F(0)
  if p!=j:M[j],M[p]=M[p],M[j];d=-d
  v=M[j][j];d*=v
  for i in range(j+1,n):
   z=M[i][j]/v
   for k in range(j+1,n):M[i][k]-=z*M[j][k]
 return d

def vertices(A,b):
 n=len(A[0]);out={}
 for I in combinations(range(len(A)),n):
  x=solve([A[i] for i in I],[b[i] for i in I])
  if x is not None and all(dot(a,x)<=v for a,v in zip(A,b)):out[x]=I
 return sorted(out)
def lp(A,b,c):
 V=vertices(A,b)
 if not V:return {'status':'infeasible','vertices':[]} # callers explicitly use bounded polyhedra
 z=max(dot(c,x) for x in V);x=min(x for x in V if dot(c,x)==z);tight=[i for i,(a,v) in enumerate(zip(A,b)) if dot(a,x)==v];dual=None
 for I in combinations(tight,len(c)):
  y=solve([[A[i][j] for i in I] for j in range(len(c))],c)
  if y is not None and all(v>=0 for v in y):
   dual=[F(0)]*len(A)
   for i,v in zip(I,y):dual[i]=v
   break
 assert dual is not None
 assert all(sum(dual[i]*A[i][j] for i in range(len(A)))==c[j] for j in range(len(c))) and dot(dual,b)==z
 return {'status':'optimal','point':x,'value':z,'dual':dual,'vertex_count':len(V),'tight_rows':tight}

def binary_model():
 n=5;A=[];b=[];names=[]
 for i in range(n):
  v=[0]*n;v[i]=-1;A.append(v);b.append(0);names.append(f'-x{i}<=0')
 for i in range(n):
  v=[0]*n;v[i]=1;v[(i+1)%n]=1;A.append(v);b.append(1);names.append(f'edge{i}')
 # edge inequalities plus nonnegativity imply every xi <=1.
 A.append([-2,0,2,0,0]);b.append(1);names.append('capacity')
 return A,b,names

def perfect_matching(S,support=None):
 # Hopcroft-Karp on positive support; deterministic increasing neighbor order.
 from collections import deque
 n=len(S);adj=([sorted(row) for row in support] if support is not None else [[j for j in range(n) if S[i][j]>0] for i in range(n)]);pu=[None]*n;pv=[None]*n;INF=10**9
 while True:
  d=[0 if pu[u] is None else INF for u in range(n)];Q=deque(u for u in range(n) if pu[u] is None);short=INF
  while Q:
   u=Q.popleft()
   if d[u]>=short:continue
   for v in adj[u]:
    if pv[v] is None:short=d[u]+1
    elif d[pv[v]]==INF:d[pv[v]]=d[u]+1;Q.append(pv[v])
  if short==INF:break
  def dfs(u):
   for v in adj[u]:
    w=pv[v]
    if (w is None and d[u]+1==short) or (w is not None and d[w]==d[u]+1 and dfs(w)):
     pu[u]=v;pv[v]=u;return True
   d[u]=INF;return False
  for u in range(n):
   if pu[u] is None:dfs(u)
 return tuple(pu) if all(v is not None for v in pu) else None

def clean(v):
 if isinstance(v,F):return v.numerator if v.denominator==1 else str(v)
 if isinstance(v,dict):return {str(k):clean(x) for k,x in v.items()}
 if isinstance(v,(list,tuple)):return [clean(x) for x in v]
 return v

def run():
 R={}
 A=[[-1,0],[0,-1],[2,1],[1,2]];b=[0,0,3,3];V=vertices(A,b);Z=[x for x in product(range(4),repeat=2) if all(dot(a,x)<=v for a,v in zip(A,b))];assert set(Z)=={(0,0),(0,1),(1,0),(1,1)}
 R['integer_hull']={'A':A,'b':b,'vertices':V,'integer_points':Z,'objective_x_LP':lp(A,b,[1,0]),'integer_optimum':1}
 # Directed incidence of a four-cycle with a diagonal.
 arcs=[(0,1),(1,2),(2,3),(3,0),(0,2)];B=[[(-1 if u==v else 1 if w==v else 0) for u,w in arcs] for v in range(4)];mins=[]
 for k in range(1,5):
  for rows in combinations(range(4),k):
   for cols in combinations(range(5),k):
    d=determinant([[B[i][j] for j in cols] for i in rows]);assert d in(-1,0,1);mins.append(d)
 triangle=[[1,0,1],[1,1,0],[0,1,1]];assert determinant(triangle)==2
 R['TU']={'incidence':B,'minors_checked':len(mins),'minor_values':sorted(set(mins)),'triangle':triangle,'triangle_determinant':2,'unimodular_not_TU':[[1,2],[0,1]]}
 # Birkhoff: exact residual mass, deterministic Hopcroft-Karp supported matching.
 perms=[(0,1,2,3),(1,2,3,0),(2,3,0,1)];weights=[F(1,2),F(1,3),F(1,6)];M=[[sum(w for p,w in zip(perms,weights) if p[i]==j) for j in range(4)] for i in range(4)];S=[r[:] for r in M];mass=F(1);steps=[];support=[{j for j in range(4) if S[i][j]>0} for i in range(4)]
 while mass:
  p=perfect_matching(S,support);assert p is not None
  alpha=min(S[i][p[i]] for i in range(4));before=[r[:] for r in S]
  for i in range(4):
   S[i][p[i]]-=alpha
   if S[i][p[i]]==0:support[i].remove(p[i])
  mass-=alpha;assert all(sum(row)==mass for row in S) and all(sum(S[i][j] for i in range(4))==mass for j in range(4))
  steps.append({'permutation':p,'weight':alpha,'before':before,'residual': [r[:] for r in S],'remaining_mass':mass})
 assert sum(x['weight'] for x in steps)==1
 assert all(M[i][j]==sum(x['weight'] for x in steps if x['permutation'][i]==j) for i in range(4) for j in range(4))
 R['birkhoff']={'matrix':M,'steps':steps,'support_size':sum(x>0 for row in M for x in row),'matching_implementation':'Hopcroft-Karp'}
 for bits in product((0,1),repeat=9):
  S=[bits[i*3:(i+1)*3] for i in range(3)];p=perfect_matching(S);exists=any(all(S[i][p[i]] for i in range(3)) for p in permutations(range(3)));assert (p is not None)==exists
  if p is not None:assert len(set(p))==3 and all(S[i][p[i]] for i in range(3))
 R['matching_support_crosscheck']={'all_3_by_3_boolean_supports':512}
 # Degree-only and full odd-set systems are checked on small exact polytopes.
 checked=[]
 for n,edges,odd in [(4,[(0,1),(1,2),(2,3),(3,0)],False),(3,[(0,1),(1,2),(2,0)],False),(3,[(0,1),(1,2),(2,0)],True),(5,[(0,1),(1,2),(2,3),(3,4),(4,0)],True)]:
  E=len(edges);AA=[[-int(i==j) for j in range(E)] for i in range(E)];bb=[0]*E
  for v in range(n):AA.append([int(v in e) for e in edges]);bb.append(1)
  if odd:
   for k in range(3,n+1,2):
    for U in combinations(range(n),k):AA.append([int(a in U and b in U) for a,b in edges]);bb.append((k-1)//2)
  VV=vertices(AA,bb);frac=[v for v in VV if any(x.denominator!=1 for x in v)]
  if odd or n==4:assert not frac
  else:assert frac==[(F(1,2),)*3]
  checked.append({'n':n,'edges':edges,'odd_sets':odd,'vertices':VV,'fractional_vertices':frac})
 R['matching_polytopes']=checked
 weights=[[5,4],[4,1]];prices=[4,3,1,0];assert all(prices[i]+prices[2+j]>=weights[i][j] for i in range(2) for j in range(2)) and sum(prices)==weights[0][1]+weights[1][0]==8
 R['TDI_example']={'edge_weights':weights,'vertex_prices':prices,'matching':[[0,1],[1,0]],'value':8,'scaled_halfline_dual_for_c1':F(1,2)}
 prism_edges=[(0,1),(1,2),(2,0),(3,4),(4,5),(5,3),(0,3),(1,4),(2,5)]
 matchings=[[(0,3),(1,2),(4,5)],[(1,4),(0,2),(3,5)],[(2,5),(0,1),(3,4)]]
 for M in matchings:assert len({v for e in M for v in e})==6
 for e in prism_edges:assert sum(F(1,3) for M in matchings if any(set(e)==set(f) for f in M))==F(1,3)
 R['prism_tight_odd_cut']={'edges':prism_edges,'edge_values':[F(1,3)]*9,'decomposition':matchings,'weights':[F(1,3)]*3,'odd_set':[0,1,2],'cut_value':1}
 assert determinant([[1,2],[0,1]])==1
 # Pure integer tableau row and a negative coefficient boundary.
 solutions=[]
 for x1,x2 in product(range(7),repeat=2):
  xb=F(3,2)-F(1,2)*x1-F(1,4)*x2
  if xb>=0 and xb.denominator==1:assert F(1,2)*x1+F(1,4)*x2>=F(1,2);solutions.append((xb,x1,x2))
 assert F(3,2)-F(3,2)*F(1,3)==1 and F(1,2)*F(1,3)<F(1,2)
 R['gomory']={'continuous_counterexample':{'z':F(1,3),'xB':1,'wrong_cut_lhs':F(1,6),'wrong_cut_rhs':F(1,2)},'basic_rhs':F(3,2),'row_coefficients':[F(1,2),F(1,4)],'cut_coefficients':[F(1,2),F(1,4)],'cut_rhs':F(1,2),'nonnegative_integer_solutions':solutions,'negative_fractional_part':F(-1,4)-math.floor(F(-1,4))}
 # Rank-4 example: finite checks support an all-integer-normal proof in prose.
 A0=[[-1,0],[1,0],[0,-1],[-4,1],[4,1]];b0=[0,1,0,0,4];normals=[(a,b) for a in range(-12,13) for b in range(-12,13) if(a,b)!=(0,0)];lower=[]
 for h in [F(2),F(3,2),F(1),F(1,2)]:
  tri=[(F(0),F(0)),(F(1),F(0)),(F(1,2),h)];q=(F(1,2),h-F(1,2))
  for c in normals:assert dot(c,q)<=math.floor(max(dot(c,p) for p in tri))
  lower.append({'h':h,'next_apex':q,'normals_checked':len(normals)})
 A=[r[:] for r in A0];b=b0[:];rounds=[]
 for cuts in [[([1,1],2),([-1,1],1)],[([0,1],1)],[([1,1],1),([-1,1],0)],[([0,1],0)]]:
  old=vertices(A,b);rows=[]
  for c,rhs in cuts:
   support=max(dot(c,x) for x in old);assert math.floor(support)==rhs;rows.append({'normal':c,'support':support,'rounded_rhs':rhs})
  for c,rhs in cuts:A.append(c);b.append(rhs)
  rounds.append({'cuts':rows,'outer_approximation_vertices':vertices(A,b)})
 assert vertices(A,b)==[(F(0),F(0)),(F(1),F(0))]
 R['CG_rank4']={'original_vertices':vertices(A0,b0),'round_certificates':rounds,'lower_witness_checks':lower,'rank':4,'warning':'These finite cut outer approximations are not asserted equal to entire intermediate closures. Exact rank lower bound uses the all-integer-normal triangle lemma in prose.'}
 # Branch and bound example, bounded two-variable LPs.
 A=[[-1,0],[0,-1],[1,0],[0,1],[3,2]];b=[0,0,2,3,7];c=[8,5]
 nodes=[]
 for name,extra,rhs in [('root',[],[]),('y<=0',[[0,1]],[0]),('y>=1',[[0,-1]],[-1]),('y>=1,x<=1',[[0,-1],[1,0]],[-1,1]),('y>=1,x>=2',[[0,-1],[-1,0]],[-1,-2])]:
  ans=lp(A+extra,b+rhs,c);nodes.append({'name':name,'A':A+extra,'b':b+rhs,**ans})
 assert nodes[0]['value']==F(37,2) and nodes[2]['value']==F(55,3) and nodes[3]['value']==18 and nodes[4]['status']=='infeasible'
 # Infeasibility certificate: 3x+2y<=7, -3x<=-6, -2y<=-2 gives 0<=-1.
 cert=[0,0,0,0,1,2,3];bad=nodes[-1];assert all(sum(t*a[j] for t,a in zip(cert,bad['A']))==0 for j in range(2)) and dot(cert,bad['b'])==-1
 R['branch_bound']={'nodes':nodes,'optimum':18,'incumbent':[1,2],'infeasible_certificate':{'multipliers':[0,0,0,0,1,2,3],'rhs':-1}}
 A,b,names=binary_model();c=[5,5,8,7,7];cases=[]
 specs=[('root',[],[],[]),('global_odd',[[1]*5],[2],['odd']),('left_before_local',[[1]*5,[1,0,0,0,0]],[2,0],['odd','x0<=0']),('left_after_local',[[1]*5,[1,0,0,0,0],[0,0,1,0,0]],[2,0,0],['odd','x0<=0','local_x2<=0']),('right',[[1]*5,[-1,0,0,0,0]],[2,-1],['odd','x0>=1']),('wrong_global_local',[[1]*5,[0,0,1,0,0]],[2,0],['odd','WRONG_x2<=0'])]
 for name,extra,rhs,labels in specs:
  ans=lp(A+extra,b+rhs,c);cases.append({'name':name,'constraint_names':names+labels,'A':A+extra,'b':b+rhs,**ans});print(name,ans['status'],clean(ans.get('value')),flush=True)
 assert [x['value'] for x in cases]==[16,F(57,4),F(27,2),12,13,12]
 feasible=[x for x in product((0,1),repeat=5) if all(dot(a,x)<=v for a,v in zip(A,b))];z=max(dot(c,x) for x in feasible);best=[x for x in feasible if dot(c,x)==z];assert z==13 and best==[(1,0,1,0,0)]
 R['branch_cut']={'objective':c,'nodes':cases,'integer_feasible':feasible,'optimum':z,'optimal_points':best,'local_cut_scope':'x0=0 descendants only','global_odd_cut_multipliers':[0]*5+[F(1,2)]*5+[0]}
 OUT.parent.mkdir(parents=True,exist_ok=True);OUT.write_text(json.dumps(clean(R),ensure_ascii=False,indent=2)+'\n');print('All exact certificates passed:',OUT)
if __name__=='__main__':run()