"""Chile ENPCCL2024: independent domain analysis, no participant-level export.
Run: python analysis.py path/to/enpccl_puf.dta output_directory
Dependencies: pandas numpy scipy. Public data from Chile's Ministry of Cultures.
"""
from pathlib import Path
import sys, json, hashlib
import pandas as pd
import numpy as np
from scipy.stats import t

def analyze(source, out):
    out=Path(out);out.mkdir(parents=True,exist_ok=True)
    cols=['kish_dummy','rph_edad','rph_sexo','region','lectura_auto','lectura_libros','lectura_tipo_3','libros_ocio','Estrato_ENPC','Varunit','fexp_p']
    raw=pd.read_stata(source,columns=cols,convert_categoricals=False)
    d=raw.loc[raw.kish_dummy.eq(1)].reset_index(drop=True)
    assert d.fexp_p.gt(0).all() and d[['Estrato_ENPC','Varunit']].ne('').all().all()
    adult=d.rph_edad.isin(range(3,9));identity=d.lectura_auto.isin([1,2,3,4])
    books=d.lectura_libros.isin(range(1,7))
    leisure=adult & books & d.lectura_tipo_3.eq(1) & d.libros_ocio.ge(1) & d.libros_ocio.le(500)
    base=leisure & identity
    w=d.fexp_p.to_numpy(); design=d[['Estrato_ENPC','Varunit']].copy()
    universe=design.drop_duplicates(); counts=universe.groupby('Estrato_ENPC').size()
    assert counts.min()>1
    df=int((counts-1).sum());critical=float(t.ppf(.975,df))
    def estimate(domain,event,label,kind='share_poco'):
        dom=np.asarray(domain,dtype=bool);y=np.asarray(event,dtype=bool)
        den=float(w[dom].sum());num=float(w[dom & y].sum());p=num/den
        # Ratio linearization over ALL selected respondents, zero outside domain.
        z=w*dom*(y-p)/den
        a=design.assign(z=z).groupby(['Estrato_ENPC','Varunit']).z.sum()
        variance=0.
        for h,vals in a.groupby(level=0):
            m=len(vals);variance+=m/(m-1)*float(((vals-vals.mean())**2).sum())
        se=variance**.5
        return dict(group=label,metric=kind,n_denominator=int(dom.sum()),n_numerator=int((dom&y).sum()),weighted_denominator=den,weighted_numerator=num,percent=100*p,se_pp=100*se,ci95_lower=100*max(0,p-critical*se),ci95_upper=100*min(1,p+critical*se),design_df=df,relative_se_percent=100*se/p if p else None,small_event_count_lt30=bool((dom&y).sum()<30),publication_flag='caution_low_event_count_or_high_relative_error' if (dom&y).sum()<30 or (se/p if p else 1)>.15 else ('sample_size_sufficient_not_official_certification' if dom.sum()>=60 else 'small_sample_do_not_feature'))
    rows=[estimate(base,d.lectura_auto.eq(1),'Total20+')]
    for code,label in [(1,'Hombres'),(2,'Mujeres')]:rows.append(estimate(base & d.rph_sexo.eq(code),d.lectura_auto.eq(1),label))
    for codes,label in [([3],'20–29'),([4],'30–39'),([5],'40–49'),([6],'50–59'),([7],'60–69'),([8],'Mayor de 70 (etiqueta de la fuente)')]:rows.append(estimate(base & d.rph_edad.isin(codes),d.lectura_auto.eq(1),label))
    for lo,hi,label in [(1,2,'1–2 libros'),(3,5,'3–5 libros'),(6,11,'6–11 libros'),(12,500,'12+ libros')]:rows.append(estimate(base & d.libros_ocio.between(lo,hi),d.lectura_auto.eq(1),label))
    regions={'1':'Tarapacá','2':'Antofagasta','3':'Atacama','4':'Coquimbo','5':'Valparaíso','6':"O’Higgins",'7':'Maule','8':'Biobío','9':'La Araucanía','10':'Los Lagos','11':'Aysén','12':'Magallanes','13':'Metropolitana','14':'Los Ríos','15':'Arica y Parinacota','16':'Ñuble'}
    for code,label in regions.items():rows.append(estimate(base & d.region.eq(code),d.lectura_auto.eq(1),label))
    for code,label in [(1,'Poco lector'),(2,'Lector moderado'),(3,'Lector frecuente'),(4,'Lector muy frecuente')]:rows.append(estimate(base,d.lectura_auto.eq(code),label,'identity_distribution'))
    weekly=base & d.lectura_libros.isin([1,2]);rows.append(estimate(weekly,d.lectura_auto.eq(1),'Libros al menos semanalmente','weekly_poco'))
    counts_out={'full_roster_rows':len(raw),'selected_respondents15plus':len(d),'selected_respondents20plus':int(adult.sum()),'leisure_book_readers20plus':int(leisure.sum()),'valid_identity_leisure20plus':int(base.sum()),'excluded_identity':int((leisure&~identity).sum()),'psus':len(universe),'strata':len(counts),'design_df':df,'source_sha256':hashlib.sha256(Path(source).read_bytes()).hexdigest()}
    pd.DataFrame(rows).to_csv(out/'results.csv',index=False,encoding='utf-8-sig')
    (out/'results.json').write_text(json.dumps({'counts':counts_out,'estimates':rows},ensure_ascii=False,indent=2)+'\n',encoding='utf-8')
    # Audit regrouping identities and explicit scalar reconstruction, no microdata export.
    assert abs(sum(x['percent'] for x in rows if x['metric']=='identity_distribution')-100)<1e-8
    first=rows[0];check=sum(float(row.fexp_p) for row in d.loc[base & d.lectura_auto.eq(1)].itertuples())/sum(float(row.fexp_p) for row in d.loc[base].itertuples())*100
    assert abs(check-first['percent'])<1e-10
    return counts_out,rows

if __name__=='__main__':
    counts,rows=analyze(sys.argv[1],sys.argv[2]);print(json.dumps({'counts':counts,'national':rows[0]},ensure_ascii=False,indent=2))