# -*- coding: utf-8 -*-
"""Wind droughts and the allocation of offshore wind capacity across European seas: Sample-average and distributionally robust conditional value-at-risk models

GOMD arastirma makalesi 07 - 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.
"""

# # Wind droughts and the allocation of offshore wind capacity across European seas: Sample-average and distributionally robust conditional value-at-risk models
#
# 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: 15–25 dakika (yaklaşık 2.500 doğrusal program).
#
# *Hızlı deneme:* ilk kod hücresinden sonra `import sys; sys.argv.append('--quick')` çalıştırılırsa önyükleme tekrar sayısı azaltılır (sonuçlar makaledekinden farklı olur).
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: NASA POWER günlük rüzgâr verisi). 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

VERI_KAYNAGI = "api"   # "api": download from the original sources; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp07", exist_ok=True); os.chdir("/content/rp07")
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')

# --- sites.py dosyasi yaziliyor (notebook'taki %%writefile hucresi) ---
open('sites.py', 'w', encoding='utf-8').write('"""Point locations used to query NASA POWER (approximate coordinates, decimal degrees)."""\n# Paper 7: offshore wind sites in European seas (existing, planned or studied offshore wind zones)\nOFFSHORE = [\n ("Dogger Bank", "GBR", "North Sea", 54.75, 1.90), ("Hornsea", "GBR", "North Sea", 53.90, 1.80), ("East Anglia", "GBR", "North Sea", 52.65, 2.60),\n ("Moray Firth", "GBR", "North Sea", 58.20, -2.80), ("Morecambe Bay", "GBR", "Irish Sea", 53.90, -3.70), ("Celtic Sea", "IRL", "Celtic Sea", 51.10, -6.30),\n ("Borssele", "NLD", "North Sea", 51.70, 3.00), ("IJmuiden Ver", "NLD", "North Sea", 52.85, 3.80), ("German Bight", "DEU", "North Sea", 54.30, 6.50),\n ("Horns Rev", "DNK", "North Sea", 55.50, 7.70), ("Sorlige Nordsjo II", "NOR", "North Sea", 56.80, 4.90), ("Utsira Nord", "NOR", "North Sea", 59.30, 4.50),\n ("Kriegers Flak", "DNK", "Baltic Sea", 55.00, 12.90), ("Baltic Power", "POL", "Baltic Sea", 54.90, 17.90), ("Gulf of Riga", "EST", "Baltic Sea", 58.30, 22.00),\n ("Gulf of Bothnia", "SWE", "Baltic Sea", 62.50, 20.50), ("South Brittany", "FRA", "Atlantic", 47.20, -3.00), ("Galicia", "ESP", "Atlantic", 43.40, -8.80),\n ("Viana do Castelo", "PRT", "Atlantic", 41.70, -9.10), ("Canary Islands", "ESP", "Atlantic", 27.90, -15.30), ("Gulf of Lion", "FRA", "Mediterranean", 43.00, 3.80),\n ("Sicily Channel", "ITA", "Mediterranean", 37.30, 12.20), ("Northern Adriatic", "ITA", "Mediterranean", 44.70, 12.90), ("Saros Gulf", "TUR", "Aegean", 40.50, 26.40),\n ("Gokceada-Bozcaada", "TUR", "Aegean", 39.90, 25.70), ("Cyclades", "GRC", "Aegean", 36.80, 25.40), ("Crete", "GRC", "Mediterranean", 35.60, 24.00),\n ("Kiyikoy", "TUR", "Black Sea", 41.75, 28.40), ("Romanian Black Sea", "ROU", "Black Sea", 44.30, 29.40), ("Cyprus South", "CYP", "Mediterranean", 34.50, 33.50),\n]\n# Paper 8: tourism destinations (two to four points per country, spanning main tourism regions)\nTOURISM = {\n "AT": [(48.21, 16.37), (47.26, 11.39), (47.80, 13.04)], "BE": [(50.85, 4.35), (51.23, 2.92)], "BG": [(43.21, 27.91), (42.70, 23.32), (42.50, 27.47)],\n "HR": [(43.51, 16.44), (42.65, 18.09), (45.81, 15.98)], "CY": [(34.99, 34.00), (34.77, 32.42)], "CZ": [(50.08, 14.44), (49.20, 16.61)],\n "DK": [(55.68, 12.57), (56.16, 10.20)], "EE": [(59.44, 24.75), (58.39, 24.50)], "FI": [(60.17, 24.94), (66.50, 25.73)],\n "FR": [(48.86, 2.35), (43.70, 7.27), (44.84, -0.58)], "DE": [(52.52, 13.40), (48.14, 11.58), (53.55, 9.99)], "EL": [(37.98, 23.73), (36.43, 28.22), (35.34, 25.14)],\n "HU": [(47.50, 19.04), (46.90, 18.05)], "IE": [(53.35, -6.26), (53.27, -9.05)], "IT": [(41.90, 12.50), (45.44, 12.32), (40.85, 14.27)],\n "LV": [(56.95, 24.11), (56.97, 23.77)], "LT": [(54.69, 25.28), (55.71, 21.13)], "LU": [(49.61, 6.13)], "MT": [(35.90, 14.51)],\n "NL": [(52.37, 4.90), (51.50, 3.60)], "PL": [(50.06, 19.94), (54.35, 18.65), (49.30, 19.95)], "PT": [(38.72, -9.14), (37.02, -7.93), (41.15, -8.61)],\n "RO": [(44.43, 26.10), (44.17, 28.64), (45.65, 25.60)], "SK": [(48.15, 17.11), (49.13, 20.22)], "SI": [(46.06, 14.51), (45.51, 13.59)],\n "ES": [(40.42, -3.70), (41.39, 2.17), (36.72, -4.42), (39.57, 2.65)], "SE": [(59.33, 18.07), (57.71, 11.97)], "IS": [(64.15, -21.94)],\n "NO": [(59.91, 10.75), (60.39, 5.32)], "CH": [(47.37, 8.54), (46.20, 6.14), (46.69, 7.86)], "TR": [(36.90, 30.70), (41.01, 28.98), (37.04, 27.43)],\n "ME": [(42.29, 18.84)], "MK": [(41.12, 20.80), (42.00, 21.43)], "AL": [(41.33, 19.82), (41.32, 19.45), (39.87, 20.00)], "RS": [(44.79, 20.45)],\n}\n# Paper 9: major container and Mediterranean ports\nPORTS = [\n ("Shanghai", 31.23, 121.47), ("Singapore", 1.26, 103.82), ("Ningbo-Zhoushan", 29.87, 121.55), ("Shenzhen", 22.54, 114.05), ("Guangzhou Nansha", 22.75, 113.60),\n ("Qingdao", 36.07, 120.38), ("Busan", 35.10, 129.04), ("Tianjin", 38.99, 117.70), ("Hong Kong", 22.30, 114.17), ("Port Klang", 3.00, 101.39),\n ("Rotterdam", 51.95, 4.14), ("Jebel Ali", 25.01, 55.06), ("Antwerp", 51.27, 4.35), ("Xiamen", 24.48, 118.08), ("Tanjung Pelepas", 1.36, 103.55),\n ("Kaohsiung", 22.61, 120.28), ("Los Angeles", 33.74, -118.27), ("Hamburg", 53.54, 9.97), ("New York-New Jersey", 40.67, -74.04), ("Laem Chabang", 13.08, 100.88),\n ("Colombo", 6.95, 79.84), ("Cai Mep", 10.53, 107.03), ("Tanjung Priok", -6.10, 106.88), ("Jeddah", 21.48, 39.17), ("Valencia", 39.45, -0.32),\n ("Piraeus", 37.94, 23.63), ("Algeciras", 36.13, -5.44), ("Nhava Sheva", 18.95, 72.95), ("Santos", -23.98, -46.30), ("Tanger Med", 35.89, -5.50),\n ("Savannah", 32.08, -81.09), ("Felixstowe", 51.95, 1.32), ("Bremerhaven", 53.55, 8.58), ("Le Havre", 49.48, 0.11), ("Barcelona", 41.35, 2.17),\n ("Gioia Tauro", 38.44, 15.90), ("Marsaxlokk", 35.83, 14.54), ("Ambarli", 40.97, 28.68), ("Mersin", 36.79, 34.64), ("Izmir", 38.44, 27.14),\n ("Port Said", 31.26, 32.30), ("Durban", -29.87, 31.03), ("Lagos Apapa", 6.44, 3.37), ("Mombasa", -4.06, 39.67), ("Djibouti", 11.60, 43.13),\n ("Salalah", 16.94, 54.00), ("Houston", 29.73, -95.27), ("Vancouver", 49.29, -123.11), ("Manzanillo", 19.05, -104.31), ("Cartagena", 10.39, -75.52),\n ("Callao", -12.05, -77.15), ("Melbourne", -37.84, 144.92), ("Sydney Botany", -33.97, 151.22), ("Tokyo", 35.62, 139.78), ("Novorossiysk", 44.72, 37.78),\n ("Constanta", 44.17, 28.66), ("Gdansk", 54.40, 18.67), ("Gothenburg", 57.69, 11.86), ("St Petersburg", 59.88, 30.21), ("Balboa", 8.95, -79.57),\n]\n')

