# -*- coding: utf-8 -*-
"""Optimistic, pessimistic and stochastic: Robust ranking of national energy-transition performance over a common weight polytope

GOMD arastirma makalesi 10 - makaledeki tum sayilari, tablolari ve sekilleri ureten kod (Google Colab surumu).
Colab: kodu tek bir hucreye yapistirip calistirin ya da .ipynb dosyasini Dosya > Not defteri yukle ile acin.
Kaynak: GOMD acik veri platformu (https://gomd.pirireis.edu.tr) - Piri Reis Universitesi ve
T.C. Cevre, Sehircilik ve Iklim Degisikligi Bakanligi Cevre Yonetimi Genel Mudurlugu.
"""

# # Optimistic, pessimistic and stochastic: Robust ranking of national energy-transition performance over a common weight polytope
#
# Bu not defteri makaledeki **tüm sayıları, tabloları ve şekilleri** üreten kodun Colab sürümüdür (replikasyon paketindeki betiklerle aynı kod, hücrelere bölünmüş hâlde).
#
# **Çalıştırma:** Çalışma zamanı → Tümünü çalıştır (*Runtime → Run all*). Tahmini süre: 2–4 dakika.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: Our World in Data/Ember, EDGAR 2026, Dünya Bankası, Dünya Karbon Fiyatlama Veritabanı (GOMD API)). Makaledeki sayıları birebir yeniden üretmek için replikasyon paketi ZIP'ini yükleyin (`VERI_KAYNAGI = "zip"`); kaynak kurumlar verilerini revize ederse canlı indirmede küçük farklar olabilir. Verilerin çekim tarihi 8 Kasım 2026'dır.
#
# Şekiller hem hücre altında görünür hem `fig/` klasörüne, sayısal sonuçlar `out/` klasörüne yazılır. Son hücre hepsini ZIP olarak indirir.

# ## 0. Kurulum

# Installs the packages used in the paper. If Colab asks to restart the session, restart and run all cells again.
# [Colab Magic] !pip install -q scipy openpyxl
# [Colab Magic] !pip install -q --upgrade gomd

VERI_KAYNAGI = "api"   # "api": download from the original sources; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp10", exist_ok=True); os.chdir("/content/rp10")
for d in ("data", "out", "fig"): os.makedirs(d, exist_ok=True)
if VERI_KAYNAGI == "zip":
    import zipfile
    from google.colab import files
    up = files.upload()                      # choose *_replication_package.zip
    with zipfile.ZipFile(next(iter(up))) as z:
        for n in z.namelist():
            if n.startswith("data/") or n == "entities.json": z.extract(n, ".")
    print(sorted(os.listdir("data")))
print("Calisma klasoru:", os.getcwd())

# ## 1. Veri

