"""Same means, different sisters, opposite extinction decisions."""
# ------------------------- EDITABLE INPUTS -------------------------
ACQUISITION_SCALE='1/10'
ACQUISITION_PROFILE=('1/10','1/10','1/10','2/5','1/10','2/5')
ACTION_A_ERASURE='3/10'
ACTION_A_EXTRA_DEATH='0'
ACTION_B_ERASURE='1/100'
ACTION_B_EXTRA_DEATH='120553/1000000'
PULSE_DURATION='10'
RESISTANT_DIVISION='1/10'
RESISTANT_DEATH='1/50'
INITIAL_COUNTS=(0,0,0,0,0,1,0) # UU,UR,RR,AU,AR,AA,M
CORRECTION_ERROR_ALLOWANCE='4/1000000'
OUTPUT_DIRECTORY='outputs'
# Synthetic model rates and per-division probabilities, not measured dosing.
# The fixed reference certificate is kept separate from custom scenarios.
# ------------------------------------------------------------------
import argparse,csv,hashlib,json,platform
from fractions import Fraction as F
from pathlib import Path
import numpy as np
import sympy as sp
from branching import *
from validated import certify
from consequence_kernel import scalar_floor,mean_matrix,linear_flow

MANUSCRIPT_SHA256='96c8eb87625ee95044617d81f9a11ed08d44cf309c63d8036f84ea0775fbf3ff'


def subtract(a,b):return (a[0]-b[1],a[1]-b[0])
def bounds(run):return tuple(map(F,run['bounds'][5]))
def winner(gap):return 'AB' if gap[0]>0 else 'BA' if gap[1]<0 else 'unresolved'


def reference_certificates():
    A=(F(3,10),F(0),F(10));B=(F(1,100),F(120553,10**6),F(10))
    runs={law+order:certify([A,B] if order=='AB' else [B,A],law) for law in ['J','I'] for order in ['AB','BA']}
    q={k:bounds(v) for k,v in runs.items()};gaps={law:subtract(q[law+'AB'],q[law+'BA']) for law in ['J','I']}
    assert gaps['J'][1]<-F(16,10**7) and gaps['I'][0]>F(16,10**7)
    left,center,right=map(F,['.120545','.120553','.120561']);K=F(6912);boundary={}
    for law,end in [('J',left),('I',right)]:
        BB=(F(1,100),end,F(10));rr={order:certify([A,BB] if order=='AB' else [BB,A],law) for order in ['AB','BA']}
        gap=subtract(bounds(rr['AB']),bounds(rr['BA']))
        num=subtract(gaps[law],gap) if law=='J' else subtract(gap,gaps[law]);distance=abs(center-end)
        secant=tuple(v/distance for v in num);derivative=(secant[0]-K*(right-left),secant[1]+K*(right-left))
        assert derivative[1]<0 and (gap[1]<0 if law=='J' else gap[0]>0)
        boundary[law]=dict(endpoint=end,gap=gap,secant=secant,derivative=derivative,runs=rr)
    return dict(runs=runs,AA_extinction=q,order_gaps=gaps,preferred={law:winner(gap) for law,gap in gaps.items()},
        continuum=dict(interval=(left,right),curvature_bound=K,endpoints=boundary,scope='Four fresh endpoint flows plus the analytic curvature and mean-value argument prove the whole interval. Not a maximal reversal set.'),
        evidence='Fresh exact outward dyadic Taylor runs; no saved certificate values loaded. Conventional soundness/probability proofs, not Lean.'),q,gaps


