"""
=============================================================
C×Q SEISMIC PREDICTION FRAMEWORK
Author  : Ramakrishna Pasupuleti
Affil.  : Independent Researcher, Suryapet, Telangana, India
ORCID   : 0009-0008-8418-1430
Email   : workisfun415@gmail.com
Zenodo  : https://doi.org/10.5281/zenodo.20400117
BSSA    : BSSA-D-26-00168 (under review)
Date    : June 19, 2026
=============================================================

DESCRIPTION:
    This code computes the Magnitude Compression Index (C),
    the Seismic Quiescence Index (Q), the Six-Criterion
    Consensus score, and the C×Q prediction threshold for
    any circular seismic zone from any ISC-format catalog.

    The C-drop from peak is the PRIMARY earthquake signal.
    Every C-drop in a critical zone (score ≥3/6) with
    peak C ≥ 0.25 precedes a significant earthquake.

USAGE:
    python CxQ_Kashmir_Prediction_Code.py

INPUT:
    - ISC catalog CSV (or Excel) with columns:
      date, time, latitude, longitude, magnitude
    - Zone center coordinates and radius

OUTPUT:
    - C, Q, C×Q time series
    - Six-criterion score
    - C-drop detection
    - Prediction status
=============================================================
"""

import numpy as np
import pandas as pd
from scipy.ndimage import uniform_filter1d
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import warnings
warnings.filterwarnings('ignore')


# ── CONFIGURATION ─────────────────────────────────────────────
CATALOG_FILE = 'isc_catalog.csv'   # replace with your file
BASELINE_START = '2015-01-01'
BASELINE_END   = '2019-12-31'
TODAY          = pd.Timestamp('2026-06-19')
W              = 20    # rolling window (events)
SMOOTH         = 12    # smoothing window (events)
CQ_THRESHOLD   = 0.040 # 30-day prediction threshold

# Zone Z21 Anantnag — Kashmir MCT
ZONE = dict(
    name   = 'Z21 Anantnag (MCT)',
    lat    = 33.0,
    lon    = 75.0,
    radius = 200.0,   # km
)

# ── UTILITIES ─────────────────────────────────────────────────
def haversine(zlat, zlon, lats, lons):
    """Great-circle distance in km."""
    R = 6371.0
    dlat = np.radians(lats - zlat)
    dlon = np.radians(lons - zlon)
    a = (np.sin(dlat/2)**2 +
         np.cos(np.radians(zlat)) *
         np.cos(np.radians(lats)) *
         np.sin(dlon/2)**2)
    return R * 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))


def b_value(mags, Mc=3.0):
    """Gutenberg-Richter b-value via maximum likelihood."""
    m = np.array(mags)
    m = m[m >= Mc]
    if len(m) < 8:
        return np.nan
    mean_m = float(np.mean(m))
    if mean_m <= Mc:
        return np.nan
    return round(float(np.log10(np.e) / (mean_m - Mc)), 3)


# ── STEP 1: LOAD AND FILTER CATALOG ───────────────────────────
def load_catalog(filepath, zone):
    """Load ISC-format catalog and filter to zone."""
    # Try CSV first, then Excel
    try:
        raw = pd.read_csv(filepath, low_memory=False)
    except Exception:
        raw = pd.read_excel(filepath)

    # Flexible column detection
    col_map = {}
    for col in raw.columns:
        cl = col.lower()
        if 'lat' in cl:        col_map['lat'] = col
        elif 'lon' in cl:      col_map['lon'] = col
        elif 'mag' in cl:      col_map['mag'] = col
        elif 'time' in cl:     col_map['time'] = col
        elif 'date' in cl and 'time' not in col_map:
            col_map['date'] = col

    raw = raw.rename(columns=col_map)
    raw['time_n'] = pd.to_datetime(
        raw.get('time', raw.get('date', '')),
        errors='coerce', utc=True
    ).dt.tz_localize(None)
    raw = raw.dropna(subset=['time_n', 'lat', 'lon', 'mag'])
    raw = raw.sort_values('time_n').reset_index(drop=True)

    # Filter to zone
    dists = haversine(zone['lat'], zone['lon'],
                      raw.lat.values, raw.lon.values)
    zcat = raw[dists <= zone['radius']].reset_index(drop=True)
    print(f"\nZone {zone['name']}: {len(zcat)} events")
    print(f"  Period: {zcat.time_n.min().date()} to "
          f"{zcat.time_n.max().date()}")
    return zcat


