import sys, warnings; warnings.filterwarnings("ignore")
import pandas as pd, numpy as np
from scipy import stats
import statsmodels.formula.api as smf
d = pd.read_csv(sys.argv[1], sep=None, engine="python", encoding="utf-8-sig")
y = d["AÑO"]; d = d.apply(lambda s: pd.to_numeric(s.astype(str).str.strip().str.replace(",", "."), errors="coerce")); d["YEAR"] = y.astype(int)
def alpha(df):
    df = df.dropna(); k = df.shape[1]; return k/(k-1)*(1-df.var().sum()/df.sum(axis=1).var()), len(df)
A = ["SOS_ACUERDO_CC","SOS_ACUERDO_CIRC","SOS_ACUERDO_PRODUCTOS","SOS_ACUERDO_COMBUSTIBLES"]
K = ["SOS_ACUERDO_SECTORES","SOS_ACUERDO_GEST_RESIDUOS","SOS_ACUERDO_BIOGAS","SOS_ACUERDO_PLANTAS"]
print("alpha ACT", alpha(d[A]), "check mean==ACT", np.corrcoef(d[A].mean(axis=1).fillna(0), d.ACT_SOST.fillna(0))[0,1])
print("alpha CONOC", alpha(d[K]), np.corrcoef(d[K].mean(axis=1).fillna(0), d.CONOC_SOST.fillna(0))[0,1])
print("ACT nonmissing", d.ACT_SOST.notna().sum(), "rows any A", d[A].notna().any(axis=1).sum(), "all A", d[A].notna().all(axis=1).sum())
Q = ["ACUERDO_BIOREF_SOSTENIBLES","ACUERDO_BIOREF_CONFIANZA","ACUERDO_BIOREF_PAGARIA_MAS","ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS"]
for combo in [Q[:3], Q, Q[:2]]:
    m = d[combo].mean(axis=1); print("ACEPT combo", combo, round(np.corrcoef(m[d.ACEPT_BIOREF.notna()], d.ACEPT_BIOREF.dropna())[0,1],4), round(m.mean(),3))
print("alpha Q17 4 pos items", alpha(d[Q]))
print("ACEPT_BIOREF mean", d.ACEPT_BIOREF.mean())
# PROB by income
d["PROB"] = d.PROB_USO_BIOCOM.where(d.PROB_USO_BIOCOM >= 1)
g = [x.dropna() for _, x in d.groupby("INGRESOS").PROB]; print("PROB by income F", stats.f_oneway(*g), d.groupby("INGRESOS").PROB.agg(["mean","count"]).round(2).to_dict())
# Q10 shares top3 per wave
R10 = [c for c in d.columns if c.startswith("COMB_RANK")]
print("Q10 range", d[R10].min().min(), d[R10].max().max())
for c in R10:
    print(c, {yr: round((s<=3).sum()/s.notna().sum()*100,1) for yr, s in d.groupby("YEAR")[c]}, "INV mean", d.groupby("YEAR")["INV_"+c].mean().round(2).to_dict() if "INV_"+c in d else "")
T3 = [c for c in d.columns if c.startswith("TOP3")]
for c in T3:
    ct = pd.crosstab(d.YEAR, d[c]); print(c, d.groupby("YEAR")[c].mean().mul(100).round(1).to_dict(), "p", round(stats.chi2_contingency(ct)[1],3))
# gender effect on ACT
dd = d[d.GENERO.isin([1,2])]
print("ACT by gender", dd.groupby("GENERO").ACT_SOST.mean().round(2).to_dict(), stats.ttest_ind(dd[dd.GENERO==1].ACT_SOST.dropna(), dd[dd.GENERO==2].ACT_SOST.dropna(), equal_var=False))
# PROB by cluster
print("PROB/ dist by cluster", d.groupby("QCL_1")[["PROB","DISTANCIA_MIN_BIOREF","ACUERDO_BIOREF_NO_CERCA","ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS"]].mean().round(2))
print("cluster gender % women", pd.crosstab(d.QCL_1, d.GENERO, normalize="index").round(3))
# Q16 aid wave
a=d[d.YEAR==2022].RAZC_RANK_NO_SIN_AYUDAS_INV.dropna(); b=d[d.YEAR==2024].RAZC_RANK_NO_SIN_AYUDAS_INV.dropna(); print("Q16 aid", stats.ttest_ind(a,b,equal_var=False), len(a), len(b))
# distance: tobit-like check: share 30 vs benefits
print("share30 by MOLESTIAS_BENEFICIOS", d.groupby("ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS").DISTANCIA_MIN_BIOREF.agg(["mean","count", lambda s: round((s==30).mean()*100,1)]).round(2))
print("dist by NINGUN", d.groupby("REQ_NINGUN_CASO").DISTANCIA_MIN_BIOREF.agg(["mean","count"]).round(2), stats.mannwhitneyu(d[d.REQ_NINGUN_CASO==1].DISTANCIA_MIN_BIOREF.dropna(), d[d.REQ_NINGUN_CASO==0].DISTANCIA_MIN_BIOREF.dropna()))
# Wave x income interaction on PROB? and wave x age on PROB
x = d[["PROB","YEAR","EDAD_CAT","GENERO","NIVEL_ESTUDIOS","INGRESOS"]].dropna(); x["W24"]=(x.YEAR==2024).astype(int)
m0 = smf.ols("PROB ~ W24 + C(EDAD_CAT) + C(GENERO) + C(NIVEL_ESTUDIOS) + C(INGRESOS)", x).fit()
m1 = smf.ols("PROB ~ W24*C(EDAD_CAT) + C(GENERO) + C(NIVEL_ESTUDIOS) + C(INGRESOS)", x).fit()
print("wave x age on PROB LR p", m1.compare_f_test(m0))
# age-standardised PROB: reweight 2024 to 2022 age x gender structure
w = d[d.GENERO.isin([1,2])].copy(); cell = w.groupby(["YEAR","EDAD_CAT","GENERO"]).size().unstack([1,2]).fillna(0)
p22 = cell.loc[2022]/cell.loc[2022].sum(); pool = cell.sum()/cell.sum().sum()
def stdmean(var, yr, ref):
    s = w[w.YEAR==yr].groupby(["EDAD_CAT","GENERO"])[var].mean(); return (s*ref).sum()/ref[s.index].sum()
for v in ["PROB","ACT_SOST","ACUERDO_BIOREF_PAGARIA_MAS","DISTANCIA_MIN_BIOREF"]:
    print("std(pooled age x gender)", v, round(stdmean(v,2022,pool),3), round(stdmean(v,2024,pool),3))
print("missing pattern completeness Q9-Q19 by wave", d.groupby("YEAR").apply(lambda g: g[["ACT_SOST","PROB","DISTANCIA_MIN_BIOREF","ACEPT_BIOREF"]].notna().all(axis=1).mean()).round(3).to_dict())