"""Reproduce the IC86_XI energy–uncertainty analysis.
Usage: python icetracks_analysis.py path/to/attached.json
Requires numpy, pandas, scipy, matplotlib. Outputs are written in the current directory.
Seed 20261008; 1000 bootstrap replicates; 200 rounding perturbations.
"""
import sys, json, io, hashlib
from pathlib import Path
import numpy as np
import pandas as pd
from scipy.stats import spearmanr, rankdata
p=Path(sys.argv[1])
payload=json.loads(p.read_text())
records={r['path']:r for r in payload['records']}
for r in records.values():
 assert hashlib.sha256(r['content'].encode()).hexdigest()==r['sha256']
s=records['IC86_XI_events.csv']['content']
cols=s.splitlines()[0].lstrip('#').split()
df=pd.read_csv(io.StringIO(s),sep=r'\s+',comment='#',names=cols)
E='log10(E/GeV)'; A='AngErr[deg]'
assert len(df)==127695 and np.isfinite(df[cols]).all().all()
assert not df.duplicated(['run','event','subevent']).any()
df['utc']=pd.to_datetime(df['MJD[days]'],unit='D',origin='1858-11-17',utc=True)
df['year']=df.utc.dt.year
df['dec_band']=pd.cut(df['Dec[deg]'],[-90,-30,0,30,90],right=False)
df['ra_band']=pd.cut(df['RA[deg]'],[0,90,180,270,360],right=False)
u=pd.read_csv(io.StringIO(records['IC86_XI_uptime.csv']['content']),sep=r'\s+',comment='#',names=['start','stop'])
ix=np.searchsorted(u.start.to_numpy(),df['MJD[days]'],side='right')-1
inside=(ix>=0)&(df['MJD[days]'].to_numpy()<=u.stop.to_numpy()[np.maximum(ix,0)])
assert inside.all()


SEED=20261008
B=1000
rng=np.random.default_rng(SEED)
def boot_rho(g,B=B):
 x,xi=np.unique(g[E],return_inverse=True); y,yi=np.unique(g[A],return_inverse=True)
 pairs,counts=np.unique(np.column_stack([xi,yi]),axis=0,return_counts=True)
 i,j=pairs.T; n=counts.sum()
 def calc(w):
  mx=np.bincount(i,weights=w,minlength=len(x)); my=np.bincount(j,weights=w,minlength=len(y))
  rx=np.cumsum(mx)-mx/2; ry=np.cumsum(my)-my/2
  dx=rx[i]-w.sum()/2; dy=ry[j]-w.sum()/2
  return np.sum(w*dx*dy)/np.sqrt(np.sum(w*dx*dx)*np.sum(w*dy*dy))
 point=calc(counts)
 assert abs(point-spearmanr(g[E],g[A]).statistic)<1e-12
 bs=np.array([calc(rng.multinomial(n,counts/n)) for _ in range(B)])
 return point,*np.quantile(bs,[.025,.975])
groups=[('overall','all',df)]
for field in ['year','dec_band','ra_band']:
 groups += [(field,str(k),g) for k,g in df.groupby(field,observed=True)]
for (d,y),g in df.groupby(['dec_band','year'],observed=True): groups.append(('declination × year',str(d)+' / '+str(y),g))
groups += [('sensitivity','above floor (>0.20°)',df[df[A]>.2]),('sensitivity','above 0.21°',df[df[A]>.21])]
results=[]
for kind,label,g in groups:
 rho,lo,hi=boot_rho(g)
 results.append(dict(group=kind,label=label,n=len(g),rho=rho,ci_low=lo,ci_high=hi,floor_fraction=np.mean(g[A]==.2)))
corr=pd.DataFrame(results)
print(corr.round(5).to_string(index=False))
qrows=[]
for kind,label,g in groups:
 if kind=='sensitivity':continue
 bins=pd.qcut(g[E],10,duplicates='drop')
 for interval,h in g.groupby(bins,observed=True):
  qs=h[A].quantile([.1,.25,.5,.75,.9]).to_numpy()
  qrows.append(dict(group=kind,label=label,n=len(h),energy_min=h[E].min(),energy_max=h[E].max(),energy_median=h[E].median(),q10=qs[0],q25=qs[1],q50=qs[2],q75=qs[3],q90=qs[4],floor_fraction=np.mean(h[A]==.2)))
quant=pd.DataFrame(qrows)
print('overall quantiles\n',quant[quant.group=='overall'].round(4).to_string(index=False))

sensitivity_rows=[]
for k,g in df.groupby('dec_band',observed=True):
 h=g[g[A]>.2]; r,l,u_=boot_rho(h)
 sensitivity_rows.append(dict(test='declination above floor',label=str(k),n=len(h),rho=r,low=l,high=u_))