# --- common.py dosyasi yaziliyor (notebook'taki %%writefile hucresi) ---
open('common.py', 'w', encoding='utf-8').write('"""Shared download helpers (public APIs; no keys needed)."""\nimport io, json, time, zipfile\nfrom pathlib import Path\nimport requests, pandas as pd\nS = requests.Session(); S.headers["User-Agent"] = "research-replication-client"\nDATA = Path("data"); DATA.mkdir(exist_ok=True)\n\ndef get(url, tries=6, **kw):\n    for k in range(tries):\n        try:\n            r = S.get(url, timeout=300, **kw)\n            if r.status_code == 200: return r\n            if r.status_code in (429, 500, 502, 503, 504): time.sleep(5 * (k + 1)); continue\n            r.raise_for_status()\n        except requests.RequestException:\n            if k == tries - 1: raise\n            time.sleep(5 * (k + 1))\n    raise RuntimeError(url)\n\ndef power_daily(lat, lon, params, start="19810101", end="20261031"):\n    url = (f"https://power.larc.nasa.gov/api/temporal/daily/point?parameters={\',\'.join(params)}&community=RE"\n           f"&longitude={lon}&latitude={lat}&start={start}&end={end}&format=JSON")\n    d = get(url).json()["properties"]["parameter"]; time.sleep(0.3)\n    return {p: {k: v for k, v in ser.items() if v != -999} for p, ser in d.items()}\n\ndef power_monthly(lat, lon, params, start=1981, end=2026):\n    url = (f"https://power.larc.nasa.gov/api/temporal/monthly/point?parameters={\',\'.join(params)}&community=RE"\n           f"&longitude={lon}&latitude={lat}&start={start}&end={end}&format=JSON")\n    d = get(url).json()["properties"]["parameter"]; time.sleep(0.3)\n    return {p: {k: v for k, v in ser.items() if v != -999 and k[4:] != "13"} for p, ser in d.items()}\n\ndef jsonstat_frame(js):\n    dims = js["id"]; size = js["size"]\n    cats = [sorted(js["dimension"][d]["category"]["index"].items(), key=lambda kv: kv[1]) for d in dims]\n    rows = []\n    for k, v in js["value"].items():\n        k = int(k); idx = []\n        for s in reversed(size):\n            idx.append(k % s); k //= s\n        idx = list(reversed(idx)); rows.append([cats[i][idx[i]][0] for i in range(len(dims))] + [v])\n    return pd.DataFrame(rows, columns=dims + ["value"])\n\ndef eurostat(dataset, query):\n    url = f"https://ec.europa.eu/eurostat/api/dissemination/statistics/1.0/data/{dataset}?format=JSON&lang=EN&{query}"\n    return jsonstat_frame(get(url).json())\n\ndef edgar(fn):\n    """EDGAR 2026 GHG release (JRC/IEA): downloads the zip and extracts the xlsx into data/."""\n    base = "https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/EDGAR/datasets/EDGAR_2026_GHG/"\n    with zipfile.ZipFile(io.BytesIO(get(base + fn).content)) as z:\n        for n in z.namelist():\n            if n.endswith(".xlsx"):\n                (DATA / Path(n).name).write_bytes(z.read(n)); print("EDGAR", Path(n).name)\n\ndef worldbank(ind, years="1985:2025"):\n    d = get(f"https://api.worldbank.org/v2/country/all/indicator/{ind}?format=json&per_page=20000&date={years}").json()\n    rows = [(r["countryiso3code"], r["country"]["value"], int(r["date"]), r["value"]) for r in d[1] if r["value"] is not None and r["countryiso3code"]]\n    df = pd.DataFrame(rows, columns=["iso3", "name", "year", "value"]); df.to_csv(DATA / f"wb_{ind}.csv", index=False); print("WB", ind, len(df))\n\ndef gomd_panel(indicator, entities, out_csv):\n    """Series retrieved through the public GOMD open data API (https://gomd.pirireis.edu.tr) with the gomd Python client."""\n    from gomd import GOMDClient\n    client = GOMDClient("https://gomd.pirireis.edu.tr"); recs = []\n    for code, name, etype in entities:\n        for k in range(4):\n            try:\n                rows = client.get_series(indicator=indicator, entity=code, limit=50000, as_df=False); break\n            except Exception:\n                if k == 3: raise\n                time.sleep(2 ** k)\n        recs += [(code, name, etype, str(r["period"]), float(r["value"])) for r in rows if str(r.get("entity_code") or "") == code and r.get("value") is not None]\n    df = pd.DataFrame(recs, columns=["entity", "name", "etype", "period", "value"]).drop_duplicates(["entity", "period"])\n    df.to_csv(out_csv, index=False); print("GOMD", indicator, len(df)); return df\n')

