#!/usr/bin/env python3
# intervals.py: the fits and 95% intervals the page quotes, computed independently of the page's JS
# engine (numpy + scipy), from digitized.json. Writes intervals.json, which the verifier embeds.
import json, numpy as np
from scipy import stats
d = json.load(open('digitized.json'))
pre = [(p['depth'], p['osm']) for p in d['pre2014']['points']]
k = [(7000, 991)]
lat = [(p['depth'], p['osm']) for p in d['later']['points']]
oth = [(p['depth'], p['osm']) for p in d['others']['points']]
sets = dict(before=pre, kermadec=pre + k, snailfish=pre + k + lat, snailonly=k + lat, all=pre + k + lat + oth)
out = dict(t975={str(df): float(stats.t.ppf(0.975, df)) for df in range(1, 41)}, sets={})
for name, pts in sets.items():
    x = np.array([p[0] for p in pts], float); y = np.array([p[1] for p in pts], float); n = len(x)
    res = stats.linregress(x, y)
    a, b = res.intercept, res.slope
    s = np.sqrt(((y - a - b * x) ** 2).sum() / (n - 2)); xb = x.mean(); sxx = ((x - xb) ** 2).sum(); t = stats.t.ppf(.975, n - 2)
    # Fieller by a fine scan, not the quadratic: every depth where 1,100 lies inside the 95% band
    grid = np.arange(5000, 12000.0001, 0.25)
    inside = np.abs(a + b * grid - 1100) <= t * s * np.sqrt(1 / n + (grid - xb) ** 2 / sxx)
    out['sets'][name] = dict(n=n, a=a, b=b, crossing=(1100 - a) / b, lo=float(grid[inside].min()), hi=float(grid[inside].max()))
F = out['sets']['before']
x = np.array([p[0] for p in pre], float); y = np.array([p[1] for p in pre], float); n = len(x)
s = np.sqrt(((y - F['a'] - F['b'] * x) ** 2).sum() / (n - 2)); sxx = ((x - x.mean()) ** 2).sum(); t = stats.t.ppf(.975, n - 2)
pi = lambda z: t * s * np.sqrt(1 + 1 / n + (z - x.mean()) ** 2 / sxx)
out['predict_before'] = [dict(depth=z, osm=F['a'] + F['b'] * z, pm=float(pi(z))) for z in [7000] + [p[0] for p in lat]]
json.dump(out, open('intervals.json', 'w'), indent=1)
for kname, v in out['sets'].items(): print(kname, v['n'], round(v['a'], 1), round(v['b'], 4), round(v['crossing']), round(v['lo']), round(v['hi']))
print(out['predict_before'][0])