#!/usr/bin/env python3
"""
=============================================================================
SUMATRA LOG-GAUSSIAN COX PROCESS (LGCP)
=============================================================================

Fits a parametric LGCP with two geological covariates as the mean function:
  μ(s) = β₀ + β₁·φ(s) + β₂·D_fault(s)

where
  φ(s)       = Sunda megathrust interseismic coupling (Chlieh et al. 2008)
  D_fault(s) = Distance to the Great Sumatran Fault trace (GEM/EOS SE Asia)

PIPELINE
  Step 1 — Load G-K mainshock catalogue (M≥4.5)
  Step 2 — Build covariate grids (coupling map + fault distance)
  Step 3 — Poisson log-linear MLE via Newton-Raphson on quadrature grid
  Step 4 — Pearson residuals on NQ×NQ grid + empirical variogram
  Step 5 — Fit exponential covariance model by weighted least squares
  Step 6 — Figures 11-13 and lgcp_parameters.csv

MODEL
  log λ(s) = μ(s) + Z(s)
  μ(s) = Xβ  (fixed trend from covariates)
  Z(s) ~ GRF(0, C(h))  with C(h) = σ² exp(−h/φ_range)  [nugget C₀]

COVARIATE SOURCES
  Chlieh et al. (2008) — Table 1 segment coupling ratios φ (digitised)
  GEM/EOS Southeast Asia / Sieh & Natawidjaja (2000) — Sumatran Fault trace

References:
  Møller et al. (1998) Scand. J. Stat.
  Diggle et al. (2013) Statistical Science
  Chlieh et al. (2008) J. Geophys. Res.
  GEM Global Active Faults / EOS Southeast Asia active-fault compilation.

Dependencies: numpy, pandas, matplotlib (no scipy)
=============================================================================
"""

import numpy as np
import pandas as pd
import matplotlib
import os, warnings, time
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes as _inset_axes

def _cartopy_cbar(fig, mappable, ax, label):
    """Colorbar that exactly matches Cartopy GeoAxes height."""
    cax = _inset_axes(ax, width="4%", height="100%", loc="lower left",
                      bbox_to_anchor=(1.03, 0., 1, 1),
                      bbox_transform=ax.transAxes, borderpad=0)
    cb = fig.colorbar(mappable, cax=cax, label=label)
    return cb
import matplotlib.gridspec as gridspec
warnings.filterwarnings('ignore')

t_start = time.time()

# ─────────────────────────────────────────────────────────────────────────────
#  CONFIGURATION
# ─────────────────────────────────────────────────────────────────────────────
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
GK_CSV   = os.path.join(SCRIPT_DIR, "declustered_gk.csv")
OUT_DIR  = SCRIPT_DIR
os.makedirs(OUT_DIR, exist_ok=True)
cartopy.config['data_dir'] = os.environ.get('CARTOPY_DATA_DIR',
                                            '/private/tmp/cartopy_data')

LAT_MIN, LAT_MAX = -6.0,  6.0
LON_MIN, LON_MAX = 95.0, 108.0
lat0, lon0 = 0.0, 101.5
KM_LON = 111.320 * np.cos(np.radians(lat0))
KM_LAT = 111.320

X_MIN = (LON_MIN - lon0) * KM_LON
X_MAX = (LON_MAX - lon0) * KM_LON
Y_MIN = (LAT_MIN - lat0) * KM_LAT
Y_MAX = (LAT_MAX - lat0) * KM_LAT
W_AREA = (X_MAX - X_MIN) * (Y_MAX - Y_MIN)   # km²

MC_LGCP  = 4.5    # magnitude completeness for LGCP (matches G-K declustering threshold)
NQ       = 60     # quadrature grid (NQ × NQ cells)
MAX_ITER = 40     # Newton-Raphson max iterations
STEP_DAMP = 0.5  # step damping for NR stability

# Variogram
H_MAX_VARIO = 700.0   # km
N_LAGS      = 28
MAX_VARIO_PTS = 2000  # subsample for speed

# ─────────────────────────────────────────────────────────────────────────────
#  COVARIATE DATA (digitised from published sources)
# ─────────────────────────────────────────────────────────────────────────────

