import sys, warnings
warnings.filterwarnings("ignore")
import pandas as pd, numpy as np
from scipy import stats
import statsmodels.formula.api as smf
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score, silhouette_score

path = sys.argv[1]
d = pd.read_csv(path, 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)
d["W24"] = (d.YEAR == 2024).astype(int)

def p(*a): print(*a)

p("== N per wave", d.YEAR.value_counts().to_dict())

# demographics
for c in ["GENERO", "EDAD_CAT", "NIVEL_ESTUDIOS", "SITUACION_LABORAL", "INGRESOS"]:
    ct = pd.crosstab(d.YEAR, d[c])
    chi = stats.chi2_contingency(ct)
    p(f"== {c} chi2={chi[0]:.3f} df={chi[2]} p={chi[1]:.4f} n={ct.values.sum()}")
p(pd.crosstab(d.YEAR, d.EDAD_CAT, normalize="index").round(3))
p(pd.crosstab(d.YEAR, d.INGRESOS, normalize="index").round(3))

# covariates for adjustment
cov = "C(GENERO_B) + C(EDAD_CAT) + C(EDU) + C(INC)"
d["GENERO_B"] = d.GENERO.where(d.GENERO.isin([1, 2]))
d["EDU"] = d.NIVEL_ESTUDIOS
d["INC"] = d.INGRESOS

def cohen(a, b):
    a, b = a.dropna(), b.dropna()
    sp = np.sqrt(((len(a)-1)*a.var()+(len(b)-1)*b.var())/(len(a)+len(b)-2))
    return (a.mean()-b.mean())/sp

def wave(var, label=None):
    a = d.loc[d.YEAR == 2022, var].dropna(); b = d.loc[d.YEAR == 2024, var].dropna()
    t = stats.ttest_ind(a, b, equal_var=False)
    df = d[[var, "W24", "GENERO_B", "EDAD_CAT", "EDU", "INC"]].dropna()
    m = smf.ols(f"{var} ~ W24 + {cov}", data=df).fit(cov_type="HC3")
    ci = m.conf_int().loc["W24"]
    p(f"{label or var}: 2022 M={a.mean():.2f} SD={a.std():.2f} n={len(a)} | 2024 M={b.mean():.2f} SD={b.std():.2f} n={len(b)} | Welch t={t.statistic:.2f} p={t.pvalue:.4f} d={cohen(a,b):.3f} | adj b={m.params['W24']:.3f} [{ci[0]:.3f},{ci[1]:.3f}] p={m.pvalues['W24']:.4f} n={int(m.nobs)}")

p("\n== WAVE COMPARISONS (raw and adjusted for gender, age, education, income)")
d["PROB"] = d.PROB_USO_BIOCOM.where(d.PROB_USO_BIOCOM >= 1)  # 0 = invalid
p("PROB zeros:", int((d.PROB_USO_BIOCOM == 0).sum()))
for v in ["ACT_SOST", "CONOC_SOST", "SOS_ACUERDO_CC", "SOS_ACUERDO_CIRC", "PROB", "PROB_USO_BIOCOM",
          "ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS", "ACUERDO_BIOREF_NO_CERCA", "ACUERDO_BIOREF_SOSTENIBLES",
          "ACUERDO_BIOREF_CONFIANZA", "ACUERDO_BIOREF_PAGARIA_MAS", "ACEPT_BIOREF", "DISTANCIA_MIN_BIOREF",
          "ACUERDO_SEPARAR", "ACUERDO_COMPOSTAJE", "ACUERDO_VERTIDO", "ACUERDO_INCINERACION", "ACUERDO_PIROLISIS",
          "TRAT_RANK_COMPOSTAJE_INV", "TRAT_RANK_PIROLISIS_INV", "TRAT_RANK_VERTIDO_INV", "TRAT_RANK_INCINERACION_INV",
          "PRED_BIOCOM_AMB", "PRED_BIOCOM_ECON", "PRED_BIOCOM_EFIC", "PRED_BIOCOM_SEG"]:
    wave(v)

