"""Selective interventions, catalytic seeds and ambient closedness on literal data.

Run: python example.py --output outputs. Presence/absence structure, not kinetics.
"""
# EDITABLE INPUTS. The source medium has 68 entries; additions change that medium.
TARGET = 'C00183'                 # L-valine
PRESERVATION_TARGET = 'C00073'     # L-methionine
DELETED_ACTIONS = ('R01209','R01210','R04441')  # removes all directions of each row
EXTRA_FOOD = ()
RESCUE_FOOD = 'C00141'             # demonstrated precursor rescue, not any barrier member
SMALL_DIRECTIONS = ('R00014::f','R04672::r','R10916::r','R10985::f','R00585::f',
                    'R01215::r','R04441::f','R00221::r','R00589::r','R_NADs_2::f')
POOLING_DIRECTION = 'R_NADs_1::f'
MAX_SUBSETS = 65536                # exact enumeration is only for small ambient sets
MANUSCRIPT_SHA256 = 'd271244c721f70339c8e26094d51d691ef1915e556e742e3dc0648a64c0aa50c'
PINNED_SOURCE_SHA256 = 'f6834863b2d272d074af5befcb6149c11edf5a64b1aa0ba716ebeb8f10b3700b'

import argparse
import csv
from dataclasses import dataclass, replace
import hashlib
import itertools
import json
from pathlib import Path
import platform
import re


@dataclass(frozen=True)
class CatalystFormula:
    """Positive formula; conjunction is never flattened into alternative catalysts."""
    kind: str
    children: tuple = ()
    atom: str = ''

    def __post_init__(self):
        if self.kind not in ('true','atom','and','or'): raise ValueError('Positive formulas only.')
        if self.kind=='atom' and (not self.atom or self.children): raise ValueError('Malformed atom.')
        if self.kind in ('and','or') and len(self.children)!=2: raise ValueError('Binary connective required.')
        if self.kind=='true' and (self.children or self.atom): raise ValueError('Malformed true constant.')
        if self.kind in ('and','or') and self.atom: raise ValueError('Connectives have no atom label.')

    def holds(self,pool):
        if self.kind=='true': return True
        if self.kind=='atom': return self.atom in pool
        if self.kind=='and': return all(c.holds(pool) for c in self.children)
        return any(c.holds(pool) for c in self.children)

    def mandatory_atoms(self):
        if self.kind=='true': return frozenset()
        if self.kind=='atom': return frozenset((self.atom,))
        a,b=(c.mandatory_atoms() for c in self.children)
        return a|b if self.kind=='and' else a&b

    def to_data(self):
        return [self.kind,self.atom] if self.kind=='atom' else [self.kind,*[c.to_data() for c in self.children]]

    @classmethod
    def parse(cls,text):
        text=text.strip().replace('|',',').replace('*','&')
        text=re.sub(r'\s*([(),&])\s*',r'\1',text)
        text=re.sub(r'\s+',',',text)
        if not text: return cls('true')
        tokens=re.findall(r'[(),&]|[^(),&]+',text); pos=0
        def term():
            nonlocal pos
            if pos>=len(tokens): raise ValueError('Incomplete catalyst formula.')
            token=tokens[pos]; pos+=1
            if token=='(':
                result=disjunction()
                if pos>=len(tokens) or tokens[pos]!=')': raise ValueError('Missing closing parenthesis.')
                pos+=1; return result
            if token in ',&)': raise ValueError('Expected catalyst atom.')
            return cls('atom',atom=token)
        def conjunction():
            nonlocal pos
            result=term()
            while pos<len(tokens) and tokens[pos]=='&':
                pos+=1; result=cls('and',(result,term()))
            return result
        def disjunction():
            nonlocal pos
            result=conjunction()
            while pos<len(tokens) and tokens[pos]==',':
                pos+=1; result=cls('or',(result,conjunction()))
            return result
        result=disjunction()
        if pos!=len(tokens): raise ValueError('Unconsumed catalyst tokens.')
        return result


@dataclass(frozen=True)
class Reaction:
    id: str
    action: str
    reactants: tuple             # (species, integer coefficient), retained for reuse
    products: tuple
    catalyst: CatalystFormula
    source_line: int = 0

    def __post_init__(self):
        if not self.id or not self.action: raise ValueError('Direction and action identifiers required.')
        if any(not s or not isinstance(n,int) or n<=0 for s,n in self.reactants+self.products):
            raise ValueError('Positive integer stoichiometric coefficients required.')

    @property
    def inputs(self): return frozenset(s for s,n in self.reactants)
    @property
    def outputs(self): return frozenset(s for s,n in self.products)