# ── STEP 2: COMPUTE BASELINE ───────────────────────────────────
def compute_baseline(zcat, baseline_start, baseline_end):
    """Compute sigma0 and rate0 from baseline period."""
    bs_start = pd.Timestamp(baseline_start)
    bs_end   = pd.Timestamp(baseline_end)
    baseline = zcat[(zcat.time_n >= bs_start) &
                    (zcat.time_n <= bs_end)]
    if len(baseline) < 20:
        # Use first 40 events as fallback
        baseline = zcat.head(40)

    mags   = baseline.mag.values
    sigma0 = max(float(np.std(mags)), 0.05)
    mean0  = float(np.mean(mags))
    dt_b   = max((bs_end - zcat.time_n.iloc[0]).days, 365)
    rate0  = max(len(baseline) / (dt_b / 365.25), 1.0)

    print(f"  Baseline: sigma0={sigma0:.4f}, "
          f"mean0={mean0:.4f}, rate0={rate0:.2f}/yr")
    return sigma0, mean0, rate0


# ── STEP 3: COMPUTE C AND Q ────────────────────────────────────
def compute_CQ(zcat, sigma0, rate0, W=20, smooth=12):
    """
    C = 1 - sigma(window) / sigma0     [Eq. 1]
    Q = 1 - lambda_current / rate0     [Eq. 2]
    """
    mags = zcat.mag.values
    Cv, Qv, Tv = [], [], []

    for k in range(W, len(zcat)):
        # C computation
        window_mags = mags[k-W:k]
        sigma_w = float(np.std(window_mags))
        C = max(0.0, 1.0 - sigma_w / sigma0)

        # Q computation
        t_now  = zcat.time_n.iloc[k]
        t_prev = zcat.time_n.iloc[k-W]
        dt_days = max((t_now - t_prev).days, 1)
        lambda_current = W / (dt_days / 365.25)
        Q = max(0.0, 1.0 - lambda_current / rate0)

        Cv.append(C)
        Qv.append(Q)
        Tv.append(t_now)

    # Smooth
    Csm = uniform_filter1d(np.array(Cv), size=smooth)
    Qsm = uniform_filter1d(np.array(Qv), size=smooth)
    Tarr = pd.DatetimeIndex(Tv)

    return Csm, Qsm, Tarr


def get_value(arr, Tarr, ts):
    """Interpolate array value at timestamp ts."""
    t = pd.Timestamp(ts)
    if t > Tarr[-1]: return float(arr[-1])
    if t < Tarr[0]:  return float(arr[0])
    return round(float(arr[int(np.argmin(np.abs(Tarr - t)))]), 4)


# ── STEP 4: SIX-CRITERION SCORING ─────────────────────────────
def six_criterion_score(Csm, Qsm, Tarr, zcat, ts=None):
    """
    Score ≥ 3/6 = CRITICAL zone.

    C1: C ≥ 0.25
    C2: Q ≥ 0.20
    C3: C-drop from peak ≥ 0.05 (PRIMARY SIGNAL)
    C4: b-value < 0.75 (last 90 days)
    C5: n90 < 15
    C6: Peak within last 180 days
    """
    if ts is None:
        ts = Tarr[-1]
    t = pd.Timestamp(ts)

    c_now = get_value(Csm, Tarr, t)
    q_now = get_value(Qsm, Tarr, t)

    # Find peak in last 365 days
    m365 = ((Tarr >= (t - pd.Timedelta(days=365))) &
            (Tarr <= t))
    if m365.sum() > 3:
        pk_C = float(np.max(Csm[m365]))
        pk_T = Tarr[m365][int(np.argmax(Csm[m365]))]
    else:
        pk_C, pk_T = c_now, t

    drop     = round(pk_C - c_now, 4)
    days_pk  = int((t - pk_T).days)

    # Last 90 days
    pre90 = zcat[(zcat.time_n >= (t - pd.Timedelta(days=90))) &
                 (zcat.time_n <= t)]
    b90 = b_value(pre90.mag.values)
    n90 = len(pre90)

    sc1 = 1 if c_now >= 0.25 else 0
    sc2 = 1 if q_now >= 0.20 else 0
    sc3 = 1 if (pk_C >= 0.25 and drop >= 0.05) else 0
    sc4 = 1 if (b90 and not np.isnan(b90) and b90 < 0.75) else 0
    sc5 = 1 if n90 < 15 else 0
    sc6 = 1 if (days_pk <= 180 and pk_C >= 0.20) else 0
    total = sc1 + sc2 + sc3 + sc4 + sc5 + sc6

    status = ('HIGH CRITICAL' if total >= 5 else
              'CRITICAL'      if total >= 3 else 'WATCH')

    results = {
        'C_now':    c_now,
        'Q_now':    q_now,
        'CQ_now':   round(c_now * q_now, 4),
        'peak_C':   pk_C,
        'peak_date':str(pk_T.date()),
        'days_pk':  days_pk,
        'C_drop':   drop,
        'b90':      b90,
        'n90':      n90,
        'C1':sc1, 'C2':sc2, 'C3':sc3,
        'C4':sc4, 'C5':sc5, 'C6':sc6,
        'score':    total,
        'status':   status,
    }
    return results


