"""Elementary cores force abundance, but arbitrary RAF abundance remains open.

Run python example.py --output outputs. Set-based RAF semantics have no kinetics.
"""
# EDITABLE INPUTS: the paper's changing-witness system, Table 1.
FOOD = ('f',)
REACTIONS = {
    'a': {'reactants':('f',), 'products':('x',), 'catalysts':('y',)},
    'b': {'reactants':('f',), 'products':('y',), 'catalysts':('x',)},
    'g': {'reactants':('f',), 'products':('y',), 'catalysts':('f',)},
    'h': {'reactants':('f',), 'products':('x',), 'catalysts':('f',)},
}
CORE = ('a','b')
MODULE = ('a','b')
# The same nonnegative weight applies to every compatible row in a context.
CONTEXT_WEIGHTS = {(): '1', ('g',): '3', ('h',): '1', ('g','h'): '1'}
UNION_GROUND = ('a','b','c','dead')
UNION_MEMBERS = ((), ('a','b'), ('a','c'), ('b','c'), ('a','b','c'))
ENUMERATION_LIMIT = 262144  # explicit rejection, never a silently partial count
MANUSCRIPT_SHA256 = '5d7518980a46f1080e39ca4631539a82ba2c8d91b4362c4b778535ef5fd1fd50'

import argparse
import csv
import hashlib
import json
import platform
from collections import deque
from dataclasses import dataclass
from fractions import Fraction as Q
from itertools import combinations
from pathlib import Path


def subsets(items, limit=ENUMERATION_LIMIT):
    items=tuple(sorted(items))
    if 2**len(items)>limit:
        raise ValueError(f'Enumeration needs {2**len(items)} subsets; limit is {limit}.')
    for size in range(len(items)+1):
        for row in combinations(items,size):yield frozenset(row)


def rows_json(rows):return [sorted(s) for s in sorted(rows,key=lambda s:(len(s),tuple(sorted(s))))]


@dataclass(frozen=True)
class Reaction:
    reactants: frozenset[str]
    products: frozenset[str]
    catalysts: frozenset[str]
    def __post_init__(self):
        for key in ('reactants','products','catalysts'):
            object.__setattr__(self,key,frozenset(getattr(self,key)))


class CatalyticSystem:
    """Literal closure first; food generation and catalysis checked afterwards."""
    def __init__(self,food,reactions):
        self.food=frozenset(food)
        self.reactions={k:v if isinstance(v,Reaction) else Reaction(**v) for k,v in reactions.items()}
        self.names=frozenset(self.reactions)

    def selected(self,selected):
        selected=frozenset(selected)
        if not selected<=self.names:raise ValueError('Unknown reaction selected.')
        return selected

    def closure_stages(self,selected):
        selected=self.selected(selected);stages=[self.food]
        while True:
            current=stages[-1]
            following=current.union(*(self.reactions[r].products for r in selected
                                      if self.reactions[r].reactants<=current))
            if following==current:return tuple(stages)
            stages.append(following)

    def closure(self,selected):return self.closure_stages(selected)[-1]

    def food_generated(self,selected):
        selected=self.selected(selected);available=self.closure(selected)
        return all(self.reactions[r].reactants<=available for r in selected)

    def fixed(self,selected):
        selected=self.selected(selected);available=self.closure(selected)
        return all(self.reactions[r].reactants<=available and bool(self.reactions[r].catalysts&available)
                   for r in selected)

    def is_raf(self,selected):return bool(selected) and self.fixed(selected)

    def max_raf(self,availability=None):
        current=self.names if availability is None else self.selected(availability)
        while True:
            available=self.closure(current)
            following=frozenset(r for r in current if self.reactions[r].reactants<=available
                                 and self.reactions[r].catalysts&available)
            if following==current:return current
            current=following

    def family(self):return tuple(s for s in subsets(self.names) if self.fixed(s))

    @property
    def food_ready(self):return frozenset(r for r in self.names if self.reactions[r].reactants<=self.food)

    def product_union(self,selected):
        return frozenset().union(*(self.reactions[r].products for r in selected))

    def predecessors(self,head,module):
        return frozenset(r for r in module if self.reactions[r].products&self.reactions[head].catalysts)

    def as_json(self):
        return {'food':sorted(self.food),'reactions':{r:{k:sorted(getattr(self.reactions[r],k))
                for k in ('reactants','products','catalysts')} for r in sorted(self.names)}}


