"""Independent Numpy/SciPy baseline; does not call or translate the browser engine."""
from pathlib import Path
import json,math,platform
import numpy as np
from scipy import stats
R=Path(__file__).resolve().parents[1]
B=json.loads((R/'spatial-statistics/data/noksan.json').read_text())['buildings'];SEED=20261113

def aggregate(size=250,offset=0,cutoff=2000,min_count=5,postal=False):
 groups={}
 for i,b in enumerate(B):
  col=math.floor((b['center'][0]-offset)/size);row=math.floor((b['center'][1]-offset)/size);key=('Z'+str(b['z'])) if postal else (row,col)
  groups.setdefault(key,[]).append(i)
 out=[]
 for key,ids in sorted(groups.items()):
  if len(ids)<min_count:continue
  years=[B[i]['y'] for i in ids if B[i]['y'] is not None and B[i]['y']>0]
  u={'id':key if postal else f'{key[1]},{key[0]}','indices':ids,'count':len(ids),'known':len(years),'share':sum(y<=cutoff for y in years)/len(years) if years else None,'height':float(np.mean([B[i]['h'] for i in ids]))}
  if not postal:u.update(col=key[1],row=key[0],center=[(key[1]+.5)*size+offset,(key[0]+.5)*size+offset])
  out.append(u)
 return out

def weight(u,mode):
 n=len(u);A=np.zeros((n,n))
 for i,a in enumerate(u):
  for j,b in enumerate(u):
   if i==j:continue
   dx=abs(a['col']-b['col']);dy=abs(a['row']-b['row'])
   A[i,j]= math.dist(a['center'],b['center'])<=500 if mode=='distance' else (dx+dy==1 if mode=='rook' else max(dx,dy)==1)
 row=A.sum(axis=1);W=np.divide(A,row[:,None],out=np.zeros_like(A),where=row[:,None]!=0)
 return A,W

def evaluate(x,A,W):
 n=len(x);z=x-x.mean();m2=float(np.mean(z*z));lag=W@z;S0=float(W.sum());den=float(z@z)
 local=np.divide(z*lag,m2,out=np.full(n,np.nan),where=m2>0);local[A.sum(axis=1)==0]=np.nan
 I=n/S0*float(z@W@z)/den if S0 and den else None
 star=A+np.eye(n);k=star.sum(axis=1);gden=np.sqrt(m2)*np.sqrt((n*(star*star).sum(axis=1)-k*k)/(n-1))
 gi=np.divide(star@x-x.mean()*k,gden,out=np.full(n,np.nan),where=gden>0);gi[A.sum(axis=1)==0]=np.nan
 return {'n':n,'mean':float(x.mean()),'m2':m2,'lag':lag.tolist(),'z':z.tolist(),'S0':S0,'expected':-1/(n-1),'I':I,'local':local.tolist(),'gi':gi.tolist()}

def uniforms():
 s=SEED
 while True:
  s=(s^((s<<13)&0xffffffff))&0xffffffff;s=(s^(s>>17))&0xffffffff;s=(s^((s<<5)&0xffffffff))&0xffffffff
  yield s/4294967296

def baseline(x,A,W,Rp=199):
 o=evaluate(x,A,W);rng=uniforms();lc=np.zeros(len(x));gc=np.zeros(len(x));larger=0;sim=[];ex=o['expected']
 for b in range(Rp):
  a=x.copy()
  for i in range(len(a)-1,0,-1):
   j=int(next(rng)*(i+1));a[i],a[j]=a[j],a[i]
  s=evaluate(a,A,W);sim.append(s['I']);larger+=abs(s['I']-ex)>=abs(o['I']-ex)-1e-12
  lc+=np.abs(np.array(s['local'])-ex)>=np.abs(np.array(o['local'])-ex)-1e-12
  gc+=np.abs(s['gi'])>=np.abs(o['gi'])-1e-12
 o.update(p=(int(larger)+1)/(Rp+1),localP=np.where(np.isfinite(o['local']),(lc+1)/(Rp+1),np.nan).tolist(),giP=np.where(np.isfinite(o['gi']),(gc+1)/(Rp+1),np.nan).tolist(),simulation=sim,permutations=Rp,seed=SEED)
 for kind in ['localP','giP']:
  ps=np.array(o[kind]);mask=np.isfinite(ps)
  for method in ['bh','by']:
   out=np.full(len(x),np.nan);out[mask]=stats.false_discovery_control(ps[mask],method=method);o[kind+'_'+method]=out.tolist()
 return o

cases=[]
for size in [250,500]:
 for mode in ['rook','queen','distance']:
  u=aggregate(size);A,W=weight(u,mode)
  for metric in ['count','share']:
   x=np.array([v[metric] for v in u],dtype=float)
   cases.append({'size':size,'mode':mode,'metric':metric,'ids':[v['id'] for v in u],'result':baseline(x,A,W)})
u=aggregate(250);A,W=weight(u,'queen');cases.append({'size':250,'mode':'queen','metric':'count','ids':[v['id'] for v in u],'result':baseline(np.array([v['count'] for v in u],float),A,W,999)})
# MAUP: rank unit counts then attach each unit rank back to its observed buildings.
mr=[];attachments=[]
for id,size,offset,postal in [('postal',250,0,True),('250',250,0,False),('500',500,0,False),('1000',1000,0,False),('500+250',500,250,False)]:
 u=aggregate(size,offset,min_count=1,postal=postal);counts=np.array([v['count'] for v in u]);att=np.empty(len(B));rank=(stats.rankdata(counts,method='average')-.5)/len(u)
 for j,v in enumerate(u):att[v['indices']]=rank[j]
 attachments.append(att)
 rh=stats.pearsonr(counts,[v['height'] for v in u]).statistic;valid=[v for v in u if v['share'] is not None];ra=stats.pearsonr([v['count'] for v in valid],[v['share'] for v in valid]).statistic
 mr.append({'id':id,'units':len(u),'total':int(counts.sum()),'min':int(counts.min()),'max':int(counts.max()),'rHeight':float(rh),'rAge':float(ra)})
for row,att in zip(mr,attachments):row['rankRho']=float(stats.spearmanr(att,attachments[1]).statistic)
def clean(v):
 if isinstance(v,dict):return {k:clean(x) for k,x in v.items()}
 if isinstance(v,list):return [clean(x) for x in v]
 if isinstance(v,float) and not math.isfinite(v):return None
 return v
out=clean({'method':'Independent NumPy dense matrices / SciPy FDR, Pearson and Spearman; deterministic total random-label permutations','seed':SEED,'cases':cases,'maup':mr,'versions':{'python':platform.python_version(),'numpy':np.__version__},'constant1996':all(v['share']==0 for v in aggregate(cutoff=1996))})
(R/'spatial-statistics/data/reference.json').write_text(json.dumps(out,separators=(',',':')))
summary={'status':'reference-generated','independentCases':len(cases),'n250':len(aggregate()),'n500':len(aggregate(500)),'maup':mr,'seed':SEED,'year30BoundaryAllZero':out['constant1996']}
(R/'evidence/G4-validation.json').write_text(json.dumps(summary,indent=2));print(json.dumps(summary,indent=2))