# More accurate Sunda Megathrust trench positions (Andaman to Java)
# Based on recent seismic rupture zone studies
TRENCH_LATS = np.array([13.5, 12.5, 11.5, 10.5, 9.5, 8.5, 8.0, 6.5, 5.5, 4.5,
                        4.0, 3.0, 2.0, 1.0, 0.0, -0.5, -1.0, -1.5, -2.0, -2.5,
                        -2.8, -3.2, -3.5, -3.8, -4.2, -4.5, -5.0, -5.5, -6.0, -6.5])
TRENCH_LONS = np.array([92.2, 92.3, 92.4, 92.5, 92.6, 92.7, 92.8, 93.0, 93.2, 93.4,
                        93.5, 93.6, 93.8, 94.2, 95.0, 95.5, 96.0, 96.5, 97.2, 98.0,
                        98.5, 99.0, 99.5, 100.0, 100.5, 101.0, 101.5, 102.0, 102.5, 103.0])

# Interpolated coupling ratios φ for each node (from Chlieh et al. 2008 segments)
# Andaman (φ≈0.08), Nicobar (φ≈0.15), N.Sumatra (φ≈0.27), C.Sumatra (φ≈0.38→0.73),
# Mentawai (φ≈0.88→0.65), S.Sumatra (φ≈0.54)
COUPLING = np.array([0.08, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, 0.17, 0.20, 0.25,
                     0.28, 0.32, 0.38, 0.50, 0.65, 0.70, 0.75, 0.80, 0.85, 0.88,
                     0.87, 0.85, 0.80, 0.75, 0.68, 0.62, 0.60, 0.57, 0.55, 0.54])

# --- Great Sumatran Fault trace ---
# Densified active-fault polyline from the GEM/EOS Southeast Asia active-fault
# compilation, consistent with Sieh & Natawidjaja (2000) segment geometry.
FAULT_LATS = np.array([
     5.912,  5.737,  5.021,  4.414,  4.293,  4.117,  3.560,
     3.218,  2.180,  2.017,  1.844,  1.507,  0.704,  0.520,
     0.251,  0.205,  0.080, -0.026, -0.072, -0.555, -0.728,
    -0.903, -1.148, -1.596, -1.797, -2.193, -2.334, -2.731,
    -2.777, -3.262, -3.441, -3.577, -4.304, -4.350, -5.177,
    -5.313, -5.834, -6.154, -6.604, -6.697
])
FAULT_LONS = np.array([
     95.326,  95.094,  95.951,  96.593,  96.813,  97.064,  97.849,
     97.987,  98.893,  98.985,  99.115,  99.331, 100.008, 100.133,
    100.054, 100.155, 100.247, 100.155, 100.247, 100.561, 100.607,
    100.727, 100.852, 101.168, 101.214, 101.606, 101.560, 101.866,
    102.034, 102.370, 102.511, 102.557, 103.210, 103.421, 104.293,
    104.339, 104.764, 104.883, 104.837, 104.764
])

# ─────────────────────────────────────────────────────────────────────────────
#  HELPER FUNCTIONS
# ─────────────────────────────────────────────────────────────────────────────

def to_km(lats, lons):
    """Convert lat/lon arrays to projected km coordinates."""
    return (lons - lon0) * KM_LON, (lats - lat0) * KM_LAT


def coupling_idw(x_km, y_km, tx, ty, phi, reg=50.0**2):
    """
    IDW interpolation of coupling ratio φ at arbitrary (x, y) points.
    reg: regularisation to avoid division by zero (50 km offset).
    """
    dx = x_km[:, None] - tx[None, :]   # (n_pts, n_nodes)
    dy = y_km[:, None] - ty[None, :]
    w  = 1.0 / (dx**2 + dy**2 + reg)
    return (w * phi[None, :]).sum(axis=1) / w.sum(axis=1)


def dist_polyline(x_km, y_km, fx, fy):
    """
    Minimum distance from each point (x_km, y_km) to the polyline (fx, fy).
    Uses exact perpendicular distance per segment, clamped to endpoints.
    """
    min_d = np.full(len(x_km), np.inf)
    for k in range(len(fx) - 1):
        ax, ay = fx[k],   fy[k]
        bx, by = fx[k+1], fy[k+1]
        abx, aby = bx - ax, by - ay
        ab2 = abx**2 + aby**2
        if ab2 < 1e-9:
            d = np.sqrt((x_km - ax)**2 + (y_km - ay)**2)
        else:
            t  = np.clip(((x_km - ax)*abx + (y_km - ay)*aby) / ab2, 0.0, 1.0)
            px = ax + t*abx - x_km
            py = ay + t*aby - y_km
            d  = np.sqrt(px**2 + py**2)
        min_d = np.minimum(min_d, d)
    return min_d


