#!/usr/bin/env python3
"""Ideal 1D cavity benchmark. Python 3 standard library only; no empirical data supplied."""
import argparse
import csv
import json
import math
from pathlib import Path
import random
import statistics

FIELDS = ['dataset_kind','sample_id','length_m','mode','frequency_hz','sigma_frequency_hz','sigma_length_m']
REFERENCE_SPEED = 343.2


def generate_demo(path):
    rng = random.Random(20261018)
    true_lengths = [0.15,0.18,0.22,0.26,0.31,0.37,0.44,0.53,0.64,0.76,0.90,1.05]
    with path.open('w', newline='') as f:
        w = csv.DictWriter(f, fieldnames=FIELDS); w.writeheader()
        for i, true_length in enumerate(true_lengths, 1):
            measured_length = true_length + rng.gauss(0, 0.00012)
            frequency = REFERENCE_SPEED/(2*true_length) + rng.gauss(0, 0.35)
            w.writerow(dict(dataset_kind='synthetic_demo',sample_id=f'G{i:02d}',
                length_m=f'{measured_length:.10f}',mode=1,frequency_hz=f'{frequency:.10f}',
                sigma_frequency_hz=0.35,sigma_length_m=0.00012))


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


def analyze(path):
    with path.open(newline='') as f: rows=list(csv.DictReader(f))
    if len(rows)<4: raise ValueError('At least four independent observations are required; the blank template contains no measurements.')
    if any(not row.get('dataset_kind') or not row.get('sample_id') for row in rows): raise ValueError('Every row needs dataset_kind and sample_id.')
    if len({r['sample_id'] for r in rows}) != len(rows): raise ValueError('sample_id values must be unique')
    x=[]; y=[]; sf=[]; sl=[]; lengths=[]; modes=[]
    for row in rows:
        length=number(row,'length_m',positive=True); mode=number(row,'mode',positive=True)
        if mode != int(mode): raise ValueError('mode must be a positive integer')
        lengths.append(length); modes.append(mode); x.append(mode/(2*length))
        y.append(number(row,'frequency_hz',positive=True))
        sf.append(number(row,'sigma_frequency_hz',positive=True)); sl.append(number(row,'sigma_length_m',nonnegative=True))
    if len(set(x))<3: raise ValueError('At least three distinct n/(2L) values are required')
    # Effective-variance approximation: small independent Gaussian length/frequency errors.
    # No claim of a full errors-in-variables likelihood; correlated calibration is excluded.
    speed=REFERENCE_SPEED
    for _ in range(20):
        variances=[a*a+(n*speed/(2*l*l)*b)**2 for a,b,n,l in zip(sf,sl,modes,lengths)]
        weights=[1/v for v in variances]
        speed=sum(w*a*b for w,a,b in zip(weights,x,y))/sum(w*a*a for w,a in zip(weights,x))
    variance_speed=1/sum(w*a*a for w,a in zip(weights,x))
    residuals=[b-speed*a for a,b in zip(x,y)]
    chi2=sum(w*r*r for w,r in zip(weights,residuals))
    sw=sum(weights); sx=sum(w*a for w,a in zip(weights,x)); sy=sum(w*b for w,b in zip(weights,y))
    sxx=sum(w*a*a for w,a in zip(weights,x)); sxy=sum(w*a*b for w,a,b in zip(weights,x,y))
    det=sw*sxx-sx*sx
    if det<=0: raise ValueError('Degenerate design')
    free_slope=(sw*sxy-sx*sy)/det; intercept=(sxx*sy-sx*sxy)/det
    kinds=sorted({r['dataset_kind'] for r in rows})
    return dict(project='geometry-resonance',dataset_kind=kinds,synthetic_demo=kinds==['synthetic_demo'],
        interpretation='Synthetic benchmark only; no empirical finding.' if kinds==['synthetic_demo'] else 'Analysis only; provenance, controls and calibration require review.',
        sample_count=len(rows),model='f_n = n v / (2 L); fitted intercept fixed at zero',
        fitted_speed_m_s=round(speed,9),conditional_standard_uncertainty_speed_m_s=round(math.sqrt(variance_speed),9),
        frequency_rmse_hz=round(math.sqrt(statistics.mean(r*r for r in residuals)),9),
        reduced_chi_squared=round(chi2/(len(rows)-1),9),
        reference_speed_m_s=REFERENCE_SPEED,
        fixed_reference_rmse_hz=round(math.sqrt(statistics.mean((b-REFERENCE_SPEED*a)**2 for a,b in zip(x,y))),9),
        free_intercept_diagnostic_hz=round(intercept,9),
        free_intercept_conditional_standard_uncertainty_hz=round(math.sqrt(sxx/det),9),
        free_intercept_slope_m_s=round(free_slope,9),
        max_abs_standardized_residual=round(max(abs(r)/math.sqrt(v) for r,v in zip(residuals,variances)),9),
        uncertainty_assumption='Small independent Gaussian errors; effective variance; supplied sigmas known; no correlated calibration, boundary correction, or temperature uncertainty.')


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('csv',nargs='?',type=Path,help='Input CSV; defaults to bundled synthetic-demo.csv')
    parser.add_argument('--generate-demo',action='store_true',help='Deterministically overwrite bundled synthetic-demo.csv, then analyze the selected input')
    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) 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()