# Downloads the inputs from the original sources (only when VERI_KAYNAGI == 'api').
if VERI_KAYNAGI == 'api':
    """Re-downloads daily wind speeds at 10 m and 50 m (NASA POWER, MERRA-2 based) for the 30 offshore sites, 1 Jan 1981 - 31 Oct 2026, into data/."""
    import pandas as pd
    from common import power_daily, DATA
    from sites import OFFSHORE
    w50, w10 = {}, {}
    for name, iso, sea, lat, lon in OFFSHORE:
        d = power_daily(lat, lon, ["WS50M", "WS10M"]); w50[name] = pd.Series(d["WS50M"]); w10[name] = pd.Series(d["WS10M"]); print(name)
    for nm, w in (("ws50m", w50), ("ws10m", w10)):
        df = pd.DataFrame(w); df.index.name = "date"; df = df[sorted(df.columns)]; df.to_csv(DATA / f"{nm}_daily.csv.gz", compression="gzip")
    pd.DataFrame(OFFSHORE, columns=["site", "iso", "sea", "lat", "lon"]).to_csv(DATA / "offshore_sites.csv", index=False)
print(sorted(os.listdir('data')))

# ## 2. Kapasite faktörleri, örneklem-ortalaması ve Wasserstein-dayanıklı CVaR modelleri, kayan değerlendirme

