"""Independent numerical audit from raw CSV/TSV using Python + NumPy only."""
import csv, json, math, pathlib, hashlib
from collections import defaultdict
import numpy as np
ROOT=pathlib.Path(__file__).resolve().parents[1]
def rows(path,delimiter=','):
 with open(ROOT/path,newline='') as f:return list(csv.DictReader(f,delimiter=delimiter))
def num(x):
 try:return float(x)
 except (ValueError,TypeError):return float('nan')
def finitepos(x):return math.isfinite(num(x)) and num(x)>0
ar=rows('data/raw/anage_anage_data.txt','\t')
a=[r for r in ar if r['Class']=='Mammalia' and all(finitepos(r[c]) for c in ['Female maturity (days)','Adult weight (g)','Maximum longevity (yrs)']) and r['Order'] and r['Family'] and r['Data quality'] in ['acceptable','high']]
orders=sorted({r['Order'] for r in a})
X=np.array([[1,math.log(num(r['Female maturity (days)'])),math.log(num(r['Adult weight (g)']))]+[float(r['Order']==o) for o in orders[1:]] for r in a]);y=np.log([num(r['Maximum longevity (yrs)']) for r in a])
beta=np.linalg.lstsq(X,y,rcond=None)[0];res=y-X@beta;groups=defaultdict(list)
for i,r in enumerate(a):groups[r['Family']].append(i)
scores=np.array([(X[idx]*res[idx,None]).sum(axis=0) for idx in groups.values()]);n,k=X.shape;g=len(groups);inv=np.linalg.inv(X.T@X);V=inv@scores.T@scores@inv*g/(g-1)*(n-1)/(n-k);se=float(np.sqrt(V[1,1]))
rm=next(r for r in rows('results/anage_models.csv') if r['model']=='primary')
assert n==int(rm['n']) and abs(beta[1]-float(rm['beta']))<1e-10 and abs(se-float(rm['se']))<1e-10
raw=rows('data/raw/drugage_drugage.csv');seen=set();d=[]
for row in raw:
 key=tuple(row.items())
 if key not in seen:d.append(row);seen.add(key)
v=[r for r in d if math.isfinite(num(r['avg_lifespan_change_percent'])) and math.isfinite(num(r['max_lifespan_change_percent'])) and r['pubmed_id']]
pg=defaultdict(list)
for r in v:pg[r['pubmed_id']].append(num(r['max_lifespan_change_percent'])-num(r['avg_lifespan_change_percent']))
u=np.array([sum(x)/len(x) for x in pg.values()]);dr=rows('results/drugage_estimates.csv')[0]
assert len(u)==int(dr['units']) and abs(u.mean()-float(dr['estimate']))<1e-12
assert abs(u.std(ddof=1)/math.sqrt(len(u))-float(dr['se']))<1e-12
rng=np.random.default_rng(20260905);boot=np.array([rng.choice(u,len(u),replace=True).mean() for _ in range(10000)]);ci=np.quantile(boot,[.025,.975]);assert max(abs(ci-np.array([float(dr['boot_lo']),float(dr['boot_hi'])])))<.15
manifest=json.loads((ROOT/'data/manifest.json').read_text())
for record in manifest:
 assert hashlib.sha256((ROOT/'data/raw'/f"{record['dataset']}.zip").read_bytes()).hexdigest()==record['sha256']
 for f in record['files']:assert hashlib.sha256((ROOT/f['path']).read_bytes()).hexdigest()==f['sha256']
report=dict(status='PASS',verification='Independent raw-data parsing, NumPy least-squares and cluster sandwich; publication aggregation; independent bootstrap RNG; all frozen data SHA-256 hashes.',anage_raw=len(ar),anage_n=n,anage_beta=float(beta[1]),anage_cluster_se=se,drugage_raw=len(raw),drugage_unique=len(d),drugage_paired=len(v),drugage_publications=len(u),drugage_estimate=float(u.mean()),drugage_bootstrap_independent_ci=ci.tolist(),numpy=np.__version__)
(ROOT/'results/independent_verification.json').write_text(json.dumps(report,indent=2)+'\n');print(json.dumps(report,indent=2))