class SupportedFamily:
    """Executable horizontal/vertical cell injection from the counting proof."""
    def __init__(self,ground,bodies,upward_predicate):
        self.ground=frozenset(ground);self.bodies={r:frozenset(bodies[r]) for r in self.ground}
        if any(not b or not b<=self.ground for b in self.bodies.values()):
            raise ValueError('Every support body must be a nonempty subset of the ground.')
        self.K=upward_predicate

    def contains(self,S):
        return self.K(S) and all(S&self.bodies[r] for r in S)

    def audit_injection(self):
        cube=tuple(subsets(self.ground));truth={s:bool(self.K(s)) for s in cube}
        if any(truth[s] and not truth[s|{r}] for s in cube for r in self.ground-s):
            raise ValueError('The additional constraint is not upward.')
        bad={D for D in cube if not self.contains(self.ground-D)}
        choices={r:min(self.bodies[r]) for r in self.ground}
        records=[];targets=set()
        for D in cube:
            if D not in bad:continue
            failed={r for r in self.ground-D if self.bodies[r]<=D}
            for r in sorted(self.ground-D):
                horizontal=not truth[self.ground-D] or bool(failed-{r})
                target=(D|{r},r) if horizontal else (D,choices[r])
                assert target[0] in bad and target[1] in target[0] and target not in targets
                targets.add(target)
                records.append({'source_deleted':sorted(D),'source_column':r,
                    'move':'horizontal' if horizontal else 'vertical',
                    'target_deleted':sorted(target[0]),'target_column':target[1]})
        family=tuple(S for S in cube if self.contains(S))
        zeros=sum(len(self.ground-D) for D in bad);ones=sum(len(D) for D in bad)
        slack=sum(2*len(S)-len(self.ground) for S in family)
        assert len(records)==zeros<=ones and slack==ones-zeros>=0
        return {'family':rows_json(family),'slack':slack,'bad_zero_cells':zeros,
                'bad_one_cells':ones,'injection':records}


class CoreAnalysis:
    def __init__(self,system,core):
        self.system=system;self.core=system.selected(core)
        if not self.core<=system.food_ready or not system.is_raf(self.core):
            raise ValueError('The designated core must itself be a nonempty food-ready RAF.')

    def fibre(self,exterior):
        q=self.system;U=self.core;T=q.selected(exterior)
        if T&U:raise ValueError('Exterior context must be disjoint from the core.')
        automatic=q.food|q.product_union(T)
        bodies={r:frozenset({r}) if q.reactions[r].catalysts&automatic else q.predecessors(r,U) for r in U}
        def constraint(S):
            W=S|T
            # Exterior products are NOT usable without this food-generation gate.
            return q.food_generated(W) and all(q.reactions[t].catalysts&(q.food|q.product_union(W)) for t in T)
        supported=SupportedFamily(U,bodies,constraint)
        audit=supported.audit_injection()
        direct=tuple(S for S in subsets(U) if q.fixed(S|T))
        assert audit['family']==rows_json(direct)
        return {'exterior':sorted(T),'bodies':{r:sorted(bodies[r]) for r in sorted(U)},
                **audit,'counts':{r:sum(r in S for S in direct) for r in sorted(U)}}

    def report(self,weights=None):
        U=self.core;records=[self.fibre(T) for T in subsets(self.system.names-U)]
        family=self.system.family();N=len(family)
        counts={r:sum(r in W for W in family) for r in sorted(self.system.names)}
        assert 2*sum(counts[r] for r in U)-len(U)*N==sum(r['slack'] for r in records)
        supplied={} if weights is None else {frozenset(k):Q(v) for k,v in weights.items()}
        if any(not k<=self.system.names-U or v<0 for k,v in supplied.items()):
            raise ValueError('Context weights must be nonnegative and refer only to exterior reactions.')
        Z=Q(0);weighted={r:Q(0) for r in sorted(U)}
        for record in records:
            h=Q(1) if weights is None else supplied.get(frozenset(record['exterior']),Q(0))
            Z+=h*len(record['family'])
            for r in U:weighted[r]+=h*record['counts'][r]
        assert 2*sum(weighted.values())>=len(U)*Z
        return {'family':rows_json(family),'N_including_empty':N,'counts':counts,
                'abundant':[r for r in sorted(counts) if 2*counts[r]>=N], 'fibres':records,
                'weighted_partition':str(Z),'weighted_core_counts':{r:str(v) for r,v in weighted.items()},
                'weighted_core_frequencies':{r:str(v/Z) for r,v in weighted.items()} if Z else None,
                'weighted_abundant':[r for r,v in weighted.items() if 2*v>=Z] if Z else []}