p("\n== Q12 rank check (raw rank 1 = top priority); share ranking each first")
for c in ["TRAT_RANK_COMPOSTAJE", "TRAT_RANK_PIROLISIS", "TRAT_RANK_VERTIDO", "TRAT_RANK_ICINERACION"]:
    p(c, d.groupby("YEAR")[c].apply(lambda s: round((s == 1).mean()*100, 1)).to_dict(), "mean raw", d.groupby("YEAR")[c].mean().round(2).to_dict())

p("\n== Q14 distribution"); p(pd.crosstab(d.YEAR, d.PROB_USO_BIOCOM))
p("\n== Q18 distance"); p(d.DISTANCIA_MIN_BIOREF.describe().round(2).to_dict(), "share=30:", round((d.DISTANCIA_MIN_BIOREF == 30).mean()*100, 1), "share<=5:", round((d.DISTANCIA_MIN_BIOREF <= 5).mean()*100, 1))
p("median by wave", d.groupby("YEAR").DISTANCIA_MIN_BIOREF.median().to_dict(), "share30 by wave", d.groupby("YEAR").DISTANCIA_MIN_BIOREF.apply(lambda s: round((s == 30).mean()*100, 1)).to_dict())
mw = stats.mannwhitneyu(d.loc[d.YEAR == 2022, "DISTANCIA_MIN_BIOREF"].dropna(), d.loc[d.YEAR == 2024, "DISTANCIA_MIN_BIOREF"].dropna())
p("MW distance p=", round(mw.pvalue, 4))

p("\n== Q19 requirements % by wave (chi2)")
for c in ["REQ_OLORDESIN", "REQ_EMPLEO_LOCAL", "REQ_COLAB_INDUSTRIA", "REQ_IMPACTO_VISUAL", "REQ_DEMANDA_LOCAL", "REQ_CIRCULACION_CAM", "REQ_SEGURIDAD", "REQ_NO_AFECTAR_PRECIO", "REQ_NINGUN_CASO"]:
    ct = pd.crosstab(d.YEAR, d[c]);
    pct = d.groupby("YEAR")[c].mean().mul(100).round(1).to_dict(); tot = round(d[c].mean()*100, 1)
    pv = stats.chi2_contingency(ct)[1] if ct.shape[1] > 1 and c != "REQ_NO_AFECTAR_PRECIO" else np.nan
    p(c, pct, "pooled", tot, "p", round(pv, 4) if pv == pv else "n/a")

p("\n== Q20 concerns")
P5 = ["PREOC_OLORES", "PREOC_IMPACTO_VISUAL", "PREOC_CAMIONES", "PREOC_INSEGURIDAD", "PREOC_OBRAS"]
P8 = P5 + ["PREOC_CONTAMINACION", "PREOC_SALUD", "PREOC_VALOR_VIV"]
s22 = d.loc[d.YEAR == 2022, P5].sum(axis=1); s24 = d.loc[d.YEAR == 2024, P8].sum(axis=1)
p("2022 row sums", s22.value_counts().to_dict()); p("2024 row sums", s24.value_counts().to_dict())
all9_22 = (d.loc[d.YEAR == 2022, P5] == 9).all(axis=1); all9_24 = (d.loc[d.YEAR == 2024, P8] == 9).all(axis=1)
p("all-9 rows 2022", int(all9_22.sum()), "2024", int(all9_24.sum()))
valid22 = (d.YEAR == 2022) & (d[P5].sum(axis=1) == 30)
valid24 = (d.YEAR == 2024) & (d[P8].sum(axis=1) == 36)
d["Q20VALID"] = valid22 | valid24
p("valid complete rankings 2022", int(valid22.sum()), "2024", int(valid24.sum()))
# convert to rank positions (1 = most concerning)
R = pd.DataFrame(index=d.index)
for c in P5: R.loc[valid22, c] = 9 - d.loc[valid22, c]
for c in P8: R.loc[valid24, c] = 9 - d.loc[valid24, c]
p("2022 rank ranges", {c: (R.loc[valid22, c].min(), R.loc[valid22, c].max()) for c in P5})
p("2024 rank ranges", {c: (R.loc[valid24, c].min(), R.loc[valid24, c].max()) for c in P8})
p("2022 mean rank (1-5)", R.loc[valid22, P5].mean().round(2).to_dict())
p("2022 % first", (R.loc[valid22, P5] == 1).mean().mul(100).round(1).to_dict())
p("2024 mean rank (1-8)", R.loc[valid24, P8].mean().round(2).sort_values().to_dict())
p("2024 % first", (R.loc[valid24, P8] == 1).mean().mul(100).round(1).to_dict())
p("2024 % top3", (R.loc[valid24, P8] <= 3).mean().mul(100).round(1).to_dict())
# relative order of 5 common items in 2024 re-ranked among themselves
RR = R.loc[valid24, P5].rank(axis=1)
p("2024 common-5 re-ranked mean (1-5)", RR.mean().round(2).to_dict())
p("2024 common-5 re-ranked % first", (RR == 1).mean().mul(100).round(1).to_dict())
for c in P5:
    a = R.loc[valid22, c]; b = RR[c]
    mwu = stats.mannwhitneyu(a, b)
    p(f"  {c}: 2022 {a.mean():.2f} vs 2024(rel) {b.mean():.2f}  MW p={mwu.pvalue:.4f}; %first {((a==1).mean()*100):.1f} vs {((b==1).mean()*100):.1f}")