print("=" * 66)
print("  SUMATRA LOG-GAUSSIAN COX PROCESS (LGCP)")
print("  Covariates: Chlieh (2008) coupling + GEM/EOS fault")
print("=" * 66)

# ─────────────────────────────────────────────────────────────────────────────
#  STEP 1 — LOAD DATA
# ─────────────────────────────────────────────────────────────────────────────
print("\n[1] Loading catalogue …")
df_raw = pd.read_csv(GK_CSV)

# Spatial filter
mask_sp = ((df_raw['Latitude']  >= LAT_MIN) & (df_raw['Latitude']  <= LAT_MAX) &
           (df_raw['Longitude'] >= LON_MIN) & (df_raw['Longitude'] <= LON_MAX))

# Use G-K mainshocks M≥4.5
mask_ms = (df_raw['is_mainshock'] == 1) & (df_raw['Magnitude'] >= MC_LGCP)
ev = df_raw[mask_sp & mask_ms].copy().reset_index(drop=True)

ev['x'], ev['y'] = to_km(ev['Latitude'].values, ev['Longitude'].values)
ev_x = ev['x'].values
ev_y = ev['y'].values
n_ev = len(ev)
lam_csr = n_ev / W_AREA

print(f"    Mainshocks M≥{MC_LGCP} (G-K declustered): {n_ev:,}")
print(f"    Study area:     {W_AREA/1e6:.2f}×10⁶ km²")
print(f"    CSR intensity:  {lam_csr:.4e} events km⁻²")

# ─────────────────────────────────────────────────────────────────────────────
#  STEP 2 — COVARIATE SURFACES
# ─────────────────────────────────────────────────────────────────────────────
print("\n[2] Building covariate surfaces …")

# Project trench and fault nodes to km
TRENCH_X, TRENCH_Y = to_km(TRENCH_LATS, TRENCH_LONS)
FAULT_X,  FAULT_Y  = to_km(FAULT_LATS,  FAULT_LONS)

# Quadrature grid
xq_1d = np.linspace(X_MIN, X_MAX, NQ)
yq_1d = np.linspace(Y_MIN, Y_MAX, NQ)
XQ, YQ = np.meshgrid(xq_1d, yq_1d)   # (NQ, NQ)
xq_f   = XQ.ravel()
yq_f   = YQ.ravel()
cell_area = W_AREA / (NQ * NQ)        # km²

# Build covariates on quadrature grid
phi_q    = coupling_idw(xq_f, yq_f, TRENCH_X, TRENCH_Y, COUPLING).reshape(NQ, NQ)
dfault_q = dist_polyline(xq_f, yq_f, FAULT_X, FAULT_Y).reshape(NQ, NQ)

# Standardise: z = (x − mean) / std  for numerical stability in NR
phi_mu, phi_sd    = phi_q.mean(),    phi_q.std()
df_mu,  df_sd     = dfault_q.mean(), dfault_q.std()

phi_q_z    = (phi_q    - phi_mu) / phi_sd
dfault_q_z = (dfault_q - df_mu)  / df_sd

# Covariates at event locations
phi_ev    = coupling_idw(ev_x, ev_y, TRENCH_X, TRENCH_Y, COUPLING)
dfault_ev = dist_polyline(ev_x, ev_y, FAULT_X, FAULT_Y)
phi_ev_z    = (phi_ev    - phi_mu) / phi_sd
dfault_ev_z = (dfault_ev - df_mu)  / df_sd

print(f"    Coupling φ:      [{phi_q.min():.3f}, {phi_q.max():.3f}]")
print(f"    Fault distance:  [{dfault_q.min():.0f}, {dfault_q.max():.0f}] km")

# ─────────────────────────────────────────────────────────────────────────────
#  STEP 3 — POISSON LOG-LINEAR MLE via NEWTON-RAPHSON
# ─────────────────────────────────────────────────────────────────────────────
print("\n[3] Poisson log-linear MLE (Newton-Raphson) …")

# Design matrices
# X_ev : (n_ev, 3)  — [intercept, φ_z, D_z] at event locations
# X_q  : (NQ², 3)   — same at quadrature nodes
X_ev = np.column_stack([np.ones(n_ev),    phi_ev_z,       dfault_ev_z])
X_q  = np.column_stack([np.ones(NQ*NQ),   phi_q_z.ravel(), dfault_q_z.ravel()])

