"""Questionnaire-routed KIHBS employment and household poverty; aggregates only."""
from pathlib import Path
import json, hashlib, sqlite3
from datetime import datetime,timezone
import numpy as np
import pandas as pd
import pyreadstat
R=Path(__file__).resolve().parents[1];B=R/'data/raw/kihbs2015_16'
m,meta=pyreadstat.read_dta(str(B/'HH_Members_Information.dta'),encoding='latin1')
h,_=pyreadstat.read_dta(str(B/'HH_Information.dta'))
w,_=pyreadstat.read_dta(str(B/'Consumption_aggregate.dta'))
x=m.merge(h[['clid','hhid','weight','strat','resid']],on=['clid','hhid'],validate='many_to_one').merge(w[['clid','hhid','padqexp','pl_abs']],on=['clid','hhid'],validate='many_to_one')
assert len(x)==len(m)
assert x[['weight','strat','clid','padqexp','pl_abs']].notna().all().all()
age=x.b05_yy.between(15,64)
screens=[f'd02_{i}' for i in range(1,7)]
assert x.loc[age,screens].stack().dropna().isin([1,2]).all()
classification_observed=x[screens].isin([1,2]).all(axis=1)
work=x[screens].eq(1).any(axis=1)
held=x[['d04_1','d04_2','d04_3']].isin(list('ABCDEF')).any(axis=1)
# Questionnaire pp11-12: work -> D08; held activity -> return timing;
# return timing4/5 -> job-search section. D06 has no printed exclusion rule.
x['employed']=work | (held & x.d07.isin([1,2,3]))
x['poor']=x.padqexp.lt(x.pl_abs)
bench=[]
def compare(label,mask,pub_rate,pub_total):
 den=x.loc[mask,'weight'].sum();total=x.loc[mask&x.employed,'weight'].sum()/1000;rate=100*total*1000/den
 bench.append(dict(group=label,employed_thousands=total,published_thousands=pub_total,rate=rate,published_rate=pub_rate,rate_match=round(rate,1)==pub_rate,total_match=round(total,1)==pub_total))
rates=[27.6,57.1,81.7,89.5,92.1,93.6,94.1,93.5,91.6,85.0]
totals=[1384.1,2321.6,2997.2,2741.8,2305.0,1878.1,1388.4,1160.3,1002.6,696.7]
for lo,rate,total in zip(range(15,65,5),rates,totals):compare(f'{lo}-{lo+4}',x.b05_yy.between(lo,lo+4),rate,total)
compare('All15-64',age,71.6,17875.7)
compare('Rural15-64',age&x.resid.eq(1),73.2,10826.0)
compare('Urban15-64',age&x.resid.eq(2),69.4,7049.6)
assert all(v['rate_match'] and v['total_match'] for v in bench)
# Ultimate-cluster, with-replacement Taylor approximation. Keep every sampled PSU,
# including zero contributions outside each domain. No invented FPC.
psus=x[['strat','clid']].drop_duplicates();counts=psus.groupby('strat').size()
assert (counts>1).all() and x.groupby('clid').strat.nunique().max()==1
rows=[]
def estimate(label,domain,definition='Native questionnaire route'):
 valid=domain & x.poor.notna();n=int(valid.sum());den=x.loc[valid,'weight'].sum();num=x.loc[valid&x.poor,'weight'].sum();p=num/den
 residual=x.weight*valid.astype(float)*(x.poor.astype(float)-p)/den
 z=psus.merge(x[['strat','clid']].assign(z=residual).groupby(['strat','clid'],as_index=False).z.sum(),on=['strat','clid'],validate='one_to_one')
 variance=sum(len(g)/(len(g)-1)*((g.z-g.z.mean())**2).sum() for _,g in z.groupby('strat'))
 se=float(np.sqrt(variance));row=dict(group=label,definition=definition,valid_n=n,weighted_people=float(den),weighted_poor_people=float(num),poverty_percent=float(p*100),se_pp=se*100,ci95_low=max(0,p-1.96*se)*100,ci95_high=min(1,p+1.96*se)*100,uncertainty='Approximate95% normal interval; stratified PSU Taylor WR; no FPC',period='2015/16',universe='Age15-64; household consumption poverty')
 assert n>=50
 rows.append(row)
 return valid
masks={}
for label,mask in [('Employed',age&x.employed),('Not employed, observed screens',age&~x.employed&classification_observed),('Employed rural',age&x.employed&x.resid.eq(1)),('Employed urban',age&x.employed&x.resid.eq(2)),('Employed male',age&x.employed&x.b04.eq(1)),('Employed female',age&x.employed&x.b04.eq(2))]:masks[label]=estimate(label,mask)
for label,codes in [('Paid employee',[1,2]),('Own-account worker',[4]),('Contributing family worker',[6])]:masks[label]=estimate(label,age&x.employed&x.d10_p.isin(codes))
masks['Employed excluding apprentice volunteer other']=estimate('Employed excluding apprentice volunteer other',age&x.employed&x.d10_p.isin([1,2,3,4,5,6]),'Sensitivity only; not harmonised ILO employment')
# Independent SQL weighted proportion recomputation (in memory, no respondent export).
with sqlite3.connect(':memory:') as db:
 audit=x[['weight','padqexp','pl_abs']].copy()
 for i,(label,mask) in enumerate(masks.items()):audit[f'd{i}']=mask.astype(int)
 audit.to_sql('audit',db,index=False)
 for i,row in enumerate(rows):
  sql=db.execute(f'SELECT COUNT(*),SUM(weight),SUM(CASE WHEN padqexp<pl_abs THEN weight ELSE 0 END) FROM audit WHERE d{i}=1').fetchone()
  assert sql[0]==row['valid_n'] and abs(100*sql[2]/sql[1]-row['poverty_percent'])<1e-9
pd.DataFrame(rows).to_csv(R/'data/processed/kihbs_work_poverty.csv',index=False)
with sqlite3.connect(R/'data/processed/beyond_employment.sqlite') as db:
 pd.DataFrame(rows).to_sql('kihbs_work_poverty',db,if_exists='replace',index=False)
 pd.DataFrame(bench).to_sql('kihbs_employment_benchmarks',db,if_exists='replace',index=False)
report={'checked_at_utc':datetime.now(timezone.utc).isoformat(),'benchmarks':bench,'benchmark_reference':'Labour Force Basic Report Table3.3, PDF page35','routing_reference':'KIHBS member questionnaire printed pp11-12/PDF12-13','unclassified_missing_work_screens_age15_64':int((age&~classification_observed&~x.employed).sum()),'proxy_disagreements_age15_64':int((age&(x.employed!=x.d10_p.notna())).sum()),'design':{'strata':len(counts),'psus':len(psus),'minimum_psus_per_stratum':int(counts.min()),'method':'Stratified ultimate-cluster Taylor WR normal95%; all PSUs retained for domains; no FPC; approximation not official KNBS interval'},'script_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),'point_estimates_independently_checked_SQL':len(rows),'classification_limit':'Native2015/16 route includes farming, unpaid family work, apprenticeship and volunteering; not modern harmonised employment. Not employed is not synonymous with unemployed.'}
(R/'outputs/kihbs_work_poverty_validation.json').write_text(json.dumps(report,indent=2,default=lambda v:v.item())+'\n')
print(pd.DataFrame(rows)[['group','valid_n','poverty_percent','ci95_low','ci95_high']].round(2).to_string(index=False))
print('PASS: 26 employment benchmark checks and 10 SQL cross-checks.')