import json
if not os.path.exists('entities.json'):
    json.dump([["TUR", "Turkiye", "country"], ["AFG", "Afghanistan", "country"], ["ALB", "Albania", "country"], ["DZA", "Algeria", "country"], ["DEU", "Germany", "country"], ["AND", "Andorra", "country"], ["AGO", "Angola", "country"], ["ATG", "Antigua and Barbuda", "country"], ["ARG", "Argentina", "country"], ["ARM", "Armenia", "country"], ["AUS", "Australia", "country"], ["AUT", "Austria", "country"], ["AZE", "Azerbaijan", "country"], ["BHS", "Bahamas, The", "country"], ["BHR", "Bahrain", "country"], ["BGD", "Bangladesh", "country"], ["BRB", "Barbados", "country"], ["BLR", "Belarus", "country"], ["BEL", "Belgium", "country"], ["BLZ", "Belize", "country"], ["BEN", "Benin", "country"], ["BTN", "Bhutan", "country"], ["BOL", "Bolivia", "country"], ["BIH", "Bosnia and Herzegovina", "country"], ["BWA", "Botswana", "country"], ["BRA", "Brazil", "country"], ["BRN", "Brunei Darussalam", "country"], ["BGR", "Bulgaria", "country"], ["BFA", "Burkina Faso", "country"], ["BDI", "Burundi", "country"], ["CPV", "Cabo Verde", "country"], ["KHM", "Cambodia", "country"], ["CMR", "Cameroon", "country"], ["CAN", "Canada", "country"], ["CAF", "Central African Republic", "country"], ["TCD", "Chad", "country"], ["CHL", "Chile", "country"], ["TWN", "Chinese Taipei", "country"], ["CHN", "China", "country"], ["COL", "Colombia", "country"], ["COM", "Comoros", "country"], ["COG", "Congo, Rep.", "country"], ["CRI", "Costa Rica", "country"], ["CIV", "Cote d'Ivoire", "country"], ["HRV", "Croatia", "country"], ["CUB", "Cuba", "country"], ["CYP", "Greek Cypriot Administration of Southern Cyprus (GCA)", "country"], ["CZE", "Czechia", "country"], ["COD", "Democratic Republic of Congo", "country"], ["DNK", "Denmark", "country"], ["DJI", "Djibouti", "country"], ["DMA", "Dominica", "country"], ["DOM", "Dominican Republic", "country"], ["ECU", "Ecuador", "country"], ["EGY", "Egypt", "country"], ["SLV", "El Salvador", "country"], ["GNQ", "Equatorial Guinea", "country"], ["ERI", "Eritrea", "country"], ["EST", "Estonia", "country"], ["SWZ", "Eswatini", "country"], ["ETH", "Ethiopia", "country"], ["FJI", "Fiji", "country"], ["FIN", "Finland", "country"], ["FRA", "France", "country"], ["GAB", "Gabon", "country"], ["GMB", "Gambia, The", "country"], ["GEO", "Georgia", "country"], ["GHA", "Ghana", "country"], ["GRD", "Grenada", "country"], ["GTM", "Guatemala", "country"], ["GIN", "Guinea", "country"], ["GNB", "Guinea-Bissau", "country"], ["GUY", "Guyana", "country"], ["HTI", "Haiti", "country"], ["HND", "Honduras", "country"], ["HKG", "Hong Kong SAR, China", "country"], ["HUN", "Hungary", "country"], ["ISL", "Iceland", "country"], ["IND", "India", "country"], ["IDN", "Indonesia", "country"], ["IRN", "Iran, Islamic Rep.", "country"], ["IRQ", "Iraq", "country"], ["IRL", "Ireland", "country"], ["ESP", "Spain", "country"], ["ISR", "Israel", "country"], ["ITA", "Italy", "country"], ["JAM", "Jamaica", "country"], ["JPN", "Japan", "country"], ["JOR", "Jordan", "country"], ["KAZ", "Kazakhstan", "country"], ["KEN", "Kenya", "country"], ["KIR", "Kiribati", "country"], ["KOR", "Korea Rep", "country"], ["PRK", "Korea, Dem. People's Rep.", "country"], ["KWT", "Kuwait", "country"], ["KGZ", "Kyrgyzstan", "country"], ["LAO", "Laos", "country"], ["LVA", "Latvia", "country"], ["LBN", "Lebanon", "country"], ["LSO", "Lesotho", "country"], ["LBR", "Liberia", "country"], ["LBY", "Libya", "country"], ["LIE", "Liechtenstein", "country"], ["LTU", "Lithuania", "country"], ["LUX", "Luxembourg", "country"], ["MAC", "Macao SAR, China", "country"], ["MDG", "Madagascar", "country"], ["MWI", "Malawi", "country"], ["MYS", "Malaysia", "country"], ["MDV", "Maldives", "country"], ["MLI", "Mali", "country"], ["MLT", "Malta", "country"], ["MHL", "Marshall Islands", "country"], ["MRT", "Mauritania", "country"], ["MUS", "Mauritius", "country"], ["MEX", "Mexico", "country"], ["FSM", "Micronesia, Fed. Sts.", "country"], ["MDA", "Moldova", "country"], ["MCO", "Monaco", "country"], ["MNG", "Mongolia", "country"], ["MNE", "Montenegro", "country"], ["MAR", "Morocco", "country"], ["MOZ", "Mozambique", "country"], ["MMR", "Myanmar", "country"], ["NAM", "Namibia", "country"], ["NRU", "Nauru", "country"], ["NPL", "Nepal", "country"], ["NLD", "Netherlands", "country"], ["NZL", "New Zealand", "country"], ["NIC", "Nicaragua", "country"], ["NER", "Niger", "country"], ["NGA", "Nigeria", "country"], ["MKD", "North Macedonia", "country"], ["NOR", "Norway", "country"], ["OMN", "Oman", "country"], ["PAK", "Pakistan", "country"], ["PLW", "Palau", "country"], ["PAN", "Panama", "country"], ["PNG", "Papua New Guinea", "country"], ["PRY", "Paraguay", "country"], ["PER", "Peru", "country"], ["PHL", "Philippines", "country"], ["POL", "Poland", "country"], ["PRT", "Portugal", "country"], ["PRI", "Puerto Rico", "country"], ["QAT", "Qatar", "country"], ["ROU", "Romania", "country"], ["RUS", "Russian Federation", "country"], ["RWA", "Rwanda", "country"], ["WSM", "Samoa", "country"], ["SMR", "San Marino", "country"], ["STP", "Sao Tome and Principe", "country"], ["SAU", "Saudi Arabia", "country"], ["SEN", "Senegal", "country"], ["SRB", "Serbia", "country"], ["SYC", "Seychelles", "country"], ["SLE", "Sierra Leone", "country"], ["SGP", "Singapore", "country"], ["SVK", "Slovakia", "country"], ["SVN", "Slovenia", "country"], ["SLB", "Solomon Islands", "country"], ["SOM", "Somalia", "country"], ["ZAF", "South Africa", "country"], ["SSD", "South Sudan", "country"], ["LKA", "Sri Lanka", "country"], ["KNA", "St. Kitts and Nevis", "country"], ["LCA", "St. Lucia", "country"], ["VCT", "St. Vincent and the Grenadines", "country"], ["SDN", "Sudan", "country"], ["SUR", "Suriname", "country"], ["SWE", "Sweden", "country"], ["CHE", "Switzerland", "country"], ["SYR", "Syrian Arab Republic", "country"], ["TJK", "Tajikistan", "country"], ["TZA", "Tanzania", "country"], ["THA", "Thailand", "country"], ["TLS", "Timor-Leste", "country"], ["TGO", "Togo", "country"], ["TON", "Tonga", "country"], ["TTO", "Trinidad and Tobago", "country"], ["TUN", "Tunisia", "country"], ["TKM", "Turkmenistan", "country"], ["TUV", "Tuvalu", "country"], ["UGA", "Uganda", "country"], ["GBR", "UK", "country"], ["UKR", "Ukraine", "country"], ["ARE", "United Arab Emirates", "country"], ["URY", "Uruguay", "country"], ["USA", "USA", "country"], ["UZB", "Uzbekistan", "country"], ["VUT", "Vanuatu", "country"], ["VAT", "Vatican", "country"], ["VEN", "Venezuela, RB", "country"], ["VNM", "Vietnam", "country"], ["PSE", "West Bank and Gaza", "country"], ["ESH", "Western Sahara", "country"], ["YEM", "Yemen", "country"], ["GRC", "Greece", "country"], ["ZMB", "Zambia", "country"], ["ZWE", "Zimbabwe", "country"]], open('entities.json', 'w', encoding='utf-8'), ensure_ascii=False)