# Perturb within the displayed rounding cells; floor ties stay tied.
jitter=[]
for b in range(200):
 xx=df[E].to_numpy()+rng.uniform(-.005,.005,len(df))
 yy=df[A].to_numpy().copy(); free=yy>.2
 yy[free]+=rng.uniform(-.005,.005,free.sum())
 jitter.append(spearmanr(xx,yy).statistic)
print('rounding jitter rho range',min(jitter),max(jitter))
# Run-cluster bootstrap, recalculating weighted midranks in every replicate.
_,xi=np.unique(df[E],return_inverse=True); _,yi=np.unique(df[A],return_inverse=True)
_,ri=np.unique(df.run,return_inverse=True); nr=ri.max()+1
cluster=[]
for b in range(1000):
 w=rng.multinomial(nr,np.ones(nr)/nr)[ri]
 mx=np.bincount(xi,weights=w); my=np.bincount(yi,weights=w)
 rx=(np.cumsum(mx)-mx/2)[xi]-w.sum()/2; ry=(np.cumsum(my)-my/2)[yi]-w.sum()/2
 cluster.append(np.sum(w*rx*ry)/np.sqrt(np.sum(w*rx*rx)*np.sum(w*ry*ry)))
print('run clusters',nr,'CI',np.quantile(cluster,[.025,.975]))
print(pd.DataFrame(sensitivity_rows).round(5).to_string(index=False))
fine=[]
for k,g in df.groupby(pd.cut(df['Dec[deg]'],np.arange(-90,91,10),right=False),observed=True):
 fine.append({'declination':str(k),'n':len(g),'rho':spearmanr(g[E],g[A]).statistic})
print('fine declination',pd.DataFrame(fine).round(3).to_string(index=False))
lowcut,highcut=df[E].quantile([.1,.9])
# Transparent extreme examples, not a new event-quality classification.
high=df[df[E]>=highcut].nlargest(3,A).copy(); high['example_type']='top energy decile, largest uncertainty'
low=df[df[E]<=lowcut].nsmallest(3,A).copy(); low['example_type']='bottom energy decile, smallest uncertainty'
examples=pd.concat([high,low]); examples['reco_muon_energy_GeV']=10**examples[E]
print('counterexample thresholds',lowcut,highcut)
print(examples[['example_type','run','event','subevent','MJD[days]',E,A,'RA[deg]','Dec[deg]']].to_string(index=False))
print('high-energy uncertainty >1deg',np.mean(df.loc[df[E]>=highcut,A]>1))

fine_results=[]
for k,g in df.groupby(pd.cut(df['Dec[deg]'],np.arange(-90,91,10),right=False),observed=True):
 r,l,h=boot_rho(g)
 fine_results.append(dict(declination=str(k),n=len(g),rho=r,ci_low=l,ci_high=h))
fine_table=pd.DataFrame(fine_results)


import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
from matplotlib.ticker import NullFormatter
plt.rcParams.update({'axes.spines.top':False,'axes.spines.right':False})
plt.rcParams.update({'font.size':10,'axes.labelsize':10,'axes.titlesize':10,'xtick.labelsize':9,'ytick.labelsize':9,'legend.fontsize':9})
figs=[]
def qpanel(ax,kind,label,title):
 t=quant[(quant.group==kind)&(quant.label==label)]
 x=t.energy_median.to_numpy()
 ax.fill_between(x,t.q10,t.q90,color='#0072B2',alpha=.13)
 ax.fill_between(x,t.q25,t.q75,color='#0072B2',alpha=.3)
 ax.plot(x,t.q50,'D-',color='#0072B2',markersize=4,lw=1.5)
 ax.axhline(.2,color='#555555',ls='--',lw=1)
 ax.set_title(title,loc='left',pad=12)
 ax.set_yscale('log'); ax.set_ylim(.17,8); ax.set_yticks([.2,.5,1,2,5]); ax.set_yticklabels(['0.2','0.5','1','2','5'])
 ax.set_xlim(2.45,5.95); ax.set_xticks([2.5,3,4,5])
 ax.grid(axis='y',alpha=.18)
 row=corr[(corr.group==kind)&(corr.label==label)].iloc[0]
 ax.text(.96,.94,f'n = {int(row.n):,} events\nρ = {row.rho:.3f}',ha='right',va='top',transform=ax.transAxes,fontsize=9)