def decision_consequences(q,gaps):
    mJ=(-gaps['J'][1],-gaps['J'][0]);mI=gaps['I']
    errors={s:subtract(q['I'+s],q['J'+s]) for s in ['AB','BA']};paired=subtract(errors['AB'],errors['BA']);corrected=subtract(gaps['I'],paired)
    regret=tuple(mJ[i]*mI[i]/(mJ[i]+mI[i]) for i in range(2))
    founders={}
    for law,h,l in [('J',q['JBA'],q['JAB']),('I',q['IAB'],q['IBA'])]:
        increasing=l[0]**31*(1-l[1])-h[1]**31*(1-h[0]);decreasing=l[1]**32*(1-l[0])-h[0]**32*(1-h[1])
        assert increasing>0 and decreasing<0
        founders[law]=dict(unique_maximum=32,increment31_lower=increasing,increment32_upper=decreasing,advantage32=(h[0]**32-l[1]**32,h[1]**32-l[0]**32))
    symmetric=subtract(gaps['I'],(-F(CORRECTION_ERROR_ALLOWANCE),F(CORRECTION_ERROR_ALLOWANCE)))
    return dict(direct_misspecification_cost=mJ,independent_margin=mI,dependence_errors=errors,paired_correction=paired,
        corrected_gap=corrected,corrected_decision=winner(corrected),symmetric_corrected_gap=symmetric,symmetric_decision=winner(symmetric),
        minimax_regret=regret,founders=founders,
        scope='AA founders, deterministic schedules and this continuation. No arbitrary-mixture ranking or global optimization. Minimax randomizes a complete course under a two-model information restriction; it is not the direct mistake cost.')


def cone_check():
    states=source.STATES;pre=lambda i,j:states[i][0]<=states[j][0] and states[i][1]>=states[j][1]
    covers=[(i,j) for i in range(6) for j in range(6) if i!=j and pre(i,j) and not any(k not in [i,j] and pre(i,k) and pre(k,j) for k in range(6))]
    upsets=[[F((mask>>i)&1) for i in range(6)] for mask in range(1,64) if all(not ((mask>>i)&1) or not pre(i,j) or (mask>>j)&1 for i in range(6) for j in range(6))]
    assert len(covers)==6 and len(upsets)==7 and all(source.KAPPA[i]<=source.KAPPA[j] for i,j in covers)
    values={}
    for label,(e,c) in [('A',(F(3,10),F(0))),('B',(F(1,100),F(120553,10**6)))]:
        A=mean_matrix(e,c,F(0));vv=[sum((A[j][k]-A[i][k])*u[k] for k in range(6)) for i,j in covers for u in upsets if u[i]==u[j]]
        assert len(vv)==32 and min(vv)>=0;values[label]=vv
    return dict(covers=covers,upsets=upsets,checks=values,scope='Exact finite cone-invariance checks. The Harris covariance and positive transport conclusions use the conventional argument; an asymptotic sign alone does not prove the finite reversal.')


def continuation_bounds(decisions):
    A=sp.Matrix(mean_matrix(F(3,10),F(0),F(0)));f=sp.Matrix([F(1,10)*k for k in source.KAPPA]);v=[F(x) for x in -A.inv()*f]
    assert min(v)>0 and all(sum(F(A[i,j])*v[j] for j in range(6))==-F(f[i]) for i in range(6))
    a=(F(3,10),F(0),F(10));b=(F(1,100),F(120553,10**6),F(10));w=list(map(F,[14,11,10,84,43,107]));U={};C={}
    for label,phases in [('AB',[a,b]),('BA',[b,a])]:
        U[label]=tuple(F(4,5)*x for x in linear_flow(phases,v));C[label]=tuple(x/10 for x in linear_flow(phases,w))
    mJ=decisions['direct_misspecification_cost'][0];mI=decisions['independent_margin'][0];eta=F(3,10**5)
    assert eta*U['BA'][1]<mJ and eta*U['AB'][1]<mI
    def exp_lower(x):
        s=term=F(1)
        for k in range(1,200):term*=x/k;s+=term
        return s
    finite={}
    for H in [151,152]:
        lower=exp_lower(F(9,100)*H);errJ=C['BA'][1]/lower;errI=C['AB'][1]/lower
        finite[str(H)]=dict(J_error_upper=errJ,I_error_upper=errI,both_preserved=errJ<mJ and errI<mI,common_added_exposure=F(29,100)*H)
    assert finite['152']['both_preserved'] and not finite['151']['both_preserved']
    return dict(occupation_reward=v,U=U,C=C,residual_acquisition_scale=eta,residual_errors=(eta*U['BA'][1],eta*U['AB'][1]),finite_clearing=finite,
        exposure_ratio=F(29,100)*152/F(29,10),scope='Separate one-sided perturbation scenarios. Do not combine them without a new bound. Finite clearing requires autonomous existing M lineages and no immigration; H=152 is sufficient, not optimal. The baseline indefinite clearing is not drug-free and has infinite charged exposure.')