# Downloads the inputs from the original sources (only when VERI_KAYNAGI == 'api').
if VERI_KAYNAGI == 'api':
    """Re-downloads the inputs of the energy-transition ranking paper into data/: Our World in Data energy dataset (Ember electricity data),
    EDGAR 2026 CO2 release, World Bank WDI and the carbon tax rate of the World Carbon Pricing Database (via the GOMD API)."""
    import json, io
    import pandas as pd
    from common import get, edgar, worldbank, gomd_panel, DATA
    e = pd.read_csv(io.StringIO(get("https://raw.githubusercontent.com/owid/energy-data/master/owid-energy-data.csv").text))
    e[["iso_code", "country", "year", "wind_share_elec", "solar_share_elec", "carbon_intensity_elec", "low_carbon_share_elec"]].to_csv(DATA / "owid_energy.csv", index=False)
    edgar("IEA_EDGAR_CO2_1970_2025.zip")
    for ind in ("SP.POP.TOTL", "NY.GDP.PCAP.PP.KD"): worldbank(ind)
    gomd_panel("carbon_tax_price_weighted", json.load(open("entities.json", encoding="utf-8")), "data/ctax.csv")
print(sorted(os.listdir('data')))

# ## 2. Kriterler, SMAA (hit-and-run), iyimser ve kötümser fayda-şüphe (BoD) modelleri, senaryolar ve şekiller

# Paper 10: optimistic, pessimistic and stochastic evaluation of national energy-transition performance over a common weight polytope. Benefit-of-the-doubt composite indicators (optimistic and pessimistic, with assurance regions) and SMAA-2 rank acceptability analysis with hit-and-run weight sampling, for the countries with 2025 electricity data. Data: Ember electricity data via Our World in Data (2015, 2025), EDGAR 2026 release (CO2, 2015 and 2025), World Bank WDI (population, GDP PPP), World Carbon Pricing Database (emissions-weighted carbon tax rate, 2025).
import json, sys, time, warnings
import numpy as np, pandas as pd
from scipy.optimize import linprog
from scipy.stats import spearmanr, kendalltau
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
plt.rcParams.update({"font.family": "DejaVu Serif", "font.size": 9, "axes.titlesize": 9.5, "legend.fontsize": 7.5,
                     "savefig.dpi": 300, "axes.spines.top": False, "axes.spines.right": False})
QUICK = "--quick" in sys.argv; NS = 2000 if QUICK else 100000
T0 = time.time(); R = {}; rng = np.random.default_rng(2025)
def log(*a): print(f"[{time.time() - T0:6.0f}s]", *a, flush=True)
LO, HI = 0.05, 0.35            # assurance region: every criterion carries between 5% and 35% of the weight