# Paper 7: wind-drought risk and the allocation of offshore wind capacity across European seas. Daily capacity factors at 30 offshore sites (NASA POWER / MERRA-2 winds at 10 m and 50 m, 1 Jan 1981 - 31 Oct 2026) for a 15 MW reference turbine; weekly lower-tail risk (CVaR) minimised by sample-average and Wasserstein distributionally robust linear programs (exact reformulation on bounded support, Mohajerin Esfahani & Kuhn, 2018); rolling out-of-sample evaluation 2001-2026 with training records of 3, 5, 10 and 20 years.
import json, sys, time, warnings
import numpy as np, pandas as pd
from math import gamma
from scipy import sparse
from scipy.optimize import linprog, minimize
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
T0 = time.time(); R = {}; rng = np.random.default_rng(7)
def log(*a): print(f"[{time.time() - T0:6.0f}s]", *a, flush=True)
ALPHA = 0.05                 # CVaR tail: worst 5% of weeks (about 2.6 weeks per year)
HUB, VIN, VR, VOUT = 150.0, 3.0, 10.59, 25.0   # IEA 15 MW reference turbine (Gaertner et al., 2020)
K_DAY = 4.5                  # Weibull shape of intra-day wind speeds around the daily mean (coefficient of variation about 0.25)
WINDOWS = [3, 5, 10, 20]     # training record lengths (years)
EPS_GRID = [0.0, 0.005, 0.01, 0.02, 0.05]
TEST_YEARS = list(range(2001, 2027)) if not QUICK else [2003, 2010, 2021, 2026]
SEA_ORDER = ["North Sea", "Irish Sea", "Celtic Sea", "Baltic Sea", "Atlantic", "Mediterranean", "Aegean", "Black Sea"]

# =============================== capacity factors ===============================
meta = pd.read_csv("data/offshore_sites.csv").set_index("site")[["iso", "sea", "lat", "lon"]]
w50 = pd.read_csv("data/ws50m_daily.csv.gz", index_col="date"); w10 = pd.read_csv("data/ws10m_daily.csv.gz", index_col="date")
assert (w50.index == w10.index).all()
w50.index = pd.to_datetime(w50.index.astype(str), format="%Y%m%d"); w10.index = w50.index
order = meta.assign(o=meta.sea.map({s: i for i, s in enumerate(SEA_ORDER)})).sort_values(["o", "lat"], ascending=[True, False]).index
w50, w10, meta = w50[order], w10[order], meta.loc[order]; SITES = list(order); NS = len(SITES)
clim = (w50.index.year >= 1991) & (w50.index.year <= 2020)
shear = (np.log(w50[clim] / w10[clim]) / np.log(5.0)).groupby(w50.index[clim].month).mean().clip(0.04, 0.30)
hub = w50 * (HUB / 50.0) ** shear.loc[w50.index.month].values
R["shear_alpha"] = {"mean": round(float(shear.values.mean()), 3), "min": round(float(shear.values.min()), 3), "max": round(float(shear.values.max()), 3)}
VG = np.linspace(0, 45, 4501)
def power_curve(v):
    p = np.where((v >= VIN) & (v < VR), (v ** 3 - VIN ** 3) / (VR ** 3 - VIN ** 3), 0.0)
    return np.where((v >= VR) & (v <= VOUT), 1.0, p)
def cf_from_mean(ws, k):
    means = np.linspace(0.05, 35, 700); tab = []
    for m in means:
        c = m / gamma(1 + 1 / k); f = (k / c) * (VG / c) ** (k - 1) * np.exp(-(VG / c) ** k); tab.append(np.trapezoid(power_curve(VG) * f, VG))
    return pd.DataFrame(np.interp(ws.values, means, np.array(tab)), index=ws.index, columns=ws.columns)
CFD = cf_from_mean(hub, K_DAY)
def weekly(cfd):
    wk = cfd.resample("W-SUN").mean(); n = cfd.resample("W-SUN").size(); return wk[n == 7]
CFW = weekly(CFD); YW = CFW.index.year.values; XW = CFW.values
CFW.to_csv("out/cf_weekly.csv")
R["data"] = {"days": int(len(CFD)), "first_day": str(CFD.index[0].date()), "last_day": str(CFD.index[-1].date()), "weeks": int(len(CFW)),
             "cf_mean_all": round(float(CFD.values.mean()), 3), "cf_best": [CFD.mean().idxmax(), round(float(CFD.mean().max()), 3)],
             "cf_worst": [CFD.mean().idxmin(), round(float(CFD.mean().min()), 3)]}
log("data", R["data"])