# Initialise β at log(CSR)
beta = np.array([np.log(lam_csr), 0.0, 0.0])
tol  = 1e-7

for it in range(MAX_ITER):
    eta_q  = X_q @ beta
    # Clip for numerical safety
    eta_q  = np.clip(eta_q, -20, 20)
    lam_q  = np.exp(eta_q)               # (NQ²,)

    # Score: ∂L/∂β = Σᵢ x(sᵢ) − Σₖ λ(sₖ)·x(sₖ)·Δs
    score  = X_ev.sum(axis=0) - cell_area * (X_q * lam_q[:, None]).sum(axis=0)

    # Fisher information (expected Hessian, negated): Σₖ λ(sₖ)·xᵢxⱼᵀ·Δs
    W_lam  = lam_q * cell_area
    FI     = (X_q.T * W_lam[None, :]) @ X_q   # (3, 3)

    # Newton step  β += FI⁻¹ · score  (damped for stability)
    try:
        step = np.linalg.solve(FI + 1e-8*np.eye(3), score)
    except np.linalg.LinAlgError:
        print(f"    Singular Fisher info at iter {it} — stopping")
        break

    beta += STEP_DAMP * step

    grad_norm = np.max(np.abs(score))
    if grad_norm < tol:
        print(f"    Converged (‖score‖∞={grad_norm:.2e}) at iteration {it+1}")
        break
else:
    print(f"    Max iterations reached (‖score‖∞={grad_norm:.2e})")

# Fitted intensity on quadrature grid and at events
eta_q_fit  = np.clip(X_q @ beta, -20, 20)
lam_q_fit  = np.exp(eta_q_fit).reshape(NQ, NQ)
Lambda_hat = lam_q_fit.sum() * cell_area       # integrated intensity

eta_ev_fit = X_ev @ beta
lam_ev_fit = np.exp(eta_ev_fit)

# Approximate log-likelihood and null deviance
LL_model = np.sum(eta_ev_fit) - Lambda_hat
LL_null  = n_ev * np.log(lam_csr) - lam_csr * W_AREA
D_stat   = -2 * (LL_null - LL_model)

print(f"    β₀ (intercept)  = {beta[0]:.4f}")
print(f"    β₁ (coupling φ) = {beta[1]:.4f}  →  exp(β₁) = {np.exp(beta[1]):.3f}")
print(f"    β₂ (fault dist) = {beta[2]:.4f}  →  exp(β₂) = {np.exp(beta[2]):.3f}")
print(f"    Λ̂ = {Lambda_hat:.1f}  (n_obs = {n_ev})")
print(f"    Likelihood-ratio statistic: {D_stat:.1f}  (df=2)")

# ─────────────────────────────────────────────────────────────────────────────
#  STEP 4 — PEARSON RESIDUALS + EMPIRICAL VARIOGRAM
# ─────────────────────────────────────────────────────────────────────────────
print("\n[4] Computing Pearson residuals and variogram …")

# Count events per quadrature cell
xq_edges = np.linspace(X_MIN, X_MAX, NQ+1)
yq_edges = np.linspace(Y_MIN, Y_MAX, NQ+1)
xi = np.clip(np.searchsorted(xq_edges, ev_x, side='right') - 1, 0, NQ-1)
yi = np.clip(np.searchsorted(yq_edges, ev_y, side='right') - 1, 0, NQ-1)
counts = np.zeros((NQ, NQ))
np.add.at(counts, (yi, xi), 1)

# Expected counts and Pearson residuals
expected = lam_q_fit * cell_area           # (NQ, NQ)
exp_safe  = np.maximum(expected, 1e-9)
pearson   = (counts - expected) / np.sqrt(exp_safe)   # (NQ, NQ)

# Cells with non-trivial expected count
mask_v = (exp_safe.ravel() > 0.02)
xv  = xq_f[mask_v]
yv  = yq_f[mask_v]
rv  = pearson.ravel()[mask_v]

# Sub-sample for tractable pairwise computation
rng = np.random.default_rng(42)
if len(xv) > MAX_VARIO_PTS:
    idx = rng.choice(len(xv), MAX_VARIO_PTS, replace=False)
    xv, yv, rv = xv[idx], yv[idx], rv[idx]