# =============================== data ===============================
e = pd.read_csv("data/owid_energy.csv", usecols=["iso_code", "country", "year", "wind_share_elec", "solar_share_elec", "carbon_intensity_elec", "low_carbon_share_elec"])
e = e[e.iso_code.notna() & ~e.iso_code.str.startswith("OWID", na=False)]
ws = e.assign(ws=e.wind_share_elec + e.solar_share_elec).pivot(index="iso_code", columns="year", values="ws")
ci = e.pivot(index="iso_code", columns="year", values="carbon_intensity_elec")
names = e.drop_duplicates("iso_code").set_index("iso_code").country
co2 = pd.read_excel("data/IEA_EDGAR_CO2_1970_2025.xlsx", sheet_name="TOTALS BY COUNTRY", header=9).groupby("Country_code_A3")[["Y_2015", "Y_2025"]].sum(min_count=1) / 1000.0
pop = pd.read_csv("data/wb_SP.POP.TOTL.csv").pivot(index="iso3", columns="year", values="value")
gdp = pd.read_csv("data/wb_NY.GDP.PCAP.PP.KD.csv").pivot(index="iso3", columns="year", values="value")
tx = pd.read_csv("data/ctax.csv", dtype={"period": str}); tx = tx[tx.etype == "country"]; tx = tx[tx.period == "2025"].set_index("entity").value
C = pd.DataFrame({
    "co2_pc": co2.Y_2025 * 1e6 / pop[2025] if 2025 in pop.columns else np.nan,
    "co2_change": np.log(co2.Y_2025 / co2.Y_2015) * 100,
    "co2_gdp": co2.Y_2025 * 1e9 / (gdp[2025] * pop[2025]),          # kg CO2 per 2021 PPP dollar
    "ci_elec": ci[2025],
    "ws_share": ws[2025],
    "ws_change": ws[2025] - ws[2015],
    "ctax": np.log1p(tx.reindex(ws.index).fillna(np.nan)),
})
C = C[C.ws_share.notna()]
C = C.dropna()
C.index.name = "iso"; C["name"] = [{"TUR": "Türkiye", "CZE": "Czechia", "KOR": "South Korea"}.get(i, names.get(i, i)) for i in C.index]
CRIT = ["co2_pc", "co2_change", "co2_gdp", "ci_elec", "ws_share", "ws_change", "ctax"]
LAB = {"co2_pc": "CO$_2$ per capita, 2025 (t)", "co2_change": "CO$_2$ change 2015–2025 (log %)", "co2_gdp": "CO$_2$ per GDP, 2025 (kg/$PPP)",
       "ci_elec": "Electricity carbon intensity, 2025 (g/kWh)", "ws_share": "Wind and solar share, 2025 (%)", "ws_change": "Wind and solar change 2015–2025 (pp)",
       "ctax": "Carbon tax rate, 2025 (log(1+US$))"}
SHORT = {"co2_pc": "CO2 pc", "co2_change": "CO2 change", "co2_gdp": "CO2/GDP", "ci_elec": "Elec. intensity", "ws_share": "W+S share", "ws_change": "W+S change", "ctax": "Carbon tax"}
DIR = {"co2_pc": -1, "co2_change": -1, "co2_gdp": -1, "ci_elec": -1, "ws_share": 1, "ws_change": 1, "ctax": 1}
C.to_csv("out/criteria.csv")
N, K = len(C), len(CRIT)
R["n_countries"] = N; R["countries"] = sorted(C.index.tolist())
R["describe"] = C[CRIT].describe().round(3).to_dict()
R["corr_spearman"] = C[CRIT].corr(method="spearman").round(3).to_dict()
log("countries", N)

def normalise(X):
    V = np.empty_like(X, dtype=float)
    for j, c in enumerate(CRIT):
        x = X[:, j]; lo, hi = x.min(), x.max(); v = (x - lo) / (hi - lo); V[:, j] = v if DIR[c] > 0 else 1 - v
    return V
X = C[CRIT].values.astype(float); V = normalise(X)
ISO = C.index.values; TUR = int(np.where(ISO == "TUR")[0][0])

# =============================== weight polytope and hit-and-run ===============================
def hit_and_run(n, lo=LO, hi=HI, k=K, thin=10, burn=2000):
    """Uniform sampling of {w : sum w = 1, lo <= w_j <= hi} by hit-and-run (Smith, 1984; Tervonen et al., 2013)."""
    w = np.full(k, 1.0 / k); out = np.empty((n, k)); it = 0; m = 0
    while m < n:
        d = rng.normal(size=k); d -= d.mean(); d /= np.linalg.norm(d)            # direction within the hyperplane sum w = 1
        with np.errstate(divide="ignore", invalid="ignore"):
            t1 = (lo - w) / d; t2 = (hi - w) / d
        tmin = np.max(np.where(d > 0, t1, np.where(d < 0, t2, -np.inf))); tmax = np.min(np.where(d > 0, t2, np.where(d < 0, t1, np.inf)))
        w = w + rng.uniform(tmin, tmax) * d; w = np.clip(w, lo, hi); w /= w.sum(); it += 1
        if it > burn and it % thin == 0: out[m] = w; m += 1
    return out
W = hit_and_run(NS)
R["weights_check"] = {"min": round(float(W.min()), 4), "max": round(float(W.max()), 4), "mean": W.mean(0).round(4).tolist()}

# =============================== SMAA-2 ===============================
def smaa(V, W):
    S = V @ W.T                                   # N x NS composite values
    order = np.argsort(-S, axis=0); ranks = np.empty_like(order); ranks[order, np.arange(S.shape[1])] = np.arange(S.shape[0])[:, None]
    B = np.zeros((V.shape[0], V.shape[0]))
    for r in range(V.shape[0]): B[:, r] = (ranks == r).mean(1)
    first = ranks == 0
    cw = np.array([W[first[i]].mean(0) if first[i].any() else np.full(W.shape[1], np.nan) for i in range(V.shape[0])])
    conf = np.array([np.mean((V @ cw[i] >= (V @ cw[i])[i] - 1e-12).sum() == 1) if np.isfinite(cw[i]).all() else np.nan for i in range(V.shape[0])])
    return S, ranks, B, cw, conf
