"""Reproduce the exploratory summaries and SVG figures from retrieved series.
Run: python analyze.py (requires numpy and matplotlib).
No interpolation or missing-to-zero replacement is performed here.
"""
from pathlib import Path
import json, datetime, html
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

ROOT=Path(__file__).resolve().parent
read=lambda name: json.loads((ROOT/'data'/name).read_text())
T=read('temperature_series.json'); F=read('fish_series.json'); E=read('early_series.json')
for item in E:
    T.append({'year':item['year'],'series':[s for s in item['series'] if 'TEMPERATURE' in s['name']]})
    F.append({'year':item['year'],'series':[s for s in item['series'] if 'CHINOOK' in s['name']]})
def day(t):return datetime.datetime.fromtimestamp(t/1000,datetime.timezone.utc)
def summer(s):return {t:v for t,v in s['pts'] if 6<=day(t).month<=8}
air={}; water={}; counts={}; coverage=[]
for item in T:
    for s in item['series']:
        m=summer(s)
        if s['name'].startswith('Hermiston'):air[item['year']]=m;coverage.append({'year':item['year'],'series':'Hermiston air','days':len(m)})
        elif 'TEMPERATURE' in s['name']:
            dam=s['name'].split('_')[0];water[item['year'],dam]=m;coverage.append({'year':item['year'],'series':dam+' water','days':len(m)})
for item in F:
    for s in item['series']:counts[item['year'],s['name'].split('_')[0]]=s['pts']
years=sorted(y for y,m in air.items() if len(m)==92 and all(len(water.get((y,d),{}))==92 for d in ['BON','MCN','LWG']))
summary=[];annual=[]
for y in years:
    for dam in ['BON','MCN','LWG']:
        w=water[y,dam];pts=counts.get((y,dam),[]);p=[(t,v) for t,v in pts if t in w]
        total=sum(v for t,v in p)
        share=sum(v for t,v in p if w[t]>20)/total if len(p)>=83 and total>0 else None
        fall=sorted((t,v) for t,v in pts if 8<=day(t).month<=10)
        median=None
        if len(fall)>=88 and sum(v for t,v in fall)>0 and all(v>=0 for t,v in fall):
            target=sum(v for t,v in fall)*.5;cum=0
            for t,v in fall:
                cum+=v
                if cum>=target:median=int(day(t).strftime('%j'));break
        annual.append({'year':y,'dam':dam,'air_C':float(np.mean(list(air[y].values()))),'water_C':float(np.mean(list(w.values()))),'passage_share_above20':share,'matched_passage_days':len(p),'fall_median_doy':median,'fall_recorded_days':len(fall)})
for dam in ['BON','MCN','LWG']:
    a=[r for r in annual if r['dam']==dam]
    x=np.array([r['air_C'] for r in a]);w=np.array([r['water_C'] for r in a])
    heat=[r for r in a if r['passage_share_above20'] is not None];fall=[r for r in a if r['fall_median_doy'] is not None]
    # Complete 13 x 92 balanced panel: remove calendar-date and year means.
    aa=np.array([[air[y][t] for t in sorted(air[y])] for y in years]);ww=np.array([[water[y,dam][t] for t in sorted(air[y])] for y in years])
    residual=lambda z:z-z.mean(axis=0,keepdims=True)-z.mean(axis=1,keepdims=True)+z.mean()
    ar=residual(aa);wr=residual(ww)
    full=sorted((y,m) for (y,d),m in water.items() if d==dam and len(m)==92)
    fy=np.array([y for y,m in full]);fv=np.array([np.mean(list(m.values())) for y,m in full]);trend=float(np.polyfit(fy-2000,fv,1)[0]*10)
    summary.append({'dam':dam,'complete_air_summers':len(a),'air_water_r':float(np.corrcoef(x,w)[0,1]),'air_hot_passage_r':float(np.corrcoef([r['air_C'] for r in heat],[r['passage_share_above20'] for r in heat])[0,1]),'hot_passage_years':len(heat),'fall_timing_r':float(np.corrcoef([r['air_C'] for r in fall],[r['fall_median_doy'] for r in fall])[0,1]),'fall_years':len(fall),'daily_two_way_anomaly_r':float(np.corrcoef(ar.ravel(),wr.ravel())[0,1]),'water_trend_C_per_decade':trend,'water_trend_years':[int(fy.min()),int(fy.max())],'water_trend_n':len(fy)})