@dataclass(frozen=True)
class ReactionSystem:
    reactions: tuple
    food: frozenset
    source_sha256: str = ''

    def __post_init__(self):
        if len({r.id for r in self.reactions})!=len(self.reactions): raise ValueError('Duplicate directed identifiers.')

    @property
    def by_id(self): return {r.id:r for r in self.reactions}

    @classmethod
    def read(cls,path,expected_hash=None):
        data=Path(path).read_bytes(); digest=hashlib.sha256(data).hexdigest()
        if expected_hash is not None and digest!=expected_hash: raise ValueError('Pinned source hash mismatch.')
        food=set(); reactions=[]; actions=set()
        def side(text):
            result=[]
            for item in text.strip().split('+'):
                parts=item.split()
                if not parts: continue
                if len(parts)==1: result.append((parts[0],1))
                elif len(parts)==2 and parts[0].isdigit() and int(parts[0])>0: result.append((parts[1],int(parts[0])))
                else: raise ValueError(f'Invalid stoichiometry: {item}')
            return tuple(result)
        for line_no,line in enumerate(data.decode('utf-8').splitlines(),1):
            if not line.strip() or line.lstrip().startswith('#'): continue
            if line.startswith('Food:'): food.update(line[5:].split()); continue
            match=re.fullmatch(r'([^:]+):\s*(.*?)\s*\[(.*?)\]\s*(<=>|=>|<=)\s*(.*)',line)
            if not match: raise ValueError(f'Unrecognized source row {line_no}.')
            action,left,cat,arrow,right=match.groups(); action=action.strip()
            if action in actions: raise ValueError('Duplicate source action.')
            actions.add(action); a,b=side(left),side(right); formula=CatalystFormula.parse(cat)
            for direction in (('f','r') if arrow=='<=>' else ('r',) if arrow=='<=' else ('f',)):
                inputs,outputs=(a,b) if direction=='f' else (b,a)
                reactions.append(Reaction(action+'::'+direction,action,inputs,outputs,formula,line_no))
        return cls(tuple(reactions),frozenset(food),digest)

    def subset(self,ids):
        ids=frozenset(ids)
        if not ids<=self.by_id.keys(): raise ValueError('Unknown direction in subset.')
        return replace(self,reactions=tuple(r for r in self.reactions if r.id in ids))

    def delete_actions(self,actions):
        actions=frozenset(actions)
        if not actions<={r.action for r in self.reactions}: raise ValueError('Unknown action in deletion.')
        return replace(self,reactions=tuple(r for r in self.reactions if r.action not in actions))

    def add_food(self,species): return replace(self,food=self.food|frozenset(species))


