import sys, warnings; warnings.filterwarnings("ignore")
import pandas as pd, numpy as np
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy import stats
src, out = sys.argv[1], sys.argv[2]
d = pd.read_csv(src, 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)
C22, C24 = "#7FA7C9", "#1F4E79"
plt.rcParams.update({"font.family": "Arial", "font.size": 9, "axes.spines.top": False, "axes.spines.right": False})
def save(f, name): f.tight_layout(); f.savefig(f"{out}/{name}.png", dpi=600); plt.close(f)
def agree(c): return d.groupby("YEAR")[c].apply(lambda s: (s >= 4).sum()/s.notna().sum()*100)

def paired_bar(ax, labels, v22, v24, na24=None, xlabel="% of respondents"):
    yy = np.arange(len(labels))[::-1]; h = 0.38
    ax.barh(yy + h/2, v22, h, color=C22, label="2022")
    ax.barh(yy - h/2, v24, h, color=C24, label="2024")
    for i, (a, b) in enumerate(zip(v22, v24)):
        ax.text(a + 0.8 if a == a else 0.8, yy[i] + h/2, "n/a" if a != a else f"{a:.1f}", va="center", fontsize=7)
        ax.text(b + 0.8, yy[i] - h/2, f"{b:.1f}", va="center", fontsize=7)
    ax.set_yticks(yy); ax.set_yticklabels(labels); ax.set_xlabel(xlabel); ax.legend(frameon=False, loc="lower right")

# Figure 1: Q9 items share agree
items = [("SOS_ACUERDO_CC", "Urgent action on climate change"), ("SOS_ACUERDO_CIRC", "Promote the circular economy"),
         ("SOS_ACUERDO_PRODUCTOS", "New solutions for organic waste"), ("SOS_ACUERDO_COMBUSTIBLES", "Authorities should promote sustainable fuels"),
         ("SOS_ACUERDO_SECTORES", "Aware of sectoral sustainability measures"), ("SOS_ACUERDO_GEST_RESIDUOS", "Informed about managing own waste"),
         ("SOS_ACUERDO_BIOGAS", "Informed about waste-derived fuels (e.g., biogas)"), ("SOS_ACUERDO_PLANTAS", "Informed about waste-to-fuel plants")]
f, ax = plt.subplots(figsize=(6.5, 3.8))
paired_bar(ax, [l for _, l in items], [agree(c)[2022] for c, _ in items], [agree(c)[2024] for c, _ in items], xlabel="% agreeing (4–5 on a 1–5 scale)")
ax.axhline(3.5, color="grey", lw=0.6, ls=":"); ax.set_xlim(0, 100)
save(f, "Figure1")

# Figure 2: Q14 distribution
P = d.PROB_USO_BIOCOM.where(d.PROB_USO_BIOCOM >= 1)
ct = pd.crosstab(d.YEAR, P, normalize="index") * 100
f, ax = plt.subplots(figsize=(6.5, 2.2))
cols = ["#C0504D", "#E6A49C", "#D9D9D9", "#9DC3E6", "#1F4E79"]; labs = ["1 Very unlikely", "2", "3", "4", "5 Totally likely"]
for r, yr in enumerate([2024, 2022]):
    left = 0
    for k in range(1, 6):
        v = ct.loc[yr, k]; ax.barh(r, v, left=left, color=cols[k-1], label=labs[k-1] if r == 0 else None)
        if v > 4: ax.text(left + v/2, r, f"{v:.0f}", ha="center", va="center", fontsize=7, color="white" if k in (1, 5) else "black")
        left += v
ax.set_yticks([0, 1]); ax.set_yticklabels(["2024 (n = 504)", "2022 (n = 498)"]); ax.set_xlim(0, 100); ax.set_xlabel("% of respondents")
ax.legend(ncol=5, frameon=False, bbox_to_anchor=(0.5, 1.25), loc="upper center", fontsize=7)
save(f, "Figure2")

# Figure 3: Q17
items = [("ACUERDO_BIOREF_CONFIANZA", "Would use biorefinery products because they generate trust"), ("ACUERDO_BIOREF_SOSTENIBLES", "Biorefineries contribute to a more sustainable model"),
         ("ACUERDO_BIOREF_NO_CERCA", "Would like their advantages, but not near my home"), ("ACUERDO_BIOREF_PAGARIA_MAS", "Would use their products even if somewhat more expensive"),
         ("ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS", "Sustainability benefits outweigh nuisances near my home")]
f, ax = plt.subplots(figsize=(6.5, 2.9))
paired_bar(ax, [l for _, l in items], [agree(c)[2022] for c, _ in items], [agree(c)[2024] for c, _ in items], xlabel="% agreeing (4–5 on a 1–5 scale)"); ax.set_xlim(0, 100)
save(f, "Figure3")

