# -*- coding: utf-8 -*-
"""How much inventory does chokepoint risk require? Regime-switching lead times and AIS-informed base-stock policies, 2019–2026

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

# # How much inventory does chokepoint risk require? Regime-switching lead times and AIS-informed base-stock policies, 2019–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: 5–15 dakika.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: IMF PortWatch günlük boğaz geçişleri). 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 hmmlearn searoute

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

# Downloads the inputs from the original sources (only when VERI_KAYNAGI == 'api').
if VERI_KAYNAGI == 'api':
    """Downloads the inputs of the chokepoint-inventory paper into data/: IMF PortWatch daily transit calls and capacity by vessel type
    at 28 maritime chokepoints (AIS-based), 1 Jan 2019 - 1 Nov 2026, and the chokepoint metadata."""
    import pandas as pd
    from common import get, DATA
    PW = "https://services9.arcgis.com/weJ1QsnbMYJlCHdG/arcgis/rest/services"

    def arcgis(layer, where, fields, order="ObjectId"):
        rows, off = [], 0
        while True:
            r = get(f"{PW}/{layer}/FeatureServer/0/query", params={"where": where, "outFields": fields, "orderByFields": order,
                    "resultOffset": off, "resultRecordCount": 1000, "f": "json"}).json()
            rows += [f["attributes"] for f in r["features"]]
            if not r.get("exceededTransferLimit"): break
            off += 1000
        return pd.DataFrame(rows)

    if __name__ == "__main__":
        arcgis("PortWatch_chokepoints_database", "1=1", "*").to_csv(DATA / "chokepoints.csv", index=False)
        d = arcgis("Daily_Chokepoints_Data", "1=1", "date,portid,portname,n_container,n_dry_bulk,n_general_cargo,n_roro,n_tanker,n_total,"
                   "capacity_container,capacity_dry_bulk,capacity_general_cargo,capacity_roro,capacity_tanker,capacity", order="date,portid")
        d.to_csv(DATA / "chokepoint_daily.csv.gz", index=False, compression="gzip")
        print(d.shape, d.date.min(), d.date.max())
print(sorted(os.listdir('data')))

# ## 2. Değişim noktası tespiti, gizli Markov modeli, rejime bağlı tedarik süreleri, baz stok politikaları, simülasyon ve şekiller

# Paper 13: chokepoint regimes, lead-time risk and inventory policy. Weekly AIS transits at maritime chokepoints (IMF PortWatch, 2019 - Nov 2026), hidden Markov and Bayesian online change-point models, sea-route lead times and a simulation of base-stock policies for Asia-Europe and Asia-US East Coast supply chains.
import json, sys, 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})
QUICK = "--quick" in sys.argv
R = {}
rng = np.random.default_rng(13)

# ---------------- weekly transits ----------------
d = pd.read_csv("data/chokepoint_daily.csv.gz", parse_dates=["date"])
d = d[d.date <= "2026-11-01"]
W = (d.set_index("date").groupby("portname")[["n_container", "n_tanker", "n_total", "capacity_container"]]
       .resample("W-SUN").sum().reset_index())
W = W[(W.date >= "2019-01-06") & (W.date <= "2026-11-01")]
C = W.pivot(index="date", columns="portname", values="n_container")
T = W.pivot(index="date", columns="portname", values="n_total")
TK = W.pivot(index="date", columns="portname", values="n_tanker")
R["weeks"] = int(len(C)); R["first_week"] = str(C.index.min().date()); R["last_week"] = str(C.index.max().date())
base = (C.index >= "2019-01-01") & (C.index <= "2023-10-31")
# Suez share of Asia-Europe container transits (Suez vs Cape of Good Hope)
share = C["Suez Canal"] / (C["Suez Canal"] + C["Cape of Good Hope"])
pi_suez = (share / share[base].mean()).clip(0, 1)                     # probability that a shipment takes the Suez route
pi_panama = (C["Panama Canal"] / C["Panama Canal"][base].mean()).clip(0, 1)
R["suez_share_base"] = round(float(share[base].mean()), 3)
R["suez_share_by_period"] = {k: round(float(share[(share.index >= a) & (share.index <= b)].mean()), 3) for k, (a, b) in
    {"2019-2022": ("2019-01-01", "2022-12-31"), "2024": ("2024-01-01", "2024-12-31"), "2025": ("2025-01-01", "2025-12-31"),
     "2026 Jan-Jun": ("2026-01-01", "2026-06-30"), "2026 Jul-Oct": ("2026-07-01", "2026-11-01")}.items()}

# ---------------- regime summary for six chokepoints (ratio to 2019-2023 weekly mean) ----------------
series = {"Suez Canal (container)": C["Suez Canal"], "Bab el-Mandeb (container)": C["Bab el-Mandeb Strait"],
          "Cape of Good Hope (container)": C["Cape of Good Hope"], "Panama Canal (all)": T["Panama Canal"],
          "Strait of Hormuz (tanker)": TK["Strait of Hormuz"], "Bosporus (all)": T["Bosporus Strait"]}
ratio = pd.DataFrame({k: v / v[base].mean() for k, v in series.items()})
ratio.to_csv("out/weekly_ratio_to_baseline.csv")
ctab = {}
for k, v in series.items():
    r4 = (v / v[base].mean()).rolling(4).mean()
    ctab[k] = {"baseline_weekly_mean": round(float(v[base].mean()), 1), "min_4wk_ratio": round(float(r4.min()), 2),
               "min_4wk_week": str(r4.idxmin().date()), "ratio_2024": round(float((v / v[base].mean())[(v.index >= "2024-01-01") & (v.index <= "2024-12-31")].mean()), 2),
               "ratio_2025": round(float((v / v[base].mean())[(v.index >= "2025-01-01") & (v.index <= "2025-12-31")].mean()), 2),
               "ratio_last_8_weeks": round(float((v / v[base].mean()).iloc[-8:].mean()), 2)}
R["chokepoint_table"] = ctab

# ---------------- Bayesian online change-point detection (Adams & MacKay, Normal-Gamma model) ----------------
def bocpd(x, hazard=1 / 104, mu0=0.0, kappa0=1.0, alpha0=1.0, beta0=1.0):
    T_ = len(x); P = np.zeros((T_ + 1, T_ + 1)); P[0, 0] = 1.0
    mu, kap, al, be = np.array([mu0]), np.array([kappa0]), np.array([alpha0]), np.array([beta0])
    cp_prob = np.zeros(T_); map_rl = np.zeros(T_, int)
    for t, xt in enumerate(x):
        scale = np.sqrt(be * (kap + 1) / (al * kap))
        pred = stats.t.pdf(xt, 2 * al, loc=mu, scale=scale)
        growth = P[t, :t + 1] * pred * (1 - hazard)
        cp = (P[t, :t + 1] * pred * hazard).sum()
        P[t + 1, 1:t + 2] = growth; P[t + 1, 0] = cp
        P[t + 1, :t + 2] /= P[t + 1, :t + 2].sum()
        cp_prob[t] = P[t + 1, :5].sum()                   # probability that the current run started within the last 4 weeks
        map_rl[t] = int(np.argmax(P[t + 1, :t + 2]))
        mu_n = (kap * mu + xt) / (kap + 1); kap_n = kap + 1; al_n = al + 0.5; be_n = be + kap * (xt - mu) ** 2 / (2 * (kap + 1))
        mu = np.concatenate([[mu0], mu_n]); kap = np.concatenate([[kappa0], kap_n]); al = np.concatenate([[alpha0], al_n]); be = np.concatenate([[beta0], be_n])
    return cp_prob, map_rl

def standardise(s):
    return (s - s[base].mean()) / s[base].std()
det = {}
CP = {}
EVENTS = {"Suez Canal (container)": "2023-12-15", "Panama Canal (all)": "2023-08-01", "Strait of Hormuz (tanker)": "2026-02-28", "Bosporus (all)": "2026-08-01"}
for k in series:
    z = standardise(np.log1p(series[k])).values
    cpp, rl = bocpd(z)
    CP[k] = pd.Series(cpp, index=C.index)
    # alarm: change-point probability above 0.5 within the last 4 weeks and a 4-week mean at least 2 baseline sd below the baseline
    lvl = pd.Series(z, index=C.index).rolling(4).mean()
    alarm = CP[k][(CP[k].rolling(4).max() > 0.5) & (lvl < -2)]
    first = alarm[alarm.index >= "2023-01-01"]
    det[k] = str(first.index.min().date()) if len(first) else None
R["bocpd_first_alarm_since_2023"] = det
pd.DataFrame(CP).to_csv("out/bocpd_changepoint_probability.csv")

# ---------------- hidden Markov model for the Suez routing share ----------------
from hmmlearn.hmm import GaussianHMM
y = np.log(share / (1 - share)).values.reshape(-1, 1)
best = None
for seed in range(10):
    m = GaussianHMM(n_components=3, covariance_type="diag", n_iter=500, random_state=seed).fit(y)
    ll = m.score(y)
    if best is None or ll > best[0]: best = (ll, m)
hmm = best[1]
order = np.argsort(hmm.means_.ravel())[::-1]                       # state 0 = normal (highest share), 2 = most disrupted
post = hmm.predict_proba(y)[:, order]
states = np.array([np.where(order == s)[0][0] for s in hmm.predict(y)])
A = hmm.transmat_[np.ix_(order, order)]
R["hmm"] = {"state_share_mean": [round(float(1 / (1 + np.exp(-hmm.means_.ravel()[o]))), 3) for o in order],
            "transition_matrix": np.round(A, 4).tolist(),
            "expected_duration_weeks": [round(float(1 / (1 - A[i, i])), 1) for i in range(3)],
            "weeks_in_state": [int((states == i).sum()) for i in range(3)], "loglik": round(float(best[0]), 1)}
# observed spells
spells = []; cur = states[0]; start = C.index[0]
for t_, s_ in zip(C.index[1:], states[1:]):
    if s_ != cur:
        spells.append((int(cur), str(start.date()), str((t_ - pd.Timedelta(days=7)).date()))); cur = s_; start = t_
spells.append((int(cur), str(start.date()), str(C.index[-1].date())))
R["hmm_spells"] = spells
# real-time filtered probabilities (expanding-window refits every 4 weeks; forward filtering with data to t only)
filt = np.full((len(y), 3), np.nan); pdis = np.full(len(y), np.nan); means_rt = np.full((len(y), 3), np.nan)
step = 8 if QUICK else 4
m_rt = None
for t in range(104, len(y)):
    if (t - 104) % step == 0:
        cand = None
        for seed in range(4):
            mm = GaussianHMM(n_components=3, covariance_type="diag", n_iter=300, random_state=seed)
            try:
                mm.fit(y[:t + 1]); sc = mm.score(y[:t + 1])
                if cand is None or sc > cand[0]: cand = (sc, mm)
            except Exception:
                pass
        m_rt = cand[1]; o_rt = np.argsort(m_rt.means_.ravel())[::-1]
    fp = m_rt.predict_proba(y[:t + 1])[-1]                          # forward-backward on data to t: the last row equals the filtered probability
    filt[t] = fp[o_rt]; means_rt[t] = 1 / (1 + np.exp(-m_rt.means_.ravel()[o_rt]))
    pdis[t] = fp[1 / (1 + np.exp(-m_rt.means_.ravel())) < 0.6].sum()      # states whose mean Suez share is below 0.6
filt_arr = filt.copy()
filt = pd.DataFrame(filt, index=C.index, columns=["normal", "partial", "diverted"]); filt["p_disrupted"] = pdis
filt.to_csv("out/hmm_filtered.csv")
onset = filt[(filt.index >= "2023-11-01") & (filt.p_disrupted > 0.5)].index.min()
R["hmm_realtime_alarm_weeks_before_nov2023"] = [str(x.date()) for x in filt.index[(filt.p_disrupted > 0.5) & (filt.index < "2023-11-01")]]
last_div = filt[(filt.index >= "2026-01-01") & (filt.p_disrupted > 0.5)].index.max()
R["hmm_realtime_last_week_disrupted"] = str(last_div.date()) if pd.notna(last_div) else None
R["hmm_realtime_first_week_disrupted"] = str(onset.date()) if pd.notna(onset) else None

# ---------------- sea-route lead times ----------------
import searoute as sr
B = ["northwest", "northeast"]
def nm(o, d_, block):
    return float(sr.searoute(o, d_, units="naut", restrictions=B + block)["properties"]["length"])
SHA, RTM, PIR, NYC = [121.64, 31.19], [4.21, 51.92], [23.61, 37.94], [-74.05, 40.67]
dist = {"Asia-North Europe": {"suez": nm(SHA, RTM, ["panama"]), "cape": nm(SHA, RTM, ["suez", "babalmandab", "panama"])},
        "Asia-Mediterranean": {"suez": nm(SHA, PIR, ["panama"]), "cape": nm(SHA, PIR, ["suez", "babalmandab", "panama"])},
        "Asia-US East Coast": {"panama": nm(SHA, NYC, ["suez", "babalmandab"]), "suez": nm(SHA, NYC, ["panama"]),
                               "cape": nm(SHA, NYC, ["panama", "suez", "babalmandab"])}}
R["distance_nm"] = {k: {r: round(v) for r, v in x.items()} for k, x in dist.items()}
SPEED = 16.0; FIXED = 7.0; SD_L = 2.0          # knots; days at origin and destination; schedule noise (days)
def sail_days(nmiles, speed=None): return nmiles / ((speed or SPEED) * 24)
R["lead_time_days"] = {k: {r: round(sail_days(v) + FIXED, 1) for r, v in x.items()} for k, x in dist.items()}

def route_probs(lane, t):
    if lane == "Asia-US East Coast":
        pp = pi_panama.iloc[t]; ps = pi_suez.iloc[t]
        return {"panama": pp, "suez": (1 - pp) * ps, "cape": (1 - pp) * (1 - ps)}
    ps = pi_suez.iloc[t]
    return {"suez": ps, "cape": 1 - ps}

# ---------------- inventory model: weekly review, daily demand, backorders, base-stock with stochastic lead times ----------------
MU, SIG = 100.0, 30.0           # weekly demand mean and standard deviation (units)
H, BCOST = 1.0, 19.0            # holding and backorder cost per unit-week (critical ratio 0.95)
NREP = 60 if QUICK else 1000
TW = len(C); TD = TW * 7
GH_x, GH_w = np.polynomial.hermite_e.hermegauss(15); GH_w = GH_w / GH_w.sum()
from scipy.optimize import brentq
_cache = {}
def base_stock(lane, probs):
    """Order-up-to level S solving P(D[L + R] <= S) = b/(b+h): daily demand N(mu/7, sig^2/7), R = 7 days, L a mixture of normal route times."""
    key = (lane, tuple(round(probs[r], 3) for r in probs), SPEED, FIXED, SD_L, SIG, BCOST)
    if key in _cache: return _cache[key]
    md, vd = MU / 7, SIG ** 2 / 7; cr = BCOST / (H + BCOST)
    comps = [(probs[r], sail_days(dist[lane][r]) + FIXED) for r in probs if probs[r] > 1e-6]
    def F(S_):
        tot = 0.0
        for w_, m_ in comps:
            l = np.maximum(m_ + SD_L * GH_x, 1.0) + 7
            tot += w_ * np.sum(GH_w * stats.norm.cdf((S_ - md * l) / np.sqrt(vd * l)))
        return tot / sum(w_ for w_, _ in comps) - cr
    val = brentq(F, 0, MU * 20); _cache[key] = val
    return val

def lt_mean(lane, probs):
    return sum(probs[r] * (sail_days(dist[lane][r]) + FIXED) for r in probs) / sum(probs.values())

def policy_levels(lane):
    pre_probs = {r: np.mean([route_probs(lane, t)[r] for t in range(TW) if base[t]]) for r in route_probs(lane, 0)}
    S_static = base_stock(lane, pre_probs)
    S = {"Static (2019-2023 calibration)": np.full(TW, S_static)}
    # reactive: lead-time distribution of shipments that have arrived in the last 13 weeks (information arrives with the goods)
    arr_t = np.array([t + lt_mean(lane, route_probs(lane, t)) / 7 for t in range(TW)])
    S_reac = np.zeros(TW)
    for t in range(TW):
        seen = np.where((arr_t <= t) & (arr_t > t - 13))[0]
        if len(seen):
            pr = {r: np.mean([route_probs(lane, u)[r] for u in seen]) for r in route_probs(lane, 0)}
            S_reac[t] = base_stock(lane, pr)
        else:
            S_reac[t] = S_static
    S["Reactive (realised lead times)"] = S_reac
    # AIS-informed: route mix of current departures, observed in real time
    S["AIS-informed (current route mix)"] = np.array([base_stock(lane, route_probs(lane, t)) for t in range(TW)])
    # regime-filtered: Suez routing probability implied by the real-time HMM state probabilities and state means (static before the first estimate)
    S_reg = np.full(TW, S_static)
    for t in range(TW):
        if np.isnan(filt_arr[t]).any(): continue
        pis = float(np.sum(filt_arr[t] * np.minimum(1, means_rt[t] / share[base].mean())))
        if lane == "Asia-US East Coast":
            pp = pi_panama.iloc[t]; pr = {"panama": pp, "suez": (1 - pp) * pis, "cape": (1 - pp) * (1 - pis)}
        else:
            pr = {"suez": pis, "cape": 1 - pis}
        S_reg[t] = base_stock(lane, {k: round(v, 3) for k, v in pr.items()})
    S["Regime-filtered (HMM probabilities)"] = S_reg
    return S

def route_matrix(lane):
    rts = list(route_probs(lane, 0))
    P = np.array([[route_probs(lane, t)[r] for r in rts] for t in range(TW)]); P = P / P.sum(1, keepdims=True)
    days = np.array([sail_days(dist[lane][r]) + FIXED for r in rts])
    return P, days

def simulate(lane, S):
    """Daily simulation, all replications vectorised; common random numbers across policies."""
    P, days = route_matrix(lane); cum = P.cumsum(1)
    rs = np.random.default_rng(2026)
    dem = np.maximum(rs.normal(MU / 7, SIG / np.sqrt(7), (NREP, TD)), 0)
    idx = (rs.random((NREP, TW))[:, :, None] > cum[None]).sum(2)
    lead = np.maximum(np.rint(days[idx] + rs.normal(0, SD_L, (NREP, TW))), 1).astype(int)     # days, order placed at week t
    res = {}
    rows = np.arange(NREP)
    for name, lev in S.items():
        net = np.full(NREP, lev[0]); pos = net.copy(); arr = np.zeros((NREP, TD + 120))
        served = np.zeros(NREP); cost_d = np.zeros((NREP, TD)); oh = np.zeros(NREP)
        for t in range(TD):
            net += arr[:, t]
            avail = np.maximum(net, 0); served += np.minimum(dem[:, t], avail)
            net -= dem[:, t]; pos -= dem[:, t]
            cost_d[:, t] = (H * np.maximum(net, 0) + BCOST * np.maximum(-net, 0)) / 7
            oh += np.maximum(net, 0)
            if t % 7 == 0:
                w = t // 7; q = np.maximum(lev[w] - pos, 0); pos += q
                np.add.at(arr, (rows, t + lead[:, w]), q)
        per = cost_d.reshape(NREP, TW, 7).sum(2)
        res[name] = {"cost": per.mean(1), "fill": served / dem.sum(1), "onhand": oh / TD, "per": per.mean(0), "per_rep": per}
    return res

PER = {"2019-2023": ("2019-03-01", "2023-11-30"), "Red Sea onset (Dec 2023 - Mar 2024)": ("2023-12-01", "2024-03-31"),
       "Disrupted (Apr 2024 - Jun 2026)": ("2024-04-01", "2026-06-30"), "Partial return (Jul - Oct 2026)": ("2026-07-01", "2026-11-01")}
SIM = {}
LANES = ["Asia-North Europe", "Asia-Mediterranean", "Asia-US East Coast"]
for lane in LANES:
    S = policy_levels(lane)
    res = simulate(lane, S)
    # oracle: AIS-informed level shifted by the realised mean lead time of the orders in the pipeline is approximated by perfect foresight of the route mix 4 weeks ahead
    SIM[lane] = {}
    for name, r in res.items():
        SIM[lane][name] = {"mean_cost": round(float(r["cost"].mean()), 1), "fill_rate": round(float(r["fill"].mean()), 4),
                           "mean_on_hand": round(float(r["onhand"].mean()), 1),
                           "cost_by_period": {p: round(float(r["per"][(C.index >= a) & (C.index <= b)].mean()), 1) for p, (a, b) in PER.items()},
                           "cost_se_by_period": {p: round(float(r["per_rep"][:, (C.index >= a) & (C.index <= b)].mean(1).std(ddof=1) / np.sqrt(NREP)), 2) for p, (a, b) in PER.items()}}
    SIM[lane]["base_stock_levels"] = {name: {p: round(float(lev[(C.index >= a) & (C.index <= b)].mean()), 1) for p, (a, b) in PER.items()} for name, lev in S.items()}
    pd.DataFrame({**{f"S_{k}": v for k, v in S.items()}, **{f"cost_{k}": res[k]["per"] for k in res}}, index=C.index).to_csv(f"out/sim_{lane.replace(' ', '_')}.csv")
R["simulation"] = SIM
R["simulation_settings"] = {"weekly_demand_mean": MU, "weekly_demand_sd": SIG, "holding_cost": H, "backorder_cost": BCOST, "replications": NREP,
                            "speed_knots": SPEED, "fixed_days": FIXED, "schedule_sd_days": SD_L}
# pipeline (in-transit) inventory: extra days of goods at sea per unit of annual imports
DIS = (C.index >= "2024-04-01") & (C.index <= "2026-06-30")
pl = {}
for lane in dist:
    m_pre = np.mean([lt_mean(lane, route_probs(lane, t)) for t in range(TW) if base[t]])
    m_dis = np.mean([lt_mean(lane, route_probs(lane, t)) for t in range(TW) if DIS[t]])
    pl[lane] = {"mean_lead_time_pre_days": round(m_pre, 1), "mean_lead_time_disrupted_days": round(m_dis, 1),
                "extra_in_transit_share_of_annual_flow_pct": round((m_dis - m_pre) / 365 * 100, 2)}
R["pipeline"] = pl
# decomposition of the change in the AIS-informed base-stock level (disrupted vs 2019-2023): mean shift vs added dispersion
dec = {}
for lane in ["Asia-North Europe", "Asia-Mediterranean", "Asia-US East Coast"]:
    rts = list(route_probs(lane, 0))
    pre = {r: np.mean([route_probs(lane, t)[r] for t in range(TW) if base[t]]) for r in rts}
    dis = {r: np.mean([route_probs(lane, t)[r] for t in range(TW) if DIS[t]]) for r in rts}
    S0 = base_stock(lane, pre); S1 = base_stock(lane, dis)
    # single-route equivalents with the same mean lead time isolate the effect of the mean shift from that of route dispersion
    def single(m):
        nm_ = f"_single_{m:.4f}"; dist[nm_] = {"x": (m - FIXED) * SPEED * 24}; v = base_stock(nm_, {"x": 1.0}); del dist[nm_]; return v
    M0, M1 = single(lt_mean(lane, pre)), single(lt_mean(lane, dis))
    dec[lane] = {"S_pre": round(S0, 1), "S_disrupted": round(S1, 1), "increase_pct": round(100 * (S1 / S0 - 1), 1),
                 "from_mean_lead_time": round(M1 - M0, 1), "from_route_dispersion": round((S1 - M1) - (S0 - M0), 1)}
R["base_stock_decomposition"] = dec
# base-stock level as a function of the Suez routing probability (two-route lanes)
curve = {}
for lane in ["Asia-North Europe", "Asia-Mediterranean"]:
    grid = np.round(np.linspace(0, 1, 21), 2)
    Sg = [base_stock(lane, {"suez": float(g), "cape": float(1 - g)}) for g in grid]
    curve[lane] = {"pi": grid.tolist(), "S": [round(x, 1) for x in Sg], "argmax_pi": float(grid[int(np.argmax(Sg))]),
                   "S_full_diversion": round(Sg[0], 1), "S_max": round(max(Sg), 1), "S_normal": round(Sg[-1], 1)}
R["base_stock_curve"] = curve
# sensitivity of the static policy's excess cost in the disrupted period (Asia-North Europe), relative to the AIS-informed policy
sens = []
grid = [("baseline", {}), ("speed 14 knots", {"SPEED": 14.0}), ("speed 18 knots", {"SPEED": 18.0}), ("fixed time 4 days", {"FIXED": 4.0}),
        ("fixed time 10 days", {"FIXED": 10.0}), ("demand CV 0.1", {"SIG": 10.0}), ("demand CV 0.5", {"SIG": 50.0}),
        ("critical ratio 0.90", {"BCOST": 9.0}), ("critical ratio 0.98", {"BCOST": 49.0}), ("schedule sd 4 days", {"SD_L": 4.0})]
keep = dict(SPEED=SPEED, FIXED=FIXED, SIG=SIG, BCOST=BCOST, SD_L=SD_L, NREP=NREP)
for name, ch in grid:
    globals().update(keep); globals().update(ch); globals()["NREP"] = 40 if QUICK else 400
    _cache.clear(); S_ = policy_levels("Asia-North Europe"); r_ = simulate("Asia-North Europe", S_)
    msk = DIS
    cs = {k: float(v["per"][msk].mean()) for k, v in r_.items()}
    msk_on = (C.index >= "2023-12-01") & (C.index <= "2024-03-31")
    sens.append({"case": name, "static_cost": round(cs["Static (2019-2023 calibration)"], 1), "reactive_cost": round(cs["Reactive (realised lead times)"], 1),
                 "ais_cost": round(cs["AIS-informed (current route mix)"], 1), "regime_cost": round(cs["Regime-filtered (HMM probabilities)"], 1),
                 "onset_reactive_cost": round(float(r_["Reactive (realised lead times)"]["per"][msk_on].mean()), 1),
                 "onset_ais_cost": round(float(r_["AIS-informed (current route mix)"]["per"][msk_on].mean()), 1),
                 "static_excess_pct": round(100 * (cs["Static (2019-2023 calibration)"] / cs["AIS-informed (current route mix)"] - 1), 1),
                 "S_increase_pct": round(100 * (S_["AIS-informed (current route mix)"][msk].mean() / S_["Static (2019-2023 calibration)"][0] - 1), 1)})
globals().update(keep); _cache.clear()
R["sensitivity"] = sens
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# ================= FIGURES =================
# Fig 1: weekly transits relative to baseline for six chokepoints
fig, axs = plt.subplots(3, 2, figsize=(9, 6.6), sharex=True)
for a, (k, s) in zip(axs.flat, ratio.items()):
    a.plot(s.index, s.values, color="k", lw=0.8); a.axhline(1, color="grey", lw=0.6, ls=":")
    a.fill_between(s.index, 0, 1, where=(CP[k].values > 0.5) & (s.index >= "2019-03-01"), color="#fdae6b", alpha=0.6, lw=0, transform=a.get_xaxis_transform())
    a.set_title(k); a.set_ylim(0, max(2.0, float(s.max()) * 1.05))
for a in axs[:, 0]: a.set_ylabel("Ratio to 2019-2023 mean")
fig.autofmt_xdate(); fig.tight_layout(); fig.savefig("fig/fig1_chokepoints.png"); plt.show()
# Fig 2: Suez routing share, HMM states and real-time filtered probability
fig, axs = plt.subplots(2, 1, figsize=(9, 5.2), sharex=True, gridspec_kw={"height_ratios": [1.3, 1]})
axs[0].plot(share.index, share.values, color="k", lw=0.9, label="Suez share of Suez + Cape container transits")
cols = ["#2ca02c", "#ff7f0e", "#d62728"]
for i in range(3):
    axs[0].fill_between(share.index, 0, 1, where=states == i, color=cols[i], alpha=0.15, lw=0, transform=axs[0].get_xaxis_transform())
axs[0].set_ylim(0, 1); axs[0].set_ylabel("Share"); axs[0].legend(frameon=False, loc="lower left"); axs[0].set_title("(a) Routing share and smoothed HMM states (green: normal, orange: partial, red: diverted)")
bz = pd.Series(bocpd(standardise(pd.Series(y.ravel(), index=C.index)).values)[0], index=C.index)
axs[1].plot(filt.index, filt.p_disrupted, color="#d62728", lw=1.3, label="HMM: real-time probability of a diverted state")
axs[1].plot(bz.index[bz.index >= "2019-03-01"], bz[bz.index >= "2019-03-01"], color="#1f77b4", lw=0.9, alpha=0.8, label="BOCPD: probability of a change point in the last 4 weeks")
axs[1].set_ylim(0, 1.05); axs[1].set_ylabel("Probability"); axs[1].legend(frameon=False, loc="upper left"); axs[1].set_title("(b) Real-time detection (expanding-window HMM and Bayesian online change-point detection)")
fig.tight_layout(); fig.savefig("fig/fig2_hmm.png"); plt.show()
# Fig 3: base-stock levels by policy and weekly cost, Asia-North Europe
sim = pd.read_csv("out/sim_Asia-North_Europe.csv", index_col=0, parse_dates=True); sim = sim[sim.index >= "2019-03-01"]
fig, axs = plt.subplots(2, 1, figsize=(9, 5.4), sharex=True)
PC = {"Static (2019-2023 calibration)": "#7f7f7f", "Reactive (realised lead times)": "#1f77b4", "AIS-informed (current route mix)": "#d62728", "Regime-filtered (HMM probabilities)": "#2ca02c"}
for k, c_ in PC.items():
    axs[0].plot(sim.index, sim[f"S_{k}"], color=c_, lw=1.1, label=k)
    axs[1].plot(sim.index, sim[f"cost_{k}"].rolling(4).mean(), color=c_, lw=1.1, label=k)
axs[0].set_ylabel("Base-stock level (units)"); axs[0].set_title("(a) Order-up-to levels, Asia-North Europe"); axs[0].legend(frameon=False, fontsize=7)
axs[1].set_ylabel("Cost per week (4-week mean)"); axs[1].set_title("(b) Expected weekly holding and backorder cost")
fig.tight_layout(); fig.savefig("fig/fig3_policies.png"); plt.show()
# Fig 4: cost and fill rate by period and policy for both lanes
fig, axs = plt.subplots(1, 3, figsize=(10, 3.6))
for a, lane in zip(axs, LANES):
    x = np.arange(len(PER)); wdt = 0.2
    for j, (k, c_) in enumerate(PC.items()):
        a.bar(x + (j - 1.5) * wdt, [SIM[lane][k]["cost_by_period"][p] for p in PER], wdt, color=c_, label=k)
    a.set_xticks(x); a.set_xticklabels(["2019-23", "Onset", "Disrupted", "Return"], fontsize=8); a.set_title(lane); a.set_ylabel("Mean cost per week")
axs[0].legend(frameon=False, fontsize=7); fig.tight_layout(); fig.savefig("fig/fig4_costs.png"); plt.show()
fig, ax = plt.subplots(figsize=(6.2, 3.4))
for lane, c_ in [("Asia-North Europe", "#1f77b4"), ("Asia-Mediterranean", "#d62728")]:
    ax.plot(curve[lane]["pi"], curve[lane]["S"], color=c_, lw=1.5, label=lane)
    obs = [pi_suez[DIS].mean()]
    ax.axvline(float(pi_suez[DIS].mean()), color="grey", ls=":", lw=0.8)
ax.set_xlabel("Probability that a shipment takes the Suez route (π)"); ax.set_ylabel("Base-stock level (units)"); ax.legend(frameon=False)
ax.axvline(float(pi_suez.iloc[-8:].mean()), color="grey", ls="--", lw=0.8)
ax.text(float(pi_suez[DIS].mean()) - 0.01, ax.get_ylim()[0] + 5, "mean, Apr 2024 – Jun 2026", fontsize=7, color="grey", ha="right")
ax.text(float(pi_suez.iloc[-8:].mean()) + 0.01, ax.get_ylim()[0] + 5, "Sep – Oct 2026", fontsize=7, color="grey")
fig.tight_layout(); fig.savefig("fig/fig5_basestock_curve.png"); plt.show()
R["pi_suez_mean_disrupted"] = round(float(pi_suez[DIS].mean()), 3); R["pi_suez_last8"] = round(float(pi_suez.iloc[-8:].mean()), 3)
json.dump(R, open("out/results.json", "w"), indent=1, default=str)
print(json.dumps(R, indent=1, default=str)[:3000])

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

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