# -*- coding: utf-8 -*-
"""Zero-shot foundation models for forecast-based decision support: Accuracy, calibration and decision costs on port weather and CO₂ emissions to 2026

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

# # Zero-shot foundation models for forecast-based decision support: Accuracy, calibration and decision costs on port weather and CO₂ emissions to 2026
#
# 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: CPU'da yaklaşık 2 saat (en uzun adım TimesFM 2.5, yaklaşık 1,5 saat; makaledeki süreler 2 iş parçacıklı CPU içindir). Model ağırlıkları Hugging Face'ten indirilir. Yalnızca değerlendirme ve şekiller için: replikasyon paketindeki `out/forecasts.parquet` dosyası çalışma klasörüne `out/forecasts.parquet` olarak konursa bölüm 2 atlanıp doğrudan bölüm 3 çalıştırılabilir.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: NASA POWER aylık liman verisi ve Carbon Monitor). 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 statsforecast==2.1.1 chronos-forecasting==2.3.2 timesfm==3.0.2 pyarrow

VERI_KAYNAGI = "api"   # "api": download from the original sources; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp09", exist_ok=True); os.chdir("/content/rp09")
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 the inputs of the forecasting paper into data/: monthly port weather (NASA POWER, 60 ports, 4 variables, 1981 - Sep 2026)
    and Carbon Monitor daily CO2 emissions by country and sector, summed to calendar months (Jan 2019 - Aug 2026)."""
    import io
    import pandas as pd
    from common import get, power_monthly, DATA
    from sites import PORTS
    rows = []
    for name, lat, lon in PORTS:
        d = power_monthly(lat, lon, ["WS10M_MAX", "T2M", "PRECTOTCORR", "ALLSKY_SFC_SW_DWN"])
        for p, ser in d.items():
            rows += [(p, f"{k[:4]}-{k[4:]}", v, name, lat, lon) for k, v in ser.items()]
        print(name)
    pd.DataFrame(rows, columns=["param", "period", "value", "port", "lat", "lon"]).to_csv(DATA / "power_ports_monthly.csv.gz", index=False, compression="gzip")
    cm = pd.read_csv(io.StringIO(get("https://datas.carbonmonitor.org/API/downloadFullDataset.php?source=carbon_global").text), usecols=[0, 1, 2, 3])
    cm["ds"] = pd.to_datetime(cm.date, format="%d/%m/%Y").values.astype("datetime64[M]")
    cm.groupby(["country", "sector", "ds"]).value.sum().reset_index().to_csv(DATA / "carbon_monitor_monthly.csv", index=False)
print(sorted(os.listdir('data')))

# ## 2. İstatistiksel modeller ve zaman serisi temel modelleriyle kayan başlangıçlı olasılıksal tahminler

# Paper 9 (stage 1): rolling-origin probabilistic forecasts from statistical models and time series foundation models. Series: monthly port weather (NASA POWER, 60 ports x 4 variables, 1981-Sep 2026) and monthly national/sectoral CO2 (Carbon Monitor, Jan 2019-Aug 2026).
import os, time, json, warnings
import numpy as np, pandas as pd
warnings.filterwarnings("ignore")
os.environ["NIXTLA_ID_AS_COL"] = "1"
import torch
torch.set_num_threads(2)
from statsforecast import StatsForecast
from statsforecast.models import SeasonalNaive, AutoETS, AutoTheta, AutoARIMA
Q = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
H = 12
t_all = time.time()

# ---------------- data ----------------
w = pd.read_csv("data/power_ports_monthly.csv.gz")
w["unique_id"] = "W|" + w.param + "|" + w.port
w["ds"] = pd.to_datetime(w.period)
W = w[["unique_id", "ds", "value"]].rename(columns={"value": "y"})
W = W[W.ds <= "2026-09-01"]
# Carbon Monitor daily emissions (Mt CO2 per day) summed to calendar months by country and sector (see fetch_data.py)
C = pd.read_csv("data/carbon_monitor_monthly.csv", parse_dates=["ds"], float_precision="round_trip")
C = C[C.ds <= "2026-08-01"]
C["unique_id"] = "C|" + C.sector + "|" + C.country
C = C.rename(columns={"value": "y"})[["unique_id", "ds", "y"]]
ok = C.groupby("unique_id").y.apply(lambda s: (s <= 0).mean() < 0.02 and len(s) == 92 and s.std() > 0)
C = C[C.unique_id.isin(ok[ok].index)]
DATA = pd.concat([W, C]).sort_values(["unique_id", "ds"]).reset_index(drop=True)
DATA.to_csv("out/series.csv", index=False)
print("series:", DATA.unique_id.nunique(), "weather", W.unique_id.nunique(), "co2", C.unique_id.nunique(), flush=True)

