#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
===============================================================================
GOMD (v2.0.0) — End-to-End Multivariate Ocean, Climate & Maritime Data Mining
===============================================================================
Platform : Global Ocean & Maritime Data Platform (https://gomd.pirireis.edu.tr)
 Paper   : https://gomd.pirireis.edu.tr/sdk_paper.php
Author   : Piri Reis University Maritime Sustainability Research Centre (COP31)

Scientific & Data Mining Pipeline Stages:
  Stage 1 — Multi-Source Panel Acquisition & SHA-256 Provenance Verification
  Stage 2 — Oceanographic & Biogeochemical Feature Engineering (TEOS-10 GSW,
            CO2SYS Aragonite Saturation, Triple Compound MHW+OAX+LOX & SST Velocity)
  Stage 3 — Unsupervised Multivariate Data Mining: Z-Score Correlation Matrix,
            Spectral PCA (PC1/PC2 Loadings & Variance) & K-Means++ Regime
            Clustering with Rousseeuw (1987) Silhouette Validation
  Stage 4 — Non-Parametric Trend, Pettitt Regime Shift, Şen IPTA Tail Amplification
            & Runge (2023) Nonlinear Shannon Transfer Entropy Causal Discovery
  Stage 5 — Supervised Machine Learning: Expanding-Window Time-Series CV
            Multivariate Lag-1 Ridge Regression vs. Naive Persistence Baseline
  Stage 6 — Maritime Decarbonization (IMO CII, EU ETS, FuelEU Maritime),
            Actuarial Marine Insurance Pricing & Mavi Vatan (UN A/74/757) Audit
===============================================================================
"""

from __future__ import annotations

import importlib
import json
import math
import os
import sys
from pathlib import Path
from typing import Any, Dict, List, Optional, Sequence

# Allow running directly from repository root or installed `gomd` package
if hasattr(sys.stdout, "reconfigure"):
    try:
        sys.stdout.reconfigure(encoding="utf-8", errors="replace")
    except Exception:
        pass
REPO_PYTHON_SDK = Path(__file__).resolve().parent.parent / "python"
if REPO_PYTHON_SDK.exists() and str(REPO_PYTHON_SDK) not in sys.path:
    sys.path.insert(0, str(REPO_PYTHON_SDK))

import gomd


def _solve_linear_system_local(a_mat: List[List[float]], b_vec: List[float]) -> List[float]:
    p = len(b_vec)
    aug = [row[:] + [float(b_vec[i])] for i, row in enumerate(a_mat)]
    for col in range(p):
        pivot_row = max(range(col, p), key=lambda r: abs(aug[r][col]))
        aug[col], aug[pivot_row] = aug[pivot_row], aug[col]
        piv = aug[col][col]
        if abs(piv) < 1e-12:
            piv = 1e-12
        for j in range(col, p + 1):
            aug[col][j] /= piv
        for r in range(p):
            if r != col:
                factor = aug[r][col]
                for j in range(col, p + 1):
                    aug[r][j] -= factor * aug[col][j]
    return [aug[i][p] for i in range(p)]


def _fallback_multivariate_data_mining_pipeline(
    feature_matrix: Sequence[Sequence[float]],
    feature_names: Optional[Sequence[str]] = None,
    period_labels: Optional[Sequence[Any]] = None,
    target_index: int = 0,
    k_clusters: int = 3,
    ridge_alphas: Sequence[float] = (0.01, 0.1, 1.0, 10.0),
) -> Dict[str, Any]:
    """Self-healing fallback implementation if an older cached `gomd` wheel is loaded in Colab."""
    rows = [[float(v) for v in r] for r in feature_matrix if len(r) > 0]
    n = len(rows)
    p = len(rows[0])
    names = [str(feature_names[j]) for j in range(p)] if (feature_names and len(feature_names) >= p) else [f"V{j+1}" for j in range(p)]
    labels = [str(period_labels[i]) for i in range(n)] if (period_labels and len(period_labels) >= n) else [str(i + 1) for i in range(n)]
    tgt_idx = min(max(int(target_index), 0), p - 1)

    means = [sum(rows[i][j] for i in range(n)) / n for j in range(p)]
    stds = [max(math.sqrt(sum((rows[i][j] - means[j]) ** 2 for i in range(n)) / max(n - 1, 1)), 1e-9) for j in range(p)]
    z_mat = [[(rows[i][j] - means[j]) / stds[j] for j in range(p)] for i in range(n)]

    corr_mat = [[0.0] * p for _ in range(p)]
    for j1 in range(p):
        for j2 in range(p):
            c = 1.0 if j1 == j2 else sum(z_mat[i][j1] * z_mat[i][j2] for i in range(n)) / max(n - 1, 1)
            corr_mat[j1][j2] = round(max(min(c, 1.0), -1.0), 4)

    def _power_eig(mat: List[List[float]]) -> Any:
        vec = [1.0 / math.sqrt(p)] * p
        eig_val = 0.0
        for _ in range(80):
            nxt = [sum(mat[r][c] * vec[c] for c in range(p)) for r in range(p)]
            norm = math.sqrt(sum(v * v for v in nxt))
            if norm < 1e-12:
                break
            vec = [v / norm for v in nxt]
            eig_val = sum(vec[r] * sum(mat[r][c] * vec[c] for c in range(p)) for r in range(p))
        if vec[0] < 0:
            vec = [-v for v in vec]
        return max(eig_val, 0.0), vec

    cov_work = [[corr_mat[r][c] for c in range(p)] for r in range(p)]
    eig1, vec1 = _power_eig(cov_work)
    for r in range(p):
        for c in range(p):
            cov_work[r][c] -= eig1 * vec1[r] * vec1[c]
    eig2, vec2 = _power_eig(cov_work)
    evr1 = eig1 / float(p)
    evr2 = eig2 / float(p)
    pc_scores = [
        {"period": labels[i], "pc1": round(sum(z_mat[i][j] * vec1[j] for j in range(p)), 4), "pc2": round(sum(z_mat[i][j] * vec2[j] for j in range(p)), 4)}
        for i in range(n)
    ]

    k = min(max(int(k_clusters), 2), max(n // 2, 2))
    pc1_vals = [sc["pc1"] for sc in pc_scores]
    sorted_idx = sorted(range(n), key=lambda idx: pc1_vals[idx])
    centroids = [z_mat[sorted_idx[min(int(round(c_idx * (n - 1) / max(k - 1, 1))), n - 1)]][:] for c_idx in range(k)]

    assignments = [0] * n
    for _ in range(40):
        changed = False
        for i in range(n):
            best_c = min(range(k), key=lambda c_idx: sum((z_mat[i][j] - centroids[c_idx][j]) ** 2 for j in range(p)))
            if assignments[i] != best_c:
                assignments[i] = best_c
                changed = True
        for c_idx in range(k):
            members = [i for i in range(n) if assignments[i] == c_idx]
            if members:
                centroids[c_idx] = [sum(z_mat[m][j] for m in members) / len(members) for j in range(p)]
        if not changed:
            break

    cluster_order = sorted(range(k), key=lambda c_idx: centroids[c_idx][tgt_idx])
    remap = {old_c: new_c for new_c, old_c in enumerate(cluster_order)}
    assignments = [remap[c] for c in assignments]
    centroids = [centroids[old_c] for old_c in cluster_order]

    def _euclid(a: Sequence[float], b: Sequence[float]) -> float:
        return math.sqrt(sum((a[j] - b[j]) ** 2 for j in range(p)))

    sil_scores = []
    for i in range(n):
        ci = assignments[i]
        same = [m for m in range(n) if assignments[m] == ci and m != i]
        a_i = sum(_euclid(z_mat[i], z_mat[m]) for m in same) / len(same) if same else 0.0
        b_i = min(
            (sum(_euclid(z_mat[i], z_mat[m]) for m in range(n) if assignments[m] == c_other) / max(sum(1 for m in range(n) if assignments[m] == c_other), 1)
             for c_other in range(k) if c_other != ci),
            default=0.0,
        )
        sil_scores.append((b_i - a_i) / max(a_i, b_i, 1e-9))

    wcss = sum(sum((z_mat[i][j] - centroids[assignments[i]][j]) ** 2 for j in range(p)) for i in range(n))
    regime_summaries = []
    for c_idx in range(k):
        members = [i for i in range(n) if assignments[i] == c_idx]
        phys_centroid = {names[j]: round(centroids[c_idx][j] * stds[j] + means[j], 4) for j in range(p)}
        regime_summaries.append({
            "regime_id": c_idx + 1,
            "n_periods": len(members),
            "period_span": f"{labels[members[0]]}–{labels[members[-1]]}" if members else "N/A",
            "centroid_physical_units": phys_centroid,
        })

    x_lag = [z_mat[t][:] for t in range(n - 1)]
    y_next = [rows[t + 1][tgt_idx] for t in range(n - 1)]
    n_pairs = len(y_next)
    min_train = max(p + 2, n_pairs // 2)

    def _fit_ridge(x_tr: List[List[float]], y_tr: List[float], alpha: float) -> Any:
        ntr = len(y_tr)
        y_mean = sum(y_tr) / ntr
        y_c = [y_tr[i] - y_mean for i in range(ntr)]
        xtx = [[sum(x_tr[i][r] * x_tr[i][c] for i in range(ntr)) + (alpha if r == c else 0.0) for c in range(p)] for r in range(p)]
        xty = [sum(x_tr[i][r] * y_c[i] for i in range(ntr)) for r in range(p)]
        return y_mean, _solve_linear_system_local(xtx, xty)

    best_alpha = float(ridge_alphas[0])
    best_cv_rmse = float("inf")
    for alpha in ridge_alphas:
        sq_errs = []
        for split in range(min_train, n_pairs):
            b0, betas = _fit_ridge(x_lag[:split], y_next[:split], float(alpha))
            pred = b0 + sum(x_lag[split][j] * betas[j] for j in range(p))
            sq_errs.append((y_next[split] - pred) ** 2)
        rmse = math.sqrt(sum(sq_errs) / max(len(sq_errs), 1))
        if rmse < best_cv_rmse:
            best_cv_rmse = rmse
            best_alpha = float(alpha)

    oof_Sq, naive_Sq = [], []
    for split in range(min_train, n_pairs):
        b0, betas = _fit_ridge(x_lag[:split], y_next[:split], best_alpha)
        pred_phys = b0 + sum(x_lag[split][j] * betas[j] for j in range(p))
        act_phys = y_next[split]
        naive_phys = rows[split][tgt_idx]
        oof_Sq.append((act_phys - pred_phys) ** 2)
        naive_Sq.append((act_phys - naive_phys) ** 2)

    test_rmse = math.sqrt(sum(oof_Sq) / max(len(oof_Sq), 1))
    naive_rmse = math.sqrt(sum(naive_Sq) / max(len(naive_Sq), 1))
    skill_pct = (1.0 - (test_rmse / max(naive_rmse, 1e-9))) * 100.0
    _, full_betas = _fit_ridge(x_lag, y_next, best_alpha)

    return {
        "n_periods": n,
        "n_features": p,
        "feature_names": names,
        "target_variable": names[tgt_idx],
        "pca": {
            "eigenvalues": [round(eig1, 4), round(eig2, 4)],
            "explained_variance_ratio": [round(evr1, 4), round(evr2, 4)],
            "cumulative_variance_pc1_pc2_pct": round((evr1 + evr2) * 100.0, 2),
            "loadings_pc1": {names[j]: round(vec1[j], 4) for j in range(p)},
            "loadings_pc2": {names[j]: round(vec2[j], 4) for j in range(p)},
            "scores": pc_scores,
        },
        "kmeans_regimes": {
            "k": k,
            "wcss_inertia": round(wcss, 4),
            "silhouette_score": round(sum(sil_scores) / max(len(sil_scores), 1), 4),
            "regimes": regime_summaries,
        },
        "ridge_ml_forecast": {
            "optimal_alpha_l2": best_alpha,
            "cv_rmse": round(best_cv_rmse, 4),
            "test_rmse": round(test_rmse, 4),
            "naive_baseline_rmse": round(naive_rmse, 4),
            "forecast_skill_score_pct": round(skill_pct, 2),
            "standardized_coefficients": {names[j]: round(full_betas[j], 4) for j in range(p)},
        },
    }


if not hasattr(gomd, "multivariate_data_mining_pipeline"):
    setattr(gomd, "multivariate_data_mining_pipeline", _fallback_multivariate_data_mining_pipeline)


def build_multivariate_ocean_panel(start_year: int = 1995, end_year: int = 2025) -> Dict[str, Any]:
    """
    Construct a 31-year (1995–2025) multivariate annual panel combining 6 core
    ocean-climate, biogeochemical, maritime decarbonization, and blue finance variables
    fetched 100% live from the GOMD server (`https://gomd.pirireis.edu.tr`).
    """
    client = gomd.GOMDClient("https://gomd.pirireis.edu.tr")

    def _fetch_annual(indicator: str, entity: str) -> Dict[int, float]:
        rows = client.get_series(
            indicator=indicator,
            entity=entity,
            period_from=start_year,
            period_to=end_year,
            limit=5000,
            as_df=False,
        )
        out: Dict[int, float] = {}
        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
            yr_str = p_str[:4]
            if yr_str.isdigit():
                yr = int(yr_str)
                if start_year <= yr <= end_year:
                    if len(p_str) == 4 or yr not in out:
                        out[yr] = float(val)
        if not out:
            raise RuntimeError(f"Failed to fetch live GOMD series: indicator={indicator}, entity={entity}")
        return out

    sst_map = _fetch_annual("copernicus_sst_med", "TUR")
    co2_map = _fetch_annual("gml_mauna_loa_co2_ppm", "WLD")
    ph_map  = _fetch_annual("socat_surface_ph_total_scale", "TUR")
    do_map  = _fetch_annual("marmara_subhalocline_do_mg_l", "TUR")
    aer_map = _fetch_annual("imo_cii_fleet_aer", "TUR")
    ins_map = _fetch_annual("iumi_marine_insurance_total_premium_usd_mn", "WLD")

    years = sorted(set(sst_map) & set(co2_map) & set(ph_map) & set(do_map) & set(aer_map) & set(ins_map))
    n = len(years)

    sst_med_c = [round(sst_map[y], 4) for y in years]
    co2_ppm = [round(co2_map[y], 2) for y in years]
    ph_total = [round(ph_map[y], 4) for y in years]
    do_mg_l = [round(do_map[y], 4) for y in years]
    fleet_aer = [round(aer_map[y], 4) for y in years]
    marine_ins_usd_mn = [round(ins_map[y], 2) for y in years]

    # Derive Aragonite Saturation State (Omega_arag) via gomd.carbonate_system for each real observation year
    omega_arag = [
        gomd.carbonate_system(
            ph_total=ph_total[i],
            temp_c=sst_med_c[i],
            salinity_psu=38.6,
            ta_umol_kg=2610.0,
        )["omega_aragonite"]
        for i in range(n)
    ]

    feature_names = [
        "SST_Med_C",
        "CO2_ppm",
        "pH_Total",
        "Omega_Arag",
        "Fleet_AER_gCO2_dwt_nm",
        "Marine_Insurance_USD_Mn",
    ]
    matrix = [
        [
            sst_med_c[i],
            co2_ppm[i],
            ph_total[i],
            omega_arag[i],
            fleet_aer[i],
            marine_ins_usd_mn[i],
        ]
        for i in range(n)
    ]

    return {
        "years": years,
        "feature_names": feature_names,
        "matrix": matrix,
        "series": {
            "sst_med_c": sst_med_c,
            "co2_ppm": co2_ppm,
            "ph_total": ph_total,
            "do_mg_l": do_mg_l,
            "omega_arag": omega_arag,
            "fleet_aer": fleet_aer,
            "marine_ins_usd_mn": marine_ins_usd_mn,
        },
    }


def run_gomd_data_mining_showcase() -> Dict[str, Any]:
    print("=" * 84)
    print("GOMD v2.0.0 — MULTIVARIATE OCEAN, CLIMATE & MARITIME DATA MINING PIPELINE")
    print("=" * 84)

    # -------------------------------------------------------------------------
    # STAGE 1: Panel Construction & Biogeochemical Feature Engineering
    # -------------------------------------------------------------------------
    panel = build_multivariate_ocean_panel(1995, 2025)
    years = panel["years"]
    s = panel["series"]

    gsw = gomd.teos10_gsw_stratification(
        salinity_psu=38.65,
        temp_c=15.4,
        pressure_dbar=120.0,
        lat=36.2,
        surface_salinity_psu=38.30,
        surface_temp_c=s["sst_med_c"][-1],
    )
    compound = gomd.compound_mhw_oax_event(
        sst_series=s["sst_med_c"],
        ph_series=s["ph_total"],
        do_mg_l_series=s["do_mg_l"],
        labels=[str(y) for y in years],
        sst_pct=80.0,
        ph_pct=20.0,
        do_pct=20.0,
    )
    sst_vel = gomd.sst_warming_velocity_profile(s["sst_med_c"], dt_days=365.25)
    triple_days = compound.get("triple_compound_days", compound.get("triple_compound_steps", 0))
    acc_val = sst_vel.get("warming_acceleration_c_per_decade2", sst_vel.get("mean_thermal_acceleration_c_per_day2", 0.0))

    print(f"\n[STAGE 1-2] Feature Engineering & Triple-Compound Extremes ({years[0]}–{years[-1]}, n={len(years)})")
    print(f"  • TEOS-10 Stratification : SA = {gsw['absolute_salinity_sa_g_kg']} g/kg | N = {gsw['buoyancy_frequency_cph']} cph ({gsw['stratification_regime']})")
    print(f"  • SST Warming Velocity   : {sst_vel['decadal_trend_c_per_decade']:+.4f} °C/decade | Acceleration = {acc_val:+.4f}")
    print(f"  • Compound MHW+OAX+LOX   : {triple_days} triple-threat years ({compound['triple_likelihood_multiplier']}x multiplier) | Mean AGI = {compound['mean_aerobic_growth_index']}")

    # -------------------------------------------------------------------------
    # STAGE 3: Unsupervised Data Mining (Correlation + PCA + K-Means++ Regimes)
    # -------------------------------------------------------------------------
    mining = gomd.multivariate_data_mining_pipeline(
        feature_matrix=panel["matrix"],
        feature_names=panel["feature_names"],
        period_labels=years,
        target_index=0,       # Target: Predict next-year Eastern Med SST (SST_Med_C)
        k_clusters=3,         # Discover k=3 macro-oceanic regimes
        ridge_alphas=(0.01, 0.1, 0.5, 1.0, 5.0, 10.0),
    )

    pca = mining["pca"]
    km = mining["kmeans_regimes"]
    print("\n[STAGE 3] Unsupervised Data Mining: PCA & K-Means++ Regime Discovery")
    print(f"  • PCA Explained Variance : PC1 = {pca['explained_variance_ratio'][0]*100:.2f}% | PC2 = {pca['explained_variance_ratio'][1]*100:.2f}% | Cumulative = {pca['cumulative_variance_pc1_pc2_pct']:.2f}%")
    print(f"  • PC1 Loadings (Thermo)  : " + ", ".join(f"{k}={v:+.3f}" for k, v in pca["loadings_pc1"].items()))
    print(f"  • K-Means++ Clustering   : k = {km['k']} Regimes | Silhouette Score = {km['silhouette_score']:.4f} | WCSS Inertia = {km['wcss_inertia']:.2f}")
    for reg in km["regimes"]:
        c_vals = reg["centroid_physical_units"]
        print(
            f"    - Regime {reg['regime_id']} ({reg['period_span']}, n={reg['n_periods']} yrs): "
            f"SST={c_vals['SST_Med_C']:.2f}°C, CO2={c_vals['CO2_ppm']:.1f} ppm, "
            f"pH={c_vals['pH_Total']:.4f}, Omega_arag={c_vals['Omega_Arag']:.3f}, "
            f"FleetAER={c_vals['Fleet_AER_gCO2_dwt_nm']:.2f}"
        )

    # -------------------------------------------------------------------------
    # STAGE 4: Non-Parametric Trend, Breakpoint & Causal Transfer Entropy
    # -------------------------------------------------------------------------
    mk_sst = gomd.mann_kendall_sen(s["sst_med_c"], years)
    cp_sst = gomd.pettitt_changepoint(s["sst_med_c"], labels=years)
    ipta_sst = gomd.innovative_polygon_trend(s["sst_med_c"])
    te_co2_sst = gomd.transfer_entropy_causal(s["co2_ppm"], s["sst_med_c"], bins=3)
    te_sst_ins = gomd.transfer_entropy_causal(s["sst_med_c"], s["marine_ins_usd_mn"], bins=3)

    print("\n[STAGE 4] Non-Parametric Econometrics, Breakpoint & Causal Information Flow")
    print(f"  • Mann-Kendall / Sen SST : Slope = {mk_sst['theil_sen_slope']:+.4f} °C/yr (Z = {mk_sst['z_stat']}, p = {mk_sst['p_value']}) -> {mk_sst['trend_classification']}")
    print(f"  • Pettitt Regime Shift   : Structural Break at Year {cp_sst['changepoint_label']} (Pre-Mean = {cp_sst['pre_regime_mean']:.2f}°C -> Post-Mean = {cp_sst['post_regime_mean']:.2f}°C, p = {cp_sst['p_value_approx']})")
    print(f"  • Şen IPTA Quantile Shift: Low = {ipta_sst['subregime_shifts']['low_quantile_shift']:+.3f}°C | Extreme High = {ipta_sst['subregime_shifts']['extreme_high_quantile_shift']:+.3f}°C (Tail Amp = {ipta_sst['tail_amplification_ratio']}x)")
    print(f"  • Shannon Transfer Entr. : TE(CO2 -> SST) = {te_co2_sst['transfer_entropy_x_to_y_bits']:.4f} bits (Net = {te_co2_sst['net_information_flow_bits']:+.4f} bits)")
    print(f"  • Climate->Insurance TE  : TE(SST -> Insurance) = {te_sst_ins['transfer_entropy_x_to_y_bits']:.4f} bits ({te_sst_ins['dominant_causal_direction']})")

    # -------------------------------------------------------------------------
    # STAGE 5: Supervised Machine Learning (Expanding-Window CV Ridge Regression)
    # -------------------------------------------------------------------------
    ml = mining["ridge_ml_forecast"]
    print("\n[STAGE 5] Supervised ML: Expanding-Window CV Multivariate Lag-1 Ridge Forecast")
    print(f"  • Target Variable        : {mining['target_variable']} (1-Year Ahead Out-of-Sample Forecast)")
    print(f"  • Optimal L2 Penalty α   : {ml['optimal_alpha_l2']} | Historical CV RMSE = {ml['cv_rmse']:.4f} °C")
    print(f"  • Out-of-Sample Test     : Ridge RMSE = {ml['test_rmse']:.4f} °C vs. Naive RMSE = {ml['naive_baseline_rmse']:.4f} °C | Skill Improvement = {ml['forecast_skill_score_pct']:+.2f}%")
    print(f"  • Standardized Weights β : " + ", ".join(f"{k}={v:+.4f}" for k, v in ml["standardized_coefficients"].items()))

    # -------------------------------------------------------------------------
    # STAGE 6: Maritime Decarbonization, Actuarial Insurance & Mavi Vatan Audit
    # -------------------------------------------------------------------------
    fueleu = gomd.calculate_fueleu_maritime(
        fuel_mix_tonnes={"VLSFO": 3800.0, "BIO_LNG": 950.0, "E_METHANOL": 350.0},
        year=2030,
    )
    ins_pricing = gomd.marine_insurance_actuarial_pricing(
        hull_value_usd=72_000_000.0,
        gross_tonnage=48_000.0,
        vessel_age_years=7,
        cii_rating="B",
        psc_detentions_3yr=0,
        war_risk_chokepoint_zone=True,
        green_dual_fuel=True,
    )
    mv_audit = gomd.mavi_vatan_eez_audit(
        mainland_coast_km=1577.0,
        opposite_coast_km=19.0,
        un_registered_area_km2=462000.0,
        provider_polygon_area_km2=254026.0,
    )

    fueleu_wtw = fueleu.get("attained_wtw_ghg_intensity_gco2e_mj", fueleu.get("attained_wtw_gco2e_mj", 0.0))
    fueleu_status = fueleu.get("compliance_status", fueleu.get("status", ""))
    ins_total = ins_pricing.get("total_annual_premium_usd", ins_pricing.get("total_annual_marine_insurance_usd", 0.0))
    ins_ppmi = ins_pricing.get("ppmi_climate_alignment_adjustment_pct", ins_pricing.get("ppmi_climate_alignment_delta_pct", 0.0))
    mv_ratio = mv_audit.get("coastal_length_ratio_mainland_to_island", mv_audit.get("coastal_length_proportionality_ratio", 0.0))
    mv_area = mv_audit.get("un_registered_mavi_vatan_km2", mv_audit.get("un_registered_mavi_vatan_area_km2", 462000.0))

    print("\n[STAGE 6] Maritime Decarbonization, Marine Underwriting & Mavi Vatan Spatial Audit")
    print(f"  • FuelEU Maritime (2030) : WtW Intensity = {fueleu_wtw:.2f} gCO2e/MJ ({fueleu_status})")
    print(f"  • PPMI Actuarial Pricing : Total Annual Premium = ${ins_total:,.2f} (PPMI Alignment Delta = {ins_ppmi:+.1f}%)")
    print(f"  • Mavi Vatan ICJ Audit   : Coastal Ratio = {mv_ratio}:1 | Official UN A/74/757 Area = {mv_area:,.0f} km²")
    print("=" * 84)

    return {
        "sdk_version": gomd.__version__,
        "stratification": gsw,
        "compound_extremes": compound,
        "sst_velocity": sst_vel,
        "multivariate_mining": mining,
        "mann_kendall_sst": mk_sst,
        "pettitt_sst": cp_sst,
        "ipta_sst": ipta_sst,
        "transfer_entropy_co2_to_sst": te_co2_sst,
        "fueleu_2030": fueleu,
        "actuarial_insurance": ins_pricing,
        "mavi_vatan_audit": mv_audit,
    }


if __name__ == "__main__":
    results = run_gomd_data_mining_showcase()
    if "--json" in sys.argv:
        print(json.dumps(results, indent=2, ensure_ascii=False))