def encode(v):
    if isinstance(v,F):return str(v)
    if isinstance(v,np.ndarray):return v.tolist()
    if isinstance(v,np.generic):return v.item()
    raise TypeError(type(v).__name__)


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',default=OUTPUT_DIRECTORY);args=parser.parse_args();out=Path(args.output);out.mkdir(exist_ok=True)
    def dump(name,data):(out/name).write_text(json.dumps(data,indent=2,default=encode)+'\n',encoding='utf-8')
    def table(name,header,rows):
        with (out/name).open('w',newline='',encoding='utf-8') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    certificate,q,gaps=reference_certificates();dump('fresh_flow_certificates.json',certificate)
    decisions=decision_consequences(q,gaps);dump('decision_consequences.json',decisions);dump('cone_invariance.json',cone_check())
    print('Eight fresh outward flow enclosures prove opposite orders and the full killing-rate interval.',flush=True)
    dump('scalar_feedback_floor.json',scalar_floor());dump('continuation_bounds.json',continuation_bounds(decisions))
    acquisition=Acquisition(F(ACQUISITION_SCALE),tuple(map(F,ACQUISITION_PROFILE)));resistant=ResistantLineage(F(RESISTANT_DIVISION),F(RESISTANT_DEATH))
    a=Phase(Action(F(ACTION_A_ERASURE),F(ACTION_A_EXTRA_DEATH)),F(PULSE_DURATION));b=Phase(Action(F(ACTION_B_ERASURE),F(ACTION_B_EXTRA_DEATH)),F(PULSE_DURATION))
    models={'J':BranchingSource(ComplementarySisters(),acquisition,resistant),'I':BranchingSource(IndependentSisters(),acquisition,resistant)}
    if len(INITIAL_COUNTS)!=7 or any(not isinstance(n,int) or n<0 for n in INITIAL_COUNTS):raise ValueError('Seven nonnegative integer founder counts required.')
    matched={};configured={};mechanisms={}
    for order,phases in [('AB',[a,b]),('BA',[b,a])]:
        matched[order]=dict(phase_matrices=[[[str(v) for v in row] for row in models['J'].mean_matrix(p.action)] for p in phases],means={})
        for phase in phases:assert models['J'].mean_matrix(phase.action)==models['I'].mean_matrix(phase.action)
        for law,model in models.items():
            terminal,slack=model.clearing_terminal();ext=model.compose(phases,terminal)
            survival=1-np.prod([ext[i]**n for i,n in enumerate(INITIAL_COUNTS)])
            endpoints=dict(eventual_total_survival=1-ext[5],total_survival_at_deadline=1-model.compose(phases,[0]*7)[5],
                M_present_at_deadline=1-model.compose(phases,[1]*6+[0])[5],M_ever_appeared=1-model.compose(phases,[1]*6+[0],killed_acquisition=True)[5])
            configured[law+order]=dict(extinction=ext,configured_population_survival=survival,AA_endpoint_diagnostics=endpoints,clearing_slack=slack)
            matched[order]['means'][law]=model.means(phases,INITIAL_COUNTS)
    dump('configured_comparison.json',dict(initial_configuration=INITIAL_COUNTS,acquisition_scale=acquisition.scale,acquisition_profile=acquisition.profile,
        resistant_rates=vars(resistant),phases={label:dict(erasure=p.action.erasure,extra_death=p.action.extra_molecular_death,duration=p.duration) for label,p in [('A',a),('B',b)]},
        results=configured,evidence='Custom numerical scenario. The fixed interval certificates are not silently transferred after editing inputs.'))
    dump('matched_means.json',dict(comparisons=matched,scope='Law-to-law equality under each common deterministic schedule; AB and BA need not have equal means, and shared feedback does not preserve this argument.'))
    # Mechanism and smooth boundary curves deliberately keep the paper reference.
    ar=Phase(Action(F(3,10)));br=Phase(Action(F(1,100),F(120553,10**6)))
    for order,phases in [('AB',[ar,br]),('BA',[br,ar])]:
        result,t,states=mechanism(phases);mechanisms[order]=result
        assert result['identity_residual']<1e-10
        table('mechanism_'+order+'.csv',['backward_time','h_AA','epsilon2_delta_AA','finite_error_AA','response_AA','absolute_bound_AA'],zip(t,states[:,17],.01*states[:,23],states[:,11]-states[:,5],states[:,29],states[:,35]))
    dump('covariance_mechanism.json',dict(results=mechanisms,evidence='Numerical sensitivities and the exact finite-amplitude response identity; no certified asymptotic remainder is claimed.'))
    curve=[]
    for cb in np.linspace(.12054,.12057,21):
        bb=Phase(Action(F(1,100),F(str(cb))));row=dict(c=cb)
        for law,cls in [('J',ComplementarySisters),('I',IndependentSisters)]:
            m=BranchingSource(cls());row[law]=m.compose([ar,bb])[5]-m.compose([bb,ar])[5]
        curve.append(row)
    table('order_gap_curve.csv',['c_B','J_gap','I_gap'],[(r['c'],r['J'],r['I']) for r in curve])
    founderrows=[]
    for n in range(1,201):
        row=[n]
        for h,l in [(q['JBA'],q['JAB']),(q['IAB'],q['IBA'])]:row.append(float(sum(h)/2)**n-float(sum(l)/2)**n)
        founderrows.append(row)
    table('founder_advantages.csv',['AA_founders','J_advantage','I_advantage'],founderrows)
    plot(out,curve,founderrows,mechanisms,decisions)
    print(f'Certified preferred order: J={certificate["preferred"]["J"]}, I={certificate["preferred"]["I"]}; direct mistake cost {float(decisions["direct_misspecification_cost"][0]):.12g}.',flush=True)
    print('Exact mean equality, founder maximum at 32, scalar feedback floor and separately priced continuation bounds checked.',flush=True)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();here=Path(__file__).parent
    dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),python=platform.python_version(),
        module_sha256={p.name:digest(p) for p in sorted(here.glob('*.py')) if p.name not in ['example.py','test_example.py']},
        output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))


