#!/usr/bin/env python3
"""
scripts/02_homola_replication.py
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Reproduce the Homola et al. 2023 (JASTP 247, 106068) cross-correlation analysis
as literally as possible from the paper.

Method (from the paper)
-----------------------
* Data:  NMDB pressure-corrected hourly neutron counts; USGS earthquake catalogue.
* CR index: Each station is normalised to its own long-term mean count rate
  (so all stations become dimensionless relative-rate series around 1.0).
  The global index is the mean across all stations that have valid data at
  each 5-day bin.
* Seismic metric: sum of Mw for all M ≥ 4.0 events per 5-day bin (Homola
  metric; zero for bins with no events).
* Binning: non-overlapping 5-day bins anchored at study_start.
* Correlation: Pearson r(τ) for τ ∈ [-1000, +1000] days  (step = 5 days).
  Convention: τ > 0 means CR leads seismic (at lag τ the seismic series is
  advanced τ days into the future relative to the CR series).
* Significance: "naive" Pearson t-test p-value (treats bins as i.i.d. — this
  is the inflated significance that Homola reports; reproduced here as the
  baseline before phase-randomised correction).

Outputs
-------
results/figs/homola_replication.png   — r(τ) plot with peak annotated
results/homola_replication.json       — machine-readable summary
"""

from __future__ import annotations

import argparse
import json
import logging
import subprocess
import sys
from datetime import datetime, timezone
from pathlib import Path
from typing import Optional

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
import pandas as pd
import scipy.stats
import yaml

PROJECT_ROOT = Path(__file__).resolve().parent.parent
sys.path.insert(0, str(PROJECT_ROOT / "src"))

from crq.ingest.nmdb import load_station, resample_daily
from crq.ingest.usgs import load_usgs, seismic_energy_per_bin

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)-8s %(name)s — %(message)s",
    datefmt="%Y-%m-%dT%H:%M:%S",
)
logger = logging.getLogger("crq.homola")

# ---------------------------------------------------------------------------
# Constants (matching Homola et al. 2023)
# ---------------------------------------------------------------------------
STUDY_START   = "1976-01-01"
STUDY_END     = "2019-12-31"
BIN_DAYS      = 5
LAG_MIN_DAYS  = -1000
LAG_MAX_DAYS  = +1000
MIN_MAG       = 4.0
SEED          = 42

NMDB_COVERAGE_THRESHOLD = 0.60   # min valid hourly fraction per station-day
MIN_STATIONS_PER_BIN    = 3      # global index → NaN if fewer stations valid


# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------

def _get_git_sha() -> str:
    try:
        r = subprocess.run(
            ["git", "-C", str(PROJECT_ROOT), "rev-parse", "--short", "HEAD"],
            capture_output=True, text=True, timeout=5,
        )
        return r.stdout.strip() if r.returncode == 0 else "unknown"
    except Exception:
        return "unknown"


def _load_station_names(config_path: Path) -> list[str]:
    with config_path.open() as fh:
        return list(yaml.safe_load(fh)["stations"].keys())