class RAFAnalysis:
    """Catalyst-free food closure, followed by final-pool RAF support tests."""
    def __init__(self,system):
        self.system=system; self.rows=system.by_id; self.universe=frozenset(self.rows)
        self.inputs={k:r.inputs for k,r in self.rows.items()}; self.outputs={k:r.outputs for k,r in self.rows.items()}

    def validate(self,ids):
        ids=frozenset(ids)
        if not ids<=self.universe: raise ValueError('Direction outside this ambient set.')
        return ids

    def closure(self,ids):
        remaining=set(self.validate(ids)); pool=set(self.system.food); ranks={}; layer=0
        while True:
            ready=sorted(k for k in remaining if self.inputs[k]<=pool)
            if not ready: break
            for k in ready: ranks[k]=layer; pool.update(self.outputs[k])
            remaining.difference_update(ready); layer+=1
        return frozenset(pool),ranks

    def supported(self,k,pool): return self.inputs[k]<=pool and self.rows[k].catalyst.holds(pool)

    def is_raf(self,ids):
        ids=self.validate(ids); pool,_=self.closure(ids)
        return bool(ids) and all(self.supported(k,pool) for k in ids)

    def maximum(self,ids=None):
        current=self.universe if ids is None else self.validate(ids); counts=[len(current)]
        while True:
            pool,_=self.closure(current); kept=frozenset(k for k in current if self.supported(k,pool)); counts.append(len(kept))
            if kept==current: return current,counts
            current=kept

    def capability(self,target,ids=None):
        maximum,trace=self.maximum(ids); pool,_=self.closure(maximum)
        return bool(maximum) and target in pool

    def constructive(self):
        """Stricter than RAF: catalysts must already hold in the prefix pool."""
        pool=set(self.system.food); remaining=set(self.universe); layers=[]
        while True:
            ready=sorted(k for k in remaining if self.supported(k,pool))
            if not ready: break
            layers.append(ready)
            for k in ready: pool.update(self.outputs[k])
            remaining.difference_update(ready)
        return frozenset(self.universe-remaining),layers

    def is_closed(self,ids):
        ids=self.validate(ids)
        if not self.is_raf(ids): return False
        pool,_=self.closure(ids)
        return all(not self.supported(k,pool) for k in self.universe-ids)

    def closed_extension(self,seed):
        current=self.validate(seed)
        if current and not self.is_raf(current): raise ValueError('Closed extension requires an initial RAF or empty seed.')
        while True:
            pool,_=self.closure(current); enlarged=current|frozenset(k for k in self.universe if self.supported(k,pool))
            if enlarged==current: return current
            current=enlarged

    def closed_uniqueness_certificate(self):
        maximum,trace=self.maximum(); active,layers=self.constructive()
        equal=active==maximum
        return {'status':'unique_closed_raf' if equal and maximum else 'no_raf' if not maximum else 'unresolved_bounds',
            'set_equality':equal,'maximum':sorted(maximum),'activated':sorted(active),
            'pruning_counts':trace,'activation_layers':layers}

    def catalogue(self,target,budget=MAX_SUBSETS):
        ids=sorted(self.universe); required=2**len(ids)
        if required>budget: return {'status':'unresolved_budget','required_subsets':required}
        rafs=[]; closed=[]; producing=[]
        for mask in range(1,required):
            candidate=frozenset(ids[i] for i in range(len(ids)) if mask>>i&1)
            if not self.is_raf(candidate): continue
            rafs.append(candidate)
            if self.is_closed(candidate): closed.append(candidate)
            if target in self.closure(candidate)[0]: producing.append(candidate)
        minimal=lambda sets:[s for s in sets if not any(t<s for t in sets)]
        encode=lambda sets:[sorted(s) for s in sets]
        return {'status':'complete','subsets_checked':required-1,'raf_count':len(rafs),
            'irreducible':encode(minimal(rafs)),'closed':encode(closed),'target_rafs':encode(producing),
            'minimal_target_rafs':encode(minimal(producing))}

    def necessary_producers(self,target):
        """Sound unique-producer propagation inside this ambient set; may be incomplete."""
        forced=set(); pending=[target]; seen=set(); reasons=[]; impossible=[]
        while pending:
            species=pending.pop(0)
            if species in seen or species in self.system.food: continue
            seen.add(species); producers=sorted(k for k in self.universe if species in self.outputs[k])
            if not producers: impossible.append(species)
            if len(producers)!=1: continue
            k=producers[0]; reasons.append({'species':species,'only_producer':k})
            if k in forced: continue
            forced.add(k); pending.extend(sorted(self.inputs[k]|self.rows[k].catalyst.mandatory_atoms()))
        return {'forced_reactions':sorted(forced),'reasons':reasons,'impossible_species':impossible,
            'scope':'Necessary reactions for this target in this ambient set; ambiguous producers and optional catalyst branches are not resolved.'}


@dataclass(frozen=True)
class SiphonCertificate:
    species: frozenset

    def check(self,system):
        intersections=[]; violations=[]
        for r in system.reactions:
            if r.outputs&self.species:
                row={'reaction':r.id,'inputs_in_barrier':sorted(r.inputs&self.species),'outputs_in_barrier':sorted(r.outputs&self.species)}
                intersections.append(row)
                if not r.inputs&self.species: violations.append(r.id)
        food_overlap=system.food&self.species
        return {'valid_siphon':not violations,'food_disjoint':not food_overlap,
            'absence_certified':not violations and not food_overlap,'food_overlap':sorted(food_overlap),
            'violations':violations,'surviving_producers':intersections}