# =============================== descriptive statistics ===============================
desc = pd.DataFrame({"sea": meta.sea, "iso": meta.iso, "lat": meta.lat, "lon": meta.lon, "ws_hub": hub.mean().round(2), "cf_mean": CFD.mean().round(3),
                     "cf_week_p05": CFW.quantile(0.05).round(3), "low_days": (CFD < 0.10).mean().round(3),
                     "winter_summer": (CFD[CFD.index.month.isin([12, 1, 2])].mean() / CFD[CFD.index.month.isin([6, 7, 8])].mean()).round(2)})
ann = CFD.groupby(CFD.index.year).mean().loc[1981:2025]
def sen(y):
    x = np.arange(len(y)); s = [(y[j] - y[i]) / (j - i) for i in range(len(y)) for j in range(i + 1, len(y))]; return float(np.median(s))
desc["sen_per_decade"] = [round(sen(ann[s].values) * 10, 4) for s in SITES]
desc.to_csv("out/site_table.csv")
anom = CFW - CFW.groupby(CFW.index.isocalendar().week.values).transform("mean")
ca = anom.corr(); cr = CFW.corr(); seas = meta.sea.values
def pair_mean(C, same):
    return round(float(np.mean([C.iloc[i, j] for i in range(NS) for j in range(i + 1, NS) if (seas[i] == seas[j]) == same])), 3)
ns_idx = meta.index[meta.sea == "North Sea"]; ae_idx = meta.index[meta.sea == "Aegean"]
R["corr"] = {"raw_within": pair_mean(cr, True), "raw_across": pair_mean(cr, False), "anom_within": pair_mean(ca, True), "anom_across": pair_mean(ca, False),
             "raw_north_aegean": round(float(cr.loc[ns_idx, ae_idx].values.mean()), 3), "anom_north_aegean": round(float(ca.loc[ns_idx, ae_idx].values.mean()), 3)}
ew_d = CFD.mean(1); ns_d = CFD[ns_idx].mean(1)
R["low_output_days"] = {"north_sea_mean_cf_lt_0.10": round(float((ns_d < 0.10).mean()), 4), "all30_mean_cf_lt_0.10": round(float((ew_d < 0.10).mean()), 4)}
log("descriptives", R["corr"])

# =============================== optimisation models ===============================
def cvar_lp(Xi, eps=0.0, alpha=ALPHA, cap=1.0, mean_floor=None, allowed=None, fix_w=None, return_val=False):
    """min_w sup_{Q: W1(Q, P_N) <= eps, supp Q in [0,1]^n} CVaR_alpha^Q(-w'xi), ground metric ||.||_1 (dual norm ||.||_inf).
    Exact LP reformulation (Mohajerin Esfahani & Kuhn, 2018, Corollary 5.1); eps = 0 is the sample-average (SAA) problem."""
    Nn, n = Xi.shape; ia = 1.0 / alpha; G = eps > 0
    nv = n + 2 + Nn + (2 * Nn * n if G else 0)
    iw = np.arange(n); itau = n; ilam = n + 1; isv = n + 2 + np.arange(Nn)
    c = np.zeros(nv); c[ilam] = eps; c[isv] = 1.0 / Nn
    I = np.arange(Nn)
    # (1) tau - s_i <= 0
    r1 = np.repeat(I, 2); c1 = np.column_stack([np.full(Nn, itau), isv]).ravel(); v1 = np.tile([1.0, -1.0], Nn)
    # (2) tau(1 - 1/a) - (1/a) w'xi_i + sum_j [gp_ij (1 - xi_ij) + gm_ij xi_ij] - s_i <= 0
    r2 = Nn + np.repeat(I, 2 + n); c2 = np.column_stack([np.full(Nn, itau), isv, np.tile(iw, (Nn, 1))]).ravel()
    v2 = np.column_stack([np.full(Nn, 1 - ia), np.full(Nn, -1.0), -ia * Xi]).ravel()
    rows, cols, vals = [r1, r2], [c1, c2], [v1, v2]; nr = 2 * Nn
    if G:
        gp = n + 2 + Nn + np.arange(Nn * n).reshape(Nn, n); gm = gp + Nn * n
        rows += [Nn + np.repeat(I, 2 * n)]; cols += [np.column_stack([gp, gm]).ravel()]; vals += [np.column_stack([1 - Xi, Xi]).ravel()]
        # (3) |gp_ij - gm_ij + w_j / a| <= lam
        k = np.arange(Nn * n); rr = nr + k; jj = np.tile(iw, Nn)
        for sgn, off in [(1.0, 0), (-1.0, Nn * n)]:
            rows += [np.repeat(rr + off, 4)]; cols += [np.column_stack([gp.ravel(), gm.ravel(), jj, np.full(Nn * n, ilam)]).ravel()]
            vals += [np.tile([sgn, -sgn, sgn * ia, -1.0], Nn * n)]
        nr += 2 * Nn * n
    b = np.zeros(nr)
    if mean_floor is not None:
        rows += [np.full(n, nr)]; cols += [iw]; vals += [-Xi.mean(0)]; b = np.append(b, -mean_floor); nr += 1
    Aub = sparse.csr_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))), shape=(nr, nv))
    Aeq = sparse.csr_matrix((np.ones(n), (np.zeros(n, int), iw)), shape=(1, nv))
    ub = np.full(n, cap) if allowed is None else np.where(allowed, cap, 0.0)
    wb = [(0, u) for u in ub] if fix_w is None else [(float(x), float(x)) for x in fix_w]
    bounds = wb + [(None, None), (0, None if G else 0)] + [(None, None)] * Nn + ([(0, None)] * (2 * Nn * n) if G else [])
    res = linprog(c, A_ub=Aub, b_ub=b, A_eq=Aeq, b_eq=[1.0], bounds=bounds, method="highs")
    if res.status != 0: raise RuntimeError(res.message)
    w = np.clip(res.x[:n], 0, None); w = w / w.sum()
    return (w, float(res.fun)) if return_val else w