def _parse_args() -> argparse.Namespace:
    p = argparse.ArgumentParser(
        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    p.add_argument("--study-start",  default=STUDY_START)
    p.add_argument("--study-end",    default=STUDY_END)
    p.add_argument("--bin-days",     type=int,   default=BIN_DAYS)
    p.add_argument("--lag-min",      type=int,   default=LAG_MIN_DAYS)
    p.add_argument("--lag-max",      type=int,   default=LAG_MAX_DAYS)
    p.add_argument("--min-mag",      type=float, default=MIN_MAG)
    p.add_argument("--min-stations", type=int,   default=MIN_STATIONS_PER_BIN)
    p.add_argument("--nmdb-dir",     type=Path,  default=PROJECT_ROOT / "data" / "raw" / "nmdb")
    p.add_argument("--usgs-dir",     type=Path,  default=PROJECT_ROOT / "data" / "raw" / "usgs")
    p.add_argument("--output-dir",   type=Path,  default=PROJECT_ROOT / "results")
    p.add_argument("--config",       type=Path,  default=PROJECT_ROOT / "config" / "stations.yaml")
    p.add_argument("--seed",         type=int,   default=SEED)
    return p.parse_args()


def _bin_series(
    series: pd.Series,
    study_start: str,
    bin_days: int,
    agg: str = "mean",
) -> pd.Series:
    """
    Resample a daily Series to *bin_days*-day bins anchored exactly at
    *study_start*.

    Uses manual floor-division binning so the anchor date is always respected,
    regardless of pandas version or calendar effects.

    Parameters
    ----------
    agg : ``"mean"`` (default) or ``"sum"``
        Aggregation function.  Use ``"mean"`` for normalised CR rates;
        use ``"sum"`` for the seismic Mw total.
    """
    t0 = pd.Timestamp(study_start)
    days_from_origin = (series.index - t0).days
    bin_number = days_from_origin // bin_days
    bin_dates  = t0 + pd.to_timedelta(bin_number * bin_days, unit="D")
    grouped = series.groupby(bin_dates)
    return grouped.sum() if agg == "sum" else grouped.mean()


# ---------------------------------------------------------------------------
# Step 1 — Global cosmic-ray index
# ---------------------------------------------------------------------------

def build_cr_index(
    station_names:     list[str],
    study_start:       str,
    study_end:         str,
    bin_days:          int,
    nmdb_dir:          Path,
    *,
    coverage_threshold: float = NMDB_COVERAGE_THRESHOLD,
    min_stations:      int   = MIN_STATIONS_PER_BIN,
) -> tuple[pd.Series, pd.DataFrame, int]:
    """
    Load all available NMDB stations, normalise each to its mean count rate
    over the study window, and return a *bin_days*-day-binned global index.

    Returns
    -------
    cr_index          : binned global CR series (NaN if < min_stations valid)
    station_norm_daily: daily normalised series per station (for diagnostics)
    effective_min     : the min_stations threshold actually applied
    """
    t0 = pd.Timestamp(study_start)
    t1 = pd.Timestamp(study_end)
    start_year, end_year = t0.year, t1.year

    norm_cols: dict[str, pd.Series] = {}

    for station in station_names:
        hourly = load_station(station, start_year, end_year, nmdb_dir)
        if hourly.empty:
            continue

        daily_df = resample_daily(hourly, station, coverage_threshold=coverage_threshold)
        daily    = daily_df[station].loc[study_start:study_end]

        n_valid = int(daily.notna().sum())
        if n_valid < 30:
            continue

        station_mean = daily.mean()
        if not np.isfinite(station_mean) or station_mean <= 0:
            logger.warning("station %s: non-positive mean — skipping", station)
            continue

        norm_cols[station] = daily / station_mean   # dimensionless ~1.0
        n_total    = (t1 - t0).days + 1
        coverage   = 100 * n_valid / n_total
        logger.info(
            "station %-6s  valid_days=%5d  coverage=%5.1f%%  mean=%.1f",
            station, n_valid, coverage, station_mean,
        )

    n_loaded = len(norm_cols)
    logger.info("Loaded %d stations with valid data.", n_loaded)

    if n_loaded == 0:
        raise RuntimeError(
            "No valid station data found — run scripts/01_download_data.py first."
        )

    # Auto-adjust min_stations when data is sparse (e.g. subset test)
    effective_min = min_stations
    if n_loaded < min_stations:
        effective_min = max(1, n_loaded - 1) if n_loaded > 1 else 1
        logger.warning(
            "Only %d stations loaded but min_stations=%d requested; "
            "using effective_min=%d.  Results with <%d stations are unreliable.",
            n_loaded, min_stations, effective_min, min_stations,
        )

    station_norm_daily = pd.DataFrame(norm_cols)

    # Global daily index
    n_valid_per_day = station_norm_daily.notna().sum(axis=1)
    global_daily    = station_norm_daily.mean(axis=1)            # NaN-safe
    global_daily[n_valid_per_day < effective_min] = np.nan

    # Bin to bin_days
    cr_index = _bin_series(global_daily, study_start, bin_days)
    cr_index.name = "cr_index"

    n_valid_bins = int(cr_index.notna().sum())
    logger.info(
        "CR index: %d %d-day bins, %d non-NaN (%.1f%%)",
        len(cr_index), bin_days, n_valid_bins,
        100 * n_valid_bins / max(len(cr_index), 1),
    )
    return cr_index, station_norm_daily, effective_min


# ---------------------------------------------------------------------------
# Step 2 — Seismic metric
# ---------------------------------------------------------------------------

def build_seismic_metric(
    study_start: str,
    study_end:   str,
    bin_days:    int,
    usgs_dir:    Path,
    *,
    min_mag:  float = MIN_MAG,
    cr_index: Optional[pd.Series] = None,
) -> pd.Series:
    """
    Sum of Mw for all M ≥ *min_mag* events per *bin_days*-day bin.

    Aligned to *cr_index*'s index if supplied; otherwise built from scratch.
    Bins with no events → 0.0.
    """
    t0 = pd.Timestamp(study_start)
    t1 = pd.Timestamp(study_end)

    events = load_usgs(t0.year, t1.year, usgs_dir)
    if events.empty:
        raise RuntimeError("No USGS data — run scripts/01_download_data.py first.")

    events = events.loc[study_start:study_end]
    events = events[events["mag"] >= min_mag]
    logger.info(
        "Seismic: %d events M≥%.1f in %s–%s",
        len(events), min_mag, study_start, study_end,
    )

    # Physically correct: E = 10^(1.5·Mw + 4.8) J per event, summed per bin,
    # returned as log10(E_sum).  See Kanamori (1977).
    t0 = pd.Timestamp(study_start)
    seismic = seismic_energy_per_bin(
        events, study_start, study_end, bin_days, t0, min_mag=min_mag,
    )
    seismic = seismic.fillna(seismic.min())   # fill rare empty bins with floor

    # Align to CR index bins
    if cr_index is not None:
        seismic = seismic.reindex(cr_index.index, fill_value=seismic.min())

    nonnan_pct = 100 * seismic.notna().mean()
    logger.info(
        "Seismic metric (log10 E): %d bins, %.1f%% non-NaN, "
        "range [%.1f, %.1f] log10(J)",
        len(seismic), nonnan_pct, float(seismic.min()), float(seismic.max()),
    )
    return seismic


# ---------------------------------------------------------------------------
# Step 3 — Cross-correlation
# ---------------------------------------------------------------------------

def pearson_lag_correlation(
    x: pd.Series,
    y: pd.Series,
    lag_bins: np.ndarray,
) -> pd.DataFrame:
    """
    Pearson r(τ) for each lag in *lag_bins* (bin units).

    Convention: τ > 0 means x leads y.
    At lag τ > 0:  x[0..N-τ]  vs  y[τ..N]
    At lag τ < 0:  x[|τ|..N]  vs  y[0..N-|τ|]

    Returns DataFrame: lag_bins, r, p_value, n_pairs.
    """
    x_arr = np.asarray(x, dtype=float)
    y_arr = np.asarray(y, dtype=float)
    N = len(x_arr)

    rs, ps, ns = [], [], []
    for lag in lag_bins:
        if lag >= 0:
            end = N - lag
            if end <= 0:
                rs.append(np.nan); ps.append(np.nan); ns.append(0)
                continue
            xa = x_arr[:end]
            ya = y_arr[lag : lag + end]
        else:
            abslag = -lag
            end = N - abslag
            if end <= 0:
                rs.append(np.nan); ps.append(np.nan); ns.append(0)
                continue
            xa = x_arr[abslag : abslag + end]
            ya = y_arr[:end]

        valid = np.isfinite(xa) & np.isfinite(ya)
        n = int(valid.sum())
        ns.append(n)
        if n < 10:
            rs.append(np.nan); ps.append(np.nan)
            continue

        r, p = scipy.stats.pearsonr(xa[valid], ya[valid])
        rs.append(float(r)); ps.append(float(p))

    return pd.DataFrame({"lag_bins": lag_bins, "r": rs, "p_value": ps, "n_pairs": ns})


def naive_sigma(p_value: float) -> float:
    """Two-tailed p → Gaussian sigma  z = Φ⁻¹(1 − p/2)."""
    if not np.isfinite(p_value) or p_value <= 0:
        return np.inf
    return float(scipy.stats.norm.isf(p_value / 2.0))


# ---------------------------------------------------------------------------
# Step 4 — Figure
# ---------------------------------------------------------------------------

def make_figure(
    lag_df:         pd.DataFrame,
    bin_days:       int,
    peak_lag_days:  float,
    peak_r:         float,
    n_stations:     int,
    study_start:    str,
    study_end:      str,
    naive_p:        float,
    sigma:          float,
    cr_5d:          pd.Series,
    seismic_5d:     pd.Series,
    n_stations_eff: int,
    r_at_15:        float,
    sigma_at_15:    float,
) -> plt.Figure:
    lags_days = lag_df["lag_bins"].values * bin_days
    r_vals    = lag_df["r"].values

    n0 = int(lag_df.loc[lag_df["lag_bins"] == 0, "n_pairs"].iloc[0])
    naive_ci95 = 1.96 / np.sqrt(n0) if n0 > 4 else 0.0

    # Naive r thresholds for σ=3 and σ=6 (two-tailed, N=n0)
    def _r_for_sigma(s: float, n: int) -> float:
        """r value corresponding to naive σ=s for sample size n."""
        p = 2.0 * scipy.stats.norm.sf(s)
        t_crit = scipy.stats.t.isf(p / 2.0, df=n - 2)
        return float(t_crit / np.sqrt(t_crit**2 + n - 2))

    r_sigma3 = _r_for_sigma(3.0, n0) if n0 > 4 else np.nan
    r_sigma6 = _r_for_sigma(6.0, n0) if n0 > 4 else np.nan

    fig, axes = plt.subplots(
        3, 1, figsize=(14, 12),
        gridspec_kw={"height_ratios": [3, 1.4, 1.4]},
    )

    # ── Panel 1: r(τ) ──────────────────────────────────────────────────────
    ax = axes[0]
    ax.axhline(0,           color="k",   lw=0.8)
    ax.axhline( naive_ci95, color="C0",  lw=0.8, ls="--",
                label=f"±{naive_ci95:.3f} naive 95 % CI  (n={n0:,})")
    ax.axhline(-naive_ci95, color="C0",  lw=0.8, ls="--")

    if np.isfinite(r_sigma3):
        ax.axhline( r_sigma3, color="C2", lw=0.9, ls=":", alpha=0.7,
                    label=f"Naive 3σ  (r={r_sigma3:.3f})")
        ax.axhline(-r_sigma3, color="C2", lw=0.9, ls=":", alpha=0.7)
    if np.isfinite(r_sigma6):
        ax.axhline( r_sigma6, color="C1", lw=0.9, ls="-.", alpha=0.8,
                    label=f"Naive 6σ  (r={r_sigma6:.3f})  ← Homola claim")
        ax.axhline(-r_sigma6, color="C1", lw=0.9, ls="-.", alpha=0.8)

    ax.axvline(0,  color="k", lw=0.6, alpha=0.35)

    ax.plot(lags_days, r_vals, color="k", lw=1.1, zorder=3)

    # Mark the τ=+15d bin (Homola claimed signal)
    lag15_bins = 15 // bin_days
    row15 = lag_df[lag_df["lag_bins"] == lag15_bins]
    if len(row15):
        r15 = float(row15["r"].iloc[0])
        if np.isfinite(r15):
            s15_str = f"{sigma_at_15:.1f}σ" if np.isfinite(sigma_at_15) and sigma_at_15 < 100 else ">8σ"
            ax.scatter([15], [r15], color="darkorange", s=60, zorder=6,
                       label=f"τ = +15 d: r = {r15:+.4f} ({s15_str} naive)")
            ax.axvline(15, color="darkorange", lw=1.0, ls=":", alpha=0.6)

    # Mark the overall peak
    ax.axvline(peak_lag_days, color="crimson", lw=1.3, ls="--",
               label=f"Peak τ = {peak_lag_days:+.0f} d,  r = {peak_r:+.4f}")
    ax.scatter([peak_lag_days], [peak_r], color="crimson", s=50, zorder=5)

    r_clean = r_vals[np.isfinite(r_vals)]
    ylo = min(-0.15, r_clean.min() * 1.15) if len(r_clean) else -0.15
    yhi = max( 0.15, r_clean.max() * 1.15) if len(r_clean) else  0.15
    ax.set_xlim(lags_days[0], lags_days[-1])
    ax.set_ylim(ylo, yhi)
    ax.set_xlabel("Lag τ (days)  [τ > 0: CR leads seismic]")
    ax.set_ylabel("Pearson r")
    ax.set_title(
        f"Homola et al. 2023 — replication baseline  |  "
        f"{n_stations} NMDB stations  |  "
        f"{study_start[:4]}–{study_end[:4]}  |  "
        f"5-day bins  |  M≥{MIN_MAG:.1f}",
        fontsize=10,
    )
    ax.legend(fontsize=8.5, loc="upper right", ncol=2)
    ax.grid(True, alpha=0.22)

    sig_str = f"{sigma:.1f}σ" if np.isfinite(sigma) and sigma < 100 else ">8σ"
    info = (
        f"Peak: r = {peak_r:+.5f}  at  τ = {peak_lag_days:+.0f} d  ({sig_str} naive)\n"
        f"τ = +15 d: r = {r_at_15:+.5f}  "
        + (f"({sigma_at_15:.1f}σ naive)" if np.isfinite(sigma_at_15) and sigma_at_15 < 100 else "")
        + f"\nn = {n0:,} valid bin-pairs at τ=0\n"
        f"⚠ Naive — ignores autocorr., solar cycle, lag-scan penalty"
    )
    ax.text(
        0.014, 0.975, info, transform=ax.transAxes,
        va="top", ha="left", fontsize=8.5,
        bbox=dict(boxstyle="round,pad=0.35", fc="lightyellow", alpha=0.92),
    )

    # ── Panel 2: CR index ──────────────────────────────────────────────────
    ax2 = axes[1]
    cr_dev = (cr_5d - 1.0) * 100   # % deviation from station-normalised mean
    ax2.plot(cr_dev.index, cr_dev.values, lw=0.7, color="navy", alpha=0.85)
    ax2.axhline(0, color="k", lw=0.5)
    ax2.set_ylabel("CR deviation (%)")
    ax2.set_title(
        f"Global CR index  ({n_stations} stations, each normalised to their long-term mean)",
        fontsize=9,
    )
    ax2.grid(True, alpha=0.18)
    ax2.set_xlim(cr_5d.index[0], cr_5d.index[-1])

    # ── Panel 3: seismic metric ────────────────────────────────────────────
    ax3 = axes[2]
    ax3.fill_between(seismic_5d.index, seismic_5d.values, color="firebrick", alpha=0.6, lw=0)
    ax3.set_ylabel(f"log₁₀(Σ E)  [log₁₀ J]  (M≥{MIN_MAG:.1f})")
    ax3.set_title(
        f"Global seismic metric  (log₁₀ of summed energy per 5-day bin, "
        f"E=10^(1.5·Mw+4.8) J)",
        fontsize=9,
    )
    ax3.grid(True, alpha=0.18)
    ax3.set_xlim(seismic_5d.index[0], seismic_5d.index[-1])

    fig.tight_layout(pad=1.2)
    return fig


# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------

def main() -> None:
    args = _parse_args()
    np.random.seed(args.seed)
    git_sha = _get_git_sha()

    logger.info("Homola replication  |  git=%s  seed=%d", git_sha, args.seed)
    logger.info("Study window: %s – %s", args.study_start, args.study_end)
    logger.info("Bin: %d days  |  Lags: %+d … %+d days",
                args.bin_days, args.lag_min, args.lag_max)

    args.output_dir.mkdir(parents=True, exist_ok=True)
    figs_dir = args.output_dir / "figs"
    figs_dir.mkdir(parents=True, exist_ok=True)

    # ── 1. CR index
    station_names = _load_station_names(args.config)
    cr_5d, station_norm_daily, effective_min = build_cr_index(
        station_names,
        args.study_start, args.study_end,
        args.bin_days, args.nmdb_dir,
        coverage_threshold=NMDB_COVERAGE_THRESHOLD,
        min_stations=args.min_stations,
    )
    n_stations = station_norm_daily.shape[1]

    # ── 2. Seismic metric
    seismic_5d = build_seismic_metric(
        args.study_start, args.study_end,
        args.bin_days, args.usgs_dir,
        min_mag=args.min_mag,
        cr_index=cr_5d,
    )

    # ── 3. Data quality assessment
    n_cr_bins   = int(cr_5d.notna().sum())
    study_days  = (pd.Timestamp(args.study_end) - pd.Timestamp(args.study_start)).days
    expect_bins = study_days // args.bin_days
    coverage_pct = 100.0 * n_cr_bins / max(expect_bins, 1)

    if coverage_pct < 50:
        logger.warning(
            "⚠  CR data coverage = %.1f%% of study window.  "
            "Results are not meaningful with less than half the data.  "
            "Run: python scripts/01_download_data.py  (full run, no --subset).",
            coverage_pct,
        )

    # ── 4. Cross-correlation
    lag_bins = np.arange(
        args.lag_min // args.bin_days,
        args.lag_max // args.bin_days + 1,
        dtype=int,
    )
    logger.info("Computing Pearson r(τ) for %d lags …", len(lag_bins))
    lag_df = pearson_lag_correlation(cr_5d, seismic_5d, lag_bins)

    # ── 5. Peak statistics
    valid_mask = lag_df["r"].notna()
    if not valid_mask.any():
        raise RuntimeError(
            "All correlation values are NaN — data is insufficient for analysis."
        )

    peak_idx      = lag_df.loc[valid_mask, "r"].abs().idxmax()
    peak_row      = lag_df.loc[peak_idx]
    peak_lag_bins = int(peak_row["lag_bins"])
    peak_lag_days = peak_lag_bins * args.bin_days
    peak_r        = float(peak_row["r"])
    naive_p       = float(peak_row["p_value"])
    n_at_peak     = int(peak_row["n_pairs"])
    sigma         = naive_sigma(naive_p)

    # Correlation specifically at the Homola ±15-day claim
    lag15_bins   = 15 // args.bin_days   # = 3 bins
    row15        = lag_df[lag_df["lag_bins"] == lag15_bins]
    r_at_15      = float(row15["r"].iloc[0]) if len(row15) else np.nan
    p_at_15      = float(row15["p_value"].iloc[0]) if len(row15) else np.nan
    sigma_at_15  = naive_sigma(p_at_15) if np.isfinite(p_at_15) else np.nan

    logger.info("─" * 55)
    logger.info("Peak r = %+.5f  at  τ = %+d days (%+d bins)",
                peak_r, peak_lag_days, peak_lag_bins)
    logger.info("Naive p = %.3e  (%.1fσ)", naive_p,
                sigma if np.isfinite(sigma) else 999)
    logger.info("n pairs at peak = %d", n_at_peak)
    logger.info("r at τ = +15 d   = %+.5f  (naive p = %.2e, %.1fσ)",
                r_at_15, p_at_15 if np.isfinite(p_at_15) else 1,
                sigma_at_15 if np.isfinite(sigma_at_15) else 0)
    logger.info("─" * 55)
    logger.info(
        "NOTE: naive significance ignores autocorrelation, the shared solar-"
        "cycle trend, and scanning over %d lags — expect gross over-significance.",
        len(lag_bins),
    )

    # ── 6. Figure
    fig = make_figure(
        lag_df, args.bin_days, peak_lag_days, peak_r,
        n_stations, args.study_start, args.study_end,
        naive_p, sigma, cr_5d, seismic_5d,
        n_stations_eff=effective_min,
        r_at_15=r_at_15,
        sigma_at_15=sigma_at_15,
    )
    fig_path = figs_dir / "homola_replication.png"
    fig.savefig(fig_path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    logger.info("Saved: %s", fig_path)

    # ── 7. JSON log
    results = {
        "script":               "scripts/02_homola_replication.py",
        "git_sha":              git_sha,
        "seed":                 args.seed,
        "timestamp_utc":        datetime.now(timezone.utc).isoformat(),
        "study_window":         [args.study_start, args.study_end],
        "bin_days":             args.bin_days,
        "lag_range_days":       [args.lag_min, args.lag_max],
        "min_magnitude":        args.min_mag,
        "nmdb_coverage_threshold": NMDB_COVERAGE_THRESHOLD,
        "min_stations_requested": args.min_stations,
        "min_stations_effective": effective_min,
        "n_stations_used":      n_stations,
        "n_cr_bins_total":      int(len(cr_5d)),
        "n_cr_bins_valid":      n_cr_bins,
        "cr_data_coverage_pct": round(coverage_pct, 2),
        "peak_lag_days":        peak_lag_days,
        "peak_lag_bins":        peak_lag_bins,
        "peak_r":               round(peak_r, 6),
        "naive_p_value":        float(f"{naive_p:.6e}"),
        "naive_sigma":          round(sigma, 3) if np.isfinite(sigma) and sigma < 1e6 else None,
        "n_pairs_at_peak":      n_at_peak,
        "r_at_tau_15d":         round(r_at_15,  6) if np.isfinite(r_at_15)      else None,
        "p_at_tau_15d":         float(f"{p_at_15:.6e}") if np.isfinite(p_at_15) else None,
        "sigma_at_tau_15d":     round(sigma_at_15, 3) if np.isfinite(sigma_at_15) else None,
        "n_lags_scanned":       len(lag_bins),
        "note_significance": (
            "Naive p from Pearson t-test treating 5-day bins as i.i.d. "
            "Not corrected for: (1) temporal autocorrelation, "
            "(2) shared solar-cycle trend, (3) scan over multiple lags. "
            "Expected to be grossly over-significant — corrected tests in phase 2."
        ),
    }
    json_path = args.output_dir / "homola_replication.json"
    json_path.write_text(json.dumps(results, indent=2))
    logger.info("Saved: %s", json_path)

    # ── Human summary
    print()
    print("=" * 62)
    print("  HOMOLA ET AL. 2023 — REPLICATION SUMMARY")
    print("=" * 62)
    print(f"  Stations used           : {n_stations}")
    print(f"  CR data coverage        : {coverage_pct:.1f}% of study window")
    print(f"  Valid 5-day bins (CR)   : {n_cr_bins:,}")
    print()
    print(f"  Peak |r|                : {abs(peak_r):.5f}  (r = {peak_r:+.5f})")
    print(f"  Peak lag τ              : {peak_lag_days:+d} days")
    print(f"  Naive p-value at peak   : {naive_p:.3e}")
    sstr = f"{sigma:.2f}σ" if np.isfinite(sigma) and sigma < 1e6 else ">8σ"
    print(f"  Naive sigma at peak     : {sstr}")
    print()
    s15str = f"{sigma_at_15:.2f}σ" if np.isfinite(sigma_at_15) else "N/A"
    print(f"  r at τ = +15 d (claim)  : {r_at_15:+.5f}  ({s15str})")
    print()
    if coverage_pct < 50:
        print("  ⚠  WARNING: data coverage < 50% — results not meaningful.")
        print("             python scripts/01_download_data.py  (full run)")
    else:
        print("  ✓  Data coverage sufficient.")
    print("=" * 62)


if __name__ == "__main__":
    main()