class ModuleDiagnostics:
    def __init__(self,system,module):
        self.system=system;self.module=system.selected(module)
        if not self.module<=system.food_ready:raise ValueError('Module diagnostics require food-ready reactions.')

    def projection_queries(self):
        q=self.system;E=self.module;failures=[]
        for r in sorted(E):
            if q.reactions[r].catalysts&q.food:continue
            W=q.max_raf(q.names-q.predecessors(r,E))
            if r in W:
                assert q.is_raf(W) and not q.fixed(W&E)
                failures.append({'head':r,'witness':sorted(W),'projection':sorted(W&E)})
        return {'preserves_projection':not failures,'failures':failures}

    def defect(self):
        q=self.system;E=self.module;U=q.max_raf(E)
        if not U:return {'status':'no_viable_elementary_core'}
        family=q.family();bad=[W for W in family if not q.fixed(W&E)]
        lhs=2*sum(len(W&U) for W in family)
        rhs=len(U)*(len(family)-len(bad))+2*sum(len(W&U) for W in bad)
        assert lhs>=rhs
        return {'core':sorted(U),'bad_rows':rows_json(bad),'lhs':lhs,'rhs':rhs}


class FoodReadyGraph:
    def __init__(self,system):
        self.system=system;self.vertices=system.food_ready
        self.edges={p:frozenset(r for r in self.vertices if system.reactions[p].products&system.reactions[r].catalysts
                    or (p==r and system.food&system.reactions[r].catalysts)) for p in self.vertices}

    def shortest_cycle(self,available=None):
        available=self.vertices if available is None else frozenset(available)&self.vertices
        candidates=[]
        for start in sorted(available):
            queue=deque([(start,(start,))]);seen={start}
            while queue:
                vertex,path=queue.popleft()
                if start in self.edges[vertex]:
                    candidates.append(path);break
                for following in sorted(self.edges[vertex]&available-seen):
                    seen.add(following);queue.append((following,path+(following,)))
        return min(candidates,key=lambda c:(len(c),c)) if candidates else None

    def greedy_packing(self):
        remaining=self.vertices;cycles=[]
        while cycle:=self.shortest_cycle(remaining):
            assert self.system.is_raf(cycle)
            cycles.append(cycle);remaining=remaining-set(cycle)
        return cycles

    def abundance_candidates(self):
        """Polynomial certificate locating a witness without counting the family."""
        cycle=self.shortest_cycle()
        if cycle:return {'reason':'viable_elementary_cycle','candidates':list(cycle)}
        if not self.system.max_raf():return {'reason':'no_nonempty_raf','candidates':[]}
        exceptional=self.system.names-self.vertices
        if len(exceptional)<=2:
            return {'reason':'at_most_two_non_food_ready_reactions','candidates':sorted(exceptional)}
        return {'reason':'outside_proved_structural_cases','candidates':None}


class UnionClosedFamily:
    def __init__(self,ground,members):
        self.ground=frozenset(ground);self.members=frozenset(frozenset(s) for s in members)
        if frozenset() not in self.members or any(not s<=self.ground for s in self.members):
            raise ValueError('Members must lie in the ground set and include the empty set.')
        if any(a|b not in self.members for a in self.members for b in self.members):
            raise ValueError('The supplied family is not union-closed.')

    def interior(self,A):
        A=frozenset(A)
        if not A<=self.ground:raise ValueError('Unknown coordinate.')
        return frozenset().union(*(s for s in self.members if s<=A))

    def blockers(self,v):
        return tuple(B for B in subsets(self.ground-{v}) if v not in self.interior(self.ground-B))

    @property
    def accessible(self):
        return all(not S or any(S-{v} in self.members for v in S) for S in self.members)