# All pairwise distances and squared differences (upper triangle)
n_v  = len(xv)
dx2  = xv[:, None] - xv[None, :]
dy2  = yv[:, None] - yv[None, :]
dr   = np.sqrt(dx2**2 + dy2**2)
dz2  = (rv[:, None] - rv[None, :])**2

iu  = np.triu_indices(n_v, k=1)
dr_f = dr[iu]
dz_f = dz2[iu]

H_BINS  = np.linspace(0, H_MAX_VARIO, N_LAGS + 1)
h_mid   = 0.5 * (H_BINS[:-1] + H_BINS[1:])
gamma_h = np.full(N_LAGS, np.nan)
n_pairs = np.zeros(N_LAGS)

for k in range(N_LAGS):
    m = (dr_f >= H_BINS[k]) & (dr_f < H_BINS[k+1])
    if m.sum() >= 10:
        gamma_h[k] = 0.5 * dz_f[m].mean()
        n_pairs[k] = m.sum()

ok = ~np.isnan(gamma_h)
print(f"    Variogram: {ok.sum()} lag classes,  {len(dr_f):,} pairs")
print(f"    γ range: [{np.nanmin(gamma_h):.4f}, {np.nanmax(gamma_h):.4f}]")

# ─────────────────────────────────────────────────────────────────────────────
#  STEP 5 — FIT EXPONENTIAL COVARIANCE  C(h) = C₀ + σ²[1 − exp(−h/φ)]
# ─────────────────────────────────────────────────────────────────────────────
print("\n[5] Fitting exponential covariance by weighted LS …")

h_ok = h_mid[ok]
g_ok = gamma_h[ok]
w_ok = n_pairs[ok] / n_pairs[ok].max()   # normalised weights

g_max = g_ok.max()
NUGGET_GRID = np.linspace(0, g_max * 0.4, 8)
SILL_GRID   = np.linspace(g_max * 0.2, g_max * 1.5, 18)
RANGE_GRID  = np.logspace(1.0, 2.9, 24)   # 10–800 km

best_wss = np.inf
nugget_hat = sill_hat = range_hat = np.nan

for c0 in NUGGET_GRID:
    for s2 in SILL_GRID:
        for phi_r in RANGE_GRID:
            g_th = c0 + s2 * (1.0 - np.exp(-h_ok / phi_r))
            wss  = np.sum(w_ok * (g_ok - g_th)**2)
            if wss < best_wss:
                best_wss  = wss
                nugget_hat = c0
                sill_hat   = s2
                range_hat  = phi_r

print(f"    Nugget C₀  = {nugget_hat:.4f}")
print(f"    Sill   σ²  = {sill_hat:.4f}")
print(f"    Range  φ   = {range_hat:.1f} km")
print(f"    Practical range (3φ) = {3*range_hat:.0f} km")

# ─────────────────────────────────────────────────────────────────────────────
#  STEP 6 — FIGURES
# ─────────────────────────────────────────────────────────────────────────────
print("\n[6] Generating figures …")

# Back-project quadrature grid to lon/lat for plotting
LON_Q = lon0 + XQ / KM_LON
LAT_Q = lat0 + YQ / KM_LAT
LON_TRENCH = lon0 + TRENCH_X / KM_LON
LAT_TRENCH = lat0 + TRENCH_Y / KM_LAT
LON_FAULT  = lon0 + FAULT_X  / KM_LON
LAT_FAULT  = lat0 + FAULT_Y  / KM_LAT

# ── Figure 11: Covariate maps (Cartopy) ─────────────────────────────────────
_PC = ccrs.PlateCarree()
_COAST_10M = cfeature.NaturalEarthFeature(
    'physical', 'coastline', '10m',
    edgecolor='#111111', facecolor='none', linewidth=1.15)
_BORDER_10M = cfeature.NaturalEarthFeature(
    'cultural', 'admin_0_boundary_lines_land', '10m',
    edgecolor='#444444', facecolor='none', linewidth=0.35, alpha=0.65)
fig, axes = plt.subplots(1, 2, figsize=(16, 7),
                         subplot_kw={'projection': _PC})

def _add_coast(ax):
    ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=_PC)
    ax.set_aspect('auto')
    ax.set_facecolor('white')
    ax.add_feature(_COAST_10M, zorder=8)
    ax.add_feature(_BORDER_10M, zorder=8)
    gl = ax.gridlines(crs=_PC, draw_labels=True, linewidth=0.3, alpha=0.4,
                      color='gray', x_inline=False, y_inline=False)
    gl.top_labels = False; gl.right_labels = False
    gl.xlabel_style = {'size': 8}; gl.ylabel_style = {'size': 8}