def lower_tail(x, alpha=ALPHA):
    """Lower-tail mean: average capacity factor over the worst alpha share of periods (= -CVaR of the loss -CF)."""
    x = np.sort(np.asarray(x)); k = max(1, int(np.ceil(alpha * len(x)))); return float(x[:k].mean())

def minvar(Xi):
    S = np.cov(Xi.T); n = S.shape[0]
    r = minimize(lambda w: w @ S @ w, np.full(n, 1 / n), jac=lambda w: 2 * S @ w, bounds=[(0, 1)] * n,
                 constraints=[{"type": "eq", "fun": lambda w: w.sum() - 1}], method="SLSQP", options={"maxiter": 1000, "ftol": 1e-14})
    w = np.clip(r.x, 0, None); return w / w.sum()

def best_site(Xi):
    w = np.zeros(Xi.shape[1]); w[np.argmax(Xi.mean(0))] = 1.0; return w

def cv_radius(Xi, yrs, grid=EPS_GRID):
    """Radius by leave-one-block-out cross-validation over contiguous years (3 blocks): eps maximising the held-out lower-tail mean."""
    uy = np.unique(yrs); blocks = np.array_split(uy, min(3, len(uy))); sc = {e: [] for e in grid}
    for bl in blocks:
        te = np.isin(yrs, bl)
        for e in grid: sc[e].append(lower_tail(Xi[te] @ cvar_lp(Xi[~te], e)))
    sc = {e: float(np.mean(v)) for e, v in sc.items()}; best = max(sc.values())
    return min(e for e, v in sc.items() if v >= best - 1e-9), sc        # smallest radius among ties

# =============================== proposition check: soft cap ===============================
# For radii small enough that no tail scenario reaches the support boundary, the worst-case CVaR equals the SAA CVaR
# plus (eps/alpha) * max_j w_j (Mohajerin Esfahani & Kuhn, 2018; Kuhn et al., 2019). We verify this numerically at a fixed w.
Xchk = XW[(YW >= 2016) & (YW <= 2025)]; chk = {}
for lab, wv in [("equal weights", np.full(NS, 1 / NS)), ("SAA optimum", None)]:
    if wv is None: wv = cvar_lp(Xchk, 0.0)
    saa = cvar_lp(Xchk, 0.0, fix_w=wv, return_val=True)[1]
    chk[lab] = {}
    for e in (0.001, 0.005, 0.01, 0.05, 0.2):
        lp = cvar_lp(Xchk, e, fix_w=wv, return_val=True)[1]; cf = saa + e / ALPHA * wv.max()
        chk[lab][str(e)] = {"lp": round(lp, 5), "closed_form": round(cf, 5)}
R["softcap_check"] = chk
log("soft-cap check", chk)

# =============================== in-sample frontier (full record) ===============================
FR = {}
groups = {"All seas": None, "Northern seas only (North, Irish, Celtic, Baltic)": meta.sea.isin(SEA_ORDER[:4]).values,
          "Southern seas only (Atlantic, Mediterranean, Aegean, Black)": meta.sea.isin(SEA_ORDER[4:]).values}
for lab, allowed in groups.items():
    Xg = XW if allowed is None else XW[:, allowed]
    w0 = cvar_lp(XW, 0.0, allowed=allowed); m0 = float(XW.mean(0) @ w0); mmax = float(Xg.mean(0).max())
    pts = [(m0, lower_tail(XW @ w0))]
    for m in np.linspace(m0, mmax, 8)[1:]:
        w = cvar_lp(XW, 0.0, mean_floor=m - 1e-6, allowed=allowed); pts.append((float(XW.mean(0) @ w), lower_tail(XW @ w)))
    FR[lab] = pts
R["frontier"] = FR
w_full = cvar_lp(XW, 0.0)
R["full_record_saa"] = {s: round(float(v), 3) for s, v in zip(SITES, w_full) if v > 1e-3}
R["full_record_saa_sea"] = {s: round(float(w_full[meta.sea.values == s].sum()), 3) for s in SEA_ORDER}
log("frontier", R["full_record_saa"])

