"""Short-lived named storms ('shorties', Landsea et al. 2010): named-status duration <= 48 h.
Counts per year; trend; and how much of the named-storm under-forecast they account for."""
import csv, numpy as np, datetime, collections
lines = open('../data/hurdat2-1851-2025-092326.txt').read().splitlines(); i = 0
dur = collections.defaultdict(list)
while i < len(lines):
    hdr = [x.strip() for x in lines[i].split(',')]; n = int(hdr[2]); y = int(hdr[0][4:8])
    fx = [[x.strip() for x in l.split(',')] for l in lines[i+1:i+1+n]]; i += 1 + n
    ts = [datetime.datetime.strptime(f[0] + f[1], '%Y%m%d%H%M') for f in fx if f[3] in ('TS', 'SS', 'HU')]
    if ts: dur[y].append((max(ts) - min(ts)).total_seconds() / 3600)
rows = []
for y in range(1951, 2026):
    d = np.array(dur[y]); rows.append((y, len(d), int((d <= 48).sum())))
w = csv.writer(open('../data/shorties.csv', 'w', newline='')); w.writerow(['year', 'ns', 'shorties_le48h']); w.writerows(rows)
for a, b in ((1951, 1970), (1971, 1994), (1995, 2007), (2008, 2025)):
    sel = [r for r in rows if a <= r[0] <= b]
    print(f'{a}-{b}: mean NS {np.mean([r[1] for r in sel]):5.1f}, shorties {np.mean([r[2] for r in sel]):4.1f}, longer-lived {np.mean([r[1]-r[2] for r in sel]):5.1f}')
# NOAA May midpoint bias, with and without shorties
noaa = list(csv.DictReader(open('../data/noaa_outlooks.csv'))); S = {r[0]: r for r in rows}
for m in ('May', 'Aug'):
    b = [(int(r['ns_lo']) + int(r['ns_hi'])) / 2 - S[int(r['year'])][1] for r in noaa if r['month'] == m]
    print(m, f'NOAA NS midpoint bias {np.mean(b):+.2f}; mean shorties 2001-2025 {np.mean([S[y][2] for y in range(2001, 2026)]):.2f}')
# Excess shorties over the 1951-1994 level vs NOAA May named-storm shortfall
long_ = {y: n - s for y, n, s in rows}; alln = {y: n for y, n, s in rows}
sh_base = np.mean([alln[y] - long_[y] for y in range(1951, 1995)])
err = np.array([alln[int(r['year'])] - (int(r['ns_lo']) + int(r['ns_hi'])) / 2 for r in noaa if r['month'] == 'May'])
exc = np.array([alln[int(r['year'])] - long_[int(r['year'])] - sh_base for r in noaa if r['month'] == 'May'])
print('May NOAA NS shortfall %.2f; excess shorties %.2f; corr %.2f (n=%d)' % (err.mean(), exc.mean(), np.corrcoef(err, exc)[0, 1], len(err)))