ORIGINS = {"W": pd.date_range("2023-09-01", "2025-09-01", freq="3MS"), "C": pd.date_range("2024-08-01", "2025-08-01", freq="3MS")}
CTX_STAT = 360          # statistical models: last 30 years
CTX_FM = 512            # foundation models: last 512 months

def origin_sets():
    for g, origins in ORIGINS.items():
        sub = DATA[DATA.unique_id.str.startswith(g + "|")]
        for o in origins:
            tr = sub[sub.ds <= o]
            te = sub[(sub.ds > o) & (sub.ds <= o + pd.DateOffset(months=H))]
            yield g, o, tr, te

def interval_quantiles(fc, name):
    out = {}
    out[0.5] = fc[name].values
    for lv, (lo, hi) in {80: (0.1, 0.9), 60: (0.2, 0.8), 40: (0.3, 0.7), 20: (0.4, 0.6)}.items():
        out[lo] = fc[f"{name}-lo-{lv}"].values; out[hi] = fc[f"{name}-hi-{lv}"].values
    return np.column_stack([out[q] for q in Q])

rows, timing = [], {}
def add(model, g, o, uid_order, ds_order, qmat):
    rows.append(pd.DataFrame({"model": model, "group": g, "origin": o, "unique_id": uid_order, "ds": ds_order,
                              **{f"q{int(q * 100)}": qmat[:, k] for k, q in enumerate(Q)}}))

# ---------------- statistical models ----------------
for g, o, tr, te in origin_sets():
    trs = tr.groupby("unique_id").tail(CTX_STAT)
    models = [SeasonalNaive(season_length=12), AutoETS(season_length=12), AutoTheta(season_length=12)]
    if g == "C": models.append(AutoARIMA(season_length=12))
    t = time.time()
    sf = StatsForecast(models=models, freq="MS", n_jobs=-1)
    fc = sf.forecast(df=trs, h=H, level=[20, 40, 60, 80]).reset_index(drop=True)
    timing.setdefault("statistical", 0.0); timing["statistical"] += time.time() - t
    for m in ["SeasonalNaive", "AutoETS", "AutoTheta"] + (["AutoARIMA"] if g == "C" else []):
        add(m, g, o, fc.unique_id.values, fc.ds.values, interval_quantiles(fc, m))
    # climatology: empirical calendar-month quantiles of the last 30 years (weather) or of all data (CO2)
    t = time.time()
    trc = trs.copy(); trc["m"] = trc.ds.dt.month
    qs = trc.groupby(["unique_id", "m"]).y.quantile(Q).unstack()
    fd = fc[["unique_id", "ds"]].copy(); fd["m"] = pd.to_datetime(fd.ds).dt.month
    qm = qs.loc[list(zip(fd.unique_id, fd.m))].values
    add("Climatology", g, o, fd.unique_id.values, fd.ds.values, qm)
    timing.setdefault("climatology", 0.0); timing["climatology"] += time.time() - t
    print("stat", g, o.date(), round(time.time() - t_all), "s", flush=True)

# ---------------- foundation models ----------------
from chronos import BaseChronosPipeline
import timesfm

def contexts(tr):
    ids, ctx = [], []
    for uid, s in tr.groupby("unique_id"):
        ids.append(uid); ctx.append(s.y.values[-CTX_FM:].astype("float32"))
    return ids, ctx

def fdates(o):
    return pd.date_range(o + pd.DateOffset(months=1), periods=H, freq="MS")