(ROOT/'results.json').write_text(json.dumps({'complete_air_years':years,'summary':summary,'annual_matched':annual,'coverage':coverage},indent=2))
plt.rcParams.update({'font.size':11,'svg.fonttype':'none','axes.spines.top':False,'axes.spines.right':False})
colors={'BON':'#276a9c','MCN':'#d47c2c','LWG':'#38846c'};names={'BON':'Bonneville','MCN':'McNary','LWG':'Lower Granite'}
fig,axes=plt.subplots(1,2,figsize=(12,4.8),layout='constrained')
a=[r for r in annual if r['dam']=='MCN'];x=np.array([r['air_C'] for r in a]);w=np.array([r['water_C'] for r in a]);share=np.array([r['passage_share_above20']*100 for r in a])
for ax,z,label,title in [(axes[0],w,'McNary mean water temperature (°C)','Warmer air, warmer river'),(axes[1],share,'Recorded summer passage on days >20°C (%)','Warmer summers, more warm-water passage')]:
    ax.scatter(x,z,s=58,c=[r['year'] for r in a],cmap='viridis',edgecolors='white',linewidth=.6)
    xx=np.linspace(x.min()-.15,x.max()+.15,100);b=np.polyfit(x,z,1);ax.plot(xx,np.polyval(b,xx),color='#557080',lw=1.3)
    for r,yy in zip(a,z):
        offsets={2007:(5,-12),2013:(6,-12),2014:(-24,11),2017:(8,5)}
        ax.annotate(str(r['year']),(r['air_C'],yy),xytext=offsets.get(r['year'],(4,4)),textcoords='offset points',fontsize=8)
    ax.set(xlabel='Hermiston mean air temperature (°C)',ylabel=label,title=title);ax.grid(alpha=.17)
fig.suptitle('June–August | 13 complete summers, 2007–2024',fontsize=15,fontweight='bold')
fig.savefig(ROOT/'air_water_and_passage.svg');fig.savefig(ROOT/'air_water_and_passage.png',dpi=150);plt.close(fig)
fig,ax=plt.subplots(figsize=(11,5),layout='constrained')
for dam in ['BON','MCN','LWG']:
    a=sorted((y,m) for (y,d),m in water.items() if d==dam and len(m)==92)
    ax.plot([y for y,m in a],[np.mean(list(m.values())) for y,m in a],marker='o',ms=4,color=colors[dam],label=names[dam])
