"""Exact finite teaching interfaces: multilinear extension, Euler continuous
 greedy and matroid swap rounding. No claim of general matroid/submodularity
 recognition or exact ODE integration. Value/oracle promises remain obligations.
 Standard library; exact rational arithmetic; stdout only; checks survive -O.
"""
from fractions import Fraction as Q
from itertools import combinations,product
from collections import Counter,defaultdict
from math import ceil,log
import random,json

def check(ok,message):
 if not ok:raise ValueError(message)
def rational(x):
 check(type(x) in (int,Q),'exact integer or Fraction');return Q(x)
def natural(x):return type(x) is int and x>=0

def subsets(items):
 items=tuple(items)
 for k in range(len(items)+1):
  for s in combinations(items,k):yield frozenset(s)

def point(x):
 x=tuple(rational(v) for v in x);check(all(0<=v<=1 for v in x),'unit cube');return x

def value(f,s):return rational(f(frozenset(s)))

def distribution(x):
 x=point(x)
 for s in subsets(range(len(x))):
  p=Q(1)
  for i,v in enumerate(x):p*=v if i in s else 1-v
  if p:yield s,p

def extension(f,x):return sum((p*value(f,s) for s,p in distribution(x)),Q(0))

def exact_gradient(f,x):
 x=point(x);n=len(x);out=[]
 for i in range(n):
  ids=[j for j in range(n) if j!=i];g=Q(0)
  for local,p in distribution([x[j] for j in ids]):
   s=frozenset(ids[j] for j in local);g+=p*(value(f,s|{i})-value(f,s))
  out.append(g)
 return tuple(out)

def exact_hessian(f,x):
 x=point(x);n=len(x);out=[[Q(0)]*n for _ in range(n)]
 for i in range(n):
  for j in range(i+1,n):
   ids=[k for k in range(n) if k not in (i,j)];h=Q(0)
   for local,p in distribution([x[k] for k in ids]):
    s=frozenset(ids[k] for k in local)
    h+=p*(value(f,s|{i,j})-value(f,s|{i})-value(f,s|{j})+value(f,s))
   out[i][j]=out[j][i]=h
 return out

def bernoulli(p,rng):
 p=rational(p);check(0<=p<=1,'probability');return rng.randrange(p.denominator)<p.numerator

def sampled_gradient(f,x,samples,rng):
 x=point(x);check(type(samples) is int and samples>0,'positive sample count');out=[];queries=0;draws=0
 for i in range(len(x)):
  total=Q(0)
  for _ in range(samples):
   s=frozenset(j for j in range(len(x)) if j!=i and bernoulli(x[j],rng));draws+=len(x)-1
   total+=value(f,s|{i})-value(f,s);queries+=2
  out.append(total/samples)
 return tuple(out),dict(value_calls=queries,bernoulli_draws=draws)

class Matroid:
 """The supplied predicate promises the finite independence axioms. A few
 tests do not certify this promise. Query sets are immutable snapshots."""
 def __init__(self,n,independent):
  check(natural(n),'ground set size');self.n=n;self.independent=independent;self.calls=0
  check(self.test(frozenset()),'empty set independent')
  self.active=tuple(i for i in range(n) if self.test({i}));active_set=set(self.active);self.loops=tuple(i for i in range(n) if i not in active_set)
  self.rank=len(self.maximum_base([Q(0)]*n))
 def test(self,s):
  s=frozenset(s);check(all(natural(i) and i<self.n for i in s),'ground set labels');self.calls+=1
  answer=self.independent(s);check(type(answer) is bool,'boolean independence oracle');return answer
 def maximum_base(self,weights):
  check(len(weights)==self.n,'weight dimension');weights=tuple(rational(x) for x in weights);chosen=set()
  for i in sorted(self.active,key=lambda i:(-weights[i],i)):
   chosen.add(i)
   if not self.test(chosen):chosen.remove(i)
  return frozenset(chosen)
 def base(self,s):
  items=tuple(s);check(len(items)==len(set(items)),'duplicate base labels');s=frozenset(items)
  check(len(s)==self.rank and self.test(s),'valid base');return s