# Figure 4: distance histogram + by benefits agreement
f, (a1, a2) = plt.subplots(1, 2, figsize=(6.5, 2.6), gridspec_kw={"width_ratios": [1.4, 1]})
bins = np.arange(0.5, 31.5, 1)
a1.hist(d.loc[d.YEAR == 2022, "DISTANCIA_MIN_BIOREF"].dropna(), bins=bins, color=C22, alpha=0.85, label="2022")
a1.hist(d.loc[d.YEAR == 2024, "DISTANCIA_MIN_BIOREF"].dropna(), bins=bins, color=C24, histtype="step", linewidth=1.4, label="2024")
a1.set_xlabel("Minimum acceptable distance (km; response range 1–30)"); a1.set_ylabel("Respondents"); a1.legend(frameon=False); a1.set_title("(a)", loc="left", fontsize=9)
g = d.groupby("ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS").DISTANCIA_MIN_BIOREF
m = g.mean(); se = g.sem()
a2.errorbar(m.index, m.values, yerr=1.96*se.values, fmt="o-", color=C24, capsize=3)
a2.set_xticks([1, 2, 3, 4, 5]); a2.set_xlabel("'Benefits outweigh nuisances' (1–5)"); a2.set_ylabel("Mean distance (km, 95% CI)"); a2.set_title("(b)", loc="left", fontsize=9)
save(f, "Figure4")

# Figure 5: Q19
items = [("REQ_OLORDESIN", "Odour and pest control"), ("REQ_EMPLEO_LOCAL", "Local employment"), ("REQ_SEGURIDAD", "All measures to make the plant safe"),
         ("REQ_DEMANDA_LOCAL", "Bioenergy supplied to local demand"), ("REQ_NO_AFECTAR_PRECIO", "No negative effect on housing/business value*"),
         ("REQ_IMPACTO_VISUAL", "Sites minimising visual impact"), ("REQ_CIRCULACION_CAM", "Mitigation of truck traffic"),
         ("REQ_COLAB_INDUSTRIA", "Collaboration with regional industry"), ("REQ_NINGUN_CASO", "Not acceptable under any circumstances")]
v22 = [np.nan if c == "REQ_NO_AFECTAR_PRECIO" else d[d.YEAR == 2022][c].mean()*100 for c, _ in items]
v24 = [d[d.YEAR == 2024][c].mean()*100 for c, _ in items]
f, ax = plt.subplots(figsize=(6.5, 3.8))
paired_bar(ax, [l for _, l in items], v22, v24, xlabel="% selecting the requirement as essential (multiple response)"); ax.set_xlim(0, 80)
save(f, "Figure5")

# Figure 6: Q20 corrected
P5 = ["PREOC_OLORES", "PREOC_IMPACTO_VISUAL", "PREOC_CAMIONES", "PREOC_INSEGURIDAD", "PREOC_OBRAS"]
P8 = P5 + ["PREOC_CONTAMINACION", "PREOC_SALUD", "PREOC_VALOR_VIV"]
L = {"PREOC_OLORES": "Odours, pests, dirt", "PREOC_IMPACTO_VISUAL": "Visual impact", "PREOC_CAMIONES": "Truck traffic", "PREOC_INSEGURIDAD": "Plant may be unsafe",
     "PREOC_OBRAS": "Construction nuisance", "PREOC_CONTAMINACION": "Air or water pollution", "PREOC_SALUD": "Health risks", "PREOC_VALOR_VIV": "Loss of property value"}
v22m = (d.YEAR == 2022) & (d[P5].sum(axis=1) == 30); v24m = (d.YEAR == 2024) & (d[P8].sum(axis=1) == 36)
R22 = 9 - d.loc[v22m, P5]; R24 = 9 - d.loc[v24m, P8]; RR = R24[P5].rank(axis=1)
f, (a1, a2) = plt.subplots(1, 2, figsize=(6.8, 2.9))
order = ["PREOC_INSEGURIDAD", "PREOC_OLORES", "PREOC_OBRAS", "PREOC_IMPACTO_VISUAL", "PREOC_CAMIONES"]
yy = np.arange(5)[::-1]; h = 0.38
a1.barh(yy + h/2, [(R22[c] == 1).mean()*100 for c in order], h, color=C22, label=f"2022 (n = {v22m.sum()})")
a1.barh(yy - h/2, [(RR[c] == 1).mean()*100 for c in order], h, color=C24, label=f"2024 (n = {v24m.sum()})")
a1.set_yticks(yy); a1.set_yticklabels([L[c] for c in order]); a1.set_xlabel("% ranking first among the five common items"); a1.legend(frameon=False, fontsize=7); a1.set_title("(a) Five items common to both waves", loc="left", fontsize=8)
o8 = R24.mean().sort_values().index.tolist(); yy = np.arange(8)[::-1]
a2.barh(yy, [(R24[c] <= 3).mean()*100 for c in o8], color=[C24 if c in P5 else "#C55A11" for c in o8])
for i, c in enumerate(o8): a2.text((R24[c] <= 3).mean()*100 + 1, yy[i], f"{(R24[c] <= 3).mean()*100:.0f}", va="center", fontsize=7)
a2.set_yticks(yy); a2.set_yticklabels([L[c] for c in o8]); a2.set_xlabel("% ranking the item in the top three"); a2.set_title("(b) 2024, eight items", loc="left", fontsize=8); a2.set_xlim(0, 65)
save(f, "Figure6")

