#!/usr/bin/env python3
"""Reproduce the Moon audit: its two tests, the drawer-evenness check and, with --sensitivity,
the smallest lunar swing the count test could have heard. Requires Python 3 and numpy.
Usage: python reproduce.py [moon-days.json] [--sensitivity]
If no path is given, downloads the museum's current public daily series.
The dark-sky baseline needs exact coordinates and cannot be reproduced from public data.
"""
import json,sys,urllib.request
from collections import defaultdict
import numpy as np
args=[a for a in sys.argv[1:] if not a.startswith('--')]; want_sensitivity='--sensitivity' in sys.argv
if args:
    with open(args[0]) as f: data=json.load(f)
else:
    with urllib.request.urlopen('https://8museum.com/data/moon-days.json') as r: data=json.load(r)
y=np.array([d['count'] for d in data],dtype=float)
deg=np.array([d['phase'] for d in data]); a=np.radians(deg); x=np.column_stack((np.cos(a),np.sin(a)))
months=defaultdict(list)
for i,d in enumerate(data): months[d['date'][:7]].append(i)
groups=[np.array(ix) for ix in months.values()]
for ix in groups: x[ix]-=x[ix].mean(axis=0)
pinv=np.linalg.pinv(x)
def stat(v):
    z=v.copy()
    for ix in groups: z[ix]-=z[ix].mean()
    fit=x@(pinv@z)
    return float(fit@fit/(z@z)) if z@z else 0.0
bins=((deg+22.5)%360//45).astype(int); masks=[bins==i for i in range(8)]
def drawer_ratio(v):
    r=np.array([v[m].mean() for m in masks]); return float(r.max()/r.min()) if r.min()>0 else float('inf')
obs=np.array([stat(y),stat((y>0).astype(float))]); ratio_obs=drawer_ratio(y)
rng=np.random.default_rng(888); exceed=np.zeros(2); ratio_null=np.empty(9999)
for k in range(9999):
    v=y.copy()
    for ix in groups: v[ix]=np.roll(y[ix],rng.integers(len(ix)))
    exceed+=np.array([stat(v),stat((v>0).astype(float))])>=obs-1e-14; ratio_null[k]=drawer_ratio(v)
p=(1+exceed)/10000; adjusted=np.zeros(2); running=0
for rank,j in enumerate(np.argsort(p)):
    running=max(running,min(1,p[j]*(2-rank)));adjusted[j]=running
for i,name in enumerate(['daily encounter count','day with any encounter']):
    print(f'{name}: partial R2={obs[i]:.6f}, p={p[i]:.4f}, Holm p={adjusted[i]:.4f}')
finite=ratio_null[np.isfinite(ratio_null)]
print(f'drawer ratio (fullest/emptiest per calendar day): observed={ratio_obs:.3f}, null median={np.median(finite):.3f}, null 95th={np.percentile(finite,95):.3f}, p={(1+(ratio_null>=ratio_obs-1e-12).sum())/10000:.4f}')
if want_sensitivity:
    # Multiply each day's count by 1 + A cos(phase - theta), theta uniform, rescale to the archive total,
    # and apply the identical shift test to the count endpoint. Same seed and draw order as the museum's build.
    srng=np.random.default_rng(8888); INJ=80; SHIFTS=200; grid=[round(float(g),1) for g in np.arange(0.1,1.01,0.1)]
    def shift_p(v,n):
        o=stat(v); hits=0
        for _ in range(n):
            s=v.copy()
            for ix in groups: s[ix]=np.roll(v[ix],srng.integers(len(ix)))
            hits+=stat(s)>=o-1e-14
        return (1+hits)/(n+1)
    for A in grid:
        detected=0
        for _ in range(INJ):
            theta=srng.uniform(0,2*np.pi); v=y*(1+A*np.cos(a-theta)); v*=y.sum()/v.sum()
            detected+=shift_p(v,SHIFTS)<.05
        print(f'swing ±{A*100:.0f}%: detected {detected/INJ:.0%} of the time')