ax = axes[0]
_add_coast(ax)
levels_phi = np.linspace(COUPLING.min(), COUPLING.max(), 22)
cf = ax.contourf(LON_Q, LAT_Q, phi_q, levels=levels_phi,
                 cmap='RdYlBu_r', transform=_PC, zorder=2)
_cartopy_cbar(fig, cf, ax, 'Coupling ratio φ')
ax.contour(LON_Q, LAT_Q, phi_q, levels=[0.5, 0.7, 0.85],
           colors='k', linewidths=[0.8, 1.2, 1.5], linestyles='-',
           transform=_PC, zorder=3)
# Historical-rupture and Mentawai-gap overlay
ax.plot(LON_FAULT, LAT_FAULT, 'k--', lw=2.0, transform=_PC, zorder=4,
        label='Sumatran Fault\n(GEM/EOS)')
ax.scatter(ev['Longitude'], ev['Latitude'], s=1.5, c='k',
           alpha=0.20, zorder=5, transform=_PC, label=f'M≥{MC_LGCP} mainshocks')
ax.text(0.97, 0.97, "(A)", transform=ax.transAxes,
    fontsize=13, fontweight="bold", ha="right", va="top", zorder=10)
ax.legend(fontsize=8, loc='lower right')

ax = axes[1]
_add_coast(ax)
dfault_plot = dist_polyline(xq_f, yq_f, FAULT_X, FAULT_Y).reshape(NQ, NQ)
cf2 = ax.contourf(LON_Q, LAT_Q, dfault_plot, levels=22,
                  cmap='YlOrRd_r', transform=_PC, zorder=2)
_cartopy_cbar(fig, cf2, ax, 'Distance to fault (km)')
ax.contour(LON_Q, LAT_Q, dfault_plot, levels=[50, 100, 200, 300],
           colors='k', linewidths=0.7, linestyles='--',
           transform=_PC, zorder=3)
# Same rupture-polygon overlay on panel B (no extra credit text — it is
# already shown in panel A)
# rupture overlay removed
ax.plot(LON_FAULT, LAT_FAULT, 'k--', lw=2.5, transform=_PC,
        zorder=4, label='Sumatran Fault')
ax.scatter(ev['Longitude'], ev['Latitude'], s=1.5, c='w',
           alpha=0.30, zorder=5, transform=_PC)
ax.text(0.97, 0.97, "(B)", transform=ax.transAxes,
    fontsize=13, fontweight="bold", ha="right", va="top", zorder=10)
ax.legend(fontsize=8, loc='lower right')

# suptitle removed
fig.savefig(f'{OUT_DIR}/fig1_covariates.png', dpi=300, bbox_inches='tight')
plt.close(fig);  print("    Saved fig1_covariates.png")

# ── Figure 12: LGCP intensity decomposition (Cartopy) ────────────────────────
fig = plt.figure(figsize=(18, 6))
gs  = gridspec.GridSpec(1, 3, wspace=0.35)

def _add_coast12(ax, left_labels=False, bottom_labels=True):
    ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=_PC)
    ax.set_aspect('auto')
    ax.set_facecolor('white')
    ax.add_feature(_COAST_10M, zorder=8)
    ax.add_feature(_BORDER_10M, zorder=8)
    gl = ax.gridlines(crs=_PC, draw_labels=True, linewidth=0.3,
                      alpha=0.45, color='#666', linestyle='--',
                      x_inline=False, y_inline=False)
    gl.top_labels   = False
    gl.right_labels = False
    gl.left_labels  = left_labels
    gl.bottom_labels = bottom_labels
    gl.xlabel_style = {'size': 8}
    gl.ylabel_style = {'size': 8}

# Panel A: kernel-smoothed (adaptive KDE, h=40 km)
ax = fig.add_subplot(gs[0], projection=_PC)
_add_coast12(ax, left_labels=True, bottom_labels=True)
h_kde = 40.0
kde_v = np.zeros(NQ*NQ)
chunk = 300
for s in range(0, n_ev, chunk):
    xe = ev_x[s:s+chunk, None]
    ye = ev_y[s:s+chunk, None]
    kde_v += np.sum(np.exp(-((xq_f[None,:]-xe)**2 +
                              (yq_f[None,:]-ye)**2) / (2*h_kde**2)), axis=0)
