"""Exposure, amplitude and duration are distinct limits in an inherited population."""
# ------------------------- EDITABLE INPUTS -------------------------
SITE_COUNT=2
BASELINE_ERASURE='1/100'
BASAL_WRITING='1/100'
RECRUITED_WRITING='1'
RECRUITED_ERASURE='1/2'
DIVISION_RATE='1/10'
PROTECTED_DIVISION_RATE='1/10'
PROTECTED_DEATH_RATE='1/100'
UNPROTECTED_DEATH_RATE='3/10'
FOUNDERS=1
INITIAL_MARK_COUNTS=(2,0)
TARGET_RISK=.01
ERASER_AMPLITUDE='29/100'
EXPOSURE_BUDGET='58/5'
HORIZON='100'
ADDED_DEATH_AMPLITUDE='29/100'
RUN_FEEDBACK_DEMONSTRATION=False
RANDOM_SEED=20260923
OUTPUT_DIRECTORY='outputs'
# These are synthetic reaction/demographic rates, not a treatment calibration.
# Fixed paper certificate replays below remain separate from edited scenarios.
# ------------------------------------------------------------------
import argparse,csv,hashlib,json,platform
from dataclasses import replace
from fractions import Fraction as F
from pathlib import Path
import numpy as np
from population import MolecularRates,InheritedPopulation,Action,Phase,population_survival,FeedbackSimulator
from certificates import ExposureCertificate,AmplitudeCertificate,WeightedDrift,collatz_bounds,m_matrix_test,mv,population_barrier_generator_identity,finite_time_floor,concentration_lower_bound

MANUSCRIPT_SHA256='1be8ca10a66dcac8fbc1175361298ce9ac2b729529221e7c6a94a8e28687b77c'


def configured_source():
    return InheritedPopulation(MolecularRates(SITE_COUNT,F(BASAL_WRITING),F(RECRUITED_WRITING),F(RECRUITED_ERASURE),F(BASELINE_ERASURE),F(DIVISION_RATE),F(PROTECTED_DIVISION_RATE),F(PROTECTED_DEATH_RATE),F(UNPROTECTED_DEATH_RATE)))


def phases_for(exposure,amplitude,horizon,placement):
    B,v,T=map(F,(exposure,amplitude,horizon))
    if B<0 or v<=0 or T<=0 or B>v*T:raise ValueError('Schedule exceeds the amplitude/horizon exposure capacity.')
    if B==0:return [Phase(T,Action())]
    pulse=B/v
    if placement=='spread':return [Phase(T,Action(B/T))]
    on=Phase(pulse,Action(v))
    if pulse==T:return [on]
    off=Phase(T-pulse,Action())
    if placement=='early':return [on,off]
    if placement=='late':return [off,on]
    raise ValueError('Unknown schedule placement.')