handles=[Line2D([0],[0],color='#0072B2',marker='D',label='Median'),Patch(color='#0072B2',alpha=.3,label='25–75% of events'),Patch(color='#0072B2',alpha=.13,label='10–90% of events'),Line2D([0],[0],color='#555555',ls='--',label='0.2° reporting floor')]
fig,axs=plt.subplots(1,3,figsize=(12,4.5),sharex=True,sharey=True)
for ax,(kind,label,title) in zip(axs,[('overall','all','All events'),('year','2021','2021: May–December'),('year','2022','2022: January–May')]):qpanel(ax,kind,label,title)
axs[0].set_ylabel('Estimated angular uncertainty (degrees)')
fig.supxlabel('Reconstructed muon energy indicator, log₁₀(E / GeV)',y=.12,fontsize=10)
fig.suptitle('The overall decrease persists in both calendar years',y=.97,fontsize=10)
fig.text(.08,.89,'Smaller uncertainty = tighter estimated localization',fontsize=9)
fig.legend(handles=handles,loc='lower center',ncol=4,bbox_to_anchor=(.5,.005),frameon=False)
fig.subplots_adjust(left=.08,right=.99,bottom=.25,top=.80,wspace=.05)
figs.append((fig,'energy_year_quantiles.png'))
fig2,axes=plt.subplots(2,4,figsize=(14,8),sharex=True,sharey=True)
for ax,label in zip(axes[0],corr[corr.group=='dec_band'].label):qpanel(ax,'dec_band',label,'Declination '+label+'°')
for ax,label in zip(axes[1],corr[corr.group=='ra_band'].label):qpanel(ax,'ra_band',label,'Right ascension '+label+'°')
for ax in axes[:,0]:ax.set_ylabel('Estimated angular uncertainty (degrees)')
fig2.supxlabel('Reconstructed muon energy indicator, log₁₀(E / GeV)',y=.075,fontsize=10)
fig2.suptitle('The relationship varies with declination; right-ascension sectors look similar',y=.98,fontsize=10)
fig2.text(.075,.925,'Smaller uncertainty = tighter estimated localization',fontsize=9)
fig2.legend(handles=handles,loc='lower center',ncol=4,bbox_to_anchor=(.5,.005),frameon=False)
fig2.subplots_adjust(left=.075,right=.99,bottom=.15,top=.87,hspace=.30,wspace=.05)
figs.append((fig2,'sky_quantiles.png'))
# Save, then verify label geometry.
for f,name in figs:
 for ax in f.axes: ax.yaxis.set_minor_formatter(NullFormatter())
 f.savefig(name,dpi=200)
 f.canvas.draw(); renderer=f.canvas.get_renderer()
 texts=[(t,t.get_window_extent(renderer)) for t in f.findobj(mpl.text.Text) if t.get_text().strip() and t.get_visible()]
 overlaps=[(a.get_text(),b.get_text()) for i,(a,ba) in enumerate(texts) for b,bb in texts[i+1:] if ba.overlaps(bb)]
 outside=[t.get_text() for t,b in texts if not f.bbox.contains(b.x0,b.y0) or not f.bbox.contains(b.x1,b.y1)]
 print(name,'text overlaps',overlaps,'outside',outside)
plt.close('all')

import zipfile, platform, scipy
summary={'dataset_doi':payload['dataset_provenance']['dataset_doi'],'version':'3.1','season':'IC86_XI','input_sha256':hashlib.sha256(p.read_bytes()).hexdigest(),'n_events':len(df),'start_utc':str(df.utc.min()),'end_utc':str(df.utc.max()),'n_runs':int(nr),'n_uptime_intervals':len(u),'livetime_days':float((u.stop-u.start).sum()),'floor_count':int((df[A]==.2).sum()),'floor_fraction':float((df[A]==.2).mean()),'rounding_jitter_rho_range':[min(jitter),max(jitter)],'cluster_bootstrap_ci':np.quantile(cluster,[.025,.975]).tolist(),'bootstrap_replicates':1000,'jitter_replicates':200,'seed':SEED,'software':{'python':platform.python_version(),'numpy':np.__version__,'pandas':pd.__version__,'scipy':scipy.__version__,'matplotlib':mpl.__version__},'counterexample_low_threshold':float(lowcut),'counterexample_high_threshold':float(highcut),'high_energy_above_1_degree_fraction':float(np.mean(df.loc[df[E]>=highcut,A]>1))}
with zipfile.ZipFile('icetracks_results.zip','w',zipfile.ZIP_DEFLATED) as z:
 for name,table in [('correlations.csv',corr),('quantiles.csv',quant),('floor_sensitivity.csv',pd.DataFrame(sensitivity_rows)),('fine_declination.csv',fine_table),('counterexamples.csv',examples)]:z.writestr(name,table.to_csv(index=False))
 z.writestr('analysis_metadata.json',json.dumps(summary,indent=2))
print(json.dumps(summary,indent=2))