FMS = {}
FMS["Chronos-Bolt-Small"] = ("chronos", "amazon/chronos-bolt-small")
FMS["Chronos-Bolt-Base"] = ("chronos", "amazon/chronos-bolt-base")
FMS["Chronos-2"] = ("chronos2", "amazon/chronos-2")
FMS["TimesFM-2.5"] = ("timesfm", "google/timesfm-2.5-200m-pytorch")
for name, (kind, path) in FMS.items():
    t0 = time.time()
    if kind == "timesfm":
        model = timesfm.TimesFM_2p5_200M_torch.from_pretrained(path)
        model.compile(timesfm.ForecastConfig(max_context=1024, max_horizon=256, normalize_inputs=True, use_continuous_quantile_head=True,
                                             force_flip_invariance=True, infer_is_positive=False, fix_quantile_crossing=True))
    else:
        model = BaseChronosPipeline.from_pretrained(path, device_map="cpu", torch_dtype=torch.float32)
    load = time.time() - t0; tinf = 0.0
    for g, o, tr, te in origin_sets():
        ids, ctx = contexts(tr); t = time.time()
        if kind == "timesfm":
            qq = []
            for i in range(0, len(ctx), 64):
                _, q = model.forecast(horizon=H, inputs=ctx[i:i + 64]); qq.append(q[:, :, 1:10])
            qarr = np.concatenate(qq)
        else:
            qq = []
            for i in range(0, len(ctx), 64):
                batch = [torch.tensor(c) for c in ctx[i:i + 64]]
                q, _ = model.predict_quantiles(batch, prediction_length=H, quantile_levels=Q)
                q = torch.stack(q) if isinstance(q, list) else q
                q = q.squeeze(1) if q.ndim == 4 else q
                qq.append(q.float().numpy())
            qarr = np.concatenate(qq)
        tinf += time.time() - t
        qarr = np.sort(qarr, axis=2)          # enforce non-crossing quantiles
        uid = np.repeat(ids, H); ds = np.tile(fdates(o), len(ids))
        add(name, g, o, uid, ds, qarr.reshape(-1, len(Q)))
    timing[name] = tinf; timing[name + "_load"] = load
    print(name, "done", round(tinf), "s", flush=True)
    del model

F = pd.concat(rows, ignore_index=True)
F["ds"] = pd.to_datetime(F.ds)
F = F.merge(DATA, on=["unique_id", "ds"], how="left")
F.to_parquet("out/forecasts.parquet", index=False)
json.dump(timing, open("out/timing.json", "w"), indent=1)
print("forecasts:", F.shape, "missing actuals:", int(F.y.isna().sum()), "total", round(time.time() - t_all), "s")

# ## 3. Doğruluk, sıralar, testler, uyumlu (conformal) kalibrasyon, karar maliyeti ve şekiller

# Paper 9 (stage 2): accuracy, calibration, conformal recalibration, decision costs, tests and figures.
import json, warnings
import numpy as np, pandas as pd
from scipy import stats
import matplotlib
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
plt.rcParams.update({"font.family": "DejaVu Serif", "font.size": 9, "axes.titlesize": 10, "legend.fontsize": 8,
                     "savefig.dpi": 300, "axes.spines.top": False, "axes.spines.right": False})
Q = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]; QC = [f"q{int(q * 100)}" for q in Q]
F = pd.read_parquet("out/forecasts.parquet"); D = pd.read_csv("out/series.csv", parse_dates=["ds"])
timing = json.load(open("out/timing.json"))
R = {}

# ---------------- data-quality screen (uses only the observed series, before any evaluation) ----------------
# mean standardised anomaly of Jan 2023 - Sep 2026 against the 1991-2020 calendar-month climatology; |z| > 5 is physically implausible
Dw = D[D.unique_id.str.startswith("W|")].copy(); Dw["m"] = Dw.ds.dt.month
base = Dw[(Dw.ds >= "1991-01-01") & (Dw.ds <= "2020-12-01")].groupby(["unique_id", "m"]).y.agg(["mean", "std"])
Z = Dw.join(base, on=["unique_id", "m"]); Z["z"] = (Z.y - Z["mean"]) / Z["std"]
zr = Z[Z.ds >= "2023-01-01"].groupby("unique_id").z.mean()
BAD = sorted(zr[zr.abs() > 5].index)
R["excluded_series"] = {u: round(float(zr[u]), 2) for u in BAD}
R["screen_z_quantiles"] = {k: round(float(v), 2) for k, v in zr.abs().quantile([0.5, 0.9, 0.99]).items()}
F = F[~F.unique_id.isin(BAD)].copy(); D = D[~D.unique_id.isin(BAD)].copy()
print("excluded by the data-quality screen:", R["excluded_series"])
VAR = {"W|WS10M_MAX": "Maximum wind speed", "W|T2M": "Air temperature", "W|PRECTOTCORR": "Precipitation", "W|ALLSKY_SFC_SW_DWN": "Solar irradiance"}
def family(uid):
    p = uid.split("|")
    return VAR.get(p[0] + "|" + p[1], "CO2 emissions") if p[0] == "W" else "CO2 emissions"