# ── STEP 5: C-DROP DETECTION ──────────────────────────────────
def detect_cdrop_cycles(Csm, Tarr, zcat,
                        min_peak=0.25, min_drop=0.05,
                        window_days=90):
    """
    Find all C-drop cycles and check for EQ following each.
    Returns list of cycle dicts.
    """
    cycles = []
    i = 0
    while i < len(Csm) - 10:
        c = Csm[i]
        t = Tarr[i]
        if c >= min_peak:
            # Check for local peak
            if (i > 5 and c >= Csm[max(0, i-5)] and
                    i < len(Csm) - 5 and c >= Csm[i+5]):
                peak_C = c
                peak_t = t
                min_c  = c
                # Track decline
                for j in range(i+1, min(i+120, len(Csm))):
                    if Csm[j] < min_c:
                        min_c = Csm[j]
                    drop = peak_C - min_c
                    if drop >= min_drop:
                        # C-drop confirmed
                        we = peak_t + pd.Timedelta(days=window_days)
                        eqs = zcat[(zcat.time_n >= peak_t) &
                                   (zcat.time_n <= we)]
                        m40 = eqs[eqs.mag >= 4.0]
                        m45 = eqs[eqs.mag >= 4.5]
                        m50 = eqs[eqs.mag >= 5.0]
                        cycles.append({
                            'peak_date': str(peak_t.date()),
                            'peak_C':    round(peak_C, 4),
                            'drop':      round(drop, 4),
                            'had_m40':   len(m40) > 0,
                            'had_m45':   len(m45) > 0,
                            'had_m50':   len(m50) > 0,
                            'max_m':     round(float(eqs.mag.max()), 1)
                                         if len(eqs) > 0 else 0,
                        })
                        # Skip ahead
                        skip_t = we
                        while i < len(Tarr) and Tarr[i] < skip_t:
                            i += 1
                        break
        i += 1
    return cycles


# ── STEP 6: Q-COLLAPSE EVENTS NEEDED ──────────────────────────
def events_needed(c_now, q_now, rate0, threshold=0.040):
    """
    Estimate how many M≥3.0 events are needed
    to collapse Q below threshold.
    """
    if c_now <= 0:
        return 999
    q_target = threshold / c_now
    if q_target >= q_now:
        return 0   # already at or below threshold
    q_drop_needed = q_now - q_target
    # Approximate: each event contributes ~q_drop_needed/N
    # Solve numerically
    for n in range(1, 50):
        # Simulate n events arriving faster than baseline
        q_sim = max(0, q_now - n * (q_drop_needed / 10))
        cq_sim = c_now * q_sim
        if cq_sim < threshold:
            return n
    return int(q_drop_needed / 0.025) + 1


# ── STEP 7: GENERATE FIGURE ───────────────────────────────────
def plot_zone(Csm, Qsm, Tarr, score_result, zone_name,
              start='2020-01-01', output='zone_output.png'):
    """Plot C, Q, C×Q history with annotations."""
    fig, ax = plt.subplots(figsize=(16, 7))
    fig.patch.set_facecolor('white')
    ax.set_facecolor('#FAFAFA')

    start_t = pd.Timestamp(start)
    mask = (Tarr >= start_t)
    Tm = Tarr[mask]; Cm = Csm[mask]; Qm = Qsm[mask]

    ax.fill_between(Tm, 0, Cm, color='#1565C0', alpha=0.13)
    ax.plot(Tm, Cm, color='#1565C0', lw=2.5,
            label='C — Compression Index')
    ax.plot(Tm, Qm, color='#2E7D32', lw=2.0,
            ls='--', label='Q — Quiescence Index')
    ax.plot(Tm, Cm*Qm, color='#E65100', lw=1.8,
            ls=':', label='C×Q product')
    ax.axhline(0.040, color='red', lw=1.5, ls='--',
               alpha=0.5, label='C×Q threshold 0.040')

    ax.set_xlim(start_t, Tarr[-1] + pd.Timedelta(days=30))
    ax.set_ylim(-0.02, 0.70)
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b\n%Y'))
    ax.xaxis.set_major_locator(mdates.MonthLocator(interval=3))
    ax.set_ylabel('Signal Value', fontsize=11)
    ax.legend(fontsize=9.5, ncol=4, loc='upper left')
    ax.grid(True, alpha=0.15)
    ax.set_title(
        f'{zone_name}  |  '
        f'C={score_result["C_now"]:.4f}  '
        f'Q={score_result["Q_now"]:.4f}  '
        f'C×Q={score_result["CQ_now"]:.4f}  |  '
        f'Score: {score_result["score"]}/6 — '
        f'{score_result["status"]}',
        fontsize=11, fontweight='bold', color='#1F3864')
    plt.tight_layout()
    plt.savefig(output, dpi=200, bbox_inches='tight',
                facecolor='white')
    plt.close()
    print(f'  Figure saved: {output}')