# =============================== rolling out-of-sample evaluation ===============================
STRATS = ["Equal weights", "Best single site", "Minimum variance", "SAA-CVaR", "DRO-CVaR"]
OOS = {}; WH = {}; EPS = {}; CVS = {}
for win in WINDOWS:
    OOS[win] = {s: [] for s in STRATS}; WH[win] = {s: {} for s in STRATS}; EPS[win] = {}; CVS[win] = {}
    for ty in TEST_YEARS:
        trm = (YW >= ty - win) & (YW <= ty - 1); tem = YW == ty
        Xtr = XW[trm]; e_star, sc = cv_radius(Xtr, YW[trm]); EPS[win][ty] = e_star; CVS[win][ty] = sc
        ws = {"Equal weights": np.full(NS, 1 / NS), "Best single site": best_site(Xtr), "Minimum variance": minvar(Xtr),
              "SAA-CVaR": cvar_lp(Xtr, 0.0)}
        ws["DRO-CVaR"] = ws["SAA-CVaR"] if e_star == 0 else cvar_lp(Xtr, e_star)
        for s in STRATS: WH[win][s][ty] = ws[s]; OOS[win][s].append((ty, ws[s]))
    log("window", win, "eps selected", EPS[win])

def apply_weights(pairs, X, Y, idx):
    """Out-of-sample series: weights chosen with data up to year ty-1 are applied in year ty."""
    out = []
    for ty, w in pairs:
        m = Y == ty; out.append(pd.Series(X[m] @ w, index=idx[m]))
    return pd.concat(out)
DY = CFD.index.year.values; XD = CFD.values
def drought_spell(x, thr=0.20):
    best = cur = 0
    for v in x:
        cur = cur + 1 if v < thr else 0; best = max(best, cur)
    return best
def metrics(wk, dy):
    return {"mean": round(float(wk.mean()), 4), "lt05_week": round(lower_tail(wk.values, 0.05), 4), "lt10_week": round(lower_tail(wk.values, 0.10), 4),
            "min_week": round(float(wk.min()), 4), "lt05_day": round(lower_tail(dy.values, 0.05), 4), "share_days_lt_0.15": round(float((dy < 0.15).mean()), 4),
            "longest_spell_lt_0.20_days": int(drought_spell(dy.values)), "mean_spell_days_per_year_lt_0.20": round(float((dy < 0.20).groupby(dy.index.year).sum().mean()), 2)}
MET = {}; SER = {}
for win in WINDOWS:
    MET[win] = {}; SER[win] = {}
    for s in STRATS:
        wk = apply_weights(OOS[win][s], XW, YW, CFW.index); dy = apply_weights(OOS[win][s], XD, DY, CFD.index)
        MET[win][s] = metrics(wk, dy); SER[win][s] = wk
    log("window", win, {s: (MET[win][s]["mean"], MET[win][s]["lt05_week"]) for s in STRATS})
R["oos"] = {str(w): MET[w] for w in WINDOWS}; R["eps_selected"] = {str(w): EPS[w] for w in WINDOWS}
R["eps_selected_share_positive"] = {str(w): round(float(np.mean([e > 0 for e in EPS[w].values()])), 3) for w in WINDOWS}
pd.DataFrame({f"{w}_{s}": SER[w][s] for w in WINDOWS for s in STRATS}).to_csv("out/oos_weekly.csv")

def block_boot(a, b, stat, B=(2000 if not QUICK else 50), L=26):
    n = len(a); nb = int(np.ceil(n / L)); out = np.empty(B)
    for k in range(B):
        st = rng.integers(0, n - L + 1, nb); idx = np.concatenate([np.arange(s, s + L) for s in st])[:n]; out[k] = stat(a[idx]) - stat(b[idx])
    return {"diff": round(float(stat(a) - stat(b)), 4), "lo": round(float(np.quantile(out, 0.025)), 4), "hi": round(float(np.quantile(out, 0.975)), 4),
            "p": round(float(min(1.0, 2 * min((out <= 0).mean(), (out >= 0).mean()))), 4)}
TESTS = {}
for win in WINDOWS:
    TESTS[str(win)] = {}
    for a_, b_ in [("DRO-CVaR", "SAA-CVaR"), ("SAA-CVaR", "Equal weights"), ("SAA-CVaR", "Minimum variance"), ("SAA-CVaR", "Best single site"), ("DRO-CVaR", "Equal weights")]:
        A_, B_ = SER[win][a_].values, SER[win][b_].values
        TESTS[str(win)][f"{a_} vs {b_}"] = {"lt05": block_boot(A_, B_, lambda x: lower_tail(x, 0.05)), "mean": block_boot(A_, B_, np.mean)}
R["tests"] = TESTS
def sea_share(w): return {s: round(float(w[meta.sea.values == s].sum()), 3) for s in SEA_ORDER}
R["sea_share_avg"] = {str(w): {s: sea_share(np.mean([WH[w][s][t] for t in TEST_YEARS], axis=0)) for s in ["Minimum variance", "SAA-CVaR", "DRO-CVaR"]} for w in WINDOWS}
R["n_sites"] = {str(w): {s: round(float(np.mean([(WH[w][s][t] > 1e-3).sum() for t in TEST_YEARS])), 1) for s in STRATS} for w in WINDOWS}
R["max_share"] = {str(w): {s: round(float(np.mean([WH[w][s][t].max() for t in TEST_YEARS])), 3) for s in STRATS} for w in WINDOWS}
R["turnover"] = {str(w): {s: round(float(np.mean([np.abs(WH[w][s][TEST_YEARS[k]] - WH[w][s][TEST_YEARS[k - 1]]).sum() / 2 for k in range(1, len(TEST_YEARS))])), 3)
                          for s in STRATS} for w in WINDOWS}
