#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Doğu Akdeniz ve Türkiye Havzalarında Çok Değişkenli Okyanus-İklim Veri Madenciliği: Spektral PCA, K-Means++, Walk-Forward Ridge CV ve Shannon Transfer Entropisi (1990–2026)
Eastern Mediterranean & Türkiye Multivariate Ocean-Climate Data Mining: Spectral PCA, K-Means++, Walk-Forward Ridge CV & Shannon Transfer Entropy (1990–2026)
==============================================================================
GOMD Python SDK v2.0.0 — Official 100% Real-Data Reproducible Case Study
Google Colab Install: !pip install -q --force-reinstall --no-cache-dir https://gomd.pirireis.edu.tr/sdk/dist/gomd-2.0.0-py3-none-any.whl
"""

# ==============================================================================
# ADIM 0: GOMD Python SDK v2.0.0 Kurulumu ve Canlı Veri Doğrulama Yardımcıları
# STEP 0: Install GOMD Python SDK v2.0.0 & Live GOMD Data Verification Helpers
# ==============================================================================
import importlib
import os
import subprocess
import sys
import math
import json
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

if hasattr(sys.stdout, "reconfigure"):
    try:
        sys.stdout.reconfigure(encoding="utf-8", errors="replace")
    except Exception:
        pass

# Yerel geliştirme ortamında sdk/python varsa önceliklendir; Colab'da resmi wheel kur
_local_sdk = Path(__file__).resolve().parent.parent / "python" if "__file__" in globals() else None
if _local_sdk and (_local_sdk / "gomd" / "__init__.py").is_file():
    if str(_local_sdk) not in sys.path:
        sys.path.insert(0, str(_local_sdk))

_IN_COLAB = "google.colab" in sys.modules
_GOMD_WHEEL_URL = "https://gomd.pirireis.edu.tr/sdk/dist/gomd-2.0.0-py3-none-any.whl"
try:
    import gomd  # noqa: F401
    if _IN_COLAB:
        subprocess.check_call([
            sys.executable, "-m", "pip", "install", "-q",
            "--force-reinstall", "--no-cache-dir", _GOMD_WHEEL_URL
        ])
except ImportError:
    subprocess.check_call([
        sys.executable, "-m", "pip", "install", "-q",
        "--force-reinstall", "--no-cache-dir", _GOMD_WHEEL_URL
    ])

import gomd
for _mod_name in list(sys.modules.keys()):
    if _mod_name == "gomd" or _mod_name.startswith("gomd."):
        importlib.reload(sys.modules[_mod_name])
import gomd
from gomd import GOMDClient

client = GOMDClient("https://gomd.pirireis.edu.tr")


def fetch_gomd_annual_series(
    indicator: str,
    entity: str = "TUR",
    period_from: int = 1950,
    period_to: int = 2026,
) -> dict:
    """
    GOMD sunucusundan (https://gomd.pirireis.edu.tr) gerçek gözlem serisini çeker
    ve {yıl: float(değer)} sözlüğü olarak döndürür. Sentetik veri üretmez.
    """
    rows = client.get_series(
        indicator=indicator,
        entity=entity,
        period_from=period_from,
        period_to=period_to,
        limit=5000,
        as_df=False,
    )
    annual = {}
    for r in rows:
        p_str = str(r.get("period") or "").strip()
        val = r.get("value")
        if val is None or not p_str:
            continue
        # Yıllık kayıtları doğrudan al; aylık kayıt varsa yıl bazında ortalamasını/değerini tut
        yr_str = p_str[:4]
        if yr_str.isdigit():
            yr = int(yr_str)
            if period_from <= yr <= period_to:
                if len(p_str) == 4 or yr not in annual:
                    annual[yr] = float(val)
    if not annual:
        raise RuntimeError(f"GOMD API'den gerçek veri alınamadı: indicator={indicator}, entity={entity}")
    print(f"📡 [GERÇEK GOMD VERİSİ] {indicator:<44} | Entity={entity:<18} | N={len(annual):<3} yıl ({min(annual)}–{max(annual)})")
    return dict(sorted(annual.items()))


print(f"✅ GOMD Python SDK Sürümü / Version : v{gomd.__version__}")
print(f"🌐 Canlı Veri Sunucusu / Server     : {client.base_url}")
print(f"🏛️ Kurum / Institution              : Pîrî Reis Üniversitesi — GOMD Bilimsel Veri Konsorsiyumu (COP31 Türkiye)")
print(f"📚 Akademik Atıf / Citation          : {gomd.cite('apa')[:115]}...")

# ==============================================================================
# ADIM 1: Canlı GOMD Sunucusundan 8 Oşinografik Göstergenin Çekimi (1990–2026)
# STEP 1: Live Retrieval of 8 Oceanographic Series from GOMD Server (1990–2026)
# ==============================================================================
from gomd.ocean import (
    seawater_density_eos80,
    oxygen_solubility,
    carbonate_system,
    marine_heatwave_detect,
    compound_mhw_oax_event,
    sst_warming_velocity_profile,
)
from gomd.stats import (
    mann_kendall_sen,
    pettitt_changepoint,
    transfer_entropy_causal,
    multivariate_data_mining_pipeline,
)

sst_map   = fetch_gomd_annual_series("copernicus_sst_med", entity="TUR", period_from=1990, period_to=2026)
sal_map   = fetch_gomd_annual_series("teos10_med_salinity_psu", entity="TUR", period_from=1990, period_to=2026)
fco2_map  = fetch_gomd_annual_series("socat_fco2", entity="TUR", period_from=1990, period_to=2026)
ph_map    = fetch_gomd_annual_series("socat_surface_ph_total_scale", entity="TUR", period_from=1990, period_to=2026)
omega_map = fetch_gomd_annual_series("socat_aragonite_saturation_omega", entity="TUR", period_from=1990, period_to=2026)
ohc_map   = fetch_gomd_annual_series("cmems_med_ohc_0_700m_zj", entity="TUR", period_from=1990, period_to=2026)
mhw_map   = fetch_gomd_annual_series("cmems_med_mhw_cumulative_intensity", entity="TUR", period_from=1990, period_to=2026)
do_map    = fetch_gomd_annual_series("marmara_subhalocline_do_mg_l", entity="TUR", period_from=1990, period_to=2026)

common_years = sorted(
    set(sst_map) & set(sal_map) & set(fco2_map) & set(ph_map)
    & set(omega_map) & set(ohc_map) & set(mhw_map) & set(do_map)
)
baseline_sst_1990_2000 = float(np.mean([sst_map[y] for y in common_years if y <= 2000]))

records = []
for yr in common_years:
    sst_degc = sst_map[yr]
    sst_anom = round(sst_degc - baseline_sst_1990_2000, 4)
    sal_psu  = sal_map[yr]
    ph_total = ph_map[yr]
    do_mg_l  = do_map[yr]

    eos = seawater_density_eos80(salinity_psu=sal_psu, temp_c=sst_degc, pressure_dbar=15.0)
    # mg/L -> umol/kg dönüşümü (O2 molar kütle = 31.9988 g/mol)
    do_umol_kg = round((do_mg_l * 1000.0) / (31.9988 * (eos["rho_kg_m3"] / 1000.0)), 3)
    oxy  = oxygen_solubility(temp_c=sst_degc, salinity_psu=sal_psu, observed_do_umol_kg=do_umol_kg)
    carb = carbonate_system(ph_total=ph_total, temp_c=sst_degc, salinity_psu=sal_psu, ta_umol_kg=2610.0)

    records.append({
        "year": yr,
        "sst_degc": round(sst_degc, 4),
        "sst_anom_c": sst_anom,
        "sal_psu": round(sal_psu, 4),
        "sigma_t": eos["sigma_t_kg_m3"],
        "do_mg_l": round(do_mg_l, 4),
        "aou_umol_kg": oxy["aou_umol_kg"],
        "ph_total": round(ph_total, 4),
        "omega_arag_obs": round(omega_map[yr], 4),
        "omega_arag_calc": carb["omega_aragonite"],
        "fco2_uatm": round(fco2_map[yr], 2),
        "ohc_0_700m_zj": round(ohc_map[yr], 4),
        "mhw_cum_intensity": round(mhw_map[yr], 2),
    })

df = pd.DataFrame(records)
mhw = marine_heatwave_detect(df["sst_degc"].tolist(), dates_or_steps=df["year"].tolist(), pct_threshold=85.0, min_duration=1)
comp = compound_mhw_oax_event(
    sst_series=df["sst_degc"].tolist(),
    ph_series=df["ph_total"].tolist(),
    do_mg_l_series=df["do_mg_l"].tolist(),
    labels=df["year"].tolist(),
)

print(f"\n📊 Doğu Akdeniz & Türkiye Gerçek GOMD Okyanus-İklim Matrisi (N={len(df)} yıl, Son 8 Yıl):")
print(df.tail(8).to_string(index=False))

# ==============================================================================
# ADIM 2: GOMD Çok Değişkenli Veri Madenciliği Hattı (PCA + K-Means++ + Ridge CV)
# STEP 2: GOMD Multivariate Data Mining Pipeline (PCA + K-Means++ + Ridge CV)
# ==============================================================================
feature_cols = ["sst_degc", "sigma_t", "aou_umol_kg", "omega_arag_obs", "fco2_uatm", "ohc_0_700m_zj", "mhw_cum_intensity"]
feature_matrix = df[feature_cols].values.tolist()

mining = multivariate_data_mining_pipeline(
    feature_matrix=feature_matrix,
    feature_names=feature_cols,
    period_labels=df["year"].tolist(),
    target_index=0,
    k_clusters=3,
    ridge_alphas=(0.01, 0.1, 1.0, 5.0, 10.0),
)

pca_info = mining["pca"]
km_info  = mining["kmeans_regimes"]
ml_info  = mining["ridge_ml_forecast"]

df["PC1"] = [s["pc1"] for s in pca_info["scores"]]
df["PC2"] = [s["pc2"] for s in pca_info["scores"]]
df["regime_id"] = km_info["regime_labels"]
regime_names = {
    1: "Rejim 1: Tarihsel Baz Dönemi",
    2: "Rejim 2: Hızlanan Isınma & Tuzlulaşma",
    3: "Rejim 3: Bileşik MHW-Asitlenme Ekstremi",
}
df["regime_label"] = df["regime_id"].map(regime_names)

mk_sst = mann_kendall_sen(df["sst_degc"].tolist(), t=df["year"].tolist())
pt_sst = pettitt_changepoint(df["sst_degc"].tolist(), labels=df["year"].tolist())
te_sst_arag = transfer_entropy_causal(df["sst_degc"].tolist(), df["omega_arag_obs"].tolist(), bins=3)
vel_sst = sst_warming_velocity_profile(df["sst_degc"].tolist(), dt_days=df["year"].tolist())

print(f"🔬 PCA Varyans Açıklama Oranı       : PC1 = %{pca_info['explained_variance_ratio'][0]*100:.2f}, PC2 = %{pca_info['explained_variance_ratio'][1]*100:.2f} (Toplam = %{pca_info['cumulative_variance_pc1_pc2_pct']:.2f})")
print(f"🧩 K-Means++ Silhouette Skoru       : {km_info['silhouette_score']:.4f} (WCSS Inertia = {km_info['wcss_inertia']:.2f})")
print(f"🤖 Lag-1 Ridge CV Test RMSE         : {ml_info['test_rmse']:.4f} °C (Yetenek Skoru / Skill = +%{ml_info['forecast_skill_score_pct']:.2f}, α = {ml_info['optimal_alpha_l2']})")
print(f"📈 Mann-Kendall SST Trendi          : +{mk_sst['theil_sen_slope']*10:.3f} °C/on-yıl (p = {mk_sst['p_value']:.2e}, {mk_sst['trend_classification']})")
print(f"⚡ Pettitt Yapısal Kırılma Yılı     : {pt_sst['changepoint_label']} (Δ = +{pt_sst['regime_shift_delta']:.3f} °C, p = {pt_sst['p_value_approx']:.4f})")
print(f"🔗 Shannon Transfer Entropisi       : T(SST -> Ω_arag) = {te_sst_arag['transfer_entropy_x_to_y_bits']:.4f} bits ({te_sst_arag['dominant_causal_direction']})")

# ==============================================================================
# ADIM 3: 4 Panelli Yayın Kalitesinde Görselleştirme (Matplotlib / Google Colab)
# STEP 3: 4-Panel Publication-Grade Visualization in Google Colab
# ==============================================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 9.5), dpi=120)
fig.patch.set_facecolor("#071324")
for ax in axes.flat:
    ax.set_facecolor("#0b1d36")
    ax.tick_params(colors="#cbd5e1")
    for spine in ax.spines.values():
        spine.set_color("#1e3a5f")

# Panel A: Gerçek GOMD SST Serisi, Pettitt Kırılma Yılı ve Theil-Sen Trend Doğrusu
ax1 = axes[0, 0]
ax1.plot(df["year"], df["sst_degc"], color="#38bdf8", marker="o", lw=2, label="Copernicus Akdeniz SST (°C)")
trend_line = [mk_sst["theil_sen_intercept"] + mk_sst["theil_sen_slope"] * yr for yr in df["year"]]
ax1.plot(df["year"], trend_line, color="#f59e0b", ls="--", lw=2, label=f"Theil-Sen (+{mk_sst['theil_sen_slope']*10:.2f} °C/10y)")
ax1.axvline(int(pt_sst["changepoint_label"]), color="#ef4444", ls=":", lw=1.8, label=f"Pettitt Kırılma ({pt_sst['changepoint_label']})")
ax1.set_title("A) Gerçek GOMD Doğu Akdeniz SST Serisi ve Yapısal Kırılma (1990–2026)", color="#f8fafc", fontweight="bold")
ax1.set_ylabel("Deniz Yüzeyi Sıcaklığı (°C)", color="#cbd5e1")
ax1.legend(facecolor="#0f172a", edgecolor="#334155", labelcolor="#f8fafc", fontsize=8.5)
ax1.grid(alpha=0.18, color="#64748b")

# Panel B: Spektral PCA (PC1 vs PC2) & K-Means++ Okyanus-İklim Rejimleri
ax2 = axes[0, 1]
palette = {1: "#22c55e", 2: "#f59e0b", 3: "#ef4444"}
for rid, rname in regime_names.items():
    sub = df[df["regime_id"] == rid]
    ax2.scatter(sub["PC1"], sub["PC2"], s=70, c=palette[rid], label=rname, edgecolors="#ffffff", linewidth=0.6)
    for _, row in sub.iloc[::3].iterrows():
        ax2.annotate(str(int(row["year"])), (row["PC1"], row["PC2"]), xytext=(4, 4), textcoords="offset points", color="#e2e8f0", fontsize=7.5)
ax2.set_title(f"B) Spektral PCA & K-Means++ Kümeleme (Silhouette = {km_info['silhouette_score']:.3f})", color="#f8fafc", fontweight="bold")
ax2.set_xlabel(f"PC1 (%{pca_info['explained_variance_ratio'][0]*100:.1f} Varyans)", color="#cbd5e1")
ax2.set_ylabel(f"PC2 (%{pca_info['explained_variance_ratio'][1]*100:.1f} Varyans)", color="#cbd5e1")
ax2.legend(facecolor="#0f172a", edgecolor="#334155", labelcolor="#f8fafc", fontsize=8)
ax2.grid(alpha=0.18, color="#64748b")

# Panel C: SOCAT Karbonat Asitlenmesi (Omega Aragonit) vs. CMEMS Okyanus Isı İçeriği (OHC)
ax3 = axes[1, 0]
ax3.plot(df["year"], df["omega_arag_obs"], color="#14b8a6", lw=2.2, label="SOCAT Aragonit Doygunluğu (Ω_arag)")
ax3_twin = ax3.twinx()
ax3_twin.plot(df["year"], df["ohc_0_700m_zj"], color="#f43f5e", lw=2.0, ls="-.", label="CMEMS OHC 0–700m (ZJ)")
ax3_twin.tick_params(colors="#fda4af")
ax3.set_title("C) SOCAT Deniz Asitlenmesi (Ω_arag) ve CMEMS Okyanus Isı İçeriği", color="#f8fafc", fontweight="bold")
ax3.set_ylabel("Ω_arag", color="#14b8a6")
ax3_twin.set_ylabel("Okyanus Isı İçeriği 0–700m (ZJ)", color="#f43f5e")
ax3.grid(alpha=0.18, color="#64748b")

# Panel D: Walk-Forward CV Lag-1 Ridge Regresyonu Örneklem Dışı Tahmin Performansı
ax4 = axes[1, 1]
preds = pd.DataFrame(ml_info["test_predictions"])
preds["period_int"] = preds["period"].astype(int)
ax4.plot(preds["period_int"], preds["observed"], color="#38bdf8", marker="o", lw=2.2, label="Gözlenen SST (°C)")
ax4.plot(preds["period_int"], preds["ridge_predicted"], color="#22c55e", marker="s", ls="--", lw=2.0, label=f"Lag-1 Ridge Tahmini (RMSE={ml_info['test_rmse']:.3f})")
ax4.plot(preds["period_int"], preds["naive_baseline"], color="#94a3b8", marker="^", ls=":", lw=1.5, label=f"Naif Kalıcılık (RMSE={ml_info['naive_baseline_rmse']:.3f})")
ax4.set_title("D) Zaman Serisi CV Lag-1 Ridge Örneklem Dışı Tahmin (Test)", color="#f8fafc", fontweight="bold")
ax4.set_ylabel("SST (°C)", color="#cbd5e1")
ax4.legend(facecolor="#0f172a", edgecolor="#334155", labelcolor="#f8fafc", fontsize=8)
ax4.grid(alpha=0.18, color="#64748b")

fig.suptitle("GOMD Python SDK v2.0.0 — Vaka Çalışması 1: Doğu Akdeniz Gerçek Veri Okyanus-İklim Madenciliği", color="#ffffff", fontsize=13, fontweight="bold", y=0.99)
plt.tight_layout()
fig.savefig("cs01_multivariate_ocean_climate_data_mining.png", dpi=140, bbox_inches="tight", facecolor=fig.get_facecolor())
plt.show()