#!/usr/bin/env python3
"""Dimensionless scaling comparison with a prespecified scale break and held-out data. Stdlib only."""
import argparse
import csv
import json
import math
from pathlib import Path
import random
import statistics

FIELDS=['dataset_kind','sample_id','x_ratio','y_ratio','split','sigma_log_y']
BREAK_X=20.0  # Chosen before demo generation. It is not optimized on held-out data.


def generate_demo(path):
    rng=random.Random(20261020)
    with path.open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=FIELDS);w.writeheader()
        for i in range(37):
            logx=5*i/36;x=math.exp(logx)
            logy=math.log(1.7)+0.60*logx+0.55*max(0,logx-math.log(BREAK_X))+rng.gauss(0,0.04)
            w.writerow(dict(dataset_kind='synthetic_demo',sample_id=f'S{i+1:02d}',
                x_ratio=f'{x:.12f}',y_ratio=f'{math.exp(logy):.12f}',split='test' if i%4==2 else 'train',sigma_log_y=0.04))


def solve(matrix,vector):
    n=len(vector);a=[list(row)+[value] for row,value in zip(matrix,vector)]
    for j in range(n):
        pivot=max(range(j,n),key=lambda k:abs(a[k][j]))
        if abs(a[pivot][j])<1e-12:raise ValueError('Singular model design; supply variation on both sides of the declared break')
        a[j],a[pivot]=a[pivot],a[j];scale=a[j][j];a[j]=[v/scale for v in a[j]]
        for k in range(n):
            if k==j:continue
            scale=a[k][j];a[k]=[v-scale*w for v,w in zip(a[k],a[j])]
    return [a[i][-1] for i in range(n)]


def features(x,broken):
    logx=math.log(x)
    return [1,logx,max(0,logx-math.log(BREAK_X))] if broken else [1,logx]


def fit(rows,broken):
    dim=3 if broken else 2;matrix=[[0.0]*dim for _ in range(dim)];vector=[0.0]*dim
    for x,y,sigma in rows:
        z=features(x,broken);weight=1/(sigma*sigma);logy=math.log(y)
        for j in range(dim):
            vector[j]+=weight*z[j]*logy
            for k in range(dim):matrix[j][k]+=weight*z[j]*z[k]
    return solve(matrix,vector)


def mse(rows,coef,broken):
    return statistics.mean((math.log(y)-sum(a*b for a,b in zip(coef,features(x,broken))))**2 for x,y,_ in rows)


def analyze(path):
    with path.open(newline='') as f:raw=list(csv.DictReader(f))
    if not raw:raise ValueError('No measurements: blank-measurements.csv contains only a schema')
    if len({r.get('sample_id') for r in raw})!=len(raw):raise ValueError('sample_id values must be unique')
    groups={'train':[],'test':[]}
    for row in raw:
        if not row.get('dataset_kind') or not row.get('sample_id'):raise ValueError('Every row needs dataset_kind and sample_id')
        split=row.get('split')
        if split not in groups:raise ValueError('split must be train or test, assigned before fitting')
        try:values=tuple(float(row[k]) for k in ['x_ratio','y_ratio','sigma_log_y'])
        except (ValueError,KeyError,TypeError):raise ValueError('x_ratio, y_ratio and sigma_log_y must be numeric')
        if any(not math.isfinite(v) or v<=0 for v in values):raise ValueError('Ratios and log-noise standard deviations must be positive and finite')
        groups[split].append(values)
    train=groups['train'];test=groups['test']
    if len(train)<8 or len(test)<4:raise ValueError('At least eight training and four held-out observations required')
    for name,rows in groups.items():
        if not any(x<BREAK_X for x,_,_ in rows) or not any(x>BREAK_X for x,_,_ in rows):raise ValueError(f'{name} must cover both sides of prespecified break x={BREAK_X}')
    single=fit(train,False);broken=fit(train,True)
    single_mse=mse(test,single,False);broken_mse=mse(test,broken,True)
    error_differences=[]
    for x,y,_ in test:
        logy=math.log(y)
        a=(logy-sum(c*z for c,z in zip(single,features(x,False))))**2
        b=(logy-sum(c*z for c,z in zip(broken,features(x,True))))**2
        error_differences.append(a-b)
    kinds=sorted({r['dataset_kind'] for r in raw})
    return dict(project='cross-scale-relationships',dataset_kind=kinds,synthetic_demo=kinds==['synthetic_demo'],
        interpretation='Synthetic scale-break demonstration; no measured universal scaling finding.' if kinds==['synthetic_demo'] else 'Descriptive held-out comparison only; split integrity and measurement model require review.',
        training_observations=len(train),held_out_observations=len(test),prespecified_break_x_ratio=BREAK_X,
        power_law=dict(amplitude=round(math.exp(single[0]),9),exponent=round(single[1],9),parameter_count=2,
            training_log_rmse=round(math.sqrt(mse(train,single,False)),9),held_out_log_rmse=round(math.sqrt(single_mse),9)),
        scale_break_baseline=dict(amplitude_below_break=round(math.exp(broken[0]),9),exponent_below_break=round(broken[1],9),exponent_above_break=round(broken[1]+broken[2],9),parameter_count=3,
            training_log_rmse=round(math.sqrt(mse(train,broken,True)),9),held_out_log_rmse=round(math.sqrt(broken_mse),9)),
        held_out_log_mse_difference_power_minus_break=round(single_mse-broken_mse,9),
        conditional_standard_error_of_mse_difference=round(statistics.stdev(error_differences)/math.sqrt(len(test)),9),
        comparison_note='Positive MSE difference favors the prespecified break baseline on this held-out sample. No p-value, universal law, or formal model-selection verdict is inferred.',
        uncertainty_assumption='Positive dimensionless ratios; independent log-scale observation noise with supplied sigma; x treated as exact; standard error is conditional on trained fits and independent test rows.')


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('csv',nargs='?',type=Path)
    parser.add_argument('--generate-demo',action='store_true',help='Regenerate 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()