TUR = meta.iso.values == "TUR"
R["turkish_share"] = {str(w): {s: round(float(np.mean([WH[w][s][t][TUR].sum() for t in TEST_YEARS])), 3) for s in ["SAA-CVaR", "DRO-CVaR"]} for w in WINDOWS}
R["site_freq_10y"] = {s: round(float(np.mean([WH[10]["SAA-CVaR"][t][j] for t in TEST_YEARS])), 3) for j, s in enumerate(SITES)}
# fixed-radius path (10-year window): how robustness changes allocations and performance
FIXP = {}
for e in [0.0, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2]:
    pairs = [(ty, cvar_lp(XW[(YW >= ty - 10) & (YW <= ty - 1)], e)) for ty in TEST_YEARS]
    wk = apply_weights(pairs, XW, YW, CFW.index)
    FIXP[str(e)] = {"mean": round(float(wk.mean()), 4), "lt05": round(lower_tail(wk.values), 4), "max_share": round(float(np.mean([w.max() for _, w in pairs])), 3),
                    "n_sites": round(float(np.mean([(w > 1e-3).sum() for _, w in pairs])), 1)}
R["fixed_radius_10y"] = FIXP
log("fixed radius", FIXP)

# =============================== sensitivity (10-year window, SAA vs DRO vs equal weights) ===============================
SENS = {}
def rerun(Xw_, Yw_, alpha=ALPHA, cap=1.0):
    out = {"Equal weights": [], "SAA-CVaR": [], "DRO-CVaR": []}
    for ty in TEST_YEARS:
        trm = (Yw_ >= ty - 10) & (Yw_ <= ty - 1); tem = Yw_ == ty; e = EPS[10][ty]
        out["Equal weights"] += list(Xw_[tem] @ np.full(NS, 1 / NS)); ws = cvar_lp(Xw_[trm], 0.0, alpha=alpha, cap=cap)
        out["SAA-CVaR"] += list(Xw_[tem] @ ws); out["DRO-CVaR"] += list(Xw_[tem] @ (ws if e == 0 else cvar_lp(Xw_[trm], e, alpha=alpha, cap=cap)))
    return {s: {"mean": round(float(np.mean(v)), 4), "lt": round(lower_tail(v, alpha), 4)} for s, v in out.items()}
for k in (3.0, 6.0):
    cw = weekly(cf_from_mean(hub, k)); SENS[f"Intra-day Weibull shape {k:g}"] = rerun(cw.values, cw.index.year.values)
hub_fixed = w50 * (HUB / 50.0) ** 0.11; cw = weekly(cf_from_mean(hub_fixed, K_DAY)); SENS["Fixed shear exponent 0.11"] = rerun(cw.values, cw.index.year.values)
SENS["Tail 10% of weeks"] = rerun(XW, YW, alpha=0.10)
SENS["Site cap 10%"] = rerun(XW, YW, cap=0.10)
cm = CFD.resample("MS").mean(); SENS["Monthly resolution, tail 10%"] = rerun(cm.values, cm.index.year.values, alpha=0.10)
R["sensitivity"] = SENS
log("sensitivity", SENS)
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# =============================== figures ===============================
sea_col = {"North Sea": "#1f77b4", "Irish Sea": "#17becf", "Celtic Sea": "#9edae5", "Baltic Sea": "#2ca02c", "Atlantic": "#9467bd",
           "Mediterranean": "#ff7f0e", "Aegean": "#d62728", "Black Sea": "#8c564b"}
fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.6), gridspec_kw={"width_ratios": [1, 1.15]})
ax = axs[0]; doy = CFD.groupby(CFD.index.dayofyear).mean().rolling(15, center=True, min_periods=1).mean()
for s in SEA_ORDER:
    ax.plot(doy.index, doy[meta.index[meta.sea == s]].mean(1), color=sea_col[s], lw=1.2, label=s)
ax.set_xticks([1, 60, 121, 182, 244, 305]); ax.set_xticklabels(["Jan", "Mar", "May", "Jul", "Sep", "Nov"]); ax.set_ylabel("Mean daily capacity factor")
ax.set_title("(a) Seasonal cycle by sea, 1981–2026", loc="left"); ax.set_ylim(0.0, 0.78); ax.legend(frameon=False, fontsize=6.0, ncol=2, loc="lower center")
ax = axs[1]; im = ax.imshow(ca.values, cmap="RdBu_r", vmin=-1, vmax=1)
ax.set_xticks(range(NS)); ax.set_yticks(range(NS)); ax.set_xticklabels([str(i + 1) for i in range(NS)], fontsize=4.6); ax.set_yticklabels([f"{i + 1} {s}" for i, s in enumerate(SITES)], fontsize=4.6)
bd = np.cumsum([0] + [int((meta.sea == s).sum()) for s in SEA_ORDER]) - 0.5
for b in bd[1:-1]: ax.axhline(b, color="k", lw=0.4); ax.axvline(b, color="k", lw=0.4)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.02, label="Correlation of weekly anomalies"); ax.set_title("(b) Correlation net of seasonality", loc="left")
fig.tight_layout(); fig.savefig("fig/fig1_seasonality_correlation.png"); plt.show()