class FamilyRealization:
    """All-blocker construction; marker count may be exponential."""
    def __init__(self,family):
        self.family=family
        self.blockers={v:family.blockers(v) for v in sorted(family.ground)}
        self.markers={(v,B):f'marker:{v}:{i}' for v in sorted(family.ground)
                      for i,B in enumerate(self.blockers[v])}

    def requirements(self,v):return frozenset(self.markers[v,B] for B in self.blockers[v])
    def outputs(self,v):return frozenset(name for (u,B),name in self.markers.items() if v in B)

    def paired(self):
        reactions={}
        for v in sorted(self.family.ground):
            reactions['p:'+v]=Reaction({'food'},{'alpha:'+v}|self.outputs(v),{'beta:'+v})
            reactions['g:'+v]=Reaction({'food'}|self.requirements(v),{'beta:'+v},{'alpha:'+v})
        return CatalyticSystem({'food'},reactions)

    @staticmethod
    def encode(A):return frozenset(prefix+v for v in A for prefix in ('p:','g:'))

    def complete_pairs(self,T):
        return frozenset(v for v in self.family.ground if 'p:'+v in T and 'g:'+v in T)

    def exact_ground(self,predecessors):
        if not self.family.accessible:raise ValueError('Exact-ground food feasibility must be accessible.')
        if set(predecessors)!=self.family.ground or any(not set(P)<=self.family.ground for P in predecessors.values()):
            raise ValueError('Supply a predecessor subset for every coordinate.')
        return CatalyticSystem({'food'},{v:Reaction({'food'}|self.requirements(v),
            {'private:'+v}|self.outputs(v),{'private:'+p for p in predecessors[v]}) for v in sorted(self.family.ground)})

    def audit_pairs(self):
        q=self.paired();expected={self.encode(A) for A in self.family.members}
        actual=set(q.family());assert actual==expected
        queries=[];max_rounds=0
        for T in subsets(q.names):
            wanted=self.encode(self.family.interior(self.complete_pairs(T)))
            found=q.max_raf(T);assert wanted==found
            rounds=len(q.closure_stages(T))-1;assert rounds<=2;max_rounds=max(max_rounds,rounds)
            queries.append({'available':sorted(T),'complete_pairs':sorted(self.complete_pairs(T)),
                            'max_raf':sorted(found),'closure_growth_rounds':rounds})
        singletons=frozenset(v for v in self.family.ground if frozenset({v}) in self.family.members)
        assert q.max_raf(q.food_ready)==self.encode(singletons)
        for r in q.names:
            assert len(q.reactions[r].catalysts)==1 and not q.reactions[r].products&q.reactions[r].catalysts
            suppliers=q.predecessors(r,q.names)
            assert suppliers==frozenset({'g:'+r[2:] if r.startswith('p:') else 'p:'+r[2:]})
        return {'source_family':rows_json(self.family.members),'encoded_family':rows_json(actual),
                'reaction_count':len(q.names),'marker_count':len(self.markers),
                'max_closure_growth_rounds':max_rounds,'elementary_max_raf':sorted(q.max_raf(q.food_ready)),
                'frequencies':{v:sum(v in A for A in self.family.members) for v in sorted(self.family.ground)},
                'queries':queries,'system':q.as_json()}


