"""Recognize support families, check certificates, and eliminate hidden vertices."""
# EDITABLE INPUTS: members are subsets of exactly these labelled vertices.
VERTICES = ('a','b','c')
FAMILY = ((),('a',),('a','b'),('a','b','c'))
RESIDUAL_POLICY = 'largest_reverse'  # largest cardinality, then largest bit mask
NODE_BUDGET = 100000
POWERSET_BUDGET = 4096
MANUSCRIPT_SHA256 = '74a88c8c0b5389c8eb0f7ff73655f82ba75f6f79246594ae36629706e225fce0'

import argparse
import csv
from dataclasses import dataclass
import hashlib
import itertools
import json
import math
from pathlib import Path
import platform


def contained(a,b): return a & ~b == 0
def bits(mask,n): return [i for i in range(n) if mask>>i&1]


@dataclass(frozen=True)
class SetFamily:
    labels: tuple
    members: frozenset

    def __post_init__(self):
        if len(set(self.labels))!=len(self.labels): raise ValueError('Ground labels must be distinct.')
        if any(not isinstance(s,int) or s<0 or s>self.full for s in self.members): raise ValueError('Set outside the ground set.')

    @property
    def n(self): return len(self.labels)
    @property
    def full(self): return (1<<self.n)-1

    @classmethod
    def from_subsets(cls,labels,members):
        index={v:i for i,v in enumerate(labels)}
        return cls(tuple(labels),frozenset(sum(1<<index[v] for v in set(s)) for s in members))

    def names(self,s): return [self.labels[i] for i in bits(s,self.n)]

    def interior(self,s):
        result=0
        for t in self.members:
            if contained(t,s):result|=t
        return result

    def maximal_residuals(self,U):
        residual=U-self.members
        return sorted(s for s in residual if not any(s!=t and contained(s,t) for t in residual))

    def eligible(self,s,used=0): return s & ~(self.interior(s)|used)

    def powerset(self,budget=POWERSET_BUDGET):
        if 1<<self.n>budget: raise ValueError('Truth-table budget exceeded; use graph certificate verification for explicit lists.')
        return frozenset(range(1<<self.n))


@dataclass(frozen=True)
class HornTheory:
    n: int
    rules: tuple  # (body mask, head index); multiple heads allowed here for general use

    def __post_init__(self):
        if self.n<0 or any(b<0 or b>=1<<self.n or h<0 or h>=self.n for b,h in self.rules): raise ValueError('Invalid Horn rule.')

    def models(self,s): return all(not contained(b,s) or s>>h&1 for b,h in self.rules)

    def closure(self,s):
        if s<0 or s>=1<<self.n: raise ValueError('Seed outside ground set.')
        while True:
            old=s
            for b,h in self.rules:
                if contained(b,s):s|=1<<h
            if old==s:return s

    def next_closure(self):
        """Ganter lectic enumeration, without constructing the whole power set."""
        current=self.closure(0)
        while True:
            yield current
            for i in reversed(range(self.n)):
                if current>>i&1:continue
                lower=(1<<i)-1; candidate=self.closure((current&lower)|(1<<i))
                if candidate&lower==current&lower:
                    current=candidate;break
            else:return

    def normalize_completion(self,accepted,residual,budget=POWERSET_BUDGET):
        """Check Lemma 4.1 premises and replace a completion rule by its G-closure."""
        if accepted.n!=self.n or not set(accepted.rules)<=set(self.rules): raise ValueError('Accepted rules must be a subtheory.')
        if len({h for b,h in self.rules})!=len(self.rules): raise ValueError('Completion must be single-head.')
        if 1<<self.n>budget:raise ValueError('Premise audit exceeds truth-table budget.')
        if not accepted.models(residual) or self.models(residual):raise ValueError('Not a residual model.')
        if any(s!=residual and contained(s,residual) and accepted.models(s) and not self.models(s) for s in range(1<<self.n)):
            raise ValueError('Residual is not inclusion-minimal.')
        for body,head in self.rules:
            if contained(body,residual) and not residual>>head&1:
                if accepted.closure(body)!=residual:raise AssertionError('Normalization invariant failed.')
                replacement=HornTheory(self.n,tuple((residual,h) if h==head else (b,h) for b,h in self.rules))
                if list(self.next_closure())!=list(replacement.next_closure()):raise AssertionError('Models changed.')
                return replacement,{'old_body':body,'new_body':residual,'head':head,'accepted_closure_of_old_body':residual}
        raise AssertionError('Missing violated rule.')