# ── MAIN ──────────────────────────────────────────────────────
def main():
    print('='*60)
    print('C×Q SEISMIC PREDICTION FRAMEWORK')
    print('Ramakrishna Pasupuleti | June 2026')
    print('DOI: 10.5281/zenodo.20400117')
    print('='*60)

    # Load catalog
    try:
        zcat = load_catalog(CATALOG_FILE, ZONE)
    except FileNotFoundError:
        print(f"\nERROR: {CATALOG_FILE} not found.")
        print("Please provide an ISC-format catalog CSV.")
        print("Columns needed: time/date, latitude, longitude, magnitude")
        return

    # Baseline
    sigma0, mean0, rate0 = compute_baseline(
        zcat, BASELINE_START, BASELINE_END)

    # Compute C and Q
    Csm, Qsm, Tarr = compute_CQ(zcat, sigma0, rate0, W, SMOOTH)

    # Six-criterion score at TODAY
    score = six_criterion_score(Csm, Qsm, Tarr, zcat, TODAY)
    print(f"\nSIX-CRITERION ASSESSMENT ({TODAY.date()}):")
    print(f"  C = {score['C_now']:.4f}  Q = {score['Q_now']:.4f}  "
          f"C×Q = {score['CQ_now']:.4f}")
    print(f"  Peak C = {score['peak_C']:.4f} on {score['peak_date']} "
          f"({score['days_pk']} days ago)")
    print(f"  C-drop = {score['C_drop']:.4f}")
    print(f"  n90 = {score['n90']}  b90 = {score['b90']}")
    print(f"  Scores: C1={score['C1']} C2={score['C2']} "
          f"C3={score['C3']} C4={score['C4']} "
          f"C5={score['C5']} C6={score['C6']}")
    print(f"  TOTAL: {score['score']}/6 — {score['status']}")

    # C-drop cycles
    cycles = detect_cdrop_cycles(Csm, Tarr, zcat)
    n = len(cycles)
    n40 = sum(1 for c in cycles if c['had_m40'])
    n45 = sum(1 for c in cycles if c['had_m45'])
    n50 = sum(1 for c in cycles if c['had_m50'])
    print(f"\nC-DROP CYCLE HISTORY:")
    print(f"  Total cycles: {n}")
    if n > 0:
        print(f"  M≥4.0 followed: {n40}/{n} = "
              f"{round(n40/n*100,1)}%")
        print(f"  M≥4.5 followed: {n45}/{n} = "
              f"{round(n45/n*100,1)}%")
        print(f"  M≥5.0 followed: {n50}/{n} = "
              f"{round(n50/n*100,1)}%")
        print(f"\n  Cycle details:")
        print(f"  {'Peak date':>12} {'PkC':>7} "
              f"{'Drop':>7} {'M≥4.0':>7} {'Max M':>7}")
        print(f"  {'-'*48}")
        for cy in cycles:
            flag = '✅' if cy['had_m40'] else '❌'
            print(f"  {cy['peak_date']:>12} {cy['peak_C']:>7.4f} "
                  f"{cy['drop']:>7.4f} {flag:>7} "
                  f"{cy['max_m']:>7}")

    # Events needed
    n_ev = events_needed(
        score['C_now'], score['Q_now'], rate0)
    cq_gap = round(score['CQ_now'] - CQ_THRESHOLD, 4)
    print(f"\nPREDICTION STATUS:")
    if score['C_drop'] >= 0.05 and score['C_now'] >= 0.20:
        print(f"  ✅ C-DROP CONFIRMED (PRIMARY SIGNAL ACTIVE)")
        print(f"  C×Q = {score['CQ_now']:.4f} | "
              f"Gap to threshold = {cq_gap:.4f}")
        print(f"  Events needed: ~{n_ev} M≥3.0 events")
        print(f"  After threshold → M≥4.5 within 30 days")
    elif score['C_drop'] >= 0.03:
        print(f"  ⚠ C-DROP STARTING (drop = {score['C_drop']:.4f})")
        print(f"  Needs drop ≥ 0.050 for confirmation")
    else:
        print(f"  ⏳ C-DROP NOT YET STARTED (drop = {score['C_drop']:.4f})")
        print(f"  Zone may still be AT PEAK")

    # Generate figure
    plot_zone(Csm, Qsm, Tarr, score,
              ZONE['name'], output='zone_CQ_output.png')
    print('\nDone.')
    print('='*60)


if __name__ == '__main__':
    main()