fig, ax = plt.subplots(figsize=(6.2, 3.6))
for (lab, pts), c in zip(FR.items(), ["k", "#1f77b4", "#d62728"]):
    p = np.array(pts); ax.plot(p[:, 1], p[:, 0], marker="o", ms=3, color=c, label=lab)
for s, w, mk in [("Equal weights", np.full(NS, 1 / NS), "s"), ("Minimum variance", minvar(XW), "D"), ("Best single site", best_site(XW), "^")]:
    ax.plot(lower_tail(XW @ w), XW.mean(0) @ w, mk, color="#7f7f7f", ms=5); ax.annotate(s, (lower_tail(XW @ w), XW.mean(0) @ w), fontsize=6.5, xytext=(4, -3), textcoords="offset points")
ax.set_xlabel("Lower-tail mean capacity factor (worst 5% of weeks)"); ax.set_ylabel("Mean capacity factor"); ax.legend(frameon=False, fontsize=6.8, loc="lower right")
fig.tight_layout(); fig.savefig("fig/fig2_frontier.png"); plt.show()

fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.3))
cols = {"Equal weights": "#7f7f7f", "Best single site": "#9467bd", "Minimum variance": "#2ca02c", "SAA-CVaR": "#1f77b4", "DRO-CVaR": "#d62728"}
xs = np.arange(len(WINDOWS)); bw = 0.16
for j, s in enumerate(STRATS):
    axs[0].bar(xs + (j - 2) * bw, [MET[w][s]["lt05_week"] for w in WINDOWS], bw, color=cols[s], label=s)
    axs[1].bar(xs + (j - 2) * bw, [MET[w][s]["mean"] for w in WINDOWS], bw, color=cols[s])
for ax, ttl, yl in [(axs[0], "(a) Lower-tail mean (worst 5% of weeks)", (0.15, 0.40)), (axs[1], "(b) Mean capacity factor", (0.35, 0.66))]:
    ax.set_xticks(xs); ax.set_xticklabels([f"{w} yr" for w in WINDOWS]); ax.set_xlabel("Training record"); ax.set_title(ttl, loc="left"); ax.set_ylim(*yl)
axs[0].legend(frameon=False, fontsize=6.0, loc="upper left", ncol=3); axs[0].set_ylabel("Out-of-sample CF, 2001–2026")
fig.tight_layout(); fig.savefig("fig/fig3_oos.png"); plt.show()

fig, axs = plt.subplots(1, 3, figsize=(7.5, 2.9), sharey=True)
for ax, s, lab in zip(axs, ["Minimum variance", "SAA-CVaR", "DRO-CVaR"], ["(a)", "(b)", "(c)"]):
    sh = pd.DataFrame({t: sea_share(WH[10][s][t]) for t in TEST_YEARS}).T
    ax.stackplot(sh.index, *[sh[c].values for c in SEA_ORDER], colors=[sea_col[c] for c in SEA_ORDER], labels=SEA_ORDER, alpha=0.9)
    ax.set_title(f"{lab} {s}", loc="left", fontsize=8.5); ax.set_xlabel("Test year"); ax.set_ylim(0, 1)
axs[0].set_ylabel("Capacity share"); axs[2].legend(frameon=False, fontsize=5.5, loc="center left", bbox_to_anchor=(1.0, 0.5))
fig.tight_layout(); fig.savefig("fig/fig4_weights.png"); plt.show()

fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.0))
ev = list(FIXP); axs[0].plot(range(len(ev)), [FIXP[e]["lt05"] for e in ev], "o-", color="#d62728", label="Lower-tail mean")
ax2 = axs[0].twinx(); ax2.plot(range(len(ev)), [FIXP[e]["mean"] for e in ev], "s--", color="#1f77b4", label="Mean"); ax2.spines["right"].set_visible(True)
axs[0].set_xticks(range(len(ev))); axs[0].set_xticklabels(ev, fontsize=7); axs[0].set_xlabel("Wasserstein radius ε (10-yr record)")
axs[0].set_ylabel("Lower-tail mean", color="#d62728"); ax2.set_ylabel("Mean", color="#1f77b4"); axs[0].set_title("(a) Out-of-sample performance", loc="left")
axs[1].plot(range(len(ev)), [FIXP[e]["max_share"] for e in ev], "o-", color="k"); axs[1].set_xticks(range(len(ev))); axs[1].set_xticklabels(ev, fontsize=7)
axs[1].set_xlabel("Wasserstein radius ε"); axs[1].set_ylabel("Largest site share (mean)"); axs[1].set_title("(b) Robustness acts as a soft cap", loc="left")
fig.tight_layout(); fig.savefig("fig/fig5_radius.png"); plt.show()
log("done")
print(json.dumps({k: R[k] for k in ["data", "corr", "low_output_days", "full_record_saa_sea", "oos", "eps_selected_share_positive", "tests", "fixed_radius_10y", "turkish_share", "sensitivity"]}, 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/rp07_outputs", "zip", ".", ".")
files.download("/content/rp07_outputs.zip")