S, RK, BA, CW, CONF = smaa(V, W)
Nn = N; meta = np.array([(Nn - r - 1) / (Nn - 1) for r in range(Nn)])        # linear meta-weights (Lahdelma & Salminen, 2001)
HAI = BA @ meta
RES = pd.DataFrame({"name": C.name, "mean_rank": RK.mean(1) + 1, "median_rank": np.median(RK, 1) + 1, "best_rank": RK.min(1) + 1, "worst_rank": RK.max(1) + 1,
                    "p_top10": (RK < 10).mean(1), "p_bottom10": (RK >= N - 10).mean(1), "a_first": BA[:, 0], "hai": HAI}, index=ISO)
RES["conf_first"] = CONF
log("SMAA done")

# =============================== benefit of the doubt: optimistic and pessimistic ===============================
def ar_rows(k=K, lo=LO, hi=HI):
    """Homogeneous assurance region on u (cone of the weight polytope): lo * sum(u) <= u_j <= hi * sum(u)."""
    A = []
    for j in range(k):
        r = np.full(k, lo); r[j] -= 1.0; A.append(r)         # lo*sum(u) - u_j <= 0
        r = np.full(k, -hi); r[j] += 1.0; A.append(r)        # u_j - hi*sum(u) <= 0
    return np.array(A)
AR = ar_rows()
def bod_opt(V, i):
    # max u.v_i  s.t. u.v_k <= 1 for all k, assurance region, u >= 0
    res = linprog(-V[i], A_ub=np.vstack([V, AR]), b_ub=np.r_[np.ones(len(V)), np.zeros(len(AR))], bounds=[(0, None)] * V.shape[1], method="highs")
    return -res.fun, res.x
def bod_pes(V, i):
    # min u.v_i  s.t. u.v_k >= 1 for all k, assurance region, u >= 0  (Zhou, Ang & Poh, 2007)
    res = linprog(V[i], A_ub=np.vstack([-V, AR]), b_ub=np.r_[-np.ones(len(V)), np.zeros(len(AR))], bounds=[(0, None)] * V.shape[1], method="highs")
    return res.fun, res.x
Vb = V + 0.01                                     # small shift keeps every country strictly positive for the pessimistic model
opt = np.array([bod_opt(Vb, i)[0] for i in range(N)]); pes = np.array([bod_pes(Vb, i)[0] for i in range(N)])
optw = np.array([bod_opt(Vb, i)[1] for i in range(N)])
optw = optw / optw.sum(1, keepdims=True)
RES["bod_opt"] = opt; RES["bod_pes"] = pes
pes_n = (pes - pes.min()) / (pes.max() - pes.min()); opt_n = (opt - opt.min()) / (opt.max() - opt.min())
RES["bod_mix"] = 0.5 * opt_n + 0.5 * pes_n                     # Zhou et al. (2007) composite with lambda = 0.5
R["bod"] = {"n_opt_efficient": int((opt >= 1 - 1e-6).sum()), "n_pes_worst": int((pes <= 1 + 1e-6).sum()),
            "opt_efficient": [ISO[i] for i in np.where(opt >= 1 - 1e-6)[0]], "pes_worst": [ISO[i] for i in np.where(pes <= 1 + 1e-6)[0]]}
# equal-weight benchmark
RES["equal_weights"] = V.mean(1)
RES["rank_equal"] = RES.equal_weights.rank(ascending=False).astype(int)
RES["rank_opt"] = RES.bod_opt.rank(ascending=False, method="min").astype(int)
RES["rank_mix"] = RES.bod_mix.rank(ascending=False).astype(int)
RES["rank_smaa"] = RES.hai.rank(ascending=False).astype(int)
RES = RES.sort_values("hai", ascending=False); RES.to_csv("out/rankings.csv")
R["rank_agreement"] = {f"{a}_vs_{b}": {"spearman": round(float(spearmanr(RES[a], RES[b]).statistic), 3), "kendall": round(float(kendalltau(RES[a], RES[b]).statistic), 3)}
                       for a, b in [("rank_smaa", "rank_equal"), ("rank_smaa", "rank_opt"), ("rank_smaa", "rank_mix"), ("rank_opt", "rank_mix")]}
R["top10_smaa"] = RES.head(10)[["name", "mean_rank", "best_rank", "worst_rank", "p_top10", "hai", "bod_opt", "bod_pes"]].round(3).to_dict(orient="index")
R["bottom10_smaa"] = RES.tail(10)[["name", "mean_rank", "best_rank", "worst_rank", "p_bottom10", "hai"]].round(3).to_dict(orient="index")
R["robust_top10"] = RES.index[RES.p_top10 >= 0.95].tolist(); R["robust_bottom10"] = RES.index[RES.p_bottom10 >= 0.95].tolist()
R["rank_range_mean"] = round(float((RES.worst_rank - RES.best_rank).mean()), 1)
log("BoD done", R["bod"])

