#!/usr/bin/env python3 """Preflight or execute the sealed nominal-calendar storage diagnostic (stdlib only). No data are downloaded. 'run' requires an internal code-review approval tied to the exact plan and runner hashes, plus a passing preflight. Raw CSVs are read-only. """ from pathlib import Path from datetime import datetime,date,timedelta,time,timezone from zoneinfo import ZoneInfo from fractions import Fraction from decimal import Decimal from collections import Counter,defaultdict import argparse,calendar,csv,hashlib,json,math,platform,random,statistics,sys,time as clock,unicodedata PARIS=ZoneInfo('Europe/Paris') FIELDS=['date_heure','date','heure','nature','consommation','eolien','solaire'] EMPTY={'','null','none','nan','nd','na','n/a'} def sha(b):return hashlib.sha256(b).hexdigest() def write_json(p,v):p.write_text(json.dumps(v,indent=2)+'\n') def write_csv(p,rows): with p.open('w',newline='',encoding='utf-8') as f: w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows) def norm(s):return ''.join(c for c in unicodedata.normalize('NFKD',s.lower()) if not unicodedata.combining(c)) def exact_integer(s): x=Decimal(s) if not x.is_finite() or x!=x.to_integral_value():raise ValueError('not finite integral MW') return int(x) def load_nominal(directory): plan_bytes=(directory/'plan-v2.json').read_bytes();plan=json.loads(plan_bytes) if plan['version']!='2.0':raise ValueError('This runner implements only the explicitly nominal v2 plan') meta_bytes=(directory/'raw/dataset-metadata.json').read_bytes();meta=json.loads(meta_bytes) schema={x['name']:x for x in meta['fields']};errors=[];years=[];data={};inputs=[] for field in ['consommation','eolien','solaire']: if schema.get(field,{}).get('type')!='int' or '(MW)' not in schema.get(field,{}).get('label',''): errors.append({'scope':'metadata','error':'unexpected units/type','field':field}) for year in plan['data']['years']: p=directory/'raw'/f'eco2mix-{year}-selected.csv';raw=p.read_bytes();receipt=json.loads((directory/'raw'/f'eco2mix-{year}-receipt.json').read_text()) if sha(raw)!=receipt['sha256']:errors.append({'year':year,'error':'raw hash differs'}) inputs.append({'year':year,'file':str(p),'bytes':len(raw),'sha256':sha(raw)}) actual={};quarter=set();year_errors=[];signs=Counter();minimum={'wind':None,'solar':None};nature=Counter();total=0 with p.open(encoding='utf-8-sig',newline='') as f: reader=csv.DictReader(f,delimiter=';') if set(reader.fieldnames or [])!=set(FIELDS):year_errors.append({'error':'wrong selected columns'}) for line,row in enumerate(reader,2): total+=1 try: day=date.fromisoformat(row['date']);tm=time.fromisoformat(row['heure']) if day.year!=year or tm.second or tm.microsecond or tm.minute not in [0,15,30,45]:raise ValueError('invalid nominal date/hour label') except Exception as exc:year_errors.append({'line':line,'error':'invalid source label','detail':str(exc)});continue key=(day,tm.hour,tm.minute);values=[row.get(x,'').strip() for x in ['consommation','eolien','solaire']];missing=[v.lower() in EMPTY for v in values] if tm.minute in [15,45]: if not all(missing):year_errors.append({'line':line,'error':'nonempty actual at forecast-only quarter label'}) if key in quarter:year_errors.append({'line':line,'error':'duplicate forecast-only source label'}) quarter.add(key);continue if key in actual:year_errors.append({'line':line,'error':'duplicate actual source label'});continue if any(missing):year_errors.append({'line':line,'error':'missing actual field','label':str(key)});continue try:load,wind,solar=map(exact_integer,values) except Exception as exc:year_errors.append({'line':line,'error':'invalid integer value','detail':str(exc)});continue if load<=0:year_errors.append({'line':line,'error':'nonpositive consumption'}) nature[row['nature']]+=1 if norm(row['nature']).strip()!='donnees definitives':year_errors.append({'line':line,'error':'not definitive nature','nature':row['nature']}) for k,v in [('wind',wind),('solar',solar)]: signs[k+'_negative']+=v<0;minimum[k]=v if minimum[k] is None else min(minimum[k],v) actual[key]={'date':day,'hour':tm.hour,'minute':tm.minute,'load':load,'wind':wind,'solar':solar,'source_date_heure':row['date_heure']} expected=set();days=[];day=date(year,1,1) while day.year==year: days.append(day) expected.update((day,h,m) for h in range(24) for m in [0,30]);day+=timedelta(days=1) missing=expected-set(actual);extra=set(actual)-expected if missing or extra:year_errors.append({'error':'nominal actual coverage mismatch','missing':len(missing),'extra':len(extra),'examples':[str(x) for x in sorted(missing)[:10]]}) rows=[actual[key] for key in sorted(actual)] if sum(r['wind'] for r in rows)<=0 or sum(r['solar'] for r in rows)<=0:year_errors.append({'error':'nonpositive annual renewable sum'}) byday=defaultdict(list) for r in rows:byday[r['date']].append(r) if any(len(byday.get(d,[]))!=48 for d in days):year_errors.append({'error':'not48complete actual labels per source date'}) years.append({'year':year,'raw_rows':total,'actual_rows':len(rows),'expected_actual_rows':len(expected),'forecast_only_quarter_rows':len(quarter),'nominal_dates':len(days),'day_length_slots':48,'negative_counts':dict(signs),'minimum_reported_MW':minimum,'nature_labels':dict(nature),'errors':year_errors}) data[year]=[byday[d] for d in days if d in byday];errors.extend({'year':year,**e} for e in year_errors) report={'status':'PASS' if not errors else 'ABORT_NO_OUTCOMES','scope':'nominal48-slot source-labelled profiles only, NOT physical UTC chronology','physical_timeline_gate':'v1 ABORT retained; not repaired or passed','plan_sha256':sha(plan_bytes),'runner_sha256':sha(Path(__file__).read_bytes()),'metadata_sha256':sha(meta_bytes),'checked_utc':datetime.now(timezone.utc).isoformat(),'raw_inputs':inputs,'years':years,'errors':errors,'storage_outcomes_computed':False,'total_actual_rows':sum(x['actual_rows'] for x in years)} return plan,report,data def dst_transition(d): a=datetime.combine(d,time(),PARIS).astimezone(timezone.utc);b=datetime.combine(d+timedelta(days=1),time(),PARIS).astimezone(timezone.utc) return (b-a).total_seconds()!=86400 def scenario_days(raw_days,scenario): kept=[];excluded=[];clipped=Counter() for day in raw_days: d=day[0]['date'] if scenario['date_policy']=='exclude_both_Paris_DST_transition_dates' and dst_transition(d):excluded.append(d.isoformat());continue rows=[] for r in day: w,s=r['wind'],r['solar'] if scenario['negative_policy']=='max_zero_W_and_S_before_normalization': clipped['wind_values']+=w<0;clipped['solar_values']+=s<0;w=max(0,w);s=max(0,s) rows.append((r['load'],w,s)) kept.append((d,rows)) return kept,excluded,dict(clipped) def summarize_days(days,wind_share): sumL=sum(x[0] for _,rs in days for x in rs);sumW=sum(x[1] for _,rs in days for x in rs);sumS=sum(x[2] for _,rs in days for x in rs) a=wind_share;b=100;den=b*sumW*sumS if den<=0:raise AssertionError('Nonpositive denominator') summaries=[];all_values=[];level=0;lo=hi=0;trajectory=[] for d,rows in days: vals=[den*L-sumL*(a*sumS*W+(b-a)*sumW*S) for L,W,S in rows] c=0;mn=mx=0 for v in vals: c+=v;mn=min(mn,c);mx=max(mx,c);level+=v;lo=min(lo,level);hi=max(hi,level) summaries.append({'date':d.isoformat(),'month':d.month,'net':c,'minimum':mn,'maximum':mx,'values':vals}) all_values.extend(vals);trajectory.append({'date':d.isoformat(),'cumulative_nominal_MWh_equivalent':float(Fraction(level,2*den))}) if sum(all_values)!=0 or level!=0:raise AssertionError('Exact cycle energy balance failed') return summaries,all_values,den,sumL,sumW,sumS,hi-lo,trajectory def range_days(days): level=low=high=0 for d in days: low=min(low,level+d['minimum']);high=max(high,level+d['maximum']);level+=d['net'] if level!=0:raise AssertionError('Permutation cycle balance changed') return high-low def range_full(days): c=lo=hi=0 for d in days: for v in d['values']:c+=v;lo=min(lo,c);hi=max(hi,c) if c!=0:raise AssertionError('Direct cycle balance failed') return hi-lo def order_for(days,year,scenario,rep): text=f'incoprea-storage-v2|year={year}|scenario={scenario}|replicate={rep}' seed=int.from_bytes(hashlib.sha256(text.encode()).digest()[:8],'big');rng=random.Random(seed);order=[] for month in range(1,13): group=[d for d in days if d['month']==month];rng.shuffle(group);order.extend(group) return seed,order def quantile(values,numerator,denominator=1000): values=sorted(values);pos=Fraction((len(values)-1)*numerator,denominator);i=pos.numerator//pos.denominator;weight=pos-i return Fraction(values[i]) if i==len(values)-1 else (1-weight)*values[i]+weight*values[i+1] def fraction_or_none(x):return None if x is None else float(x) def main(): ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('mode',choices=['preflight','run']);ap.add_argument('--directory',type=Path,default=Path(__file__).resolve().parent);ap.add_argument('--approval',type=Path);args=ap.parse_args();directory=args.directory.resolve() plan,preflight,data=load_nominal(directory) if args.mode=='preflight': write_json(directory/'preflight-v2.json',preflight) if preflight['status']=='PASS': rows=[{'year':year,'source_date':r['date'].isoformat(),'source_hour':f"{r['hour']:02}:{r['minute']:02}",'source_date_heure':r['source_date_heure'],'consumption_MW':r['load'],'wind_MW':r['wind'],'solar_MW':r['solar']} for year,days in data.items() for day in days for r in day] write_csv(directory/'validated-nominal-actuals.csv',rows) print(json.dumps({'status':preflight['status'],'total_actual_rows':preflight['total_actual_rows'],'errors':preflight['errors'],'plan_sha256':preflight['plan_sha256'],'runner_sha256':preflight['runner_sha256']},indent=2));return if preflight['status']!='PASS':raise SystemExit('Nominal preflight failed; no outcomes') if not args.approval or not args.approval.is_file():raise SystemExit('Exact plan/runner internal approval is required before outcomes') approval=json.loads(args.approval.read_text()) if not approval.get('approved_for_execution') or approval.get('plan_sha256')!=preflight['plan_sha256'] or approval.get('runner_sha256')!=preflight['runner_sha256']:raise SystemExit('Approval is not bound to this plan and code') output=directory/'results' if output.exists():raise SystemExit('Results directory exists; no repeat or overwrite') output.mkdir();start=datetime.now(timezone.utc).isoformat();tick=clock.monotonic() lock={'plan_sha256':preflight['plan_sha256'],'runner_sha256':preflight['runner_sha256'],'raw_inputs':preflight['raw_inputs'],'approval_path':str(args.approval),'approval_sha256':sha(args.approval.read_bytes()),'started_utc':start,'scope':preflight['scope']};write_json(output/'EXECUTION-LOCK.json',lock) control=[] for label,seq in [('persistent',[1,1,-1,-1]),('alternating',[1,-1,1,-1])]: c=lo=hi=0;trace=[0] for v in seq:c+=v;lo=min(lo,c);hi=max(hi,c);trace.append(c) control.append({'label':label,'residual_MW':seq,'cumulative_nominal_MWh':[float(Fraction(x,2)) for x in trace],'sum_MW_samples':sum(seq),'mean':0,'population_variance':1,'positive_peak_MW':1,'capacity_MWh':float(Fraction(hi-lo,2))}) if [x['capacity_MWh'] for x in control]!=[1,.5]:raise AssertionError('Deterministic control failed') write_json(output/'deterministic-control.json',control) summaries=[];permutations=[];day_records=[];original_paths=[] for scenario in plan['scenarios']: sid=scenario['id'] for year in plan['data']['years']: kept,excluded,clipped=scenario_days(data[year],scenario) days,values,D,sumL,sumW,sumS,original,trajectory=summarize_days(kept,scenario['wind_share_percent']) if scenario['date_policy']=='exclude_both_Paris_DST_transition_dates' and len(excluded)!=2:raise AssertionError('Expected two DST transition dates removed') if range_days(days)!=original or range_full(days)!=original:raise AssertionError('Original daily-summary optimization differs') marginal_hash=sha(('\n'.join(map(str,sorted(values)))+'\n').encode());month_nets={str(m):sum(d['net'] for d in days if d['month']==m) for m in range(1,13)} for d in days:day_records.append({'year':year,'scenario':sid,'date':d['date'],'month':d['month'],'slots':len(d['values']),'net_numerator':str(d['net']),'min_prefix_numerator':str(d['minimum']),'max_prefix_numerator':str(d['maximum']),'residual_denominator':str(D)}) ranges=[];ids={d['date'] for d in days} for rep in range(plan['permutation_design']['replicates_per_year_scenario']): seed,order=order_for(days,year,sid,rep) if {d['date'] for d in order}!=ids or len(order)!=len(days):raise AssertionError('Day multiset changed') if any(str(m) not in month_nets or sum(d['net'] for d in order if d['month']==m)!=month_nets[str(m)] for m in range(1,13)):raise AssertionError('Monthly total changed') R=range_days(order) if rep==0 and range_full(order)!=R:raise AssertionError('Permuted daily-summary optimization differs') ranges.append(R);capacity=Fraction(R,2*D);hours=capacity/Fraction(sumL,len(values));order_hash=sha(','.join(d['date'] for d in order).encode()) permutations.append({'year':year,'scenario':sid,'replicate':rep,'seed':seed,'range_numerator':str(R),'capacity_nominal_MWh_equivalent':float(capacity),'capacity_nominal_GWh_equivalent':float(capacity/1000),'equivalent_mean_demand_nominal_hours':float(hours),'order_sha256':order_hash}) med=quantile(ranges,500);low=quantile(ranges,25);high=quantile(ranges,975);upper=Fraction(1+sum(R>=original for R in ranges),1000);lower=Fraction(1+sum(R<=original for R in ranges),1000);ratio=Fraction(original,med) if med else None C=Fraction(original,2*D);M=med/(2*D);mean_load=Fraction(sumL,len(values)) record={'year':year,'scenario':sid,'wind_share_percent':scenario['wind_share_percent'],'date_policy':scenario['date_policy'],'negative_policy':scenario['negative_policy'],'nominal_days':len(days),'nominal_slots':len(values),'excluded_dates':excluded,'clipped_value_counts':clipped,'sum_consumption_MW_samples':sumL,'sum_wind_MW_samples':sumW,'sum_solar_MW_samples':sumS,'mean_demand_MW':float(mean_load),'normalization_denominator':str(D),'exact_cycle_residual_numerator':str(sum(values)),'residual_population_variance_MW2':float(Fraction(sum(v*v for v in values),len(values)*D*D)),'residual_positive_peak_MW':float(Fraction(max(values),D)),'residual_absolute_peak_MW':float(Fraction(max(map(abs,values)),D)),'residual_multiset_sha256':marginal_hash,'month_net_numerators':{m:str(v) for m,v in month_nets.items()},'original_range_numerator':str(original),'original_nominal_GWh_equivalent':float(C/1000),'original_equivalent_mean_demand_nominal_hours':float(C/mean_load),'permutation_median_nominal_GWh_equivalent':float(M/1000),'permutation_q025_nominal_GWh_equivalent':float(low/(2*D*1000)),'permutation_q975_nominal_GWh_equivalent':float(high/(2*D*1000)),'permutation_min_nominal_GWh_equivalent':float(Fraction(min(ranges),2*D*1000)),'permutation_max_nominal_GWh_equivalent':float(Fraction(max(ranges),2*D*1000)),'original_to_permutation_median_ratio':fraction_or_none(ratio),'difference_original_minus_median_nominal_GWh_equivalent':float((C-M)/1000),'difference_original_minus_median_nominal_hours':float((C-M)/mean_load),'upper_reference_rank':float(upper),'lower_reference_rank':float(lower),'two_sided_reference_rank':float(min(Fraction(1),2*min(upper,lower))),'original_above_permutation_median':original>med,'original_below_permutation_median':originalhigh,'permutations':len(ranges),'daily_summary_original_and_first_permutation_match_direct_slots':True,'original_date_order':[d['date'] for d in days]} summaries.append(record) for r in trajectory:original_paths.append({'year':year,'scenario':sid,**r}) print(f"{sid} {year}: original/median={float(ratio):.6f}",flush=True) primary=[r for r in summaries if r['scenario']=='primary_70_30'];ranked=sorted(primary,key=lambda r:r['two_sided_reference_rank']);adjusted=0 for i,r in enumerate(ranked):adjusted=max(adjusted,min(1,(len(ranked)-i)*r['two_sided_reference_rank']));r['two_sided_reference_rank_holm_descriptive']=adjusted across=[] for scenario in plan['scenarios']: group=[r for r in summaries if r['scenario']==scenario['id']];ratios=[r['original_to_permutation_median_ratio'] for r in group if r['original_to_permutation_median_ratio'] is not None];logs=[math.log(x) for x in ratios if x>0] across.append({'scenario':scenario['id'],'years':len(group),'median_original_to_reference_ratio':statistics.median(ratios) if ratios else None,'minimum_ratio':min(ratios) if ratios else None,'maximum_ratio':max(ratios) if ratios else None,'years_above_reference_median':sum(r['original_above_permutation_median'] for r in group),'years_below_reference_median':sum(r['original_below_permutation_median'] for r in group),'years_outside_central95_reference':sum(r['original_outside_central95_reference'] for r in group),'mean_year_log_ratio':statistics.mean(logs) if logs else None,'sd_year_log_ratio':statistics.stdev(logs) if len(logs)>1 else None}) if len(permutations)!=len(plan['scenarios'])*len(plan['data']['years'])*999:raise AssertionError('Unexpected permutation count') write_json(output/'year-scenario-results.json',summaries);write_json(output/'permutation-results.json',permutations);write_csv(output/'permutation-results.csv',permutations);write_csv(output/'day-summaries.csv',day_records);write_csv(output/'original-daily-cumulative.csv',original_paths) scalar=[] for r in summaries: z={k:v for k,v in r.items() if not isinstance(v,(dict,list))};z.setdefault('two_sided_reference_rank_holm_descriptive',None);scalar.append(z) write_csv(output/'year-scenario-results.csv',scalar);write_json(output/'across-year-summaries.json',across);write_csv(output/'across-year-summaries.csv',across) metadata={**lock,'finished_utc':datetime.now(timezone.utc).isoformat(),'elapsed_seconds':clock.monotonic()-tick,'python':sys.version,'platform':platform.platform(),'scenario_year_conditions':len(summaries),'permutations':len(permutations),'physical_timeline_claim_allowed':False,'plan_modified_during_run':sha((directory/'plan-v2.json').read_bytes())!=lock['plan_sha256']} write_json(output/'run-metadata.json',metadata);write_json(output/'results.json',{'metadata':metadata,'controls':control,'across_years':across,'conditions':summaries});print(json.dumps({'finished':metadata['finished_utc'],'conditions':len(summaries),'permutations':len(permutations),'across_years':across},indent=2)) if __name__=='__main__':main()