def plot(out,curve,founders,mechanisms,decisions):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for law in ['J','I']:axs[0].plot([r['c'] for r in curve],[1e6*r[law] for r in curve],label=law+(' complementary' if law=='J' else ' independent'))
    axs[0].axvspan(.120545,.120561,color='gray',alpha=.15,label='Certified interval');axs[0].axhline(0,color='black',lw=.7)
    axs[0].set(xlabel='Additional death rate in B',ylabel='AB minus BA extinction (millionths)',title='Extinction advantage of AB over BA by\nsister law');axs[0].legend(fontsize=8)
    for i,law in enumerate(['J','I']):axs[1].plot([r[0] for r in founders],[1e6*r[i+1] for r in founders],label=law)
    axs[1].axvline(32,color='gray',ls=':');axs[1].set(xlabel='Identical AA founders',ylabel='Extinction advantage (millionths)',title='Pulse-order extinction advantage versus\nfounder count');axs[1].legend()
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'decisions.png',dpi=180);fig.savefig(out/'decisions.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for order in ['AB','BA']:
        data=np.genfromtxt(out/('mechanism_'+order+'.csv'),delimiter=',',names=True)
        axs[0].plot(data['backward_time'],1e6*data['epsilon2_delta_AA'],label=order+' leading order')
        axs[0].scatter([20],[1e6*mechanisms[order]['finite_dependence_error'][5]],marker='*',s=90)
    axs[0].axvline(10,color='gray',ls=':');axs[0].set(xlabel='Backward time from the deadline',ylabel='Dependence correction (millionths)',title='Sister-dependence correction through each\nschedule');axs[0].legend(fontsize=8)
    vals=[sum(decisions['independent_margin'])/2,sum(decisions['paired_correction'])/2,sum(decisions['corrected_gap'])/2]
    axs[1].bar(['I order gap','Paired correction','Corrected J gap'],[1e6*float(v) for v in vals]);axs[1].axhline(0,color='black',lw=.7)
    axs[1].set(ylabel='Probability difference (millionths)',title='Dependence correction and pulse-order\ndecision margin');axs[1].tick_params(axis='x',labelsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'dependence.png',dpi=180);fig.savefig(out/'dependence.svg');plt.close(fig)


if __name__=='__main__':main()