def coverage(sets,weights=None):
 sets=tuple(frozenset(s) for s in sets);universe=set().union(*sets) if sets else set()
 weights={u:Q(1) for u in universe} if weights is None else {u:rational(w) for u,w in weights.items()}
 check(universe<=weights.keys() and all(w>=0 for w in weights.values()),'nonnegative coverage weights')
 def f(s):return sum((weights[u] for u in set().union(*(sets[i] for i in s))),Q(0))
 return f

def euler(f,matroid,steps,samples=None,seed=0):
 check(type(steps) is int and steps>0,'positive Euler step count');check(value(f,[])==0,'normalized objective')
 n=matroid.n;r=matroid.rank;delta=Q(1,steps);x=[Q(0)]*n;decomposition=Counter();history=[];rng=random.Random(seed);sample_cost=Counter()
 # Restrict objective and random subsets to NONLOOP elements.
 active=matroid.active
 def restricted(s):return f(frozenset(active[j] for j in s))
 M=max((value(f,[i]) for i in active),default=Q(0));check(M>=0,'monotone nonnegative singleton promise')
 if r==0 or M==0:
  b=matroid.maximum_base([0]*n);x=[Q(i in b) for i in range(n)]
  return dict(point=x,decomposition=[(Q(1),b)],history=[],M=M,rank=r,boundary='rank zero' if r==0 else 'zero singletons',sample_cost={})
 for t in range(steps):
  local=[x[i] for i in active]
  if samples is None:g0=exact_gradient(restricted,local)
  else:
   g0,cost=sampled_gradient(restricted,local,samples,rng);sample_cost.update(cost)
  g=[Q(0)]*n
  for i,v in zip(active,g0):g[i]=v
  b=matroid.maximum_base(g);check(len(b)==r,'oracle rank promise');before=x.copy()
  for i in b:x[i]+=delta
  check(all(0<=v<=1 for v in x),'Euler cube invariant');decomposition[b]+=delta
  history.append(dict(round=t,gradient=g,base=sorted(b),before=before,after=x.copy()))
 check(sum(decomposition.values())==1,'base convex weights')
 return dict(point=x,decomposition=sorted(((p,b) for b,p in decomposition.items()),key=lambda pb:tuple(sorted(pb[1]))),history=history,M=M,rank=r,boundary=None,sample_cost=dict(sample_cost))

def normalize_decomposition(matroid,decomposition):
 check(type(decomposition) in (list,tuple) and decomposition,'explicit nonempty base decomposition');out=[]
 for p,b in decomposition:
  p=rational(p);check(p>0,'strictly positive base weight');out.append((p,matroid.base(b)))
 check(sum(p for p,b in out)==1,'weights sum to one');return tuple(out)

def fractional_point(n,decomposition):return [sum((p for p,b in decomposition if i in b),Q(0)) for i in range(n)]

def exchange(matroid,b1,b2):
 i=min(b1-b2)
 for j in sorted(b2-b1):
  if matroid.test((b1-{i})|{j}) and matroid.test((b2-{j})|{i}):return i,j
 raise ValueError('symmetric exchange absent: input violates matroid/base promise')

def swap_round(matroid,decomposition,seed=0,trace=False):
 terms=normalize_decomposition(matroid,decomposition);alpha,b1=terms[0];rng=random.Random(seed);records=[];operations=0
 for beta,b2 in terms[1:]:
  while b1!=b2:
   i,j=exchange(matroid,b1,b2);old1,old2=b1,b2;p=alpha/(alpha+beta)
   if bernoulli(p,rng):b2=(b2-{j})|{i};which=2
   else:b1=(b1-{i})|{j};which=1
   operations+=1
   if trace:records.append(dict(alpha=alpha,beta=beta,i=i,j=j,probability_change_second=p,changed=which,before=[sorted(old1),sorted(old2)],after=[sorted(b1),sorted(b2)]))
  alpha+=beta
 check(alpha==1 and len(b1)==matroid.rank,'final base');return dict(base=sorted(b1),exchanges=operations,trace=records)