# =============================== Turkiye focus ===============================
tr = RES.loc["TUR"]
R["turkiye"] = {"criteria": C.loc["TUR", CRIT].astype(float).round(3).to_dict(), "percentile_by_criterion": {c: round(float((C[c] * DIR[c] <= C.loc["TUR", c] * DIR[c]).mean()), 3) for c in CRIT},
                "smaa": {k: (round(float(v), 3) if isinstance(v, (float, np.floating)) else v) for k, v in tr.items()}}
# what-if scenarios for Turkiye (only Turkiye's values change; normalisation bounds kept fixed when inside the observed range)
def scen(mod):
    Xs = X.copy()
    for c, f in mod.items():
        j = CRIT.index(c); Xs[TUR, j] = f(Xs[TUR, j])
    Vs = normalise(Xs); Ss = Vs @ W.T; rk = (Ss > Ss[TUR]).sum(0) + 1
    return {"mean_rank": round(float(rk.mean()), 1), "best": int(rk.min()), "worst": int(rk.max()), "p_top_quartile": round(float((rk <= N / 4).mean()), 3)}
TS = {"observed": scen({}),
      "wind_solar_plus_10pp": scen({"ws_share": lambda x: x + 10, "ws_change": lambda x: x + 10}),
      "electricity_intensity_minus_25pct": scen({"ci_elec": lambda x: 0.75 * x}),
      "carbon_price_eu_ets_like_70usd": scen({"ctax": lambda x: np.log1p(70.0)}),
      "combined": scen({"ws_share": lambda x: x + 10, "ws_change": lambda x: x + 10, "ci_elec": lambda x: 0.75 * x, "ctax": lambda x: np.log1p(70.0)})}
R["turkiye_scenarios"] = TS
log("Turkiye", TS)