# share of 2024 respondents placing a new item (pollution/health/value) first
p("2024 share with a NEW item ranked first", round(((R.loc[valid24, ['PREOC_CONTAMINACION','PREOC_SALUD','PREOC_VALOR_VIV']] == 1).any(axis=1)).mean()*100, 1))

p("\n== CORRELATIONS (Spearman)")
def sp(a, b):
    x = d[[a, b]].dropna(); r = stats.spearmanr(x[a], x[b]); p(f"{a} ~ {b}: rho={r.correlation:.3f} p={r.pvalue:.4g} n={len(x)}")
sp("ACT_SOST", "CONOC_SOST")
sp("PROB", "ACT_SOST"); sp("PROB", "CONOC_SOST"); sp("PROB", "PRED_BIOCOM_AMB"); sp("PROB", "PRED_BIOCOM_ECON"); sp("PROB", "PRED_BIOCOM_EFIC"); sp("PROB", "PRED_BIOCOM_SEG")
sp("PROB", "RAZF_RANK_AYUDAS_INV")
sp("ACUERDO_SEPARAR", "ACUERDO_COMPOSTAJE"); sp("ACUERDO_SEPARAR", "TRAT_RANK_COMPOSTAJE_INV"); sp("ACUERDO_SEPARO_SIEMPRE", "ACUERDO_COMPOSTAJE"); sp("ACUERDO_SEPARO_SIEMPRE", "TRAT_RANK_COMPOSTAJE_INV")
sp("ACT_SOST", "ACUERDO_COMPOSTAJE"); sp("ACT_SOST", "ACUERDO_PIROLISIS")
Q17 = ["ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS", "ACUERDO_BIOREF_NO_CERCA", "ACUERDO_BIOREF_SOSTENIBLES", "ACUERDO_BIOREF_CONFIANZA", "ACUERDO_BIOREF_PAGARIA_MAS", "DISTANCIA_MIN_BIOREF", "ACT_SOST", "CONOC_SOST", "PROB"]
p(d[Q17].corr(method="spearman").round(3))
p("Q17 n complete", len(d[Q17].dropna()))

p("\n== H3 model: PROB ~ motivations ranks + covariates (OLS HC3)")
df = d[["PROB", "PRED_BIOCOM_AMB", "PRED_BIOCOM_ECON", "RAZF_RANK_AYUDAS_INV", "ACT_SOST", "CONOC_SOST", "W24", "GENERO_B", "EDAD_CAT", "EDU", "INC"]].dropna()
m = smf.ols(f"PROB ~ ACT_SOST + CONOC_SOST + PRED_BIOCOM_AMB + RAZF_RANK_AYUDAS_INV + W24 + {cov}", data=df).fit(cov_type="HC3")
p(m.summary().tables[1].as_text().split("\n")[0:3]);
for k in ["ACT_SOST", "CONOC_SOST", "PRED_BIOCOM_AMB", "RAZF_RANK_AYUDAS_INV", "W24"]:
    ci = m.conf_int().loc[k]; p(f"  {k}: b={m.params[k]:.3f} [{ci[0]:.3f},{ci[1]:.3f}] p={m.pvalues[k]:.4f}")