def merge_distribution(matroid,alpha,b1,beta,b2):
 # Diagnostic exact branching, NOT the polynomial-time sampling algorithm.
 pending={(b1,b2):Q(1)};finished=defaultdict(Q)
 while pending:
  next_states=defaultdict(Q)
  for (left,right),prob in pending.items():
   if left==right:finished[left]+=prob;continue
   i,j=exchange(matroid,left,right);p=alpha/(alpha+beta)
   next_states[(left,(right-{j})|{i})]+=prob*p
   next_states[((left-{i})|{j},right)]+=prob*(1-p)
  pending=next_states
 return dict(finished)

def swap_distribution(matroid,decomposition):
 terms=normalize_decomposition(matroid,decomposition);alpha,b=terms[0];states={b:Q(1)}
 for beta,right in terms[1:]:
  new=defaultdict(Q)
  for left,prob in states.items():
   for base,q in merge_distribution(matroid,alpha,left,beta,right).items():new[base]+=prob*q
  alpha+=beta;states=dict(new)
 check(sum(states.values())==1,'distribution normalization');return states

def distribution_checks(f,matroid,terms,law):
 x=fractional_point(matroid.n,terms);marginals=[sum((p for b,p in law.items() if i in b),Q(0)) for i in range(matroid.n)]
 check(marginals==x and all(len(b)==matroid.rank and matroid.test(b) for b in law),'marginals and bases')
 mean=sum((p*value(f,b) for b,p in law.items()),Q(0));F=extension(f,x);check(mean>=F,'submodular expected value')
 products=0
 for s in subsets(range(matroid.n)):
  yes=no=Q(1)
  for i in s:yes*=x[i];no*=1-x[i]
  check(sum((p for b,p in law.items() if s<=b),Q(0))<=yes,'inclusion product')
  check(sum((p for b,p in law.items() if not(s&b)),Q(0))<=no,'exclusion product');products+=2
 return dict(point=x,marginals=marginals,extension=F,expected_value=mean,product_checks=products,law=[dict(base=sorted(b),probability=p,value=value(f,b)) for b,p in sorted(law.items(),key=lambda item:tuple(sorted(item[0])))])

def budget(rank,steps,eta=Q(0),rho=Q(0)):
 check(natural(rank) and type(steps) is int and steps>0,'budget dimensions');eta=rational(eta);rho=rational(rho);check(eta>=0 and 0<=rho<1,'budget tolerances')
 ideal=1-(1-Q(1,steps))**steps;loss=2*rank*eta+Q(rank*(rank-1),2*steps)
 return dict(discrete_growth=ideal,loss=loss,good_event_factor=max(Q(0),ideal-loss),unconditional_factor=(1-rho)*max(Q(0),ideal-loss))

def gradient_budget(n,steps,eta,rho):
 check(type(n) is int and n>0 and type(steps) is int and steps>0,'nonempty sampling budget')
 eta=rational(eta);rho=rational(rho);check(eta>0 and 0<rho<1,'sampling tolerances')
 # Integer-certified upper bound: ln(q) <= ceil(log2(q)) for q>1.
 q=Q(2*n*steps)/rho;upper=(ceil(q)-1).bit_length();K=ceil(Q(upper)/(2*eta*eta))
 return dict(samples_per_coordinate=K,log_upper_integer=upper,value_calls=2*n*steps*K,bernoulli_draws=n*(n-1)*steps*K)