@dataclass(frozen=True)
class SupportGraph:
    labels: tuple
    predecessors: tuple

    def __post_init__(self):
        if len(self.labels)!=len(self.predecessors) or len(set(self.labels))!=len(self.labels):raise ValueError('One row per distinct vertex required.')
        if any(p<0 or p>=1<<len(self.labels) for p in self.predecessors):raise ValueError('Predecessor outside ground set.')

    def supports(self,s):
        if s<0 or s>=1<<len(self.labels):raise ValueError('Set outside ground set.')
        return all(s&p for r,p in enumerate(self.predecessors) if s>>r&1)

    def family(self,budget=POWERSET_BUDGET):
        n=len(self.labels)
        if 1<<n>budget:raise ValueError('Family materialization exceeds budget.')
        return SetFamily(self.labels,frozenset(s for s in range(1<<n) if self.supports(s)))

    def interior(self,s):
        if s<0 or s>=1<<len(self.labels):raise ValueError('Set outside ground set.')
        while True:
            reduced=sum(1<<r for r,p in enumerate(self.predecessors) if s>>r&1 and s&p)
            if reduced==s:return s
            s=reduced

    def horn(self):return HornTheory(len(self.labels),tuple((p,r) for r,p in enumerate(self.predecessors)))

    def verify_explicit_family(self,family):
        """Polynomial in n + listed-family size; stop at first extra supported set."""
        if family.labels!=self.labels:raise ValueError('Ground set mismatch.')
        for s in sorted(family.members):
            if not self.supports(s):return {'equal':False,'reason':'listed_set_unsupported','witness':s,'models_checked':0}
        checked=0
        for model in self.horn().next_closure():
            checked+=1; supported=family.full^model
            if supported not in family.members:
                return {'equal':False,'reason':'unlisted_supported_set','witness':supported,'models_checked':checked}
        return {'equal':True,'models_checked':checked}

    def eliminate(self,label):
        """Existential projection of one vertex; different from induced deletion."""
        v=self.labels.index(label); remaining=[r for r in range(len(self.labels)) if r!=v]; rows=[]
        for r in remaining:
            p=self.predecessors[r]
            if p>>v&1:
                p=(1<<r) if self.predecessors[v]>>v&1 else (p&~(1<<v))|self.predecessors[v]
            rows.append(sum(1<<j for j,old in enumerate(remaining) if p>>old&1))
        return SupportGraph(tuple(self.labels[r] for r in remaining),tuple(rows))

    def project(self,visible):
        visible=frozenset(visible)
        if not visible<=set(self.labels):raise ValueError('Unknown visible vertex.')
        graph=self
        for label in self.labels:
            if label not in visible:graph=graph.eliminate(label)
        return graph

    def delete(self,visible):
        visible=frozenset(visible)
        if not visible<=set(self.labels):raise ValueError('Unknown visible vertex.')
        kept=[i for i,s in enumerate(self.labels) if s in visible]
        return SupportGraph(tuple(self.labels[i] for i in kept),tuple(sum(1<<j for j,k in enumerate(kept) if self.predecessors[i]>>k&1) for i in kept))

    def independent_product(self,other):
        if set(self.labels)&set(other.labels):raise ValueError('Disjoint product requires disjoint ground labels.')
        n=len(self.labels)
        return SupportGraph(self.labels+other.labels,self.predecessors+tuple(p<<n for p in other.predecessors))


def select_residual(maximal,policy):
    if policy not in ('largest','largest_reverse'):raise ValueError('Unknown residual policy.')
    return min(maximal,key=lambda s:(-s.bit_count(),s if policy=='largest' else -s))