kde_grid = (kde_v / (n_ev * 2*np.pi*h_kde**2)).reshape(NQ, NQ)
cf = ax.contourf(LON_Q, LAT_Q, kde_grid * 1e6, levels=20,
                 cmap='hot_r', transform=_PC, zorder=2)
_cartopy_cbar(fig, cf, ax, 'λ (×10⁻⁶ km⁻²)')
ax.plot(LON_FAULT, LAT_FAULT, 'w--', lw=1.5, transform=_PC, zorder=3)
ax.text(0.97, 0.97, "(A)", transform=ax.transAxes,
    fontsize=13, fontweight="bold", ha="right", va="top", zorder=10)

# Panel B: LGCP fitted mean μ(s)
ax = fig.add_subplot(gs[1], projection=_PC)
_add_coast12(ax, left_labels=False, bottom_labels=True)
cf2 = ax.contourf(LON_Q, LAT_Q, lam_q_fit * 1e6, levels=20,
                  cmap='hot_r', transform=_PC, zorder=2)
_cartopy_cbar(fig, cf2, ax, 'λ̂(s) (×10⁻⁶ km⁻²)')
# Rupture-polygon overlay on the fitted LGCP trend (Fig 12 panel B):
# shows that the fitted intensity ridge sits where the great events ruptured,
# and that the Mentawai gap (locked) is the highest-rate predicted patch.
# rupture overlay removed
ax.plot(LON_FAULT, LAT_FAULT, 'w--', lw=1.5,
        transform=_PC, zorder=3, label='S. Fault')
b0, b1, b2 = beta
sign1 = '+' if b1 >= 0 else '−'
sign2 = '+' if b2 >= 0 else '−'
ax.text(0.97, 0.97, "(B)", transform=ax.transAxes,
    fontsize=13, fontweight="bold", ha="right", va="top", zorder=10)
ax.legend(fontsize=8, loc='lower right')

# Panel C: Pearson residuals
ax = fig.add_subplot(gs[2], projection=_PC)
_add_coast12(ax, left_labels=False, bottom_labels=True)
vmax = np.percentile(np.abs(pearson), 97)
cf3  = ax.contourf(LON_Q, LAT_Q, pearson,
                   levels=np.linspace(-vmax, vmax, 22),
                   cmap='RdBu_r', transform=_PC, zorder=2)
_cartopy_cbar(fig, cf3, ax, 'Pearson residual')
ax.contour(LON_Q, LAT_Q, pearson, levels=[0],
           colors='k', linewidths=0.8, transform=_PC, zorder=3)
ax.plot(LON_FAULT, LAT_FAULT, 'k--', lw=1.5,
        transform=_PC, zorder=3, label='S. Fault')
ax.text(0.97, 0.97, "(C)", transform=ax.transAxes,
    fontsize=13, fontweight="bold", ha="right", va="top", zorder=10)
ax.legend(fontsize=8, loc='lower right')

# suptitle removed', fontsize=13, fontweight='bold')
fig.savefig(f'{OUT_DIR}/fig3_intensity.png', dpi=300, bbox_inches='tight')
plt.close(fig);  print("    Saved fig3_intensity.png")

# ── Figure 13: Variogram + covariance ───────────────────────────────────────
h_plot = np.linspace(0, H_MAX_VARIO, 400)
g_fit  = nugget_hat + sill_hat * (1.0 - np.exp(-h_plot / range_hat))
C_fit  = sill_hat  * np.exp(-h_plot / range_hat)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

ax = axes[0]
dot_size = 10 + 100 * n_pairs[ok] / n_pairs[ok].max()
ax.scatter(h_ok, g_ok, s=dot_size, c='steelblue', zorder=5, edgecolors='navy',
           linewidths=0.5, label='Empirical γ̂(h)\n(size ∝ pair count)')
ax.plot(h_plot, g_fit, 'r-', lw=2.5,
        label=(fr'Exp. fit: $\hat C_0={nugget_hat:.3f}$, '
               fr'$\hat\tau^2={sill_hat:.3f}$, '
               fr'$\hat\phi_r={range_hat:.0f}$ km'))
ax.axhline(nugget_hat + sill_hat, color='k', ls='--', lw=1.0,
           label=f'Sill $C_0+\\tau^2 = {nugget_hat+sill_hat:.3f}$')