# Figure 7: clusters
Z = ["ZBIOCOM_ECON", "ZBIOCOM_AMB", "ZBIOCOM_SEG", "ZACT_SOST", "ZCONOC_SOST"]
ZL = ["Economic motive\n(biofuels)", "Environmental motive\n(biofuels)", "Safety motive\n(biofuels)", "Sustainability\nattitude", "Perceived\nknowledge"]
cen = d.groupby("QCL_1")[Z].mean()
names = {1: "C1 Environment- and safety-oriented", 2: "C2 Economy-oriented", 3: "C3 Low sustainability engagement"}
f, ax = plt.subplots(figsize=(6.5, 2.8)); x = np.arange(5); w = 0.27
for i, (k, col) in enumerate(zip([1, 2, 3], ["#2E86AB", "#1F4E79", "#C55A11"])):
    n = int((d.QCL_1 == k).sum()); ax.bar(x + (i-1)*w, cen.loc[k], w, color=col, label=f"{names[k]} (n = {n}; {n/d.QCL_1.notna().sum()*100:.1f}%)")
ax.axhline(0, color="black", lw=0.6); ax.set_xticks(x); ax.set_xticklabels(ZL, fontsize=8); ax.set_ylabel("Cluster centroid (z-score)")
ax.legend(frameon=False, fontsize=7, loc="lower left")
save(f, "Figure7")

# print support index stats for text
Q = ["ACUERDO_BIOREF_SOSTENIBLES", "ACUERDO_BIOREF_CONFIANZA", "ACUERDO_BIOREF_PAGARIA_MAS", "ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS"]
d["SUPPORT"] = d[Q].mean(axis=1).where(d[Q].notna().sum(axis=1) >= 3)
print("SUPPORT by wave", d.groupby("YEAR").SUPPORT.agg(["mean", "std", "count"]).round(2).to_dict())
print("SUPPORT by cluster", d.groupby("QCL_1").SUPPORT.agg(["mean", "std", "count"]).round(2).to_dict())
g = [x.dropna() for _, x in d.groupby("QCL_1").SUPPORT]; print(stats.f_oneway(*g))
print("PROB by cluster", d.assign(P=P).groupby("QCL_1").P.agg(["mean", "std"]).round(2).to_dict())
print("ideology by cluster", d.groupby("QCL_1").ESPECTRO_IZQ_DER.agg(["mean", "std"]).round(2).to_dict())
print("dist by cluster", d.groupby("QCL_1").DISTANCIA_MIN_BIOREF.agg(["mean"]).round(2).to_dict())
x = d.assign(P=P)[["SUPPORT","P"]].dropna(); print("SUPPORT~PROB", stats.spearmanr(x.SUPPORT, x.P))
x = d[["SUPPORT","DISTANCIA_MIN_BIOREF"]].dropna(); print("SUPPORT~DIST", stats.spearmanr(x.SUPPORT, x.DISTANCIA_MIN_BIOREF), len(x))
x = d[["ACUERDO_BIOREF_NO_CERCA","DISTANCIA_MIN_BIOREF"]].dropna(); print("NOCERCA~DIST", stats.spearmanr(x.iloc[:,0], x.iloc[:,1]))
x = d[["ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS","DISTANCIA_MIN_BIOREF"]].dropna(); print("BEN~DIST", stats.spearmanr(x.iloc[:,0], x.iloc[:,1]), len(x))
x = d[["ACUERDO_BIOREF_MOLESTIAS_BENEFICIOS","ACUERDO_BIOREF_SOSTENIBLES"]].dropna(); print("BEN~SOST", stats.spearmanr(x.iloc[:,0], x.iloc[:,1]), len(x))
x = d[["ACUERDO_BIOREF_CONFIANZA","ACUERDO_BIOREF_SOSTENIBLES"]].dropna(); print("CONF~SOST", stats.spearmanr(x.iloc[:,0], x.iloc[:,1]), len(x))