p("  n=", int(m.nobs), "R2adj=", round(m.rsquared_adj, 3))

p("\n== H4 model: ACCEPT (willingness to use, sustainable, trust) ~ benefits-outweigh + NO_CERCA + covariates")
df = d[["ACEPT_BIOREF", "ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS", "ACUERDO_BIOREF_NO_CERCA", "DISTANCIA_MIN_BIOREF", "W24", "GENERO_B", "EDAD_CAT", "EDU", "INC"]].dropna()
m = smf.ols(f"DISTANCIA_MIN_BIOREF ~ ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS + ACUERDO_BIOREF_SOSTENIBLES + ACUERDO_BIOREF_CONFIANZA + W24 + {cov}", data=d[["DISTANCIA_MIN_BIOREF","ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS","ACUERDO_BIOREF_SOSTENIBLES","ACUERDO_BIOREF_CONFIANZA","W24","GENERO_B","EDAD_CAT","EDU","INC"]].dropna()).fit(cov_type="HC3")
for k in ["ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS", "ACUERDO_BIOREF_SOSTENIBLES", "ACUERDO_BIOREF_CONFIANZA", "W24"]:
    ci = m.conf_int().loc[k]; p(f"  dist ~ {k}: b={m.params[k]:.3f} [{ci[0]:.3f},{ci[1]:.3f}] p={m.pvalues[k]:.4f}")
p("  n=", int(m.nobs), "R2adj=", round(m.rsquared_adj, 3))

p("\n== CLUSTERS")
p("QCL_1 by wave", pd.crosstab(d.YEAR, d.QCL_1).to_dict("index"), "n", int(d.QCL_1.notna().sum()))
ct = pd.crosstab(d.YEAR, d.QCL_1); chi = stats.chi2_contingency(ct); p(f"chi2={chi[0]:.3f} df={chi[2]} p={chi[1]:.4f}")
p("cluster % ", d.QCL_1.value_counts(normalize=True).mul(100).round(1).to_dict())
cands = {
    "A: BIOCOM_ECON/AMB/SEG + ACT + CONOC": ["ZBIOCOM_ECON", "ZBIOCOM_AMB", "ZBIOCOM_SEG", "ZACT_SOST", "ZCONOC_SOST"],
    "B: A + ZACEPT_BIOREF": ["ZBIOCOM_ECON", "ZBIOCOM_AMB", "ZBIOCOM_SEG", "ZACT_SOST", "ZCONOC_SOST", "ZACEPT_BIOREF"],
    "C: B + PREOC impact/risk": ["ZBIOCOM_ECON", "ZBIOCOM_AMB", "ZBIOCOM_SEG", "ZACT_SOST", "ZCONOC_SOST", "ZACEPT_BIOREF", "ZPREOC_BIOREF_IMPACTO", "ZPREOC_BIOREF_RIESGO"],
    "D: C + PROB": ["ZBIOCOM_ECON", "ZBIOCOM_AMB", "ZBIOCOM_SEG", "ZACT_SOST", "ZCONOC_SOST", "ZACEPT_BIOREF", "ZPREOC_BIOREF_IMPACTO", "ZPREOC_BIOREF_RIESGO", "ZPROB_BIOCOM"],
}
best = None
for name, cols in cands.items():
    sub = d[cols + ["QCL_1"]].dropna()
    km = KMeans(n_clusters=3, n_init=50, random_state=1).fit(sub[cols])
    ari = adjusted_rand_score(sub.QCL_1, km.labels_)
    p(f"{name}: n={len(sub)} ARI vs QCL_1={ari:.3f}")
    if best is None or ari > best[1]: best = (name, ari, cols)
p("best", best[0], round(best[1], 3))
cols = best[2]; sub = d[cols + ["QCL_1", "YEAR"]].dropna()
p("cluster centroids (z) of reported solution:"); p(sub.groupby("QCL_1")[cols].mean().round(2))
for k in range(2, 7):
    km = KMeans(n_clusters=k, n_init=50, random_state=1).fit(sub[cols])
    p(f"k={k} silhouette={silhouette_score(sub[cols], km.labels_):.3f}")