F["family"] = F.unique_id.map(family)
F["h"] = ((F.ds.dt.year - pd.to_datetime(F.origin).dt.year) * 12 + F.ds.dt.month - pd.to_datetime(F.origin).dt.month).astype(int)
# scale: in-sample seasonal-naive MAE on the training window of each origin
sc = []
for (uid, o), _ in F.groupby(["unique_id", "origin"]):
    s = D[(D.unique_id == uid) & (D.ds <= o)].y.values[-360:]
    sc.append((uid, o, np.mean(np.abs(s[12:] - s[:-12]))))
SC = pd.DataFrame(sc, columns=["unique_id", "origin", "scale"]); F = F.merge(SC, on=["unique_id", "origin"])
F = F[F.scale > 0]
def pinball(y, q, tau): return np.maximum(tau * (y - q), (tau - 1) * (y - q))
F["ae"] = np.abs(F.y - F.q50) / F.scale
F["sql"] = np.mean([pinball(F.y.values, F[c].values, q) for c, q in zip(QC, Q)], axis=0) * 2 / F.scale     # scaled CRPS approximation
F["cov80"] = ((F.y >= F.q10) & (F.y <= F.q90)).astype(float)
F["width80"] = (F.q90 - F.q10) / F.scale
F["is80"] = (F.width80 + 2 / 0.2 * np.maximum(F.q10 - F.y, 0) / F.scale + 2 / 0.2 * np.maximum(F.y - F.q90, 0) / F.scale)   # interval score
MODELS = ["SeasonalNaive", "Climatology", "AutoETS", "AutoTheta", "AutoARIMA", "Chronos-Bolt-Small", "Chronos-Bolt-Base", "Chronos-2", "TimesFM-2.5"]
MODELS = [m for m in MODELS if m in F.model.unique()]
FM = ["Chronos-Bolt-Small", "Chronos-Bolt-Base", "Chronos-2", "TimesFM-2.5"]
R["n_series"] = {k: int(v) for k, v in F.groupby("family").unique_id.nunique().items()}
R["n_forecasts"] = int(len(F[F.model == "AutoETS"]))
def summary(df):
    return df.groupby("model").agg(MASE=("ae", "mean"), SQL=("sql", "mean"), cov80=("cov80", "mean"), IS80=("is80", "mean")).round(4)
R["overall_weather"] = summary(F[F.family != "CO2 emissions"]).to_dict("index")
R["overall_co2"] = summary(F[F.family == "CO2 emissions"]).to_dict("index")
R["by_family"] = {fam: summary(g).to_dict("index") for fam, g in F.groupby("family")}
# ranks per series (mean SQL over origins and horizons), Friedman + Nemenyi
def ranks(df, models):
    p = df[df.model.isin(models)].groupby(["unique_id", "model"]).sql.mean().unstack()[models].dropna()
    rk = p.rank(axis=1); fr = stats.friedmanchisquare(*[p[m] for m in models])
    k, n = len(models), len(p); qa = {5: 2.728, 6: 2.850, 7: 2.949, 8: 3.031, 9: 3.102}[k]
    return {"mean_rank": rk.mean().round(3).to_dict(), "friedman_chi2": round(float(fr.statistic), 1), "p": float(fr.pvalue), "N": int(n),
            "CD": round(qa * np.sqrt(k * (k + 1) / (6 * n)), 3)}, p
MW = [m for m in MODELS if m != "AutoARIMA"]
R["ranks_weather"], PW = ranks(F[F.family != "CO2 emissions"], MW)
R["ranks_co2"], PC = ranks(F[F.family == "CO2 emissions"], MODELS)
# robustness: only targets that postdate the public release of all four foundation models (Chronos-2, October 2025)
POST = F[F.ds >= "2025-11-01"]
R["post_release_n"] = {"weather": int(len(POST[(POST.family != "CO2 emissions") & (POST.model == "AutoETS")])), "co2": int(len(POST[(POST.family == "CO2 emissions") & (POST.model == "AutoETS")]))}
R["post_release_weather"] = summary(POST[POST.family != "CO2 emissions"]).to_dict("index")
R["post_release_co2"] = summary(POST[POST.family == "CO2 emissions"]).to_dict("index")
R["ranks_weather_post"], _ = ranks(POST[POST.family != "CO2 emissions"], MW)
R["ranks_co2_post"], _ = ranks(POST[POST.family == "CO2 emissions"], MODELS)
# paired comparisons on series-level SQL
def wil(p, a, b):
    d = p[a] - p[b]; w = stats.wilcoxon(p[a], p[b])
    return {"median_ratio": round(float((p[a] / p[b]).median()), 3), "share_a_better": round(float((d < 0).mean()), 3), "p": float(w.pvalue)}