def main():
 sets=[{0,1,2},{0,3,4},{0,1,2},{5}];f=coverage(sets)
 independent=lambda s:len(s&{0,1})<=1 and len(s&{2,3})<=1
 m=Matroid(4,independent);run=euler(f,m,12);terms=run['decomposition'];law=swap_distribution(m,terms);checks=distribution_checks(f,m,terms,law)
 optimum=max((value(f,b),sorted(b)) for b in subsets(range(4)) if independent(b))[0]
 greedy=set()
 while len(greedy)<2:
  feasible=[i for i in range(4) if i not in greedy and independent(greedy|{i})];i=max(feasible,key=lambda i:(value(f,greedy|{i})-value(f,greedy),-i));greedy.add(i)
 half=[(Q(1,2),frozenset({0,2})),(Q(1,2),frozenset({1,3}))];half_law=swap_distribution(m,half);half_checks=distribution_checks(f,m,half,half_law)
 half_checks['sample_original_base_mean']=sum((p*value(f,b) for p,b in half),Q(0));half_checks['independent_feasibility']=sum((p for b,p in distribution([Q(1,2)]*4) if independent(b)),Q(0))
 # Real graphic-matroid migration on square plus diagonal, no partition shortcut.
 edges=[(0,1),(1,2),(2,3),(3,0),(0,2)]
 def forest(s):
  components=[{i} for i in range(4)]
  for e in sorted(s):
   u,v=edges[e];a=next(i for i,c in enumerate(components) if u in c);b=next(i for i,c in enumerate(components) if v in c)
   if a==b:return False
   merged=components[a]|components[b];components=[c for i,c in enumerate(components) if i not in(a,b)]+[merged]
  return True
 graph=Matroid(5,forest);gf=coverage([{0,1},{1,2},{2,3},{0,3},{0,2}]);gt=[(Q(1,3),{0,1,2}),(Q(2,3),{0,3,4})];graph_checks=distribution_checks(gf,graph,gt,swap_distribution(graph,gt))
 # Loop with huge singleton is absent from every feasible solution.
 loop_f=coverage(sets+[{100}],{**{i:Q(1) for i in range(6)},100:Q(10**6)})
 lm=Matroid(5,lambda s:4 not in s and independent(s));loop_run=euler(loop_f,lm,12);check(loop_run['M']==3 and loop_run['point'][4]==0,'remove loops before M')
 sampled=euler(f,Matroid(4,independent),12,samples=64,seed=3)
 boundary=dict(empty=euler(lambda s:Q(0),Matroid(0,lambda s:True),3),all_loops=euler(lambda s:Q(len(s)),Matroid(2,lambda s:not s),3),zero=euler(lambda s:Q(0),Matroid(2,lambda s:len(s)<=1),3))
 rejected=[]
 for label,fn in [('negative probability',lambda:extension(f,[-1,0,0,0])),('zero steps',lambda:euler(f,m,0)),('empty decomposition',lambda:swap_round(m,[])),('wrong weight sum',lambda:swap_round(m,[(Q(1,2),{0,2})])),('dependent base',lambda:swap_round(m,[(Q(1),{0,1})])),('zero weight',lambda:swap_round(m,[(Q(0),{0,2}),(Q(1),{1,3})]))]:
  try:fn()
  except ValueError:rejected.append(label)
  else:raise RuntimeError('missing rejection '+label)
 check(checks['expected_value']>=budget(2,12)['good_event_factor']*optimum,'finite-step theorem')
 return dict(status='PASS',coverage=[sorted(s) for s in sets],optimum=optimum,ordinary_greedy=dict(base=sorted(greedy),value=value(f,greedy)),euler=run,euler_rounding=checks,finite_step_budget=budget(2,12),half_basis_rounding=half_checks,one_random_path=swap_round(m,half,seed=7,trace=True),graphic_migration=graph_checks,loop=dict(removed=lm.loops,M=loop_run['M'],point=loop_run['point']),certified_sample_budget=gradient_budget(4,12,Q(1,20),Q(1,10)),sampled_demo=dict(point=sampled['point'],decomposition=sampled['decomposition'],cost=sampled['sample_cost'],note='64 samples demonstrate execution; no claimed uniform accuracy event'),boundaries=boundary,rejections=rejected)

def encode(x):
 if isinstance(x,Q):return str(x)
 if isinstance(x,(set,frozenset)):return sorted(x)
 raise TypeError(type(x).__name__)
if __name__=='__main__':print(json.dumps(main(),default=encode,indent=2))