aris = []
rng = np.random.default_rng(7)
for i in range(50):
    idx = rng.choice(len(sub), size=int(0.8*len(sub)), replace=False)
    km = KMeans(n_clusters=3, n_init=10, random_state=i).fit(sub[cols].iloc[idx])
    aris.append(adjusted_rand_score(sub.QCL_1.iloc[idx], km.labels_))
p(f"bootstrap 80% subsample ARI vs reported: mean={np.mean(aris):.3f} min={np.min(aris):.3f}")
# per-wave clustering consistency
for yr in [2022, 2024]:
    s = sub[sub.YEAR == yr]; km = KMeans(n_clusters=3, n_init=50, random_state=1).fit(s[cols])
    p(f"wave {yr} separate k-means ARI vs reported={adjusted_rand_score(s.QCL_1, km.labels_):.3f}")
# profiles
prof = d.groupby("QCL_1")[["ACEPT_BIOREF", "ESPECTRO_IZQ_DER", "PROB", "DISTANCIA_MIN_BIOREF", "ACT_SOST", "CONOC_SOST"]].mean().round(2)
p(prof)
p(pd.crosstab(d.QCL_1, d.EDAD_CAT, normalize="index").round(3))
for c in ["EDAD_CAT", "GENERO_B", "EDU", "INC"]:
    ct = pd.crosstab(d.QCL_1, d[c]); chi = stats.chi2_contingency(ct)
    n = ct.values.sum(); v = np.sqrt(chi[0]/(n*(min(ct.shape)-1)))
    p(f"cluster x {c}: chi2={chi[0]:.2f} df={chi[2]} p={chi[1]:.4f} V={v:.3f}")
for c in ["ACEPT_BIOREF", "ESPECTRO_IZQ_DER", "DISTANCIA_MIN_BIOREF"]:
    g = [x.dropna() for _, x in d.groupby("QCL_1")[c]]; f = stats.f_oneway(*g); p(f"ANOVA {c} by cluster F={f.statistic:.2f} p={f.pvalue:.4g}")

p("\n== Q15 ipsative check: correlations among motivation ranks")
p(d[["PRED_BIOCOM_AMB", "PRED_BIOCOM_ECON", "PRED_BIOCOM_EFIC", "PRED_BIOCOM_SEG"]].corr(method="spearman").round(2))
p("Q15 mean importance (inv, 1-6)", d[["RAZF_RANK_MEDIOAMBIENTE_INV", "RAZF_RANK_SI_ECONOMICOS_INV", "RAZF_RANK_SI_EFICIENTES_INV", "RAZF_RANK_SI_SEGUROS_INV", "RAZF_RANK_AYUDAS_INV", "RAZF_RANK_EMPLEO_LOCAL_INV"]].mean().round(2).to_dict())
p("Q16 mean importance (inv, 1-7)", d[[c for c in d.columns if c.startswith("RAZC") and c.endswith("_INV")]].mean().round(2).to_dict())
p("Q16 by wave ayudas", d.groupby("YEAR").RAZC_RANK_NO_SIN_AYUDAS_INV.mean().round(2).to_dict())
p("Q15 % first", {c: round((d[c] == 1).mean()*100, 1) for c in ["RAZF_RANK_MEDIOAMBIENTE", "RAZF_RANK_SI_ECONOMICOS", "RAZF_RANK_SI_EFICIENTES", "RAZF_RANK_SI_SEGUROS", "RAZF_RANK_AYUDAS", "RAZF_RANK_EMPLEO_LOCAL"]})
p("Q16 % first", {c: round((d[c] == 1).mean()*100, 1) for c in [c for c in d.columns if c.startswith("RAZC") and not c.endswith("_INV")]})

p("\n== Q9 items by wave share agree (4-5)")
for c in [c for c in d.columns if c.startswith("SOS_ACUERDO")]:
    p(c, d.groupby("YEAR")[c].apply(lambda s: round((s >= 4).sum()/s.notna().sum()*100, 1)).to_dict(), "M", d.groupby("YEAR")[c].mean().round(2).to_dict())
p("Q17 share agree pooled", {c: round((d[c] >= 4).sum()/d[c].notna().sum()*100, 1) for c in Q17[:5]})
p("Q17 means pooled", d[Q17[:5]].agg(["mean", "std"]).round(2).to_dict())