ax.axhline(nugget_hat, color='gray', ls=':', lw=1.0, label=f'Nugget $C_0 = {nugget_hat:.3f}$')
ax.set_xlabel('Lag distance h (km)', fontsize=10)
ax.set_ylabel('γ̂(h)', fontsize=10)
ax.text(0.97, 0.97, "(A)", transform=ax.transAxes,
    fontsize=13, fontweight="bold", ha="right", va="top", zorder=10)
ax.legend(fontsize=8, loc='lower right');  ax.grid(True, alpha=0.25)
ax.set_xlim(0, H_MAX_VARIO)

ax = axes[1]
ax.plot(h_plot, C_fit + nugget_hat, 'r-', lw=2.5, label=r'$C(h)=C_0+\tau^2 e^{-h/\phi_r}$')
ax.fill_between(h_plot, nugget_hat, C_fit + nugget_hat,
                alpha=0.12, color='red', label=r'Structured variance $\tau^2$')
ax.axvline(range_hat,   color='navy',  ls='--', lw=1.5,
           label=f'Range $\\hat{{\\phi}}_r = {range_hat:.0f}$ km')
ax.axvline(3*range_hat, color='green', ls='--', lw=1.5,
           label=f'Practical range $3\\hat{{\\phi}}_r = {3*range_hat:.0f}$ km')
ax.axhline(nugget_hat,  color='k',     ls=':',  lw=1.0,
           label=f'Nugget $C_0 = {nugget_hat:.3f}$')
ax.set_xlabel('Lag distance h (km)', fontsize=10)
ax.set_ylabel('Covariance C(h)', fontsize=10)
ax.text(0.97, 0.97, "(B)", transform=ax.transAxes,
    fontsize=13, fontweight="bold", ha="right", va="top", zorder=10)
ax.legend(fontsize=8);  ax.grid(True, alpha=0.25)
ax.set_xlim(0, H_MAX_VARIO)

# suptitle removed
fig.savefig(f'{OUT_DIR}/fig2_variogram.png', dpi=300, bbox_inches='tight')
plt.close(fig);  print("    Saved fig2_variogram.png")

# ─────────────────────────────────────────────────────────────────────────────
#  STEP 7 — SAVE RESULTS CSV
# ─────────────────────────────────────────────────────────────────────────────
params = {
    'Parameter': ['MC_LGCP', 'n_events', 'W_area_km2', 'NQ_grid',
                  'beta0', 'beta1_coupling', 'beta2_fault_dist',
                  'exp_beta1_coupling', 'exp_beta2_fault',
                  'Lambda_hat', 'LR_statistic_df2',
                  'nugget_C0', 'sill_sigma2', 'range_phi_km',
                  'practical_range_km', 'phi_mean', 'phi_std',
                  'dfault_mean_km', 'dfault_std_km'],
    'Value':     [MC_LGCP, n_ev, round(W_AREA, 0), NQ,
                  round(beta[0],5), round(beta[1],5), round(beta[2],5),
                  round(np.exp(beta[1]),4), round(np.exp(beta[2]),4),
                  round(Lambda_hat,1), round(D_stat,1),
                  round(nugget_hat,5), round(sill_hat,5), round(range_hat,1),
                  round(3*range_hat,1), round(phi_mu,4), round(phi_sd,4),
                  round(df_mu,1), round(df_sd,1)]
}
pd.DataFrame(params).to_csv(f'{OUT_DIR}/lgcp_parameters.csv', index=False)

elapsed = time.time() - t_start
print(f"\n{'='*66}")
print(f"  LGCP RESULTS SUMMARY")
print(f"{'─'*66}")
print(f"  Events M≥{MC_LGCP}: {n_ev:,}")
print(f"  β₀ = {beta[0]:.4f},  β₁(φ) = {beta[1]:.4f},  β₂(D) = {beta[2]:.4f}")
print(f"  Coupling effect:  exp(β₁) = {np.exp(beta[1]):.3f}")
print(f"  Fault-dist effect: exp(β₂) = {np.exp(beta[2]):.3f}")
print(f"  LR statistic = {D_stat:.1f}  (χ²₂ 1% critical = 9.21)")
print(f"  GRF:  C₀={nugget_hat:.4f},  σ²={sill_hat:.4f},  φ={range_hat:.0f} km")
print(f"  Elapsed: {elapsed:.1f} s")
print(f"{'='*66}")