class IntervalRecognizer:
    def __init__(self,family,policy=RESIDUAL_POLICY,node_budget=NODE_BUDGET,powerset_budget=POWERSET_BUDGET):
        self.family=family;self.policy=policy;self.node_budget=node_budget;self.initial=family.powerset(powerset_budget)
        select_residual([0],policy)
        if node_budget<1:raise ValueError('Positive node budget required.')

    def solve(self):
        family=self.family;nodes=[];memo={}
        def visit(U,used):
            key=(U,used)
            if key in memo:return memo[key]
            if len(nodes)>=self.node_budget:
                return 'unresolved',None,[]
            node_id=len(nodes);node={'U':sorted(U),'used':used};nodes.append(node)
            if U==family.members:
                node.update(status='success');answer=('success',node_id,[]);memo[key]=answer;return answer
            s=select_residual(family.maximal_residuals(U),self.policy);roots=bits(family.eligible(s,used),family.n)
            node.update(residual=s,eligible=roots,children=[])
            for root in roots:
                remaining=frozenset(x for x in U if not (x>>root&1 and contained(x,s)))
                status,child,trace=visit(remaining,used|1<<root)
                node['children'].append({'root':root,'node':child})
                if status=='success':
                    node['status']='success';answer=('success',node_id,[(s,root),*trace]);memo[key]=answer;return answer
                if status=='unresolved':
                    node['status']='unresolved';answer=('unresolved',node_id,[]);memo[key]=answer;return answer
            node['status']='rejected';answer=('rejected',node_id,[]);memo[key]=answer;return answer
        status,root,trace=visit(self.initial,0)
        b=sum((1<<r) not in family.members for r in range(family.n))
        result={'labels':list(family.labels),'target_family':sorted(family.members),
            'status':status,'policy':self.policy,'trace':[list(step) for step in trace],'unique_states':len(nodes),
            'unmemoized_node_bound':sum(math.factorial(b)//math.factorial(b-j) for j in range(b+1)),
            'excluded_singletons':bits(sum(1<<r for r in range(family.n) if 1<<r not in family.members),family.n),
            'root':root,'nodes':nodes}
        if status=='success':
            graph=check_trace(family,trace)
            result['predecessors']=list(graph.predecessors);result['explicit_list_verification']=graph.verify_explicit_family(family)
        elif status=='rejected':check_rejection(family,result)
        return result


def check_trace(family,trace,budget=POWERSET_BUDGET):
    U=family.powerset(budget);used=0;rows=[1<<r for r in range(family.n)]
    for s,r in trace:
        if s not in family.maximal_residuals(U) or r<0 or r>=family.n or not family.eligible(s,used)>>r&1:
            raise ValueError('Illegal maximal residual or root.')
        U=frozenset(x for x in U if not (x>>r&1 and contained(x,s)));used|=1<<r;rows[r]=family.full^s
    if U!=family.members:raise ValueError('Trace does not reach target family.')
    graph=SupportGraph(family.labels,tuple(rows))
    if not graph.verify_explicit_family(family)['equal']:raise AssertionError('Reconstruction failed.')
    return graph


def check_rejection(family,certificate,budget=POWERSET_BUDGET):
    """Replay every branch of a rejection DAG; absent children are not rejection."""
    if certificate['status']!='rejected':raise ValueError('No completed rejection certificate.')
    nodes=certificate['nodes'];checked=set();active=set()
    def visit(i,U,used):
        if not isinstance(i,int) or i<0 or i>=len(nodes):raise ValueError('Missing rejection child.')
        node=nodes[i]
        if node['U']!=sorted(U) or node['used']!=used or node['status']!='rejected':raise ValueError('State mismatch.')
        if i in active:raise ValueError('Cyclic proof.')
        if i in checked:return
        if U==family.members:raise ValueError('Rejected a successful leaf.')
        active.add(i);s=select_residual(family.maximal_residuals(U),certificate['policy'])
        roots=bits(family.eligible(s,used),family.n)
        if node['residual']!=s or node['eligible']!=roots or [c['root'] for c in node['children']]!=roots:
            raise ValueError('Incomplete root branching.')
        for child in node['children']:
            r=child['root'];remaining=frozenset(x for x in U if not (x>>r&1 and contained(x,s)))
            visit(child['node'],remaining,used|1<<r)
        active.remove(i);checked.add(i)
    visit(certificate['root'],family.powerset(budget),0)
    return {'valid':True,'checked_states':len(checked)}


def hall_matching(family,U=None,used=0):
    """Necessary fresh-root matching, never treated as a sufficient recognizer."""
    U=family.powerset() if U is None else U;left=family.maximal_residuals(U)
    adjacency={s:bits(family.eligible(s,used),family.n) for s in left};owner={}
    def augment(s,seen):
        for r in adjacency[s]:
            if r in seen:continue
            seen.add(r)
            if r not in owner or augment(owner[r],seen):owner[r]=s;return True
        return False
    for s in left:augment(s,set())
    complete=len(owner)==len(left)
    reached=set(left)-set(owner.values());rights=set();pending=list(reached)
    while pending:
        s=pending.pop()
        for r in adjacency[s]:
            rights.add(r)
            if r in owner and owner[r] not in reached:reached.add(owner[r]);pending.append(owner[r])
    return {'matching_exists':complete,'matching':{str(s):r for r,s in sorted(owner.items())},
        'deficient_residual_sets':sorted(reached) if not complete else [],'their_eligible_roots':sorted(rights) if not complete else []}


@dataclass(frozen=True)
class StructuralCRS:
    """Small presence/absence model, with OR catalysis and genuine reactant conjunction."""
    labels: tuple
    food: frozenset
    reactants: tuple
    products: tuple
    catalysts: tuple

    def __post_init__(self):
        if len(set(self.labels))!=len(self.labels) or any(len(rows)!=len(self.labels) for rows in (self.reactants,self.products,self.catalysts)):
            raise ValueError('One reactant/product/catalyst set per distinct reaction required.')

    def fixed_family(self,budget=POWERSET_BUDGET):
        if 1<<len(self.labels)>budget:raise ValueError('CRS fixed-family materialization exceeds budget.')
        members={0}
        for mask in range(1,1<<len(self.labels)):
            selected=bits(mask,len(self.labels));pool=set(self.food)
            while True:
                before=set(pool)
                for i in selected:
                    if self.reactants[i]<=pool:pool.update(self.products[i])
                if before==pool:break
            if all(self.reactants[i]<=pool and self.catalysts[i]&pool for i in selected):members.add(mask)
        return SetFamily(self.labels,frozenset(members))

    @classmethod
    def elementary(cls,graph):
        n=len(graph.labels)
        return cls(graph.labels,frozenset(('food',)),tuple(frozenset(('food',)) for _ in range(n)),
            tuple(frozenset((f'product_{i}',)) for i in range(n)),
            tuple(frozenset(f'product_{j}' for j in bits(p,n)) for p in graph.predecessors))


def hall_counterexample():
    K={63}|{s for s in range(64) if (s&7).bit_count()<=1 and (s&56).bit_count()<=1}
    return SetFamily(tuple('abcdef'),frozenset(63^s for s in K))


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));args=parser.parse_args()
    out=args.output;out.mkdir(parents=True,exist_ok=True)
    custom=SetFamily.from_subsets(VERTICES,FAMILY);examples={'editable':custom,'cycle':SetFamily(tuple('abc'),frozenset((0,7))),
        'conjunction':SetFamily(tuple('abc'),frozenset((0,1,2,3,7))),'hall_counterexample':hall_counterexample()}
    results={}
    for name,family in examples.items():
        result=IntervalRecognizer(family).solve();(out/f'{name}_certificate.json').write_text(json.dumps(result,indent=2)+'\n')
        results[name]={k:v for k,v in result.items() if k!='nodes'};results[name]['hall']=hall_matching(family)
    chain=SupportGraph(tuple('abc'),(1,1,2));chain_f=chain.family()
    bad_prefix=[(6,2)]; U=chain_f.powerset();U=frozenset(s for s in U if not (s&4 and contained(s,6)))
    next_s=select_residual(chain_f.maximal_residuals(U),'largest_reverse')
    results['bad_root']={'family':sorted(chain_f.members),'legal_prefix':bad_prefix,'next_residual':next_s,'eligible_unused':bits(chain_f.eligible(next_s,4),3)}
    accepted=HornTheory(3,((1,1),)); completion=HornTheory(3,((1,1),(1,2)))
    normalized,receipt=completion.normalize_completion(accepted,3);results['normalization']=receipt
    cycle=SupportGraph(('a','h'),(2,1));projected=cycle.project(('a',));deleted=cycle.delete(('a',))
    results['projection_vs_deletion']={'original_family':sorted(cycle.family().members),'projected_family':sorted(projected.family().members),
        'deleted_family':sorted(deleted.family().members),'projected_row':projected.predecessors[0]}
    general=StructuralCRS(tuple('abc'),frozenset(('f',)),(frozenset(('f',)),frozenset(('f',)),frozenset(('x','y'))),
        (frozenset(('x',)),frozenset(('y',)),frozenset(('z',))),(frozenset(('f',)),)*3)
    results['general_crs_conjunction_family']=sorted(general.fixed_family().members)
    graph_families={SupportGraph(tuple('abc'),p).family().members for p in itertools.product(range(8),repeat=3)}
    all_families=[frozenset(s for s in range(8) if mask>>s&1) for mask in range(256)]
    union_closed=[F for F in all_families if 0 in F and all(a|b in F for a in F for b in F)]
    results['three_vertex_census']={'graphs':512,'realizable_families':len(graph_families),'union_closed_with_empty':len(union_closed)}
    (out/'results.json').write_text(json.dumps(results,indent=2)+'\n')
    with (out/'trace_steps.csv').open('w',newline='') as f:
        writer=csv.writer(f);writer.writerow(['example','step','residual','root','predecessor_row'])
        for name,family in examples.items():
            for i,(s,r) in enumerate(results[name]['trace']):writer.writerow([name,i+1,','.join(family.names(s)),family.labels[r],','.join(family.names(family.full^s))])
    lines=[f'Editable family: {results["editable"]["status"]}; trace {results["editable"]["trace"]}.',
        f'Conjunction family: {results["conjunction"]["status"]}; completed rejection DAGs are independently checked.',
        f'Six-vertex Hall example: initial matching {results["hall_counterexample"]["hall"]["matching_exists"]}, recognizer {results["hall_counterexample"]["status"]}, {results["hall_counterexample"]["unique_states"]} memoized states.',
        'Projection of a <-> h gives a supported visible singleton; deletion of h does not.',
        f'Three-vertex census: {len(graph_families)} graph-realizable families among {len(union_closed)} union-closed empty-containing families.',
        'Budget exhaustion is unresolved. These finite certificates do not rerun Lean or prove polynomial-time recognition.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    import numpy as np
    display=[chain_f,examples['cycle'],examples['conjunction']];names=['Chain: realizable','Cycle: realizable','Conjunction: rejected']
    fig,ax=plt.subplots(figsize=(10,3.8),layout='constrained');matrix=np.array([[s in family.members for s in range(8)] for family in display])
    ax.imshow(matrix,cmap=matplotlib.colors.ListedColormap(['#efefeb','#187f91']),vmin=0,vmax=1,aspect='auto')
    ax.set_yticks(range(3),names);ax.set_xticks(range(8),['empty','a','b','ab','c','ac','bc','abc']);ax.set(xlabel='Subset of the same three labelled vertices',title='Subset families and directed-graph\nrealizability')
    ax.set_xticks(np.arange(-.5,8),minor=True);ax.set_yticks(np.arange(-.5,3),minor=True);ax.grid(which='minor',color='white',linewidth=2);ax.tick_params(which='minor',bottom=False,left=False)
    fig.savefig(out/'families.png',dpi=180);fig.savefig(out/'families.svg');plt.close(fig)
    fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    for ax,steps,title in [(axes[0],[(6,1),(5,2)],'Successful roots b, then c'),(axes[1],[(6,2)],'Root c first: a dead end')]:
        current=set(range(8));hist=[len(current)];labels=['Start']
        for s,r in steps:
            current={x for x in current if not (x>>r&1 and contained(x,s))};hist.append(len(current));labels.append('Delete at '+chain_f.labels[r])
        ax.plot(range(len(hist)),hist,'o-',color='#187f91');ax.axhline(len(chain_f.members),ls='--',color='#b95024',label='Target family size')
        ax.set_xticks(range(len(hist)),labels);ax.set(ylim=(3.5,8.5),ylabel='Sets remaining in U',title=title);ax.legend(fontsize=8);ax.grid(alpha=.2)
    fig.savefig(out/'root_choices.png',dpi=180);fig.savefig(out/'root_choices.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    (out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),
        'python':platform.python_version(),'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')


if __name__=='__main__':main()