class InterventionStudy:
    def __init__(self,system,actions,target,preserve,barrier):
        self.system=system; self.actions=frozenset(actions); self.target=target; self.preserve=preserve; self.barrier=barrier

    def evaluate(self,restoration):
        cut=self.system.delete_actions(self.actions); analysis=RAFAnalysis(cut); barrier=self.barrier.check(cut)
        negative=barrier['absence_certified'] and self.target in self.barrier.species
        witnesses={}
        for action in sorted(self.actions):
            available=self.system.delete_actions(self.actions-{action}); evaluator=RAFAnalysis(available)
            ids=frozenset(restoration.get(action,()))
            valid=ids<=evaluator.universe and evaluator.is_raf(ids)
            pool,ranks=evaluator.closure(ids) if ids<=evaluator.universe else (frozenset(),{})
            witnesses[action]={'valid_target_raf':valid and self.target in pool,'size':len(ids),
                'reactions':sorted(ids),'reactant_ranks':ranks,'final_closure':sorted(pool)}
        return {'deleted_actions':sorted(self.actions),'deleted_directions':sorted(set(self.system.by_id)-set(cut.by_id)),
            'target_absence_by_siphon':negative,'target_capability_by_maxraf':analysis.capability(self.target),
            'preservation_capability':analysis.capability(self.preserve),'barrier':barrier,'restoration':witnesses,
            'inclusion_minimal_cut_certified':negative and all(w['valid_target_raf'] for w in witnesses.values()),
            'global_minimum_cost_claimed':False}


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);base=Path(__file__).resolve().parent
    system=ReactionSystem.read(base/'prokaryotic-network.txt',PINNED_SOURCE_SHA256).add_food(EXTRA_FOOD)
    inputs=json.loads((base/'certificate_inputs.json').read_text())
    if inputs['source_sha256']!=system.source_sha256: raise ValueError('Certificates use a different source.')
    barrier=SiphonCertificate(frozenset(inputs['barrier']))
    study=InterventionStudy(system,DELETED_ACTIONS,TARGET,PRESERVATION_TARGET,barrier); result=study.evaluate(inputs['restoration'])
    parent=RAFAnalysis(system); cut_system=system.delete_actions(DELETED_ACTIONS); cut=RAFAnalysis(cut_system)
    small=RAFAnalysis(system.subset(SMALL_DIRECTIONS)); augmented=RAFAnalysis(system.subset((*SMALL_DIRECTIONS,POOLING_DIRECTION)))
    small_catalogue=small.catalogue(TARGET); augmented_catalogue=augmented.catalogue(TARGET)
    parent_certificate=parent.closed_uniqueness_certificate(); cut_certificate=cut.closed_uniqueness_certificate()
    L=frozenset(SMALL_DIRECTIONS[i] for i in (0,1,4,7,8)); H=frozenset(SMALL_DIRECTIONS)
    rescue_ids=tuple(SMALL_DIRECTIONS[i] for i in (7,8,4,5)); rescue=RAFAnalysis(cut_system.add_food((RESCUE_FOOD,)))
    preserved=RAFAnalysis(cut_system.subset(('R09639::f',)))
    necessity=small.necessary_producers(TARGET)
    result.update({'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':system.source_sha256,
        'food':sorted(system.food),'source_actions':len({r.action for r in system.reactions}),'source_directions':len(system.reactions),
        'small_catalogue':small_catalogue,'augmented_catalogue':augmented_catalogue,'target_necessary_producers':necessity,
        'same_parent_closed_extensions':parent.closed_extension(L)==parent.closed_extension(H)==frozenset(parent_certificate['maximum']),
        'methionine_singleton_valid':preserved.is_raf(preserved.universe) and preserved.capability(PRESERVATION_TARGET),
        'precursor_rescue':{'food_added':RESCUE_FOOD,'reactions':list(rescue_ids),'valid_target_raf':rescue.is_raf(rescue_ids) and TARGET in rescue.closure(rescue_ids)[0]},
        'parent_bounds':{k:v for k,v in parent_certificate.items() if k not in ('maximum','activated','activation_layers')},
        'cut_bounds':{k:v for k,v in cut_certificate.items() if k not in ('maximum','activated','activation_layers')},
        'scope':'Exact finite source computations, not a Lean rerun. Structural capability is not flux, growth or concentration equilibria.'})
    for name,certificate in [('parent',parent_certificate),('postcut',cut_certificate)]:
        (out/f'{name}_closed_certificate.json').write_text(json.dumps(certificate,indent=2)+'\n')
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    def csv_file(name,headers,rows):
        with (out/name).open('w',newline='') as f:
            writer=csv.writer(f);writer.writerow(headers);writer.writerows(rows)
    csv_file('restoration_witnesses.csv',['restored_action','direction','reactant_rank','catalyst_in_final_closure'],
        [[action,k,w['reactant_ranks'].get(k,''),system.by_id[k].catalyst.holds(set(w['final_closure']))] for action,w in result['restoration'].items() for k in w['reactions']])
    csv_file('siphon_producers.csv',['direction','reactants_in_barrier','products_in_barrier'],
        [[r['reaction'],';'.join(r['inputs_in_barrier']),';'.join(r['outputs_in_barrier'])] for r in result['barrier']['surviving_producers']])
    csv_file('activation_layers.csv',['model','synchronous_layer','new_directions','cumulative_directions'],
        [[name,i+1,len(layer),sum(map(len,c['activation_layers'][:i+1]))] for name,c in [('parent',parent_certificate),('postcut',cut_certificate)] for i,layer in enumerate(c['activation_layers'])])
    csv_file('small_reactions.csv',['index','direction','action','reactants','products','catalyst_formula_tree'],
        [[i,k,system.by_id[k].action,json.dumps(system.by_id[k].reactants),json.dumps(system.by_id[k].products),json.dumps(system.by_id[k].catalyst.to_data())] for i,k in enumerate(SMALL_DIRECTIONS)])
    lines=[f'Source: {len(system.reactions)} directions, {len(system.food)} food species/markers.',
        f'Cut removes {len(result["deleted_directions"])} directions; selective inclusion-minimal certificate: {result["inclusion_minimal_cut_certified"]}.',
        f'Restoration witness sizes: {[w["size"] for w in result["restoration"].values()]}; methionine singleton: {result["methionine_singleton_valid"]}.',
        f'H: {small_catalogue["raf_count"]} RAFs, {len(small_catalogue["irreducible"])} irrRAFs, {len(small_catalogue["closed"])} closed RAFs, {len(small_catalogue["target_rafs"])} target RAF.',
        f'Adding the food-enabled pooling row: {len(augmented_catalogue["closed"])} closed RAF.',
        f'Parent bounds: {parent_certificate["pruning_counts"]}; A=M: {parent_certificate["set_equality"]}.',
        f'Postcut bounds: {cut_certificate["pruning_counts"]}; A=M: {cut_certificate["set_equality"]}.',
        'No complete parent irrRAF catalogue, global cut optimum, depletion or kinetic uniqueness is claimed.']
    (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
    fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    for name,c,color in [('Parent',parent_certificate,'#187f91'),('After cut',cut_certificate,'#b95024')]:
        cumulative=np.cumsum([len(layer) for layer in c['activation_layers']])
        axes[0].step(range(1,len(cumulative)+1),cumulative,where='post',label=name,color=color)
        axes[0].axhline(len(c['maximum']),color=color,ls='--',lw=.8)
        axes[1].plot(range(len(c['pruning_counts'])),c['pruning_counts'],'o-',label=name,color=color)
    axes[0].set(xlabel='Constructive activation layer',ylabel='Activated directions',title='Constructive lower bound reaches maxRAF')
    axes[1].set(xlabel='Pruning step (last equality retained)',ylabel='Directions in RAF upper bound',title='Catalysts checked in final reactant closure')
    for ax in axes:ax.grid(alpha=.2);ax.legend(fontsize=8)
    fig.savefig(out/'parent_bounds.png',dpi=180);fig.savefig(out/'parent_bounds.svg');plt.close(fig)
    fig,ax=plt.subplots(figsize=(10,4.2),layout='constrained')
    rows=[set([SMALL_DIRECTIONS[0]]),set([SMALL_DIRECTIONS[7]]),set(L),set(H),set((*H,POOLING_DIRECTION))]
    labels=['irrRAF seed 0','irrRAF seed 7','Closed L in H','Only valine RAF in H','Only closed RAF in H + q']
    matrix=np.array([[k in row for k in (*SMALL_DIRECTIONS,POOLING_DIRECTION)] for row in rows])
    ax.imshow(matrix,cmap=matplotlib.colors.ListedColormap(['#f0f0ed','#187f91']),vmin=0,vmax=1,aspect='auto')
    ax.set_yticks(range(5),labels);ax.set_xticks(range(11),[str(i) for i in range(10)]+['q']);ax.set(xlabel='Literal subsystem reaction index (source mapping in CSV)',title='Reaction membership in catalytic seeds and\nvaline-producing sets')
    ax.set_xticks(np.arange(-.5,11),minor=True);ax.set_yticks(np.arange(-.5,5),minor=True);ax.grid(which='minor',color='white',linewidth=2);ax.tick_params(which='minor',bottom=False,left=False)
    fig.savefig(out/'subsystem.png',dpi=180);fig.savefig(out/'subsystem.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__)),
        'input_sha256':{name:digest(base/name) for name in ('prokaryotic-network.txt','certificate_inputs.json')},'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()