"""Sharp fresh-conversion certificate with material histories and design sweeps."""
from fractions import Fraction as F
from pathlib import Path
from dataclasses import asdict
import argparse,csv,hashlib,itertools,json,platform,random
from material import Window,Recovery,History,attaining_history,amount_error
from certificate import DirectReserve,UptakeBound,UnboundedReserve,MaterialCertificate,Observation,report

# EDITABLE INPUTS: micromol product equivalents per ORIGINAL 10 mL aliquot.
ORIGINAL_VOLUME_ML='10'
INITIAL_STOCK_BOUND='10'
PREWASH_RESERVE_BOUND='2'
INITIAL_RESERVE_BOUND='2'          # distinct uptake-premise input
CUMULATIVE_UPTAKE_BOUND='0.1'      # gross uptake, not net disappearance
MOBILE_RETENTION='0.05'
RESERVE_RETENTION='0.9'
RECOVERY_INPUT_BOUND='0.2'
REPORTED_COLLECTIONS=('6','4')
AMOUNT_ERROR_BOUNDS=('0.2','0.2')
FRESH_TARGET='1.6'                # paper's principal target, not legacy 1.5
FRESH_UPPER=None                   # separately calibrated; never inferred from low output
ADDITIONAL_ACTIVITY_RETAINED='1'   # relative to already-sampled unwashed comparator
RESERVE_SLACK='0'
WINDOW_DURATIONS_HOURS=('1','1')   # illustrative only: average, not steady-state rate
SEED=62092026
RANDOM_REPLAYS=400
MANUSCRIPT_SHA256='e6e69d291896c7e93c96d47920557c1a1eb16bd3f6c15c83c5697207f75f45df'


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)
    def dump(name,data):
        (out/name).write_text(json.dumps(data,indent=2,default=lambda x:{'exact':str(x),'decimal':float(x)} if isinstance(x,F) else str(x))+'\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)
    B,J,R,T,e,s,H=map(F,(INITIAL_STOCK_BOUND,PREWASH_RESERVE_BOUND,INITIAL_RESERVE_BOUND,CUMULATIVE_UPTAKE_BOUND,MOBILE_RETENTION,RESERVE_RETENTION,RECOVERY_INPUT_BOUND))
    observation=Observation(*map(F,REPORTED_COLLECTIONS),*map(F,AMOUNT_ERROR_BOUNDS));threshold=F(FRESH_TARGET)
    cert=lambda premise,ee=e:MaterialCertificate(B,ee,s,H,premise)
    direct=cert(DirectReserve(J));uptake=cert(UptakeBound(R,T));unbounded=cert(UnboundedReserve())
    decisions={name:report(observation,[c],threshold,FRESH_UPPER) for name,c in [('direct',direct),('uptake',uptake),('no-reserve-information',unbounded)]}
    decisions['combined']=report(observation,[direct,uptake],threshold,FRESH_UPPER)
    witnesses={};q=observation.lower_collections()
    if q is not None:
        for name,c in [('direct',direct),('uptake',uptake),('no-reserve-information',unbounded)]:
            h=c.witness(*q);witnesses[name]=dict(history=asdict(h),fresh=h.fresh,valid_for_premise=c.valid_for(h),events=h.ledger())
        # Attainment from any split applies to the DIRECT reserve class.
        splits=[direct.witness(*q,v).fresh for v in [F(0),B/2,B]]
        unbounded_fresh=direct.witness(*q).extra_unobserved_formation(100)
        witnesses['same-collections-extra-fresh-loss']=dict(history=asdict(unbounded_fresh),fresh=unbounded_fresh.fresh,
            valid_for_premise=direct.valid_for(unbounded_fresh),note='Low collection does not bound fresh input above: the extra material may be lost.')
    else:splits=[]
    # Paper reference histories are fixed explicit examples, separate from edits.
    first=Window(F(8),F(2),collection=F(6));washed=History(first,Window(F(1,10),F(9,5),fresh=F(21,10),release=F(39,10),collection=F(4)),Recovery())
    unwashed=History(first,Window(F(9,5),F(9,5),fresh=F(21,10),release=F(39,10),collection=F(57,10)),Recovery(F(9,10),F(9,10)))
    stored=History(Window(F(8),F(2),uptake=F(2),collection=F(6)),Window(F(0),F(18,5),release=F(18,5),collection=F(18,5)),Recovery())
    ambiguity=History(Window(F(8),F(2),uptake=F(2),collection=F(29,5)),Window(F(1,100),F(19,5),release=F(19,5),collection=F(381,100)),Recovery(reserve_input=F(1,5)))
    paper_cert=lambda j,ee=F(1,20),hh=F(1,5):MaterialCertificate(F(10),ee,F(9,10),hh,DirectReserve(j))
    paper_obs=Observation(F(6),F(4),F(1,5),F(1,5))
    paper_w=report(paper_obs,[paper_cert(F(2))],F(8,5));paper_u=report(Observation(F(6),F(57,10),F(1,5),F(1,5)),[paper_cert(F(2),F(9,10))],F(8,5))
    counter=dict(storage=asdict(stored),true_fresh=stored.fresh,
        wrong_initial_substitution=paper_cert(F(2),hh=F(0)).evaluate(6,F(18,5))['lower'],
        wrong_premise_rejected=not paper_cert(F(2),hh=F(0)).valid_for(stored),
        correct_direct=paper_cert(F(4),hh=F(0)).evaluate(6,F(18,5))['lower'],
        correct_uptake=MaterialCertificate(F(10),F(1,20),F(9,10),F(0),UptakeBound(F(2),F(2))).evaluate(6,F(18,5))['lower'],
        same_record_zero_fresh_history=asdict(ambiguity),within_output_envelopes=abs(ambiguity.first.collection-6)<=F(1,5) and abs(ambiguity.second.collection-4)<=F(1,5))
    reserve_rows=[];activity_rows=[]
    for n in range(101):
        j=F(n,20);r=paper_cert(j).evaluate(F(29,5),F(19,5))
        reserve_rows.append([float(j),float(r['lower']),';'.join(r['active'])])
    for delta in [F(0),F(1,20),F(1,10),F(1,5)]:
        for n in range(101):
            a=F(n,100);second=F(17,10)+F(21,10)*a
            r=paper_cert(2+delta).evaluate(F(29,5),second)
            assert r['lower']==max(F(0),F(21,10)*a-F(41,100)-F(17,20)*delta)
            activity_rows.append([float(a),float(delta),float(r['lower']),r['lower']>F(38,25),r['lower']>=F(8,5)])
    a,delta=map(F,(ADDITIONAL_ACTIVITY_RETAINED,RESERVE_SLACK))
    if not 0<=a<=1 or delta<0:raise ValueError('Activity must be in [0,1]; reserve slack must be nonnegative.')
    design=dict(activity=a,slack=delta,certificate=paper_cert(2+delta).evaluate(F(29,5),F(17,10)+F(21,10)*a)['lower'],
        strict_improvement_activity_boundary=F(193,210)+F(17,42)*delta,
        target_activity_boundary=F(67,70)+F(17,42)*delta,
        margin_budget=F(17,100)-F(17,20)*delta-F(21,10)*(1-a),
        reserve_boundaries=dict(positive_strict=F(339,85),target_inclusive=F(179,85),improvement_strict=F(11,5)),
        uptake_boundaries=dict(target_inclusive=F(9,85),improvement_strict=F(1,5)))
    # Matched source: second collection changes with mobile retention.
    matched=[]
    for ee in [F(0),F(1,20),F(1,2),F(9,10)]:
        Q2=ee*2+F(9,10)*2+F(21,10)
        nominal=paper_cert(F(2),ee).evaluate(6,Q2)['lower']
        robust=paper_cert(F(2),ee).evaluate(F(29,5),Q2-F(1,5))['lower']
        matched.append(dict(e=ee,second_collection=Q2,nominal=nominal,robust=robust))
    # Exact endpoint and random schedule replay, not a replacement for the proof.
    replay_count=0
    for q1,q2,stock,reserve,input_ in itertools.product([F(0),F(1),F(3)],repeat=5):
        for ee,ss in [(F(0),F(0)),(F(0),F(1)),(F(1,2),F(1,2)),(F(1),F(1))]:
            c=MaterialCertificate(stock,ee,ss,input_,DirectReserve(reserve))
            for split in [F(0),stock/2,stock]:
                h=c.witness(q1,q2,split);h.ledger();assert c.valid_for(h);replay_count+=1
    rng=random.Random(SEED)
    for _ in range(RANDOM_REPLAYS):
        draw=lambda:F(rng.randrange(101),10)
        stock,q1,q2,rb,tb,inp=draw(),draw(),draw(),draw(),draw(),draw();ee=F(rng.randrange(11),10);ss=F(rng.randrange(int(ee*10),11),10)
        c=MaterialCertificate(stock,ee,ss,inp,UptakeBound(rb,tb));h=c.witness(q1,q2);h.ledger();assert c.valid_for(h);replay_count+=1
    durations=tuple(map(F,WINDOW_DURATIONS_HOURS))
    if min(durations)<=0:raise ValueError('Positive recorded window durations required.')
    summary=dict(configured_decisions=decisions,all_initial_splits_attain=splits,original_volume_ml=F(ORIGINAL_VOLUME_ML),
        total_window_average_lower=None if decisions['direct']['lower'] is None else decisions['direct']['lower']/sum(durations),
        paper_reference=dict(washed=paper_w,unwashed=paper_u,washed_history=asdict(washed),unwashed_history=asdict(unwashed),realized_fresh=washed.fresh),
        exact_replay_count=replay_count,matched_retention=matched,
        measurement_example=dict(volume_ml=F(2),concentration_umol_per_ml=F(3),volume_error=F(1,20),concentration_error=F(1,10),
             amount_error=amount_error(2,3,F(1,20),F(1,10))),
        interpretation=['Credits fresh entry into declared inventory, which may include intermediates; not substrate-specific atom attribution.',
        'No inference of viability, growth, abundance, future function or a steady-state rate.',
        'Initial stock must include all old convertible reserves. Challenge substrate must be assigned once, outside before fresh entry or inside old stock.',
        'Collection physically leaves the accounting boundary once; wash waste is separate.',
        'All inputs are synthetic. Lean is not rerun; exact rational computation checks the source algebra and constructed histories.'])
    dump('summary.json',summary);dump('attaining_histories.json',witnesses);dump('counterexamples.json',counter);dump('design.json',design)
    table('reserve_sweep.csv',['reserve_ceiling','lower','active_expressions'],reserve_rows)
    table('activity_sweep.csv',['activity','reserve_slack','lower','strictly_beats_unwashed','meets_target'],activity_rows)
    plot(out,reserve_rows,activity_rows,washed,ambiguity)
    print(f'Configured direct-reserve report: {decisions["direct"]["verdict"]}; lower = {decisions["direct"]["lower"]}')
    print(f'{replay_count} exact endpoint/random histories replayed. Paper reference: washed 1.69, unwashed 1.52.',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={n:digest(here/n) for n in ['material.py','certificate.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,reserve,activity,positive,zero):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    import numpy as np
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    axs[0].plot([r[0] for r in reserve],[r[1] for r in reserve]);axs[0].axhline(1.6,color='black',ls=':',label='Fresh target 1.6');axs[0].axhline(1.52,color='gray',ls='--',label='Unwashed certificate 1.52')
    axs[0].set(xlabel='Pre-wash reserve ceiling (micromol equivalents)',ylabel='Certified fresh lower amount',title='Fresh-conversion bound versus initial\nreserve allowance');axs[0].legend(fontsize=8)
    for d in [0,.05,.1,.2]:
        rows=[r for r in activity if r[1]==d];axs[1].plot([r[0] for r in rows],[r[2] for r in rows],label=f'Reserve slack {d}')
    axs[1].axhline(1.52,color='gray',ls='--');axs[1].axhline(1.6,color='black',ls=':');axs[1].set(xlabel='Additional productive activity retained',ylabel='Washed fresh lower amount',title='Fresh-conversion bound versus retained\nactivity');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'design.png',dpi=180);fig.savefig(out/'design.svg');plt.close(fig)
    fig,ax=plt.subplots(figsize=(10,4),layout='constrained');x=np.arange(4)
    for shift,h,label in [(-.18,positive,'Positive reference history'),(.18,zero,'Zero-fresh alternative')]:
        values=[h.first.collection,h.second.collection,h.first.final_reserve,h.fresh]
        ax.bar(x+shift,list(map(float,values)),width=.36,label=label)
    for i,v in [(0,6),(1,4)]:ax.errorbar(i,v,yerr=.2,fmt='none',color='black',capsize=12,lw=2)
    ax.set(xticks=x,xticklabels=['First collection','Second collection','Pre-wash reserve','Fresh total'],ylabel='Micromol equivalents per original aliquot',title='Two material histories compatible with the\ncollection data');ax.legend();ax.grid(axis='y',alpha=.2)
    fig.savefig(out/'ambiguity.png',dpi=180);fig.savefig(out/'ambiguity.svg');plt.close(fig)


if __name__=='__main__':main()
