#!/usr/bin/env python3
"""Retained-state benchmark with passive exponential relaxation. Synthetic demonstration, stdlib only."""
import argparse
import csv
import json
import math
from pathlib import Path
import random
import statistics

FIELDS=['dataset_kind','pair_id','condition','delay_s','response_au','sigma_response_au']


def generate_demo(path):
    rng=random.Random(20261019)
    delays=[0,2,4,7,10,14,20,28]
    with path.open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=FIELDS);w.writeheader()
        for di,t in enumerate(delays):
            for replicate in range(6):
                pair=f'M{di+1:02d}-{replicate+1:02d}'
                # Each pair is a new reset/retrain/probe run, not serial probing of one specimen.
                shared_offset=0.40+0.003*t+rng.gauss(0,0.02)
                for condition in ['trained','untrained']:
                    retained=1.5*math.exp(-t/12.0) if condition=='trained' else 0.0
                    response=shared_offset+retained+rng.gauss(0,0.015)
                    w.writerow(dict(dataset_kind='synthetic_demo',pair_id=pair,condition=condition,
                        delay_s=t,response_au=f'{response:.10f}',sigma_response_au=0.015))


def number(row,key,positive=False,nonnegative=False):
    try:v=float(row[key])
    except (KeyError, ValueError, TypeError):raise ValueError(f'{row.get("pair_id","row")}: {key} must be numeric')
    if not math.isfinite(v):raise ValueError(f'{key} must be finite')
    if positive and v<=0:raise ValueError(f'{key} must be positive')
    if nonnegative and v<0:raise ValueError(f'{key} must be nonnegative')
    return v


def analyze(path):
    with path.open(newline='') as f:rows=list(csv.DictReader(f))
    if not rows:raise ValueError('No measurements: fill a copy of blank-measurements.csv; do not relabel the synthetic demonstration.')
    pairs={}
    for row in rows:
        if not row.get('dataset_kind') or not row.get('pair_id'):raise ValueError('Every row needs dataset_kind and pair_id')
        condition=row.get('condition')
        if condition not in ['trained','untrained']:raise ValueError('condition must be trained or untrained')
        pair=pairs.setdefault(row['pair_id'],{})
        if condition in pair:raise ValueError('Each pair_id must have exactly one row per condition')
        pair[condition]=(number(row,'delay_s',nonnegative=True),number(row,'response_au'),number(row,'sigma_response_au',positive=True))
    groups={}
    for pair_id,pair in pairs.items():
        if set(pair)!={'trained','untrained'}:raise ValueError(f'{pair_id}: missing matched condition')
        trained=pair['trained'];untrained=pair['untrained']
        if trained[0]!=untrained[0]:raise ValueError('Matched pair delays must be equal')
        groups.setdefault(trained[0],[]).append((trained[1]-untrained[1],trained[2]**2+untrained[2]**2))
    if len(groups)<4:raise ValueError('At least four distinct delays are required')
    times=sorted(groups);differences=[];variances=[];counts=[]
    for t in times:
        values=groups[t]
        if len(values)<2:raise ValueError('At least two independent pairs per delay are required')
        # Known-noise weighted mean. Shared pair offsets cancel by construction in demo.
        total_weight=sum(1/v for _,v in values)
        mean=sum(d/v for d,v in values)/total_weight
        if mean<=0:raise ValueError('All mean paired differences must be positive for this log-linear method. Do not drop nonpositive means; use a justified original-scale model instead.')
        differences.append(mean);variances.append(1/total_weight);counts.append(len(values))
    logs=[math.log(v) for v in differences]
    # Delta-method variance of log(mean difference): approximate when noise is small.
    logvars=[v/(d*d) for v,d in zip(variances,differences)];weights=[1/v for v in logvars]
    sw=sum(weights);sx=sum(w*t for w,t in zip(weights,times));sy=sum(w*y for w,y in zip(weights,logs))
    sxx=sum(w*t*t for w,t in zip(weights,times));sxy=sum(w*t*y for w,t,y in zip(weights,times,logs))
    det=sw*sxx-sx*sx
    if det<=0:raise ValueError('Degenerate delay design')
    slope=(sw*sxy-sx*sy)/det;intercept=(sxx*sy-sx*sxy)/det
    if slope>=0:raise ValueError('The fitted slope is nonnegative; the proposed decaying exponential is not supported by this fit.')
    tau=-1/slope;amplitude=math.exp(intercept)
    predictions=[amplitude*math.exp(-t/tau) for t in times]
    residuals=[d-p for d,p in zip(differences,predictions)]
    kinds=sorted({r['dataset_kind'] for r in rows})
    return dict(project='retained-physical-memory',dataset_kind=kinds,synthetic_demo=kinds==['synthetic_demo'],
        interpretation='Synthetic passive-relaxation benchmark; no empirical learning or universal-memory finding.' if kinds==['synthetic_demo'] else 'Analysis only; paired controls and physical substrate require independent review.',
        paired_runs=len(pairs),delay_count=len(times),fit='Weighted least squares of log(mean trained-minus-untrained response) = log(A) - t/tau',
        fitted_amplitude_au=round(amplitude,9),conditional_standard_uncertainty_amplitude_au=round(amplitude*math.sqrt(sxx/det),9),
        fitted_tau_s=round(tau,9),conditional_standard_uncertainty_tau_s=round(math.sqrt(sw/det)/(slope*slope),9),
        exponential_mean_response_rmse_au=round(math.sqrt(statistics.mean(r*r for r in residuals)),9),
        zero_difference_baseline_rmse_au=round(math.sqrt(statistics.mean(d*d for d in differences)),9),
        minimum_mean_paired_difference_au=round(min(differences),9),
        reduced_chi_squared_log_fit=round(sum((y-intercept-slope*t)**2/v for y,t,v in zip(logs,times,logvars))/(len(times)-2),9),
        delay_results=[dict(delay_s=t,pairs=n,mean_difference_au=round(d,9),standard_uncertainty_au=round(math.sqrt(v),9),fitted_response_au=round(p,9)) for t,n,d,v,p in zip(times,counts,differences,variances,predictions)],
        uncertainty_assumption='Known independent response noise; matched shared offsets cancel; small-noise log delta method; no time error or condition-specific offset modeled.')


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('csv',nargs='?',type=Path)
    parser.add_argument('--generate-demo',action='store_true',help='Regenerate the bundled deterministic synthetic CSV')
    args=parser.parse_args();here=Path(__file__).resolve().parent
    if args.generate_demo:generate_demo(here/'synthetic-demo.csv')
    try:result=analyze(args.csv or here/'synthetic-demo.csv')
    except (ValueError,OSError,OverflowError) as exc:parser.exit(2,f'Error: {exc}\n')
    print(json.dumps(result,indent=2,sort_keys=True,allow_nan=False))

if __name__=='__main__':main()