def worked_system(name):
    # Literal remaining manuscript examples; bracketed molecules are catalysts.
    examples={
        'feedback':{'a':(('f',),('ca',),('cb',)), 'b':(('f',),('cb','z'),('ca',)), 'g':(('z',),('cb','d'),('f',))},
        'false_average':{'a':(('f',),('z',),('f',)), 'g':(('z',),('cg',),('cl',)), 'l':(('f',),('cl',),('cg',))},
        'failed_projection':{'a':(('f',),('u',),('f',)), 'b':(('f',),('z',),('c',)), 'g':(('z',),('c',),('f',))},
    }
    return CatalyticSystem({'f'},{r:Reaction(*data) for r,data in examples[name].items()})


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)
    system=CatalyticSystem(FOOD,REACTIONS)
    report=CoreAnalysis(system,CORE).report(CONTEXT_WEIGHTS)
    diagnostics=ModuleDiagnostics(system,MODULE)
    graph=FoodReadyGraph(system);packing=graph.greedy_packing()
    rare={r for r,count in report['counts'].items() if 2*count<report['N_including_empty']}
    assert graph.shortest_cycle(rare) is None
    source=FamilyRealization(UnionClosedFamily(UNION_GROUND,UNION_MEMBERS)).audit_pairs()
    worked={}
    for name in ('feedback','false_average','failed_projection'):
        q=worked_system(name);U={'a','b'} if name=='feedback' else {'a'}
        E={'a','b'} if name!='false_average' else {'a','l'}
        worked[name]={'analysis':CoreAnalysis(q,U).report(),
                      'projection':ModuleDiagnostics(q,E).projection_queries(),
                      'defect':ModuleDiagnostics(q,E).defect()}
    result={'input_system':system.as_json(),'core_report':report,
            'module_projection':diagnostics.projection_queries(),'module_defect':diagnostics.defect(),
            'greedy_disjoint_cycles':[list(c) for c in packing],
            'abundance_candidate_certificate':graph.abundance_candidates(),
            'abundant_reaction_lower_bound':len(packing),'rare_food_ready_subgraph_acyclic':True,
            'paired_realization':source,'worked_examples':worked,
            'scope':'finite exact audits; no kinetics, no Lean rerun, no unrestricted half-frequency proof'}
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    def csv_file(name,headers,rows):
        with (out/name).open('w',newline='') as stream:
            writer=csv.writer(stream);writer.writerow(headers);writer.writerows(rows)
    csv_file('fibres.csv',['exterior','member_count',*sorted(CORE),'slack'],
             [('|'.join(r['exterior']),len(r['family']),*(r['counts'][v] for v in sorted(CORE)),r['slack']) for r in report['fibres']])
    csv_file('frequencies.csv',['reaction','count','N_including_empty','frequency'],
             [(r,n,report['N_including_empty'],str(Q(n,report['N_including_empty']))) for r,n in report['counts'].items()])
    csv_file('availability_queries.csv',['available','complete_pairs','max_raf','closure_growth_rounds'],
             [('|'.join(r['available']),'|'.join(r['complete_pairs']),'|'.join(r['max_raf']),r['closure_growth_rounds']) for r in source['queries']])
    csv_file('cell_injections.csv',['exterior','source_deleted','source_column','move','target_deleted','target_column'],
             [('|'.join(f['exterior']),'|'.join(c['source_deleted']),c['source_column'],c['move'],'|'.join(c['target_deleted']),c['target_column'])
              for f in report['fibres'] for c in f['injection']])
    lines=[f'RAF family including empty: N={report["N_including_empty"]}; frequencies={report["counts"]}.',
           f'Fibre slacks: {[f["slack"] for f in report["fibres"]]}.',
           f'Weighted core frequencies: {report["weighted_core_frequencies"]}.',
           f'Projection-preserving module: {result["module_projection"]["preserves_projection"]}.',
           f'Greedy disjoint-cycle certificate: at least {len(packing)} distinct abundant reactions.',
           f'Producer-gate benchmark: {source["reaction_count"]} reactions; {source["marker_count"]} markers; {len(source["queries"])} availability queries checked.',
           f'Elementary maxRAF in benchmark: {source["elementary_max_raf"]}.',
           'Abundance counts structural sets including the empty set; it is not a kinetic probability.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axes=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    contexts=report['fibres'];core=sorted(CORE);width=.8/max(1,len(core))
    for i,r in enumerate(core):
        axes[0].bar([j-.4+width*(i+.5) for j in range(len(contexts))],
                    [f['counts'][r]/len(f['family']) if f['family'] else float('nan') for f in contexts],width,label=r)
    axes[0].set_xticks(range(len(contexts)),['empty' if not f['exterior'] else '+'.join(f['exterior']) for f in contexts])
    axes[0].set(xlabel='Fixed exterior context',ylabel='Fraction of compatible core subsets',ylim=(0,1),title='Core-reaction frequencies by exterior\nselection')
    axes[0].legend(title='Core reaction')
    names=sorted(report['counts'])
    axes[1].bar(names,[report['counts'][r]/report['N_including_empty'] for r in names],color=['#287a99' if r in CORE else '#999999' for r in names])
    axes[1].set(xlabel='Reaction',ylabel='Fraction of all RAF-family members',ylim=(0,1),title='Reaction frequencies across the full RAF\nfamily')
    for ax in axes:ax.axhline(.5,color='black',ls='--',lw=1);ax.grid(axis='y',alpha=.2)
    fig.savefig(out/'abundance.png',dpi=180);fig.savefig(out/'abundance.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()