# =============================== robustness ===============================
ROB = {}
for lab, (lo, hi) in {"no preference (0-1)": (0.0, 1.0), "narrow AR (0.10-0.20)": (0.10, 0.20)}.items():
    Wr = hit_and_run(NS // 2, lo=lo, hi=hi); _, rk, ba, _, _ = smaa(V, Wr); hai = ba @ meta
    ROB[lab] = {"spearman_with_main_hai": round(float(spearmanr(hai, BA @ meta).statistic), 3), "turkiye_mean_rank": round(float(rk[TUR].mean() + 1), 1),
                "top5": [ISO[i] for i in np.argsort(-hai)[:5]]}
# measurement uncertainty: multiplicative lognormal noise (CO2 and GDP 5%, electricity indicators 2%), resampled jointly with weights
noise_sd = {"co2_pc": 0.05, "co2_change": 0.0, "co2_gdp": 0.07, "ci_elec": 0.02, "ws_share": 0.02, "ws_change": 0.0, "ctax": 0.0}
M_ = 1000 if not QUICK else 100; rk_all = np.zeros((N, M_))
Wn = hit_and_run(M_)
for m_ in range(M_):
    Xn = X.copy()
    for j, c in enumerate(CRIT):
        if noise_sd[c] > 0: Xn[:, j] = X[:, j] * np.exp(rng.normal(0, noise_sd[c], N))
    co2_15 = co2.Y_2015.reindex(ISO).values; co2_25 = co2.Y_2025.reindex(ISO).values
    Xn[:, CRIT.index("co2_change")] = np.log(co2_25 * np.exp(rng.normal(0, 0.03, N)) / (co2_15 * np.exp(rng.normal(0, 0.03, N)))) * 100
    s = normalise(Xn) @ Wn[m_]; rk_all[:, m_] = (-s).argsort().argsort()
ROB["with_measurement_noise"] = {"spearman_mean_rank_vs_main": round(float(spearmanr(rk_all.mean(1), RK.mean(1)).statistic), 3),
                                 "turkiye_mean_rank": round(float(rk_all[TUR].mean() + 1), 1), "turkiye_rank_5_95": [int(np.percentile(rk_all[TUR], 5) + 1), int(np.percentile(rk_all[TUR], 95) + 1)]}
# leave-one-criterion-out
for c in CRIT:
    keep = [j for j, cc in enumerate(CRIT) if cc != c]; Vr = V[:, keep]
    k2 = len(keep); Wr = hit_and_run(NS // 4, lo=min(LO, 1 / k2), hi=max(HI, 1 / k2), k=k2); _, _, ba, _, _ = smaa(Vr, Wr)
    ROB[f"drop_{c}"] = {"spearman_with_main_hai": round(float(spearmanr(ba @ meta, BA @ meta).statistic), 3)}
R["robustness"] = ROB
log("robustness", ROB)
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# =============================== figures ===============================
order = RES.index.tolist()
fig, ax = plt.subplots(figsize=(7.5, 9.6))
P = pd.DataFrame(BA, index=ISO).loc[order]
cum = np.zeros(N); bands = [(0, 10, "#1a9850", "Ranks 1–10"), (10, 30, "#91cf60", "11–30"), (30, 50, "#fee08b", "31–50"), (50, 70, "#fc8d59", "51–70"), (70, N, "#d73027", f"71–{N}")]
for a_, b_, col, lab in bands:
    v = P.iloc[:, a_:b_].sum(1).values; ax.barh(range(N), v, left=cum, color=col, label=lab, height=0.85); cum += v
ax.set_yticks(range(N)); ax.set_yticklabels([f"{C.loc[i, 'name']}" for i in order], fontsize=5.2); ax.invert_yaxis(); ax.set_xlim(0, 1)
for t in ax.get_yticklabels():
    if t.get_text() == C.loc["TUR", "name"]: t.set_fontweight("bold"); t.set_color("#d62728")
ax.set_ylim(N - 0.5, -0.5); ax.set_xlabel("Rank acceptability (share of feasible weight vectors)"); ax.legend(frameon=False, ncol=5, fontsize=6.5, loc="lower center", bbox_to_anchor=(0.5, 1.0))
fig.tight_layout(); fig.savefig("fig/fig1_rank_acceptability.png"); plt.show()

fig, ax = plt.subplots(figsize=(6.3, 4.4))
ax.scatter(RES.bod_pes, RES.bod_opt, s=12, c=RES.hai, cmap="viridis")
for i, r in RES.iterrows():
    if i in RES.index[:8].tolist() + RES.index[-5:].tolist() + ["TUR", "DEU", "USA", "CHN", "IND", "FRA", "GBR"]:
        ax.annotate(i, (r.bod_pes, r.bod_opt), fontsize=6, xytext=(2, 2), textcoords="offset points", color="#d62728" if i == "TUR" else "k")
ax.set_xlabel("Pessimistic BoD score (distance from worst practice)"); ax.set_ylabel("Optimistic BoD score (best practice = 1)")
sm = plt.cm.ScalarMappable(cmap="viridis", norm=plt.Normalize(RES.hai.min(), RES.hai.max())); fig.colorbar(sm, ax=ax, label="SMAA holistic acceptability")
fig.tight_layout(); fig.savefig("fig/fig2_bod.png"); plt.show()

fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.4), gridspec_kw={"width_ratios": [1, 1.1]})
ax = axs[0]; top = RES.index[:12]
cwm = pd.DataFrame(CW, index=ISO, columns=CRIT).loc[[i for i in top if np.isfinite(CW[list(ISO).index(i)]).all()]]
im = ax.imshow(cwm.values, cmap="Blues", vmin=0, vmax=HI); ax.set_yticks(range(len(cwm))); ax.set_yticklabels(cwm.index, fontsize=6.5)
ax.set_xticks(range(K)); ax.set_xticklabels([SHORT[c] for c in CRIT], rotation=60, ha="right", fontsize=6.5); ax.set_title("(a) Central weights", loc="left")
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.03, label="Weight")
ax = axs[1]; sc = list(TS); vals = [TS[s]["mean_rank"] for s in sc]; lo_ = [TS[s]["best"] for s in sc]; hi_ = [TS[s]["worst"] for s in sc]
ax.errorbar(vals, range(len(sc)), xerr=[np.array(vals) - lo_, np.array(hi_) - vals], fmt="o", color="#d62728", capsize=2)
ax.set_yticks(range(len(sc))); ax.set_yticklabels(["Observed 2025", "+10 pp wind & solar", "−25% electricity CO$_2$ intensity", "Carbon price US$70/t", "All three"], fontsize=7)
ax.invert_yaxis(); ax.invert_xaxis(); ax.set_xlabel(f"Rank among {N} (mean; bar = best to worst)"); ax.set_title("(b) Türkiye: what-if scenarios", loc="left")
fig.tight_layout(); fig.savefig("fig/fig3_weights_turkiye.png"); plt.show()

fig, ax = plt.subplots(figsize=(6.0, 3.6))
ax.scatter(RES.rank_equal, RES.mean_rank, s=10, color="#1f77b4")
ax.errorbar(RES.rank_equal, RES.mean_rank, yerr=[RES.mean_rank - RES.best_rank, RES.worst_rank - RES.mean_rank], fmt="none", lw=0.4, color="#9ecae1")
ax.plot([1, N], [1, N], color="k", lw=0.5); ax.annotate("TUR", (RES.loc["TUR", "rank_equal"], RES.loc["TUR", "mean_rank"]), color="#d62728", fontsize=7)
ax.set_xlabel("Rank under equal weights"); ax.set_ylabel("SMAA mean rank (best–worst)"); fig.tight_layout(); fig.savefig("fig/fig4_rank_ranges.png"); plt.show()
log("done")
print(json.dumps({k: R[k] for k in ["n_countries", "bod", "rank_agreement", "robust_top10", "robust_bottom10", "rank_range_mean", "turkiye", "turkiye_scenarios", "robustness"]}, indent=1, default=str))

# ## Sonuçlar ve çıktılar

# Packs figures and numerical results for download.
import shutil
from google.colab import files
shutil.make_archive("/content/rp10_outputs", "zip", ".", ".")
files.download("/content/rp10_outputs.zip")