best_stat_w = min(["Climatology", "AutoETS", "AutoTheta", "SeasonalNaive"], key=lambda m: PW[m].mean())
best_stat_c = min(["Climatology", "AutoETS", "AutoTheta", "SeasonalNaive", "AutoARIMA"], key=lambda m: PC[m].mean())
R["best_stat"] = {"weather": best_stat_w, "co2": best_stat_c}
R["wilcoxon_weather"] = {f"{m}_vs_{best_stat_w}": wil(PW, m, best_stat_w) for m in FM}
R["wilcoxon_co2"] = {f"{m}_vs_{best_stat_c}": wil(PC, m, best_stat_c) for m in FM}
R["wilcoxon_weather_by_family"] = {}
for fam in VAR.values():
    p = F[F.family == fam].groupby(["unique_id", "model"]).sql.mean().unstack()
    bs = min(["Climatology", "AutoETS", "AutoTheta", "SeasonalNaive"], key=lambda m: p[m].mean())
    bf = min(FM, key=lambda m: p[m].mean())
    R["wilcoxon_weather_by_family"][fam] = {"best_stat": bs, "best_fm": bf, **wil(p, bf, bs)}
# by horizon
R["sql_by_h"] = {fam: {m: [round(float(v), 4) for v in g[g.model == m].groupby("h").sql.mean().values] for m in MODELS if m in g.model.unique()}
                 for fam, g in F.assign(fam2=np.where(F.family == "CO2 emissions", "CO2", "Weather")).groupby("fam2")}