def replay_sites(data):
    results=dict(mean_sign_endpoints={},exposure={},amplitude={},added_death={},spectral_at_point={})
    for key,row in data['sites'].items():
        N=int(key);source=InheritedPopulation(MolecularRates(sites=N));lo,hi=F(row['ec_lower']),F(row['ec_upper'])
        al,ah=Action(lo-F(1,100)),Action(hi-F(1,100));wl=tuple(map(F,row['witness_lower']));wh=tuple(map(F,row['witness_upper']))
        stored_cwl=collatz_bounds(source,al,wl);stored_cwh=collatz_bounds(source,ah,wh)
        # Most released six-decimal Perron vectors straddle zero this close to
        # the threshold. Propose finer vectors and accept ONLY exact signs.
        def sharpen(action):
            A=source.mean_matrix(action);ev,vec=np.linalg.eig(np.array(A,dtype=float))
            v=np.real(vec[:,np.argmax(ev.real)]);v*=np.sign(v.sum());v/=v.min()
            for digits in [10,12,14]:
                w=tuple(F(str(round(float(x),digits))) for x in v)
                if min(w)<=0:continue
                bracket=collatz_bounds(source,action,w)
                if bracket[0]>0 or bracket[1]<0:return w,bracket
            raise RuntimeError('Cannot certify endpoint sign with the proposed weights.')
        wl,cwl=sharpen(al);wh,cwh=sharpen(ah)
        assert cwl[0]>0 and cwh[1]<0 and hi-lo==F(1,10**6)
        lowtest=m_matrix_test(source,al);hitest=m_matrix_test(source,ah)
        assert not lowtest['subcritical'] and hitest['subcritical']
        w=tuple(map(F,row['w_half']));drift=WeightedDrift(w,F(row['gamma_half']),Action(F(49,100)));slack=drift.verify(source)
        results['mean_sign_endpoints'][key]=dict(intrinsic_erasure_bracket=(lo,hi),amplitude_bracket=(lo-F(1,100),hi-F(1,100)),
            lower_collatz=cwl,upper_collatz=cwh,lower_weights=wl,upper_weights=wh,
            published_weight_collatz=(stored_cwl,stored_cwh),lower_M_matrix=lowtest,upper_M_matrix=hitest,
            contraction_gamma=drift.gamma,contraction_slack=slack,founder_weight_ratio=w[-1]/min(w),
            scope='Fresh endpoint sign certificates. A global first-crossing/uniqueness assertion is not inferred from bisection or endpoint checks.')
    for key,row in data['exposure'].items():
        source=InheritedPopulation(MolecularRates(sites=int(key)));cert=ExposureCertificate(tuple(map(F,row['q'])),F(row['c']))
        results['exposure'][key]=dict(q=cert.q,c=cert.c,eta=1-cert.q[-1],residuals=cert.verify(source),necessary_one_founder_point01=cert.necessary_exposure(1,.01,source.size-1))
    for key,row in data['amplitude'].items():
        source=InheritedPopulation(MolecularRates(sites=row['N']));cert=AmplitudeCertificate(tuple(map(F,row['q'])),F(row['e_hi'])-F(1,100))
        z=[0]*source.size;z[-1]=1
        results['amplitude'][key]=dict(N=row['N'],maximum=cert.maximum,q=cert.q,residuals=cert.verify(source),fully_marked_floor=cert.survival_floor_exact(z,source),any_type_floor=1-max(cert.q))
    for key,row in data['spectral_at_three_tenths'].items():
        source=InheritedPopulation(MolecularRates(sites=int(key)));cw=collatz_bounds(source,Action(F(29,100)),tuple(map(F,row['witness'])))
        assert cw==tuple(map(F,[row['lower'],row['upper']]))
        results['spectral_at_point'][key]=dict(bounds=cw,scope='Mean growth at one constant action only.')
    for key,pair in data['added_death'].items():
        N=int(key.split('_')[0][1:]);source=InheritedPopulation(MolecularRates(sites=N));eraser=F(0) if key.endswith('no_eraser') else F(29,100)
        lo,hi=map(F,pair);low=m_matrix_test(source,Action(eraser,lo));high=m_matrix_test(source,Action(eraser,hi))
        assert not low['subcritical'] and high['subcritical'] and hi<=F(9,100)
        results['added_death'][key]=dict(bracket=(lo,hi),lower=low,upper=high,scope='Added death makes the Metzler matrix componentwise smaller, so this endpoint bracket has a monotone threshold interpretation.')
    return results


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);here=Path(__file__).parent
    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)
    data=json.loads((here/'site_certificates.json').read_text());sites=replay_sites(data);dump('rechecked_site_certificates.json',sites)
    source=InheritedPopulation();q=tuple(map(F,['.94','.982','.987','.27','.735','.235']));exposure=ExposureCertificate(q,F(67,73))
    weight=WeightedDrift(tuple(map(F,[1419,1105,1000,8423,4305,10765])),F(21,200),Action(F(29,100)))
    exposure_res=exposure.verify(source);weight_slack=weight.verify(source)
    assert min(-v for v in exposure_res['baseline_residual'])==F(1,5000)
    assert min(exposure_res['control_slack'])==0
    rows=[]
    for count in [1,10,100,1000,1000000]:
        z=[0]*6;z[-1]=count
        rows.append(dict(founders=count,necessary=exposure.necessary_exposure(count,TARGET_RISK,5),**weight.sufficient(source,z,TARGET_RISK)))
    table('founder_exposure_bounds.csv',['founders','necessary_exposure','sufficient_exposure','sufficient_duration'],[(r['founders'],r['necessary'],r['eraser_exposure'],r['duration']) for r in rows])
    # Directly check the exact population generator, not a product assumption.
    identity=population_barrier_generator_identity(source,q,[2,0,1,3,0,4],Action(F(13,100)))
    segment=[]
    for a in [F(0),F(1,3),F(1)]:
        qa=[1-a*(1-x) for x in q];h=[1-x for x in q];lhs=source.phi(qa);D=source.daughter(h,h)
        rhs=[a*p-a*(1-a)*b*dd for p,b,dd in zip(source.phi(q),source.birth,D)];assert lhs==rhs
        segment.append(dict(a=a,residual=lhs))
    dump('two_site_barriers.json',dict(exposure_q=q,c=exposure.c,residuals=exposure_res,weight=weight.weights,weight_slack=weight_slack,founder_bounds=rows,
        population_generator_identity=identity,segment_checks=segment,scope='Common predictable action with a pathwise rate-time budget. Exact source inequalities feed the conventional probability theorem; no independence is assumed under shared feedback.'))
    print('Exact exposure, amplitude and mean-drift witnesses replayed against reconstructed sources N=2..8.',flush=True)
    from validated_pgf import certify
    policy=certify([(F(3,10),F(40)),(F(1,100),F(60))]);assert F(policy['survival_AA_upper'])<F(1,100)
    policy['evidence']='Fresh outward-rounded dyadic Taylor enclosure of the fixed two-site source. Degree 11, step 1/20, invariant cube and one-sided Lipschitz bound; not a floating ODE tolerance.'
    policy['exposure']=F(58,5);policy['founders']=1;policy['endpoint']='Any living descendant at time 100; eraser withdrawn at 40 while original killing continues.'
    eps=F(policy['survival_AA_upper']);policy['multiple_founder_upper']={str(n):1-(1-eps)**n for n in [1,10,100]}
    dump('validated_policy_replay.json',policy)
    # All custom configurations use their own source, never the fixed witnesses.
    configured=configured_source();idx=configured.index[INITIAL_MARK_COUNTS];z=[0]*configured.size;z[idx]=FOUNDERS
    configured_records=[]
    for placement in ['early','late','spread']:
        phases=phases_for(EXPOSURE_BUDGET,ERASER_AMPLITUDE,HORIZON,placement);u=configured.schedule(phases)
        configured_records.append(dict(placement=placement,extinction_vector=u,population_survival=population_survival(u,z),
            phases=[dict(duration=p.duration,eraser=p.action.eraser) for p in phases]))
    dump('configured_schedules.json',dict(states=configured.states,rates=vars(configured.rates),initial_configuration=z,schedules=configured_records,evidence='Numerical PGF for predetermined schedules; founder independence is valid only here.'))
    sweep=[]
    for B in np.linspace(0,20,21):
        row=dict(exposure=float(B))
        for placement in ['early','late','spread']:
            u=source.schedule(phases_for(F(str(B)),F(29,100),100,placement));row[placement]=1-u[5]
        sweep.append(row)
    table('schedule_sweep.csv',['exposure','early','late','spread'],[[r[k] for k in ['exposure','early','late','spread']] for r in sweep])
    # Endpoints distinguish live-at-100 from eventual survival after withdrawal.
    # Long baseline integration from zero approximates the minimal fixed point;
    # this is a numerical diagnostic, not a certified root enclosure.
    off=source.schedule([Phase(2000,Action())]);terminal_full=np.full(6,.1)
    early=phases_for(F(58,5),F(29,100),100,'early');on40=[Phase(40,Action(F(29,100)))]
    endpoints={name:1-source.schedule(phases,terminal)[np.array([5,2])] for name,phases,terminal in [
        ('live_at_100',early,np.zeros(6)),('eventual_background_retained',early,off),('full_withdrawal_at_100',early,terminal_full),('full_withdrawal_at_40',on40,terminal_full)]}
    lo=source.taylor_coefficients(Action());hi=source.taylor_coefficients(Action(F(29,100)))
    difference=hi[4][2]-lo[4][2];assert difference==-F(49619,600000000) and all(hi[k][2]==lo[k][2] for k in range(4))
    rr_low=source.schedule([Phase(1,Action())])[2];rr_high=source.schedule([Phase(1,Action(F(29,100)))])[2]
    dump('preparation_and_withdrawal.json',dict(founder_order=['AA','RR'],endpoints=endpoints,baseline_terminal_extinction=off,
        RR_fourth_order_extinction_difference=difference,RR_t1_numeric_extinction_difference=rr_high-rr_low,
        scope='The negative exact fourth-order coefficient disproves universal erasure monotonicity near time zero. t=1 and withdrawal curves are numerical. Full withdrawal stipulates all-state death .01 and division .1, with extinction vector .1.'))
    # Protected-division uncertainty: affine endpoint checks establish each interval.
    robust=[]
    for row in json.loads((here/'source_certificates.json').read_text()):
        results=[]
        for b in map(F,row['birth_interval']):
            ss=InheritedPopulation(MolecularRates(sites=row['N'],protected_division=b))
            ec=ExposureCertificate(tuple(map(F,row['q'])),F(row['c']));wc=WeightedDrift(tuple(map(F,row['w'])),F(row['gamma']),Action(F(29,100)))
            results.append(dict(protected_division=b,exposure=ec.verify(ss),weighted_slack=wc.verify(ss)))
        robust.append(dict(N=row['N'],interval=row['birth_interval'],checks=results,scope='Same witness at both endpoints; affine dependence extends inequalities through this fixed-parameter interval.'))
    dump('demographic_sensitivity.json',robust)
    added=[]
    for N in range(2,9):
        ss=InheritedPopulation(MolecularRates(sites=N));drift=WeightedDrift((F(1),)*ss.size,F(1,5),Action(F(0),F(29,100)));zz=[0]*ss.size;zz[-1]=FOUNDERS
        added.append(dict(N=N,slack=drift.verify(ss),**drift.sufficient(ss,zz,TARGET_RISK)))
    # .025 repairs N=6 at e=.3; exact positive weight proposed numerically.
    ss=InheritedPopulation(MolecularRates(sites=6));act=Action(F(29,100),F(1,40));A=ss.mean_matrix(act)
    eig,vec=np.linalg.eig(np.array(A,dtype=float));v=np.real(vec[:,np.argmax(eig.real)]);v*=np.sign(v.sum());v/=v.min()
    w=tuple(F(str(round(float(x),8))) for x in v);bounds=collatz_bounds(ss,act,w);assert bounds[1]<0
    configured_gamma=min(d+F(ADDED_DEATH_AMPLITUDE)*p-b for d,p,b in zip(configured.death,configured.protected,configured.birth))
    configured_second=dict(amplitude=F(ADDED_DEATH_AMPLITUDE),gamma=configured_gamma,status='unit_weight_drift_not_negative')
    if configured_gamma>0:
        configured_weight=WeightedDrift((F(1),)*configured.size,configured_gamma,Action(F(0),F(ADDED_DEATH_AMPLITUDE)))
        configured_second.update(status='certified',slack=configured_weight.verify(configured),bound=configured_weight.sufficient(configured,z,TARGET_RISK))
    dump('second_actuator.json',dict(configured=configured_second,uniform_policy=added,six_site_repair=dict(action=vars(act),weights=w,collatz=bounds),
        scope='Uniform sufficient amplitude is >.09, not an exact N-independent critical threshold. At .29 the sufficient exposure is 1.45*log(n/risk). Costs of distinct actuators must be declared separately.'))
    Breq=exposure.necessary_exposure(FOUNDERS,TARGET_RISK,5)
    dump('duration_and_units.json',dict(finite_time_floor_at_100=finite_time_floor(FOUNDERS,100,.3),
        conditional_illustration=concentration_lower_bound(Breq,100,.29,1.,.1,1.),
        scope='Illustrative uncalibrated saturating engagement v(C)=vmax*C/(K+C), and one-compartment clearance. AUC/amount outputs are conditional necessary bounds, not sufficient schedules; exposure is dimensionless additional reaction-rate*time. Time and exposure lower bounds are simultaneous, not multiplied.'))
    if RUN_FEEDBACK_DEMONSTRATION:
        sim=FeedbackSimulator(configured,F(ERASER_AMPLITUDE),F(EXPOSURE_BUDGET),RANDOM_SEED)
        def feedback_policy(t,count,spent,history):return float(F(ERASER_AMPLITUDE)) if sum(count[i] for i,p in enumerate(configured.protected) if p)>0 else 0.
        dump('feedback_path.json',sim.run(z,float(F(HORIZON)),feedback_policy))
    plot(out,rows,sweep,sites,policy)
    print(f'Fresh finite-time enclosure: one AA founder, exposure 11.6, survival upper {float(eps):.15f} at time 100.',flush=True)
    print('Amplitude .29 leaves certified survival floors at N=5..8; a second actuator restores contraction. Exact and numerical endpoints remain distinct.',flush=True)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    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']},
        input_sha256={p.name:digest(p) for p in sorted(here.glob('*certificates.json'))},
        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,rows,sweep,sites,policy):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    axs[0].semilogx([r['founders'] for r in rows],[r['necessary'] for r in rows],'o-',label='Necessary, any feedback')
    axs[0].semilogx([r['founders'] for r in rows],[r['eraser_exposure'] for r in rows],'s-',label='Sufficient weighted policy')
    axs[0].set(xlabel='AA founder count',ylabel='Reaction exposure',title='Exposure bounds versus founder count');axs[0].legend(fontsize=8)
    for name in ['early','late','spread']:axs[1].semilogy([r['exposure'] for r in sweep],[r[name] for r in sweep],label=name.capitalize())
    axs[1].scatter([11.6],[float(F(policy['survival_AA_upper']))],marker='*',s=100,color='black',label='Certified upper bound')
    axs[1].axhline(.01,color='gray',ls=':');axs[1].set(xlabel='Reaction exposure',ylabel='One AA founder survival at time 100',title='Survival under equal-exposure treatment\nschedules');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'exposure_and_timing.png',dpi=180);fig.savefig(out/'exposure_and_timing.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    Ns=list(range(2,9));vals=[float(sites['mean_sign_endpoints'][str(n)]['intrinsic_erasure_bracket'][1]) for n in Ns]
    axs[0].plot(Ns,vals,'o-',label='Certified mean-sign endpoint');axs[0].axhline(.3,color='black',ls='--',label='Baseline + .29 actuator')
    axs[0].set(xlabel='Memory sites',ylabel='Intrinsic erasure rate',title='Mean-growth erasure thresholds versus\nmemory size');axs[0].legend(fontsize=8)
    selected=[v for v in sites['amplitude'].values() if v['N']>=5 and v['maximum']==F(29,100)]
    axs[1].bar([v['N'] for v in selected],[float(v['fully_marked_floor']) for v in selected])
    axs[1].set(xlabel='Memory sites',ylabel='Certified survival floor, one founder',ylim=(0,1),title='All-policy survival lower bounds by memory\nsize')
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'amplitude_floor.png',dpi=180);fig.savefig(out/'amplitude_floor.svg');plt.close(fig)


if __name__=='__main__':main()