"""Reproduce Pluraterra's exploratory diagnostics from a pinned public evidence release.
Usage: python3 reproduce.py [directory-containing-release-json-files]
With no argument, retrieves only the six required public snapshot files.
All calculations use the captured 2026-11-02-r3 release, not live source revisions.
"""
import json,math,statistics,hashlib,sys
from pathlib import Path
from urllib.request import urlopen
RELEASE='2026-11-02-r3'
BASE='https://www.pluraterra.com/data/releases/'+RELEASE+'/'
INPUTS=['atlas','households','regions','financial-strain','rankings','robustness']
def analyse(data,hashes):
 D,H,G,F,R,S=[data[k] for k in INPUTS];C=D['countries'];B={c['id']:c for c in C}
 def val(c,k,y=2023):return c['values'].get(k,{}).get(str(y))
 def name(id):return B[id]['name']
 hh={}
 for r in H['rows']:
  scenario=next(k for k,s in H['scenarios'].items() if s['key']==[r['household_type'],r['principal'],r['spouse']])
  hh.setdefault(r['country'],{})[scenario]=r
 pay=[]
 for id,h in hh.items():
  a,b=[h[k]['values'] for k in ['taxcouple100','taxcouple167']];gross=b['GEBT_XDC']-a['GEBT_XDC'];net=b['NIAT_XDC']-a['NIAT_XDC']
  pay.append(dict(country=id,name=name(id),year=2025,yield_percent=100*net/gross,extra_gross_local=gross,extra_net_local=net,one_earner_net_local=a['NIAT_XDC'],two_earner_net_local=b['NIAT_XDC'],one_earner_gross_local=a['GEBT_XDC'],two_earner_gross_local=b['GEBT_XDC']))
 pay.sort(key=lambda r:r['yield_percent'])
 strain={(r['country'],r['household'],r['income_group']):r for r in F['rows']};buffer=[]
 for id in B:
  a,b=[strain.get((id,k,'A_60')) for k in ['A1_DCH','TOTAL']]
  if a is None or b is None:continue
  buffer.append(dict(country=id,name=name(id),year=2024,parent_strain=a['value'],all_strain=b['value'],gap_pp=a['value']-b['value'],source_flags={'parent':a.get('source_status'),'all':b.get('source_status')}))
 buffer.sort(key=lambda r:-r['gap_pp'])
 lift=[]
 for c in C:
  if val(c,'life') is None or not (val(c,'income') or 0)>0:continue
  variants=[]
  for band in [.20,.25,.35]:
   peers=[o for o in C if o['id']!=c['id'] and val(o,'life') is not None and (val(o,'income') or 0)>0 and abs(math.log(val(o,'income')/val(c,'income'))) <=math.log(1+band)]
   median=statistics.median(val(p,'life') for p in peers) if len(peers)>=5 else None
   variants.append(dict(band=band,n=len(peers),peer_ids=[p['id'] for p in peers],peer_median_life=median,lift_years=val(c,'life')-median if median is not None else None))
  baseline=variants[1];defined=[v['lift_years'] for v in variants if v['lift_years'] is not None]
  stable=(all(x>0 for x in defined) or all(x<0 for x in defined)) if len(defined)==3 else None
  lift.append(dict(country=c['id'],name=c['name'],year=2023,income=val(c,'income'),life=val(c,'life'),lift_years=baseline['lift_years'],peer_median_life=baseline['peer_median_life'],peer_n=baseline['n'],same_sign_all_bands=stable,band_sensitivity=variants))
 lift.sort(key=lambda r:-(r['lift_years'] if r['lift_years'] is not None else -1000))
 floor=[]
 for c in C:
  inputs={k:val(c,k) for k in ['water','sanitation','electricity','cleancooking']}
  if any(v is None for v in inputs.values()):continue
  minimum=min(inputs.values());mean=sum(inputs.values())/4
  floor.append(dict(country=c['id'],name=c['name'],year=2023,inputs=inputs,floor_percent=minimum,weakest=[k for k,v in inputs.items() if v==minimum],mean_access=mean,average_floor_gap_pp=mean-minimum,joint_lower_percent=max(0,sum(inputs.values())-300),joint_upper_percent=minimum))
 floor.sort(key=lambda r:-r['average_floor_gap_pp'])
 drift=[dict(country=id,name=name(id),year=2023,best=r['best'],worst=r['worst'],drift_places=r['worst']-r['best'],histogram=r['histogram']) for id,r in S['countries'].items()];drift.sort(key=lambda r:-r['drift_places'])
 common=[r for r in R['rows'] if val(B[r['id']],'lifeevaluation') is not None];lived=[]
 for r in common:
  id=r['id'];v=val(B[id],'lifeevaluation');index_rank=1+sum(o['totals']['balanced']>r['totals']['balanced']+1e-10 for o in common);rating_rank=1+sum(val(B[o['id']],'lifeevaluation')>v for o in common)
  lived.append(dict(country=id,name=name(id),cohort_n=len(common),index_rank=index_rank,life_rating_rank=rating_rank,gap_places=index_rank-rating_rank,life_rating=v,source_interval=B[id].get('intervals',{}).get('lifeevaluation',{}).get('2023'),index_score=r['totals']['balanced'],index_year=2023,survey_start=2021,survey_end=2023))
 lived.sort(key=lambda r:-r['gap_places'])
 health=[]
 for y in range(2015,2025):
  us=B['USA'];pairs=[c for c in C if val(c,'life',y) is not None and (val(c,'income',y) or 0)>0]
  selected=[c for c in pairs if val(c,'income',y)<val(us,'income',y) and val(c,'life',y)>val(us,'life',y)]
  health.append(dict(year=y,cohort_n=len(pairs),count=len(selected),countries=[c['id'] for c in selected],costa_rica_life_gap=val(B['CRI'],'life',y)-val(us,'life',y)))
 health_rows=[dict(country=c['id'],name=c['name'],year=2023,income=val(c,'income'),life=val(c,'life')) for c in C if val(c,'life') is not None and (val(c,'income') or 0)>0]
 belgium=[r for r in G['rows'] if r['country']=='BEL'];brussels=next(r for r in belgium if r['id']=='BE10')
 brussels_ranks={k:1+sum(r['values'].get(k,-1)>brussels['values'][k] for r in belgium) for k in ['gdp','life','employment']}
 return dict(release=RELEASE,published_on='2026-11-02',analysis_version='1',inputs_sha256=hashes,method_status='Exploratory descriptive diagnostics; selected after inspecting available data. No causal effects, preregistered tests or universal welfare index claimed.',models={'paycheck-yield':pay,'buffer-gap':buffer,'longevity-lift':lift,'essentials-floor':floor,'rank-drift':drift,'lived-gap':lived},health={'annual_checks':health,'countries_2023':health_rows},brussels={'country_regions':belgium,'ranks_high_to_low':brussels_ranks})
def load(directory=None):
 data={};hashes={}
 for k in INPUTS:
  b=(Path(directory)/(k+'.json')).read_bytes() if directory else urlopen(BASE+k+'.json',timeout=60).read()
  data[k]=json.loads(b);hashes[k]=hashlib.sha256(b).hexdigest()
 return analyse(data,hashes)
if __name__=='__main__':print(json.dumps(load(sys.argv[1] if len(sys.argv)>1 else None),indent=2,allow_nan=False))