ax.set(title='The long record does not tell one uniform warming story',xlabel='Year',ylabel='June–August mean water temperature (°C)');ax.grid(alpha=.17);ax.legend(ncol=3)
fig.savefig(ROOT/'summer_water_record.svg');fig.savefig(ROOT/'summer_water_record.png',dpi=150);plt.close(fig)
table=''.join(f"<tr><td>{names[r['dam']]}</td><td>{r['air_water_r']:.2f}</td><td>{r['air_hot_passage_r']:.2f}</td><td>{r['fall_timing_r']:.2f} (n={r['fall_years']})</td></tr>" for r in summary)
trendtable=''.join(f"<tr><td>{names[r['dam']]}</td><td>{r['water_trend_years'][0]}–{r['water_trend_years'][1]}</td><td>{r['water_trend_n']}</td><td>{r['water_trend_C_per_decade']:+.2f}</td></tr>" for r in summary)
methods='''<p>Exploratory calculations on RainbowStats-returned daily series, retrieved October 22, 2026. This is a research memo, not a publication-ready blog.</p><ul><li>Air: AQS site 41-059-1003, outdoor temperature 62101, POC 1, duration 1 HOUR, blank standard, No Events; arithmetic_mean. Fahrenheit converted to Celsius. Daily summaries are not hourly records. Observation-hour completeness and instrument continuity remain to be audited.</li><li>Water: DART(BON/MCN/LWG,TEMPERATURE). Fish: DART(BON/MCN/LWG,CHINOOK). No interpolation or missing-to-zero replacement in this analysis. Chart-returned series are preserved in data/; these are not original EPA hourly records.</li><li>Primary air comparisons use 13 summers with all 92 June–August dates present in air and all three water series: 2007, 2008, 2010–2015, 2017, 2020, 2022–2024. Excluded from the 2007–2025 air comparison: 2009, 2016, 2018, 2019, 2021, 2025. In 2016 the selected AQS series returned no numeric observations; 2019 had 91 days and 2025 had 86. In 2021 air had 62 days and Lower Granite water had 76. With a relaxed 90% air-day rule, 2019 and 2025 enter and the annual air–water correlations are 0.88, 0.86, 0.82.</li><li>Warm-water passage share is the fraction of recorded June–August Chinook passage counts on dates where the matching reported daily water temperature exceeded 20°C. At least 83 matching days are required. It is a station-level association, not individual exposure, survival, abundance, or a universal biological threshold.</li><li>Fall timing uses the date at which cumulative recorded passage within August 1–October 31 reaches 50%, requiring at least 88 of 92 dates and nonnegative counts. This is an operational window, not a biological stock classification; incomplete dates are never called zeros. Bonneville and McNary have 13 eligible years; Lower Granite has 11.</li><li>Daily anomaly comparisons subtract both calendar-date means and year means on the balanced 13 × 92 panel. Correlations are 0.48 (Bonneville), 0.40 (McNary), 0.20 (Lower Granite). They remain descriptive; discharge, operations, serial dependence, and regional representativeness are not controlled here.</li><li>Long-record trends are simple OLS slopes of complete-summer water means, not climate attribution estimates. Bonneville begins in 2001 because 2000 has 88 summer dates; McNary temperature was not returned for 2000–2004; Lower Granite excludes incomplete 2021. The Lower Granite result is sensitive to start date: its full-record slope is negative, while the recent matched-air-period slope is positive. Verify station location, measurement definitions, source changes, and cold-water operations before interpreting this difference.</li><li>USACE documents deliberate cold-water releases from Dworshak to regulate Lower Granite tailwater. This is a reason to investigate operations, not proof that operations explain the apparent historical change.</li><li>Repeated exploration across years means this is retrospective evidence. There are 13 independent seasonal units, not 1,196 independent daily climate observations. No causal temperature effect, mortality claim, climate attribution, or general migration acceleration is established.</li></ul>'''
sources='''<ul><li><a href="https://rainbowstats.com/commandDesc/AQS.html">RainbowStats AQS commands</a></li><li><a href="https://www.epa.gov/aqs">EPA AQS</a></li><li><a href="https://www.cbr.washington.edu/dart/metadata/river">DART river metadata</a></li><li><a href="https://www.cbr.washington.edu/dart/metadata/adult">DART passage metadata</a></li><li><a href="https://www.nww.usace.army.mil/Media/News-Releases/Article/3824182/dworshak-dam-to-increase-releases-beginning-july-1st/">USACE cold-water operations</a></li><li><a href="https://rainbowstats.com/blog/Oct2026/columbia-river1/index.html">Remembering the River</a></li></ul>'''
report='''<!doctype html><html lang="en"><meta charset="utf-8"><meta name="viewport" content="width=device-width,initial-scale=1"><title>Columbia River: Air, Water, and Passage — Research Memo</title><style>body{font:18px/1.65 Georgia,serif;color:#243747;background:#f5f8fb;margin:0}main{max-width:1050px;margin:auto;background:white;padding:36px}h1,h2,table{font-family:Arial,sans-serif}h1{color:#254d69;line-height:1.2}h2{background:#e8f3fb;padding:10px 14px;font-size:23px}table{border-collapse:collapse;width:100%;font-size:16px}th,td{text-align:left;padding:10px;border-bottom:1px solid #d8e2ea}th{background:#edf5fa}svg{width:100%;height:auto}li{margin:10px 0}.note{background:#fff5dc;padding:15px}a{color:#176692}</style><main><p>RAINBOWSTATS™ · EXPLORATORY RESEARCH · OCTOBER 22, 2026</p><h1>Warmer Air, Warmer Water—But What Changes for the Fish?</h1><p>The strongest initial result is a relationship between summer warmth and the share of recorded summer Chinook passage occurring on warm-water days. Fall passage timing does not show a clear corresponding association. The long water-temperature record differs by dam and must be checked before making a broad warming claim.</p><h2>What the complete summers show</h2><table><tr><th>Dam</th><th>Air–water correlation</th><th>Air–warm-water passage correlation</th><th>Air–fall median date correlation</th></tr>'''+table+'</table>'+(ROOT/'air_water_and_passage.svg').read_text().split('<svg',1)[1].join(['<svg',''])+'''<h2>The full record changes the interpretation</h2>'''+(ROOT/'summer_water_record.svg').read_text().split('<svg',1)[1].join(['<svg',''])+'''<table><tr><th>Dam</th><th>Complete-summer span</th><th>Years</th><th>Descriptive °C/decade</th></tr>'''+trendtable+'''</table><p class="note">These different slopes prevent a simple claim that all three dams have warmed uniformly since 2000. The Lower Granite discontinuity and possible operational influences need investigation.</p><h2>Methods and limitations</h2>'''+methods+'''<h2>Reproduce and inspect</h2><p>Run <code>python analyze.py</code> with NumPy and Matplotlib. The packet includes returned daily series, computed JSON results, two SVG figures, PNG previews, and example RainbowStats retrieval equations. Inspect the data provenance before publication.</p><h2>Sources</h2>'''+sources+'</main></html>'
(ROOT/'research_memo.html').write_text(report)
print(json.dumps({'complete_years':years,'summary':summary},indent=2))