#!/usr/bin/env python3
"""Analyse preregistered paired block means without treating sensor samples as replicates.

Usage: python analyse.py measurements.jsonl run-metadata.json --output result.json
Requires numpy and scipy. Empty data returns NO_EVIDENCE. Synthetic fixtures can
exercise the pipeline but can never produce a physical-evidence classification.
"""
from pathlib import Path
import argparse,json,math,hashlib
import numpy as np
from scipy.stats import t as student
HERE=Path(__file__).resolve().parent

def analyse(data,metadata,protocol=None):
    p=protocol or json.loads((HERE/'protocol.json').read_text());rows=[json.loads(s) for s in Path(data).read_text().splitlines() if s.strip()];meta=json.loads(Path(metadata).read_text());classification=meta.get('classification','NO_MEASUREMENTS')
    report={'classification':classification,'physical_evidence':False,'status':'NO_EVIDENCE','raw_records':len(rows),'data_sha256':hashlib.sha256(Path(data).read_bytes()).hexdigest(),'metadata_sha256':hashlib.sha256(Path(metadata).read_bytes()).hexdigest(),'protocol_sha256':hashlib.sha256(json.dumps(p,sort_keys=True).encode()).hexdigest(),'analysis_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),'limitations':['Consistency and statistical checks do not authenticate operator-entered measurements.','This tests one low-pressure channel; it cannot verify whole-car performance or net energy benefit.']}
    if not rows:return report|{'reason':'No measurement records supplied.'}
    if classification not in ['PHYSICAL_MEASUREMENTS','SYNTHETIC_PIPELINE_TEST']:return report|{'reason':'Dataset classification is not physical measurements or a labelled pipeline test.'}
    required=['operator','article_id','flow_calibration_id','pressure_calibration_id','as_built_record','leak_control_record','warm_warm_control_record']
    missing=[k for k in required if not meta.get(k)]
    if len(meta.get('temperature_calibration_ids',[]))<4:missing.append('four temperature calibration identifiers')
    if missing:return report|{'status':'INVALID_METADATA','missing':missing}
    runs={};invalid=[];plan={r['run_id']:r for r in p['randomisation']['runs']};target=p['boundary_targets'];q=p['run_quality'];expected_order=[]
    for r in rows:
        rid=r.get('run_id');expected=plan.get(rid)
        if not expected:invalid.append({'run_id':rid,'reason':'Run not in frozen primary schedule.'});continue
        if r.get('pair_id')!=expected['pair_id'] or r.get('condition')!=expected['condition'] or r.get('record_order')!=expected['order']:invalid.append({'run_id':rid,'reason':'Pair, condition or order differs from frozen schedule.'});continue
        if not expected_order or expected_order[-1]!=r['record_order']:expected_order.append(r['record_order'])
        runs.setdefault(rid,[]).append(r)
    if expected_order!=sorted(expected_order) or len(expected_order)!=len(set(expected_order)):invalid.append({'reason':'Primary run records are interleaved or out of frozen order.'})
    summary={}
    for rid,rr in runs.items():
        reasons=[]
        try:
            times=np.array([r['elapsed_s'] for r in rr],float);flow=np.array([r['mass_flow_g_s'] for r in rr],float);head=np.array([r['head_Pa'] for r in rr],float);wall=np.array([r['wall_K'] for r in rr],float);inlet=np.array([r['inlet_K'] for r in rr],float);dew=np.array([r['dewpoint_K'] for r in rr],float);power=np.array([[r[k] for k in ['blower_W','cooling_W','pump_W','trim_heater_W']] for r in rr],float)
            numeric=[times,flow,head,wall,inlet,dew,power]
            if not all(np.all(np.isfinite(a)) for a in numeric):reasons.append('Non-finite numeric data.')
            if len(rr)<q['minimum_samples'] or times[-1]-times[0]<q['minimum_duration_s']:reasons.append('Insufficient measurement duration or sample count.')
            if np.any(np.diff(times)<=0) or np.max(np.diff(times))>q['maximum_gap_between_samples_s']:reasons.append('Invalid or interrupted sample timing.')
            if wall.shape!=(len(rr),4):reasons.append('Four wall probes are required.')
            if min(flow)<=0 or min(head)<=0 or min(power.ravel())<0:reasons.append('Nonphysical flow, head or input-power value.')
            if max(abs(head-target['head_Pa']))>target['maximum_head_deviation_Pa']:reasons.append('Pressure head outside target.')
            if max(abs(inlet-target['inlet_temperature_K']))>target['inlet_tolerance_K']+1e-9:reasons.append('Inlet temperature outside target.')
            tw=target[rr[0]['condition']+'_wall_K']
            if max(abs(wall.mean(axis=1)-tw))>target['wall_mean_tolerance_K']+1e-9:reasons.append('Mean wall temperature outside target.')
            if max(np.ptp(wall,axis=1))>target['maximum_wall_probe_spread_K']+1e-9:reasons.append('Wall temperature spatial spread outside target.')
            if max(dew)>target['maximum_dewpoint_K']:reasons.append('Dry-air dewpoint outside target.')
            if min(r['prestable_s'] for r in rr)<q['stabilisation_minimum_s'] or max(abs(r['prestable_wall_drift_K_min']) for r in rr)>q['wall_stability_before_measurement_K_per_min']:reasons.append('Thermal prestabilisation criterion not met.')
            if any(r.get('quality_note','').strip() for r in rr):reasons.append('Operator marked an exception; retain and repeat the entire pair.')
            split=len(flow)//2
            if abs(np.mean(flow[:split])-np.mean(flow[split:]))/np.mean(flow)>q['maximum_flow_half_window_change_fraction']:reasons.append('Flow drift exceeds the steady-window criterion.')
            if reasons:invalid.append({'run_id':rid,'reasons':reasons});continue
            # The run mean is the independent unit, not individual sensor samples.
            summary[rid]={'pair_id':rr[0]['pair_id'],'condition':rr[0]['condition'],'mass_flow_g_s':float(np.mean(flow)),'head_Pa':float(np.mean(head)),'wall_mean_K':wall.mean(axis=0).tolist(),'inlet_K':float(np.mean(inlet)),'electrical_input_W':float(np.sum(np.mean(power,axis=0))),'component_input_W':dict(zip(['blower','cooling','pump','trim_heater'],map(float,np.mean(power,axis=0)))),'sample_count':len(rr),'duration_s':float(times[-1]-times[0])}
        except (KeyError,ValueError,TypeError,IndexError) as e:invalid.append({'run_id':rid,'reason':'Malformed or incomplete instrument record: '+str(e)})
    pairs=[]
    for i in range(1,13):
        pid=f'P{i:02d}';both=[r for r in summary.values() if r['pair_id']==pid]
        if len(both)!=2:continue
        w=next(r for r in both if r['condition']=='warm');c=next(r for r in both if r['condition']=='cold');ratio=c['mass_flow_g_s']/w['mass_flow_g_s'];ew=max(meta['flow_bound_relative'],meta['flow_bound_full_scale_fraction']*meta['flow_full_scale_g_s']/w['mass_flow_g_s']);ec=max(meta['flow_bound_relative'],meta['flow_bound_full_scale_fraction']*meta['flow_full_scale_g_s']/c['mass_flow_g_s']);assert 0<=ew<1 and 0<=ec<1
        uw=meta['pressure_bound_relative']*w['head_Pa']+meta['pressure_bound_absolute_Pa'];uc=meta['pressure_bound_relative']*c['head_Pa']+meta['pressure_bound_absolute_Pa'];pressure_penalty=max(0.,math.log((c['head_Pa']+uc)/(w['head_Pa']-uw)));flow_penalty=math.log((1+ew)/(1-ec));pairs.append({'pair_id':pid,'log_ratio':math.log(ratio),'raw_gain_percent':100*(ratio-1),'flow_systematic_log_bound':flow_penalty,'pressure_systematic_log_bound':pressure_penalty})
    report|={'valid_runs':len(summary),'valid_pairs':len(pairs),'run_means':summary,'pairs':pairs,'invalid_records_or_runs':invalid,'pressure_bound_assumption':'Conservative local mass-flow pressure elasticity <=1 for this low-Mach laminar coupon. Actual head mismatch and transducer bounds are included; no favourable head normalisation is applied.'}
    if invalid or len(pairs)!=12:return report|{'status':'INCOMPLETE_OR_REJECTED','reason':'All twelve preregistered pairs must pass. Rejected records are retained in this report.'}
    logs=np.array([r['log_ratio'] for r in pairs]);mean=float(logs.mean());half=float(student.ppf(.975,len(logs)-1)*logs.std(ddof=1)/math.sqrt(len(logs)));systematic=max(r['flow_systematic_log_bound']+r['pressure_systematic_log_bound'] for r in pairs);lower=mean-half-systematic;passed=lower>0
    report|={'mean_gain_percent':100*math.expm1(mean),'paired_statistical_95_interval_percent':[100*math.expm1(mean-half),100*math.expm1(mean+half)],'additional_worst_systematic_log_penalty':systematic,'systematic_penalised_lower_bound_percent':100*math.expm1(lower),'decision':'POSITIVE_EFFECT_SUPPORTED' if passed else 'POSITIVE_EFFECT_NOT_ESTABLISHED','status':'ANALYSED','physical_evidence':bool(passed and classification=='PHYSICAL_MEASUREMENTS'),'inference':'The statistical interval concerns variation between thermal blocks. The additional instrument penalty is a deterministic stated bound, not a claimed 95% metrology coverage interval.'}
    if classification=='SYNTHETIC_PIPELINE_TEST':report['status']='SYNTHETIC_PIPELINE_TEST_ONLY';report['physical_evidence']=False
    return report

def main():
    a=argparse.ArgumentParser();a.add_argument('data',type=Path);a.add_argument('metadata',type=Path);a.add_argument('--output',type=Path,required=True);v=a.parse_args();result=analyse(v.data,v.metadata);v.output.write_text(json.dumps(result,indent=2)+'\n');print(json.dumps({k:result.get(k) for k in ['status','physical_evidence','valid_pairs','mean_gain_percent','systematic_penalised_lower_bound_percent']},indent=2))
if __name__=='__main__':main()