# ---------------- conformal recalibration of the 80% interval (CQR, calibrated on the earlier half of origins) ----------------
cq = []
for (g, m, fam, h), d in F.groupby(["group", "model", "family", "h"]):
    os_ = sorted(d.origin.unique()); cal = d[d.origin.isin(os_[:len(os_) // 2])]; test = d[d.origin.isin(os_[len(os_) // 2:])]
    s = np.maximum((cal.q10 - cal.y) / cal.scale, (cal.y - cal.q90) / cal.scale)
    n = len(s); c = np.quantile(s, min(1.0, np.ceil((n + 1) * 0.8) / n)) if n > 0 else 0.0
    lo = test.q10 - c * test.scale; hi = test.q90 + c * test.scale
    cov = ((test.y >= lo) & (test.y <= hi)).mean(); cov0 = test.cov80.mean()
    isc = ((hi - lo) / test.scale + 10 * np.maximum(lo - test.y, 0) / test.scale + 10 * np.maximum(test.y - hi, 0) / test.scale).mean()
    cq.append((g, m, fam, h, len(test), cov0, cov, test.is80.mean(), isc))
CQ = pd.DataFrame(cq, columns=["group", "model", "family", "h", "n", "cov_raw", "cov_cqr", "is_raw", "is_cqr"])
CQ.to_csv("out/conformal.csv", index=False)
agg = CQ.assign(w=CQ.n).groupby("model").apply(lambda d: pd.Series({"cov_raw": np.average(d.cov_raw, weights=d.w), "cov_cqr": np.average(d.cov_cqr, weights=d.w),
                                                                     "is_raw": np.average(d.is_raw, weights=d.w), "is_cqr": np.average(d.is_cqr, weights=d.w)}))
R["conformal"] = agg.round(4).to_dict("index")
R["conformal_by_group"] = {g: {m: {k: round(float(np.average(d[k], weights=d.n)), 4) for k in ["cov_raw", "cov_cqr", "is_raw", "is_cqr"]}
                               for m, d in gg.groupby("model")} for g, gg in CQ.groupby("group")}

# ---------------- decision costs: newsvendor-type asymmetric costs at critical ratio tau ----------------
TAUS = {0.1: "q10", 0.2: "q20", 0.5: "q50", 0.8: "q80", 0.9: "q90"}
dc = {}
for tau, c in TAUS.items():
    F[f"cost_{tau}"] = pinball(F.y.values, F[c].values, tau) / F.scale
    for fam2, g in F.assign(fam2=np.where(F.family == "CO2 emissions", "CO2", "Weather")).groupby("fam2"):
        base = g[g.model == "Climatology"][f"cost_{tau}"].mean()
        dc.setdefault(fam2, {})[str(tau)] = {m: round(float(g[g.model == m][f"cost_{tau}"].mean() / base), 4) for m in MODELS if m in g.model.unique()}
R["decision_cost_rel_climatology"] = dc

# ---------------- where do foundation models gain? forecastability (spectral entropy) ----------------
def spectral_entropy(x):
    x = (x - x.mean()) / (x.std() + 1e-12); p = np.abs(np.fft.rfft(x)) ** 2; p = p[1:] / p[1:].sum()
    return float(-(p * np.log(p + 1e-12)).sum() / np.log(len(p)))
se = D.groupby("unique_id").y.apply(lambda s: spectral_entropy(s.values[-360:]))
gain = (PW[best_stat_w] - PW[[m for m in FM]].min(axis=1)) / PW[best_stat_w]
G = pd.DataFrame({"entropy": se.reindex(gain.index), "gain": gain}).dropna()
rho = stats.spearmanr(G.entropy, G.gain); R["gain_vs_entropy_weather"] = {"rho": round(float(rho.correlation), 3), "p": float(rho.pvalue), "n": int(len(G))}
R["timing_seconds"] = {k: round(v, 1) for k, v in timing.items()}
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# ================= FIGURES =================
COL = {"SeasonalNaive": "#9e9e9e", "Climatology": "#8c564b", "AutoETS": "#1f77b4", "AutoTheta": "#17becf", "AutoARIMA": "#bcbd22",
       "Chronos-Bolt-Small": "#ff9896", "Chronos-Bolt-Base": "#d62728", "Chronos-2": "#9467bd", "TimesFM-2.5": "#2ca02c"}
LAB = {"SeasonalNaive": "Seasonal naive", "Climatology": "Climatology", "AutoETS": "ETS", "AutoTheta": "Theta", "AutoARIMA": "ARIMA",
       "Chronos-Bolt-Small": "Chronos-Bolt (small)", "Chronos-Bolt-Base": "Chronos-Bolt (base)", "Chronos-2": "Chronos-2", "TimesFM-2.5": "TimesFM 2.5"}
# Fig 1: examples from the last origin
ex = [("W|WS10M_MAX|Ambarli", "Maximum wind speed, Ambarlı (m/s)"), ("W|T2M|Piraeus", "Air temperature, Piraeus (°C)"),
      ("W|PRECTOTCORR|Nhava Sheva", "Precipitation, Nhava Sheva (mm/day)"), ("C|Power|China", "CO$_2$ from power, China (Mt per month)")]
fig, axs = plt.subplots(2, 2, figsize=(9, 5.6))
for a, (uid, ttl) in zip(axs.flat, ex):
    s = D[D.unique_id == uid]; last = F[(F.unique_id == uid)].origin.max()
    a.plot(s.ds[s.ds >= pd.Timestamp(last) - pd.DateOffset(months=48)], s.y[s.ds >= pd.Timestamp(last) - pd.DateOffset(months=48)], color="k", lw=1)
    for m in ["Climatology", "TimesFM-2.5"]:
        f = F[(F.unique_id == uid) & (F.origin == last) & (F.model == m)].sort_values("ds")
        a.plot(f.ds, f.q50, color=COL[m], lw=1.3, label=LAB[m]); a.fill_between(f.ds, f.q10, f.q90, color=COL[m], alpha=0.18, lw=0)
    a.axvline(pd.Timestamp(last), color="grey", ls=":", lw=0.8); a.set_title(ttl)
axs[0, 0].legend(frameon=False, fontsize=7); fig.autofmt_xdate(); fig.tight_layout(); fig.savefig("fig/fig1_examples.png"); plt.show()
# Fig 2: SQL by horizon
fig, axs = plt.subplots(1, 2, figsize=(9, 3.9))
for a, fam2 in zip(axs, ["Weather", "CO2"]):
    for m, v in R["sql_by_h"][fam2].items():
        a.plot(range(1, len(v) + 1), v, marker="o", ms=2.5, color=COL[m], lw=1.5 if m in FM else 0.9, label=LAB[m])
    a.set_title(f"({'ab'[['Weather', 'CO2'].index(fam2)]}) {'Port weather' if fam2 == 'Weather' else 'CO$_2$ emissions'}: scaled quantile loss"); a.set_xlabel("Horizon (months)")
h_, l_ = axs[1].get_legend_handles_labels(); fig.legend(h_, l_, loc="lower center", ncol=5, frameon=False, fontsize=7.5)
fig.tight_layout(rect=(0, 0.13, 1, 1)); fig.savefig("fig/fig2_horizon.png"); plt.show()
# Fig 3: mean ranks with critical difference
fig, axs = plt.subplots(1, 2, figsize=(9, 3.6))
for a, (key, ttl) in zip(axs, [("ranks_weather", "(a) Port weather"), ("ranks_co2", "(b) CO$_2$ emissions")]):
    r = R[key]; mr = pd.Series(r["mean_rank"]).sort_values()
    a.errorbar(mr.values, range(len(mr)), xerr=r["CD"] / 2, fmt="o", color="k", capsize=3)
    a.axvspan(mr.iloc[0] - r["CD"] / 2, mr.iloc[0] + r["CD"] / 2, color="#c6dbef", alpha=0.5)
    a.set_yticks(range(len(mr))); a.set_yticklabels([LAB[m] for m in mr.index]); a.invert_yaxis(); a.set_xlabel("Mean rank of scaled quantile loss"); a.set_title(ttl)
fig.tight_layout(); fig.savefig("fig/fig3_ranks.png"); plt.show()
# Fig 4: coverage raw vs conformal
fig, ax = plt.subplots(figsize=(6.5, 3.4))
c = pd.DataFrame(R["conformal"]).T.loc[[m for m in MODELS if m in R["conformal"]]]
x = np.arange(len(c)); ax.bar(x - 0.2, c.cov_raw * 100, 0.4, color="#9ecae1", label="Raw 80% interval"); ax.bar(x + 0.2, c.cov_cqr * 100, 0.4, color="#08519c", label="Conformalised")
ax.axhline(80, color="k", ls="--", lw=0.8); ax.set_xticks(x); ax.set_xticklabels([LAB[m] for m in c.index], rotation=30, ha="right", fontsize=7)
ax.set_ylabel("Empirical coverage (%)"); ax.set_ylim(0, 100); ax.legend(frameon=False, ncol=2, loc="lower center", bbox_to_anchor=(0.5, 1.0)); fig.tight_layout(); fig.savefig("fig/fig4_conformal.png"); plt.show()
# Fig 5: decision cost relative to climatology by critical ratio
fig, axs = plt.subplots(1, 2, figsize=(9, 3.9), sharey=False)
for a, fam2 in zip(axs, ["Weather", "CO2"]):
    d = pd.DataFrame(R["decision_cost_rel_climatology"][fam2])
    for m in d.index:
        a.plot([float(t) for t in d.columns], d.loc[m].values, marker="o", ms=3, color=COL[m], lw=1.5 if m in FM else 0.9, label=LAB[m])
    a.axhline(1, color="k", lw=0.6); a.set_xlabel("Critical ratio τ = c$_u$/(c$_u$+c$_o$)"); a.set_ylabel("Expected cost relative to climatology")
    a.set_title("(a) Port weather" if fam2 == "Weather" else "(b) CO$_2$ emissions")
h_, l_ = axs[1].get_legend_handles_labels(); fig.legend(h_, l_, loc="lower center", ncol=5, frameon=False, fontsize=7.5)
fig.tight_layout(rect=(0, 0.13, 1, 1)); fig.savefig("fig/fig5_decision.png"); plt.show()
print(json.dumps({k: R[k] for k in ["n_series", "n_forecasts", "overall_weather", "overall_co2", "ranks_weather", "ranks_co2", "best_stat", "wilcoxon_weather", "wilcoxon_co2",
                                    "wilcoxon_weather_by_family", "conformal", "decision_cost_rel_climatology", "gain_vs_entropy_weather", "timing_seconds",
                                    "post_release_n", "post_release_weather", "post_release_co2", "ranks_weather_post", "ranks_co2_post"]}, 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/rp09_outputs", "zip", ".", ".")
files.download("/content/rp09_outputs.zip")