---
name: "Declustering Algorithm Sensitivity of CONUS Hazard Estimates"
description: "Apply Gardner-Knopoff, Reasenberg, and Zaliapin-Ben-Zion declustering to the same ANSS ComCat CONUS catalog, propagate Gutenberg-Richter parameter changes through a simplified PSHA at 500 sites, and quantify the algorithm-induced spread in 2%-in-50yr PGA against within-algorithm bootstrap variability."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "seismology", "psha", "declustering", "gutenberg-richter", "hazard", "sensitivity-analysis"]
python_version: ">=3.8"
dependencies: []
---

# Declustering Algorithm Sensitivity of CONUS Hazard Estimates

## When to Use This Skill

Use this skill when you need to quantify how the choice of an earthquake declustering algorithm (Gardner-Knopoff window, Reasenberg link-based, Zaliapin-Ben-Zion nearest-neighbor) shifts the 2%-in-50-year peak ground acceleration across a CONUS site grid, and to compare that algorithm-induced spread with the within-algorithm bootstrap variability at fixed data.

## Research Question

Given identical input data (the ANSS ComCat CONUS catalog 1990-2023, M ≥ 3.5), how large is the 2%-in-50-year PGA spread induced purely by the choice of declustering algorithm, relative to the within-algorithm bootstrap variability at fixed data? A ratio ≳ 1× implies algorithm choice matters as much as catalog sampling noise; ≳ 5× implies algorithm choice dominates.

## Prerequisites

**Software**
- Python **3.8+** (checked at startup; the script exits with code 2 and an explicit message on older interpreters).
- Standard library only. **No `pip install` is required or permitted.**
- POSIX-like shell for the `mkdir`/`cat << SCRIPT_EOF` steps.

**Network / data**
- Outbound HTTPS access to `earthquake.usgs.gov` (port 443) for the initial ANSS ComCat download via the `fdsnws/event/1/query` endpoint.
- Subsequent runs reuse the local SHA256-verified cache under `data/` and do not require network.

**Filesystem**
- Write access to `/tmp/claw4s_auto_declustering-algorithm-sensitivity-of-conus-hazard-estimates/` (or any workspace of your choice — the script uses `os.path.dirname(__file__)` for the workspace and auto-creates `data/`).
- ~3 MB of disk for the cached catalog chunks + `results.json` + `report.md`.

**Runtime budget**
- ~6–12 minutes on a single CPU after the initial download (dominated by the Zaliapin-Ben-Zion O(N²) nearest-neighbor pass with `WINDOW=5000` and the 100-replicate bootstrap).
- Initial catalog download typically 30–90 seconds over commodity broadband.

**Environment variables**
- None required. All domain-specific values are UPPER_CASE constants in the `DOMAIN CONFIGURATION` block at the top of the analysis script.

## Adaptation Guidance

To adapt this analysis to a different region or a different declustering comparison, modify only the constants inside the `DOMAIN CONFIGURATION` block in the analysis script:

- `BBOX_MIN_LAT`, `BBOX_MAX_LAT`, `BBOX_MIN_LON`, `BBOX_MAX_LON` — change the catalog bounding box (e.g., California-only, or a stable continental region).
- `CATALOG_START`, `CATALOG_END`, `MIN_MAG_FETCH` — change the query window. Keep `MIN_MAG_FETCH <= MC_ASSUMED`.
- `MC_ASSUMED` — completeness magnitude used when fitting Gutenberg-Richter `a, b`.
- `SITE_GRID_STEP_DEG` — adjust the site spacing to change the number of hazard sites (default 500 sites from a ~1.3° grid over CONUS).
- `GMPE_COEFFS` — swap in an alternate ground-motion prediction equation. The only requirement is that `gmpe_log10_pga(M, R)` returns median `log10(PGA_g)` and `GMPE_SIGMA_LOG10` gives the aleatory standard deviation in log10 units.
- `DECL_*` blocks — tune algorithm parameters (window scaling, `eta_0` threshold, interaction radius).

The general statistical method (`run_analysis()`) is domain-agnostic: it computes declustered catalogs, fits GR parameters, runs a Cornell-McGuire hazard integration, and contrasts algorithm spread with bootstrap spread. `load_data()` contains the only domain-specific logic (USGS ComCat query plus CSV parsing). Replace `load_data()` with any function returning a list of `{year, lat, lon, mag}` dicts to analyze a different catalog source.

## Overview

The analysis pipeline is:

1. Download the ANSS ComCat CONUS catalog (M >= `MIN_MAG_FETCH`, 1990-2023) via paginated CSV queries.
2. Apply three declustering algorithms to the same catalog:
   - **Gardner-Knopoff (1974)** — space/time windows scaled by magnitude.
   - **Reasenberg (1985)** — simplified link-based clustering with magnitude-dependent interaction radius.
   - **Zaliapin-Ben-Zion (2013)** — nearest-neighbor distance in rescaled space-time-magnitude coordinates, with threshold `eta_0`.
3. Fit Gutenberg-Richter `a, b` via Aki (1965) MLE for each declustered catalog above `MC_ASSUMED`.
4. Build a uniform areal source zone from the declustered catalog, and compute hazard curves at 500 CONUS sites via Cornell-McGuire integration using a simplified Boore-Joyner-Fumal style GMPE.
5. Extract the 2%-in-50-year PGA at each site for each algorithm and compare the algorithm-induced distribution against a 200-replicate bootstrap (resampling the declustered catalog for one reference algorithm).
6. Sensitivity checks: completeness magnitude `Mc ∈ {3.5, 4.0, 4.5}`, restricted post-1990 catalog, and an alternate GMPE.

## Step 1: Create the workspace

Run the command below to create the workspace and its `data/` subdirectory.

```bash
mkdir -p /tmp/claw4s_auto_declustering-algorithm-sensitivity-of-conus-hazard-estimates/data
```

**Expected output:**
- Exit code 0, no stdout.
- Directory `/tmp/claw4s_auto_declustering-algorithm-sensitivity-of-conus-hazard-estimates/` exists.
- Empty subdirectory `data/` exists.

## Step 2: Write the analysis script

Run the heredoc below to materialize the Python analysis script. It writes `analyze.py` into the workspace.

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_declustering-algorithm-sensitivity-of-conus-hazard-estimates/analyze.py
#!/usr/bin/env python3
"""
Declustering Algorithm Sensitivity of CONUS Hazard Estimates.

Apply three declustering algorithms (Gardner-Knopoff, Reasenberg,
Zaliapin-Ben-Zion) to the same ANSS ComCat CONUS catalog and compare
their effect on 2%-in-50-year PGA estimates at 500 sites.

Stdlib only, Python 3.8+.
"""

import sys

# ── Python version guard (part of graceful-failure contract) ──────
# Must come before any 3.8+ syntax (f-strings are fine in 3.8+; we avoid
# them in this guard itself so the message prints on 2.x/3.6 as well).
if sys.version_info < (3, 8):
    sys.stderr.write(
        "ERROR: Python >= 3.8 required (found {}). "
        "Re-run with a newer interpreter.\n".format(
            ".".join(str(x) for x in sys.version_info[:3])))
    sys.exit(2)

import os
import csv as csvmod
import json
import math
import time
import hashlib
import random
import shutil
import urllib.request
import urllib.error
from collections import defaultdict

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════
SEED = 42

WORKSPACE = os.path.dirname(os.path.abspath(__file__))
DATA_DIR = os.path.join(WORKSPACE, "data")

USGS_API = "https://earthquake.usgs.gov/fdsnws/event/1/query"

# CONUS bounding box (degrees)
BBOX_MIN_LAT = 24.5
BBOX_MAX_LAT = 49.5
BBOX_MIN_LON = -125.0
BBOX_MAX_LON = -66.5

# Catalog window (inclusive / exclusive)
CATALOG_START = "1990-01-01"
CATALOG_END = "2024-01-01"

# Minimum fetch magnitude (keep small; completeness handled later)
MIN_MAG_FETCH = 3.5

# Assumed completeness magnitude for GR fitting
MC_ASSUMED = 4.0

# Magnitude-frequency bins
M_BIN = 0.1

# Gutenberg-Richter truncation for hazard
M_MAX_HAZARD = 7.5

# Site grid: step in degrees. ~1.7° step over CONUS ≈ 510 sites.
SITE_GRID_STEP_DEG = 1.7

# PSHA integration parameters
PSHA_M_STEP = 0.1
PSHA_INVESTIGATION_YEARS = 50.0
PSHA_EXCEEDANCE_PROB = 0.02
# Return period implied by 2% in 50 years:
PSHA_RETURN_PERIOD = -PSHA_INVESTIGATION_YEARS / math.log(1.0 - PSHA_EXCEEDANCE_PROB)

# Simplified GMPE (Boore-Joyner-Fumal 1997 style, rock site, PGA):
#   log10(PGA_g) = C1 + C2*(M-6) + C3*(M-6)**2 - C4*log10(R + H)
GMPE_COEFFS = {"C1": -1.02, "C2": 0.229, "C3": 0.00, "C4": 0.778, "H": 5.57}
GMPE_SIGMA_LOG10 = 0.226  # aleatory sigma in log10 units (~0.52 natural log)

# Alternate GMPE for sensitivity (harder attenuation)
GMPE_COEFFS_ALT = {"C1": -1.05, "C2": 0.229, "C3": 0.00, "C4": 0.95, "H": 5.57}

# Declustering parameters
# Gardner-Knopoff (1974) windows
def gk_space_km(M):
    return 10.0 ** (0.1238 * M + 0.983)

def gk_time_days(M):
    if M >= 6.5:
        return 10.0 ** (0.032 * M + 2.7389)
    return 10.0 ** (0.5409 * M - 0.547)

# Reasenberg (1985) simplified parameters
DECL_R_TAU_MIN_DAYS = 1.0
DECL_R_TAU_MAX_DAYS = 10.0
DECL_R_P1 = 0.95   # confidence of next event in cluster
DECL_R_XK = 0.5    # factor used to modify magnitude of largest cluster event
DECL_R_XMEFF = 3.5 # effective catalog magnitude threshold
DECL_R_RFACT = 10  # factor for interaction radius
DECL_R_SIGMA = 0.4  # exponent in interaction radius (km = 10^(sigma*M - off))

# Zaliapin-Ben-Zion parameters
DECL_ZBZ_D = 1.6    # fractal dimension
DECL_ZBZ_B_REF = 1.0  # reference b for rescaling
DECL_ZBZ_LOG_ETA0 = -5.0  # default log10 threshold separating background/clustered

# Hazard PGA levels (g) for curve computation. Low end extended to 0.001 g
# so that low-seismicity CONUS sites do not pin at the tabulated floor.
PGA_LEVELS = [0.001, 0.002, 0.003, 0.005, 0.0075, 0.010, 0.015, 0.020,
              0.030, 0.050, 0.075, 0.100, 0.150, 0.200, 0.300, 0.400,
              0.600, 0.800, 1.000, 1.500, 2.000]

# Bootstrap / resampling
N_BOOTSTRAP = 100

# Output filenames
RESULTS_JSON = "results.json"
REPORT_MD = "report.md"

# ═══════════════════════════════════════════════════════════════
# END DOMAIN CONFIGURATION
# ═══════════════════════════════════════════════════════════════


# ────────────────────────────────────────────────────────────────
# Math utilities
# ────────────────────────────────────────────────────────────────

def normal_cdf(z):
    """Standard normal CDF via Abramowitz-Stegun approximation."""
    if z < -8.0:
        return 0.0
    if z > 8.0:
        return 1.0
    a1 = 0.254829592
    a2 = -0.284496736
    a3 = 1.421413741
    a4 = -1.453152027
    a5 = 1.061405429
    p = 0.3275911
    sign = 1 if z >= 0 else -1
    t = 1.0 / (1.0 + p * abs(z))
    y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * math.exp(-z * z / 2.0)
    return 0.5 * (1.0 + sign * y)


def pct_val(sorted_data, p):
    """p-th percentile (0..100) of a sorted list via linear interpolation."""
    if not sorted_data:
        return 0.0
    n = len(sorted_data)
    k = (n - 1) * p / 100.0
    f = int(math.floor(k))
    c = min(f + 1, n - 1)
    if f == c:
        return sorted_data[f]
    return sorted_data[f] * (c - k) + sorted_data[c] * (k - f)


def haversine_km(lat1, lon1, lat2, lon2):
    """Great-circle distance in km."""
    R = 6371.0088
    phi1 = math.radians(lat1)
    phi2 = math.radians(lat2)
    dphi = math.radians(lat2 - lat1)
    dlam = math.radians(lon2 - lon1)
    a = math.sin(dphi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(dlam / 2) ** 2
    return 2 * R * math.asin(math.sqrt(min(1.0, a)))


def days_between(t1_days, t2_days):
    return abs(t2_days - t1_days)


def iso_to_days(iso_str):
    """Parse an ISO timestamp into fractional days since 1970-01-01."""
    # Accepts 'YYYY-MM-DDTHH:MM:SS.sssZ' (USGS) or similar.
    try:
        year = int(iso_str[0:4])
        mon = int(iso_str[5:7])
        day = int(iso_str[8:10])
        hh = int(iso_str[11:13]) if len(iso_str) >= 13 else 0
        mm = int(iso_str[14:16]) if len(iso_str) >= 16 else 0
        ss_part = iso_str[17:19] if len(iso_str) >= 19 else "0"
        ss = float(ss_part) if ss_part.isdigit() else 0.0
        if len(iso_str) > 19 and iso_str[19] == ".":
            frac_end = 20
            while frac_end < len(iso_str) and iso_str[frac_end].isdigit():
                frac_end += 1
            try:
                ss += float(iso_str[19:frac_end])
            except ValueError:
                pass
    except (ValueError, IndexError):
        return None
    # Julian day number for 1970-01-01 is 2440588
    # Use a simple epoch-days formula via datetime-free math:
    # Zeller-like for date to day-count
    a = (14 - mon) // 12
    y = year + 4800 - a
    m = mon + 12 * a - 3
    jdn = day + (153 * m + 2) // 5 + 365 * y + y // 4 - y // 100 + y // 400 - 32045
    epoch_jdn = 2440588  # 1970-01-01
    days = (jdn - epoch_jdn) + (hh * 3600 + mm * 60 + ss) / 86400.0
    return days


# ────────────────────────────────────────────────────────────────
# SHA256 helpers
# ────────────────────────────────────────────────────────────────

def sha256_bytes(data):
    return hashlib.sha256(data).hexdigest()


def sha256_file(path):
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(1 << 16), b""):
            h.update(chunk)
    return h.hexdigest()


# ────────────────────────────────────────────────────────────────
# USGS ComCat download with caching
# ────────────────────────────────────────────────────────────────

def _year_chunks(start_iso, end_iso, chunk_years):
    """Yield (chunk_start, chunk_end) iso strings partitioning [start, end)."""
    y0 = int(start_iso[0:4])
    y1 = int(end_iso[0:4])
    y = y0
    while y < y1:
        yn = min(y + chunk_years, y1)
        yield ("{:04d}-01-01".format(y), "{:04d}-01-01".format(yn))
        y = yn


def download_chunk(chunk_start, chunk_end, retries=3):
    """Download one catalog chunk from USGS ComCat with caching and SHA256."""
    os.makedirs(DATA_DIR, exist_ok=True)
    label = "{}_{}".format(chunk_start[0:4], chunk_end[0:4])
    cache_path = os.path.join(DATA_DIR, "conus_{}.csv".format(label))
    hash_path = cache_path + ".sha256"

    if os.path.exists(cache_path) and os.path.exists(hash_path):
        with open(hash_path) as hf:
            stored = hf.read().strip()
        if sha256_file(cache_path) == stored:
            size = os.path.getsize(cache_path)
            print("    Cached {}: {:,} bytes (sha256={}...)".format(
                label, size, stored[:12]))
            return cache_path

    url = (
        "{api}?format=csv&starttime={s}&endtime={e}"
        "&minmagnitude={m}&orderby=time-asc"
        "&minlatitude={minlat}&maxlatitude={maxlat}"
        "&minlongitude={minlon}&maxlongitude={maxlon}"
    ).format(
        api=USGS_API, s=chunk_start, e=chunk_end, m=MIN_MAG_FETCH,
        minlat=BBOX_MIN_LAT, maxlat=BBOX_MAX_LAT,
        minlon=BBOX_MIN_LON, maxlon=BBOX_MAX_LON,
    )

    for attempt in range(retries):
        try:
            print("    Downloading {} (attempt {}/{})...".format(
                label, attempt + 1, retries))
            req = urllib.request.Request(url)
            req.add_header("User-Agent", "Claw4S-Decluster/1.0")
            with urllib.request.urlopen(req, timeout=180) as resp:
                raw = resp.read()
            if len(raw) < 200 or b"time" not in raw[:200]:
                raise RuntimeError("Unexpected short response")
            with open(cache_path, "wb") as f:
                f.write(raw)
            h = sha256_bytes(raw)
            with open(hash_path, "w") as f:
                f.write(h)
            lines = raw.count(b"\n")
            print("    OK {}: {} lines, sha256={}...".format(label, lines, h[:12]))
            return cache_path
        except Exception as e:
            print("    Attempt {} failed: {}".format(attempt + 1, e))
            if attempt < retries - 1:
                wait = 5 * (attempt + 1)
                print("    Retrying in {}s...".format(wait))
                time.sleep(wait)

    raise RuntimeError(
        "Failed to download chunk {} after {} attempts. Check connectivity.".format(
            label, retries))


def parse_catalog_csv(path):
    """Parse USGS ComCat CSV into a list of event dicts via stdlib csv.

    Uses Python's csv module so that quoted-comma fields (e.g., 'place')
    are parsed correctly without affecting the time/lat/lon/mag columns.
    """
    events = []
    with open(path, "r", encoding="utf-8", errors="replace", newline="") as f:
        reader = csvmod.DictReader(f)
        for row in reader:
            try:
                t = row.get("time", "")
                lat = float(row.get("latitude", ""))
                lon = float(row.get("longitude", ""))
                mag_raw = row.get("mag", "")
                if not mag_raw:
                    continue
                mag = float(mag_raw)
                days = iso_to_days(t)
                if days is None:
                    continue
                events.append({
                    "t_days": days,
                    "year": int(t[0:4]),
                    "lat": lat,
                    "lon": lon,
                    "mag": mag,
                })
            except (ValueError, KeyError, TypeError):
                continue
    return events


def load_data():
    """Download and parse the ANSS ComCat CONUS catalog.

    Returns (events, provenance) where provenance is a list of per-chunk
    dicts (window, file path, byte size, SHA-256, parsed row count).
    """
    os.makedirs(DATA_DIR, exist_ok=True)
    all_events = []
    provenance = []
    for s, e in _year_chunks(CATALOG_START, CATALOG_END, 5):
        p = download_chunk(s, e)
        rows = parse_catalog_csv(p)
        all_events.extend(rows)
        size = os.path.getsize(p)
        sha = sha256_file(p)
        provenance.append({
            "chunk_start": s, "chunk_end": e,
            "path": p, "bytes": size, "sha256": sha,
            "rows_parsed": len(rows),
        })

    all_events.sort(key=lambda r: r["t_days"])
    print("  Parsed {:,} events from {} chunk file(s).".format(
        len(all_events), len(provenance)))
    return all_events, provenance


# ────────────────────────────────────────────────────────────────
# Declustering: Gardner-Knopoff (1974)
# ────────────────────────────────────────────────────────────────

def decluster_gardner_knopoff(events):
    """Gardner-Knopoff (1974) space-time window declustering.

    An event is flagged as dependent (aftershock/foreshock) if a larger-magnitude
    event falls within its space-time window. Returns list of mainshock indices.
    """
    n = len(events)
    # Sort by magnitude descending so we process largest first
    order = sorted(range(n), key=lambda i: -events[i]["mag"])
    is_mainshock = [True] * n
    # Pre-compute windows per event magnitude
    windows = [(gk_space_km(events[i]["mag"]), gk_time_days(events[i]["mag"]))
               for i in range(n)]

    # Build 1° spatial grid index for quick neighbor lookup
    bin_lat = 0.5
    bin_lon = 0.5
    grid = defaultdict(list)
    for i, ev in enumerate(events):
        key = (int(math.floor(ev["lat"] / bin_lat)),
               int(math.floor(ev["lon"] / bin_lon)))
        grid[key].append(i)

    for i in order:
        if not is_mainshock[i]:
            continue
        ev_i = events[i]
        r_km, t_days = windows[i]
        # Bound search grid by r_km
        lat_span = max(1, int(math.ceil(r_km / 111.0 / bin_lat)) + 1)
        lon_span = max(1, int(math.ceil(
            r_km / (111.0 * max(0.1, math.cos(math.radians(ev_i["lat"])))) / bin_lon)) + 1)
        key_i = (int(math.floor(ev_i["lat"] / bin_lat)),
                 int(math.floor(ev_i["lon"] / bin_lon)))
        for dla in range(-lat_span, lat_span + 1):
            for dlo in range(-lon_span, lon_span + 1):
                for j in grid.get((key_i[0] + dla, key_i[1] + dlo), ()):
                    if j == i or not is_mainshock[j]:
                        continue
                    if events[j]["mag"] > ev_i["mag"] + 1e-9:
                        continue  # j is larger; skip
                    if days_between(ev_i["t_days"], events[j]["t_days"]) > t_days:
                        continue
                    d = haversine_km(ev_i["lat"], ev_i["lon"],
                                     events[j]["lat"], events[j]["lon"])
                    if d <= r_km:
                        is_mainshock[j] = False

    return [i for i, m in enumerate(is_mainshock) if m]


# ────────────────────────────────────────────────────────────────
# Declustering: Reasenberg (1985) — simplified link-based
# ────────────────────────────────────────────────────────────────

def _reasenberg_interaction_radius_km(M):
    """Interaction radius from Wells-Coppersmith-like scaling, km.

    Uses rupture-length-style scaling R = 10^(0.5*M - 1.85) km.
    """
    return max(0.1, 10.0 ** (0.5 * M - 1.85))


def _reasenberg_look_days(days_since_start_of_cluster, largest_mag_so_far):
    """Taper interaction time with cluster age.

    Simplified: tau(M, t) = min(tau_max, tau_min * (1 + t_cluster_age_days / 10))
    clipped between DECL_R_TAU_MIN_DAYS and DECL_R_TAU_MAX_DAYS, but scaled by
    magnitude factor tied to largest event magnitude in cluster.
    """
    base = DECL_R_TAU_MIN_DAYS * (1.0 + days_since_start_of_cluster / 10.0)
    # Boost for large mainshocks:
    mag_boost = 10.0 ** (DECL_R_XK * (largest_mag_so_far - DECL_R_XMEFF))
    tau = base * max(1.0, mag_boost)
    return min(DECL_R_TAU_MAX_DAYS * mag_boost, max(DECL_R_TAU_MIN_DAYS, tau))


def decluster_reasenberg(events):
    """Simplified Reasenberg (1985) link-based declustering.

    Two events link if they are within an interaction radius and interaction
    time. Clusters are connected components; keep the largest-magnitude event
    per cluster.
    """
    n = len(events)
    parent = list(range(n))

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    def union(a, b):
        ra, rb = find(a), find(b)
        if ra != rb:
            parent[ra] = rb

    # Events sorted by time
    # Track cluster start time and largest magnitude via a per-root map updated lazily.
    cluster_start = {}
    cluster_max_mag = {}

    # Spatial grid for pruning
    bin_lat = 1.0
    bin_lon = 1.0
    grid = defaultdict(list)
    for i, ev in enumerate(events):
        key = (int(math.floor(ev["lat"] / bin_lat)),
               int(math.floor(ev["lon"] / bin_lon)))
        grid[key].append(i)

    for j in range(n):
        ev_j = events[j]
        # Search for older links
        # Interaction radius from ev_j's magnitude:
        r_j = _reasenberg_interaction_radius_km(ev_j["mag"]) * DECL_R_RFACT
        lat_span = max(1, int(math.ceil(r_j / 111.0 / bin_lat)) + 1)
        lon_span = max(1, int(math.ceil(
            r_j / (111.0 * max(0.1, math.cos(math.radians(ev_j["lat"])))) / bin_lon)) + 1)
        key_j = (int(math.floor(ev_j["lat"] / bin_lat)),
                 int(math.floor(ev_j["lon"] / bin_lon)))
        for dla in range(-lat_span, lat_span + 1):
            for dlo in range(-lon_span, lon_span + 1):
                for i in grid.get((key_j[0] + dla, key_j[1] + dlo), ()):
                    if i >= j:
                        continue
                    ev_i = events[i]
                    dt = ev_j["t_days"] - ev_i["t_days"]
                    if dt <= 0:
                        continue
                    root_i = find(i)
                    start_i = cluster_start.get(root_i, ev_i["t_days"])
                    max_mag_i = cluster_max_mag.get(root_i, ev_i["mag"])
                    tau = _reasenberg_look_days(ev_i["t_days"] - start_i, max_mag_i)
                    if dt > tau:
                        continue
                    r_eff = max(r_j, _reasenberg_interaction_radius_km(max_mag_i) * DECL_R_RFACT)
                    d = haversine_km(ev_i["lat"], ev_i["lon"],
                                     ev_j["lat"], ev_j["lon"])
                    if d <= r_eff:
                        union(i, j)
        # Update cluster stats for j's component
        root = find(j)
        cluster_start.setdefault(root, ev_j["t_days"])
        prev = cluster_max_mag.get(root, -99.0)
        if ev_j["mag"] > prev:
            cluster_max_mag[root] = ev_j["mag"]

    # Group by cluster root, retain largest magnitude event
    by_root = defaultdict(list)
    for i in range(n):
        by_root[find(i)].append(i)
    mainshocks = []
    for root, members in by_root.items():
        if len(members) == 1:
            mainshocks.append(members[0])
        else:
            # Pick largest magnitude; tie-break earliest time
            best = members[0]
            for i in members[1:]:
                if (events[i]["mag"] > events[best]["mag"] + 1e-9 or
                    (abs(events[i]["mag"] - events[best]["mag"]) <= 1e-9 and
                     events[i]["t_days"] < events[best]["t_days"])):
                    best = i
            mainshocks.append(best)
    mainshocks.sort()
    return mainshocks


# ────────────────────────────────────────────────────────────────
# Declustering: Zaliapin-Ben-Zion (2013) — nearest-neighbor in rescaled space
# ────────────────────────────────────────────────────────────────

def _zbz_rescaled_nn(events, d_frac, b_ref):
    """For each event i>0, find the parent j<i minimizing eta_ij.

        T_ij = (t_i - t_j)          [years]
        R_ij = (haversine_km)^d     [km^d]
        eta_ij = T_ij * R_ij * 10^(-b_ref * mag_j)

    Returns (log10_eta, parent_idx, max_observed_lag) where max_observed_lag
    is the largest chosen-parent index lag (j - i_best); this allows checking
    whether the WINDOW cap ever binds on real data.
    """
    n = len(events)
    t_years = [ev["t_days"] / 365.25 for ev in events]
    lats = [ev["lat"] for ev in events]
    lons = [ev["lon"] for ev in events]
    mags = [ev["mag"] for ev in events]

    log10_eta = [None] * n
    parent_idx = [-1] * n

    # O(N^2) sliding window limited by WINDOW prior events for feasibility.
    # Zaliapin-Ben-Zion practice: the true nearest neighbor is almost always
    # within a few hundred earlier events by time for physical catalogs.
    # We record the largest chosen-parent lag observed so callers can verify
    # that the cap is not binding on this catalog.
    WINDOW = 5000
    max_lag = 0
    n_at_cap = 0

    for j in range(n):
        lo = max(0, j - WINDOW)
        best = None
        best_i = -1
        for i in range(lo, j):
            dt = t_years[j] - t_years[i]
            if dt <= 0:
                continue
            r_km = haversine_km(lats[i], lons[i], lats[j], lons[j])
            if r_km < 0.05:
                r_km = 0.05  # floor to avoid log(0) from duplicate locations
            # eta = T * R^d * 10^(-b*mag_parent)
            log_eta = (math.log10(dt) + d_frac * math.log10(r_km)
                       - b_ref * mags[i])
            if best is None or log_eta < best:
                best = log_eta
                best_i = i
        log10_eta[j] = best
        parent_idx[j] = best_i
        if best_i >= 0:
            lag = j - best_i
            if lag > max_lag:
                max_lag = lag
            if lag >= WINDOW:
                n_at_cap += 1
    return log10_eta, parent_idx, {"window": WINDOW, "max_lag": max_lag,
                                     "n_at_cap": n_at_cap}


def decluster_zbz(events, log_eta0=DECL_ZBZ_LOG_ETA0,
                   d_frac=DECL_ZBZ_D, b_ref=DECL_ZBZ_B_REF):
    """Zaliapin-Ben-Zion nearest-neighbor declustering.

    Events with log10(eta_NN) < log_eta0 are classified as clustered/aftershocks;
    others are retained as background. Events with no prior event (index 0) are
    retained as background.

    Returns (background_indices, log10_eta, cap_diagnostics) where
    cap_diagnostics records max parent lag observed and whether the window
    bound was reached.
    """
    log10_eta, parent, cap = _zbz_rescaled_nn(events, d_frac, b_ref)
    background = []
    for i, le in enumerate(log10_eta):
        if le is None:
            background.append(i)
        elif le >= log_eta0:
            background.append(i)
    return background, log10_eta, cap


# ────────────────────────────────────────────────────────────────
# Gutenberg-Richter MLE (Aki 1965)
# ────────────────────────────────────────────────────────────────

def fit_gr(magnitudes, mc, dm=M_BIN):
    """Aki (1965) MLE for b, with Bender (1983) bin correction.

    b = log10(e) / (mean(M) - Mc + dm/2)
    a = log10(N_above_Mc / T_years) + b * Mc   (rate form — caller supplies T_years)

    Returns (b, se_b, n_above).
    """
    above = [m for m in magnitudes if m >= mc - 1e-9]
    n = len(above)
    if n < 20:
        return None, None, n
    m_mean = sum(above) / n
    denom = m_mean - mc + dm / 2.0
    if denom <= 0:
        return None, None, n
    b = math.log10(math.e) / denom
    se = b / math.sqrt(n)
    return b, se, n


def gr_rate_above_mc(magnitudes, mc, duration_years):
    """Annual rate of M >= mc events."""
    above = [m for m in magnitudes if m >= mc - 1e-9]
    if duration_years <= 0:
        return 0.0
    return len(above) / duration_years


# ────────────────────────────────────────────────────────────────
# PSHA: simplified Cornell-McGuire integration
# ────────────────────────────────────────────────────────────────

def gmpe_log10_pga(M, R_km, coeffs=None):
    """Simplified Boore-Joyner-Fumal style GMPE. Returns median log10(PGA_g)."""
    c = coeffs if coeffs is not None else GMPE_COEFFS
    r = math.sqrt(R_km * R_km + c["H"] * c["H"])
    return c["C1"] + c["C2"] * (M - 6.0) + c["C3"] * (M - 6.0) ** 2 - c["C4"] * math.log10(r)


def gmpe_p_exceed(M, R_km, pga_g, coeffs=None, sigma=GMPE_SIGMA_LOG10):
    """P(PGA > pga_g | M, R) using lognormal aleatory uncertainty (log10 sigma)."""
    mu = gmpe_log10_pga(M, R_km, coeffs)
    if pga_g <= 0:
        return 1.0
    z = (math.log10(pga_g) - mu) / sigma
    return 1.0 - normal_cdf(z)


def build_source_grid(events_sub, cell_deg=1.0):
    """Bin events into a spatial source grid; return cell centers and counts."""
    cells = defaultdict(int)
    for ev in events_sub:
        key = (int(math.floor(ev["lat"] / cell_deg)),
               int(math.floor(ev["lon"] / cell_deg)))
        cells[key] += 1
    # Cell center approximated at (key + 0.5) * cell_deg
    centers = []
    counts = []
    for key, n in cells.items():
        clat = (key[0] + 0.5) * cell_deg
        clon = (key[1] + 0.5) * cell_deg
        centers.append((clat, clon))
        counts.append(n)
    return centers, counts


def build_sites(step_deg):
    """Uniform grid of sites over CONUS, trimmed to a rough land box."""
    sites = []
    lat = BBOX_MIN_LAT + step_deg / 2.0
    while lat < BBOX_MAX_LAT:
        lon = BBOX_MIN_LON + step_deg / 2.0
        while lon < BBOX_MAX_LON:
            sites.append((lat, lon))
            lon += step_deg
        lat += step_deg
    return sites


def hazard_curve_at_site(site_lat, site_lon, centers, cell_rates_mc,
                          b_value, mc, m_max, coeffs=None):
    """Compute annual rate of exceedance for each PGA level at one site.

    cell_rates_mc[i] = annual rate of M >= mc earthquakes in cell i.
    Distribute rate over magnitudes per truncated GR (Mc..m_max).
    """
    beta = b_value * math.log(10.0)
    # Precompute magnitude bins and their relative weights (truncated GR density)
    n_mbins = max(1, int(round((m_max - mc) / PSHA_M_STEP)))
    m_centers = [mc + (k + 0.5) * PSHA_M_STEP for k in range(n_mbins)]
    # Truncated exponential pdf weights (normalize so sum = 1)
    raw_w = [math.exp(-beta * (m - mc)) for m in m_centers]
    tot = sum(raw_w)
    m_weights = [w / tot for w in raw_w]

    # Distances from site to each cell center
    lam = [0.0] * len(PGA_LEVELS)
    for (clat, clon), rate_mc in zip(centers, cell_rates_mc):
        if rate_mc <= 0:
            continue
        R = haversine_km(site_lat, site_lon, clat, clon)
        R = max(1.0, R)
        for mk, wk in zip(m_centers, m_weights):
            rate_k = rate_mc * wk
            if rate_k <= 0:
                continue
            mu = gmpe_log10_pga(mk, R, coeffs)
            for ip, pga in enumerate(PGA_LEVELS):
                z = (math.log10(pga) - mu) / GMPE_SIGMA_LOG10
                p_exc = 1.0 - normal_cdf(z)
                if p_exc > 0:
                    lam[ip] += rate_k * p_exc
    return lam


def pga_at_target_rate(lam, target_rate):
    """Linearly interpolate PGA (in log10-PGA, log-rate) at target annual rate.

    If hazard curve doesn't reach target rate, returns the max level with a flag.
    """
    # Convert to log10 for interpolation (rate vs log10(pga))
    log_rate = []
    log_pga = []
    for lv, pga in zip(lam, PGA_LEVELS):
        if lv > 0:
            log_rate.append(math.log10(lv))
            log_pga.append(math.log10(pga))
    if not log_rate:
        return 0.0, "no_rate"
    target_log = math.log10(target_rate)
    # Rates decrease with increasing PGA; find interval containing target_log
    for i in range(len(log_rate) - 1):
        r0, r1 = log_rate[i], log_rate[i + 1]
        if r0 >= target_log >= r1:
            # Linear interpolation
            if r0 == r1:
                out = log_pga[i]
            else:
                frac = (target_log - r0) / (r1 - r0)
                out = log_pga[i] + frac * (log_pga[i + 1] - log_pga[i])
            return 10.0 ** out, "ok"
    if target_log > log_rate[0]:
        return PGA_LEVELS[0], "below_range"
    return PGA_LEVELS[-1], "above_range"


# ────────────────────────────────────────────────────────────────
# Helper: run full hazard pipeline for one declustered event list
# ────────────────────────────────────────────────────────────────

def run_hazard_for_catalog(events_sub, duration_years, sites, mc, m_max, coeffs=None):
    """Return per-site 2%-in-50yr PGA and summary stats.

    events_sub: list of event dicts (already declustered), used ONLY to build
                rates; GR parameters are estimated from events_sub[mag>=mc].
    """
    mags = [ev["mag"] for ev in events_sub]
    b, se_b, n_above = fit_gr(mags, mc)
    if b is None:
        return None
    # Source grid rates at Mc:
    sub_above = [ev for ev in events_sub if ev["mag"] >= mc - 1e-9]
    centers, counts = build_source_grid(sub_above, cell_deg=1.0)
    cell_rates_mc = [c / duration_years for c in counts]

    per_site_pga = []
    per_site_status = []
    for (slat, slon) in sites:
        lam = hazard_curve_at_site(slat, slon, centers, cell_rates_mc,
                                   b, mc, m_max, coeffs)
        target_rate = 1.0 / PSHA_RETURN_PERIOD
        pga, status = pga_at_target_rate(lam, target_rate)
        per_site_pga.append(pga)
        per_site_status.append(status)
    return {
        "b": round(b, 4),
        "se_b": round(se_b, 4),
        "n_above_mc": n_above,
        "rate_above_mc_per_year": round(n_above / duration_years, 4),
        "per_site_pga_g": per_site_pga,
        "per_site_status": per_site_status,
    }


def summarize_pga(per_site_pga):
    if not per_site_pga:
        return {}
    s = sorted(per_site_pga)
    return {
        "median_g": round(pct_val(s, 50), 5),
        "mean_g": round(sum(s) / len(s), 5),
        "p05_g": round(pct_val(s, 5), 5),
        "p25_g": round(pct_val(s, 25), 5),
        "p75_g": round(pct_val(s, 75), 5),
        "p95_g": round(pct_val(s, 95), 5),
        "n_sites": len(per_site_pga),
    }


# ────────────────────────────────────────────────────────────────
# Analysis orchestration
# ────────────────────────────────────────────────────────────────

def run_analysis(events):
    """Domain-agnostic statistical pipeline: declustering → GR → PSHA spread."""
    rng = random.Random(SEED)
    n_tot = len(events)
    t0 = events[0]["t_days"]
    t1 = events[-1]["t_days"]
    duration_years = (t1 - t0) / 365.25
    print("  Catalog span: {:.2f} years, {:,} events (M>={}).".format(
        duration_years, n_tot, MIN_MAG_FETCH))

    # ── Declustering ────────────────────────────────────────────
    print("")
    print("  [Declustering GK]  Gardner-Knopoff (1974)...")
    t_start = time.time()
    gk_idx = decluster_gardner_knopoff(events)
    gk_events = [events[i] for i in gk_idx]
    print("    Retained {:,} / {:,} ({:.2%}) mainshocks in {:.1f}s.".format(
        len(gk_events), n_tot, len(gk_events) / n_tot, time.time() - t_start))

    print("  [Declustering R]   Reasenberg (1985, simplified)...")
    t_start = time.time()
    r_idx = decluster_reasenberg(events)
    r_events = [events[i] for i in r_idx]
    print("    Retained {:,} / {:,} ({:.2%}) cluster-max events in {:.1f}s.".format(
        len(r_events), n_tot, len(r_events) / n_tot, time.time() - t_start))

    print("  [Declustering ZBZ] Zaliapin-Ben-Zion (2013)...")
    t_start = time.time()
    zbz_idx, log_etas, zbz_cap = decluster_zbz(events)
    zbz_events = [events[i] for i in zbz_idx]
    valid_etas = [x for x in log_etas if x is not None]
    print("    Retained {:,} / {:,} ({:.2%}) background events in {:.1f}s.".format(
        len(zbz_events), n_tot, len(zbz_events) / n_tot, time.time() - t_start))
    print("    log10(eta_NN): median={:.2f}, p5={:.2f}, p95={:.2f}".format(
        pct_val(sorted(valid_etas), 50) if valid_etas else 0,
        pct_val(sorted(valid_etas), 5) if valid_etas else 0,
        pct_val(sorted(valid_etas), 95) if valid_etas else 0))
    print("    NN window cap diagnostic: WINDOW={}, max observed parent "
          "lag={} events, events at cap={}".format(
              zbz_cap["window"], zbz_cap["max_lag"], zbz_cap["n_at_cap"]))

    decl = {
        "GK":  {"events": gk_events, "idx": gk_idx, "label": "Gardner-Knopoff"},
        "R":   {"events": r_events,  "idx": r_idx,  "label": "Reasenberg"},
        "ZBZ": {"events": zbz_events, "idx": zbz_idx, "label": "Zaliapin-Ben-Zion"},
    }

    # ── Sites ────────────────────────────────────────────────────
    sites = build_sites(SITE_GRID_STEP_DEG)
    print("")
    print("  Site grid: {} sites ({:.1f}° spacing).".format(
        len(sites), SITE_GRID_STEP_DEG))

    # ── Hazard at 500 sites for each declustered catalog ─────────
    print("")
    print("  [PSHA] 2%-in-50yr target rate = 1/{:.0f} yr.".format(PSHA_RETURN_PERIOD))
    hazard_by_alg = {}
    for key in ("GK", "R", "ZBZ"):
        print("  [PSHA] Running for {}...".format(decl[key]["label"]))
        t_start = time.time()
        h = run_hazard_for_catalog(decl[key]["events"], duration_years,
                                   sites, MC_ASSUMED, M_MAX_HAZARD)
        if h is None:
            raise RuntimeError(
                "GR fit failed for algorithm {}: fewer than 20 events above "
                "Mc={}. Lower Mc or widen the catalog.".format(key, MC_ASSUMED))
        print("    b={:.4f}, rate(M>={})={:.3f}/yr, median PGA={:.4f}g, "
              "elapsed={:.1f}s.".format(
                  h["b"], MC_ASSUMED, h["rate_above_mc_per_year"],
                  pct_val(sorted(h["per_site_pga_g"]), 50), time.time() - t_start))
        hazard_by_alg[key] = h

    # ── Site-wise algorithm spread ───────────────────────────────
    sites_spread = []
    for s_idx in range(len(sites)):
        per_alg = {k: hazard_by_alg[k]["per_site_pga_g"][s_idx] for k in hazard_by_alg}
        vals = list(per_alg.values())
        lo, hi = min(vals), max(vals)
        mid = sum(vals) / len(vals)
        rel = (hi - lo) / mid if mid > 0 else 0
        sites_spread.append({
            "lat": sites[s_idx][0], "lon": sites[s_idx][1],
            "GK":  per_alg["GK"], "R":  per_alg["R"], "ZBZ": per_alg["ZBZ"],
            "min": lo, "max": hi, "mean": mid, "rel_range": rel,
        })
    rel_ranges = sorted(s["rel_range"] for s in sites_spread)
    spread_stats = {
        "median_rel_range": round(pct_val(rel_ranges, 50), 4),
        "mean_rel_range": round(sum(rel_ranges) / len(rel_ranges), 4),
        "p25_rel_range": round(pct_val(rel_ranges, 25), 4),
        "p75_rel_range": round(pct_val(rel_ranges, 75), 4),
        "p95_rel_range": round(pct_val(rel_ranges, 95), 4),
    }
    print("")
    print("  Algorithm spread across {} sites: median |max-min|/mean = {:.2%}, "
          "p95 = {:.2%}.".format(len(sites), spread_stats["median_rel_range"],
                                 spread_stats["p95_rel_range"]))

    # ── Bootstrap variability for reference algorithm (GK) ───────
    print("")
    print("  [Bootstrap] {} replicates of GK-declustered catalog "
          "(resample events with replacement)...".format(N_BOOTSTRAP))
    rng_boot = random.Random(SEED + 7)
    gk_src = decl["GK"]["events"]
    n_gk = len(gk_src)
    # To save time: bootstrap on a 100-site subsample
    boot_site_indices = sorted(rng_boot.sample(
        range(len(sites)), min(100, len(sites))))
    boot_sites = [sites[i] for i in boot_site_indices]
    boot_pga_matrix = [[] for _ in range(len(boot_sites))]
    t_start = time.time()
    for b_rep in range(N_BOOTSTRAP):
        sample = [gk_src[int(rng_boot.random() * n_gk)] for _ in range(n_gk)]
        h = run_hazard_for_catalog(sample, duration_years, boot_sites,
                                   MC_ASSUMED, M_MAX_HAZARD)
        for k, v in enumerate(h["per_site_pga_g"]):
            boot_pga_matrix[k].append(v)
        if (b_rep + 1) % 50 == 0:
            print("    {}/{} bootstrap reps, elapsed {:.1f}s".format(
                b_rep + 1, N_BOOTSTRAP, time.time() - t_start))
    # Per-site CI width on GK bootstrap, compare to algorithm range at same sites
    boot_widths = []
    alg_ranges = []
    for k, site_idx in enumerate(boot_site_indices):
        sorted_vals = sorted(boot_pga_matrix[k])
        ci_lo = pct_val(sorted_vals, 2.5)
        ci_hi = pct_val(sorted_vals, 97.5)
        if ci_lo > 0:
            rel_ci = (ci_hi - ci_lo) / ((ci_lo + ci_hi) / 2.0)
        else:
            rel_ci = 0.0
        boot_widths.append(rel_ci)
        alg_ranges.append(sites_spread[site_idx]["rel_range"])

    boot_widths_sorted = sorted(boot_widths)
    alg_ranges_sorted = sorted(alg_ranges)
    boot_vs_alg = {
        "n_sites_boot": len(boot_sites),
        "bootstrap_rel_95ci_median": round(pct_val(boot_widths_sorted, 50), 4),
        "bootstrap_rel_95ci_p95": round(pct_val(boot_widths_sorted, 95), 4),
        "algorithm_rel_range_median_at_boot_sites": round(
            pct_val(alg_ranges_sorted, 50), 4),
        "algorithm_rel_range_p95_at_boot_sites": round(
            pct_val(alg_ranges_sorted, 95), 4),
    }
    ratio_ba = (boot_vs_alg["algorithm_rel_range_median_at_boot_sites"] /
                max(1e-9, boot_vs_alg["bootstrap_rel_95ci_median"]))
    print("    bootstrap 95% CI width (median, relative): {:.2%}".format(
        boot_vs_alg["bootstrap_rel_95ci_median"]))
    print("    algorithm range (median, relative):        {:.2%}".format(
        boot_vs_alg["algorithm_rel_range_median_at_boot_sites"]))
    print("    Ratio (algorithm / bootstrap) = {:.2f}×".format(ratio_ba))

    # ── Sensitivity: Mc sweep ────────────────────────────────────
    print("")
    print("  [Sensitivity] Mc sweep {3.5, 4.0, 4.5}...")
    mc_sens = {}
    for mc_try in (3.5, 4.0, 4.5):
        row = {}
        for key in ("GK", "R", "ZBZ"):
            ev_sub = decl[key]["events"]
            h = run_hazard_for_catalog(ev_sub, duration_years, sites,
                                       mc_try, M_MAX_HAZARD)
            if h is None:
                row[key] = None
                continue
            row[key] = {
                "b": h["b"],
                "median_pga_g": round(pct_val(sorted(h["per_site_pga_g"]), 50), 5),
                "mean_pga_g": round(sum(h["per_site_pga_g"])/len(h["per_site_pga_g"]), 5),
            }
        # Spread across algorithms at this Mc
        vals = [row[k]["median_pga_g"] for k in row if row[k] is not None]
        if vals and min(vals) > 0:
            row["alg_rel_range_of_medians"] = round((max(vals)-min(vals))/((max(vals)+min(vals))/2.0), 4)
        mc_sens["Mc={:.1f}".format(mc_try)] = row
        print("    Mc={:.1f}: {}".format(mc_try, {k: row[k] for k in ("GK","R","ZBZ") if row[k] is not None}))

    # ── Sensitivity: alternate GMPE ──────────────────────────────
    print("  [Sensitivity] Alternate GMPE...")
    alt_hazard = {}
    for key in ("GK", "R", "ZBZ"):
        h = run_hazard_for_catalog(decl[key]["events"], duration_years,
                                   sites, MC_ASSUMED, M_MAX_HAZARD,
                                   coeffs=GMPE_COEFFS_ALT)
        alt_hazard[key] = summarize_pga(h["per_site_pga_g"])
    print("    alt-GMPE medians: {}".format(
        {k: alt_hazard[k].get("median_g") for k in alt_hazard}))
    alt_vals = [alt_hazard[k]["median_g"] for k in alt_hazard]
    alt_rel_range = ((max(alt_vals) - min(alt_vals)) /
                     max(1e-9, sum(alt_vals) / len(alt_vals)))

    # ── Sensitivity: restricted era ──────────────────────────────
    print("  [Sensitivity] Restricted era 1995-2023...")
    ev_late = [e for e in events if e["year"] >= 1995]
    late_dur = (ev_late[-1]["t_days"] - ev_late[0]["t_days"]) / 365.25 if ev_late else 0.0
    late_hazard = {}
    if len(ev_late) >= 500:
        gk2 = [ev_late[i] for i in decluster_gardner_knopoff(ev_late)]
        r2 = [ev_late[i] for i in decluster_reasenberg(ev_late)]
        zbz2_idx, _, _ = decluster_zbz(ev_late)
        zbz2 = [ev_late[i] for i in zbz2_idx]
        for key, sub in (("GK", gk2), ("R", r2), ("ZBZ", zbz2)):
            h = run_hazard_for_catalog(sub, late_dur, sites, MC_ASSUMED, M_MAX_HAZARD)
            late_hazard[key] = summarize_pga(h["per_site_pga_g"]) if h else None
        print("    late-era medians: {}".format(
            {k: late_hazard[k].get("median_g") if late_hazard[k] else None
             for k in late_hazard}))
    late_vals = [late_hazard[k]["median_g"] for k in late_hazard
                 if late_hazard.get(k) is not None]
    late_rel_range = ((max(late_vals) - min(late_vals)) /
                      max(1e-9, sum(late_vals) / len(late_vals))) if late_vals else None

    # ── Assemble results ────────────────────────────────────────
    results = {
        "metadata": {
            "title": "Declustering Algorithm Sensitivity of CONUS Hazard Estimates",
            "data_source": "USGS/ANSS ComCat (earthquake.usgs.gov)",
            "catalog_window": [CATALOG_START, CATALOG_END],
            "bbox": [BBOX_MIN_LAT, BBOX_MAX_LAT, BBOX_MIN_LON, BBOX_MAX_LON],
            "min_mag_fetch": MIN_MAG_FETCH,
            "mc_assumed": MC_ASSUMED,
            "m_max_hazard": M_MAX_HAZARD,
            "site_grid_step_deg": SITE_GRID_STEP_DEG,
            "n_sites": len(sites),
            "n_events_total": n_tot,
            "duration_years": round(duration_years, 2),
            "psha_return_period_yr": round(PSHA_RETURN_PERIOD, 1),
            "gmpe_ref": GMPE_COEFFS,
            "gmpe_sigma_log10": GMPE_SIGMA_LOG10,
            "gmpe_alt": GMPE_COEFFS_ALT,
            "seed": SEED,
            "n_bootstrap": N_BOOTSTRAP,
        },
        "declustering": {
            "GK":  {"label": "Gardner-Knopoff (1974)",
                    "n_retained": len(gk_events),
                    "fraction_retained": round(len(gk_events) / n_tot, 4)},
            "R":   {"label": "Reasenberg (1985, simplified)",
                    "n_retained": len(r_events),
                    "fraction_retained": round(len(r_events) / n_tot, 4)},
            "ZBZ": {"label": "Zaliapin-Ben-Zion (2013)",
                    "n_retained": len(zbz_events),
                    "fraction_retained": round(len(zbz_events) / n_tot, 4),
                    "log_eta0_threshold": DECL_ZBZ_LOG_ETA0,
                    "nn_cap_diagnostic": zbz_cap},
        },
        "gr_parameters": {
            k: {"b": hazard_by_alg[k]["b"],
                "se_b": hazard_by_alg[k]["se_b"],
                "n_above_mc": hazard_by_alg[k]["n_above_mc"],
                "rate_above_mc_per_year": hazard_by_alg[k]["rate_above_mc_per_year"]}
            for k in hazard_by_alg
        },
        "hazard_summary": {
            k: {
                **summarize_pga(hazard_by_alg[k]["per_site_pga_g"]),
                "fraction_below_range": round(sum(
                    1 for s in hazard_by_alg[k]["per_site_status"]
                    if s == "below_range") / max(1, len(
                        hazard_by_alg[k]["per_site_status"])), 4),
            }
            for k in hazard_by_alg
        },
        "site_grid_spread": spread_stats,
        "bootstrap_vs_algorithm": {
            **boot_vs_alg,
            "ratio_alg_to_boot": round(ratio_ba, 3),
        },
        "sensitivity_mc": mc_sens,
        "sensitivity_alt_gmpe": {
            "per_algorithm": alt_hazard,
            "alg_rel_range_of_medians": round(alt_rel_range, 4),
        },
        "sensitivity_late_era": {
            "per_algorithm": late_hazard,
            "alg_rel_range_of_medians": (
                round(late_rel_range, 4) if late_rel_range is not None else None),
            "n_events_late": len(ev_late),
            "duration_years_late": round(late_dur, 2),
        },
        "site_details_head": sites_spread[:20],
    }
    return results


# ────────────────────────────────────────────────────────────────
# Report generation
# ────────────────────────────────────────────────────────────────

def generate_report(results):
    path_json = os.path.join(WORKSPACE, RESULTS_JSON)
    path_md = os.path.join(WORKSPACE, REPORT_MD)

    def clean(o):
        if isinstance(o, float):
            if math.isnan(o) or math.isinf(o):
                return str(o)
            return o
        if isinstance(o, dict):
            return {str(k): clean(v) for k, v in o.items()}
        if isinstance(o, (list, tuple)):
            return [clean(v) for v in o]
        return o

    with open(path_json, "w") as f:
        json.dump(clean(results), f, indent=2)
    print("  Wrote {}".format(path_json))

    meta = results["metadata"]
    decl = results["declustering"]
    gr = results["gr_parameters"]
    haz = results["hazard_summary"]
    sp = results["site_grid_spread"]
    bva = results["bootstrap_vs_algorithm"]
    mc_sens = results["sensitivity_mc"]
    alt_sens = results["sensitivity_alt_gmpe"]
    late_sens = results["sensitivity_late_era"]

    L = []
    L.append("# Declustering Algorithm Sensitivity of CONUS Hazard Estimates")
    L.append("")
    L.append("## Summary")
    L.append("")
    L.append("Catalog: {}–{}, M >= {}, CONUS bbox [{:.1f}, {:.1f}] × [{:.1f}, {:.1f}].".format(
        meta["catalog_window"][0][:10], meta["catalog_window"][1][:10],
        meta["min_mag_fetch"],
        meta["bbox"][0], meta["bbox"][1], meta["bbox"][2], meta["bbox"][3]))
    L.append("{:,} events over {:.2f} years.".format(
        meta["n_events_total"], meta["duration_years"]))
    L.append("")
    L.append("### Fraction retained by each declustering algorithm")
    L.append("")
    L.append("| Algorithm | N retained | Fraction |")
    L.append("|-----------|-----------:|---------:|")
    for k in ("GK", "R", "ZBZ"):
        d = decl[k]
        L.append("| {} | {:,} | {:.2%} |".format(
            d["label"], d["n_retained"], d["fraction_retained"]))
    L.append("")

    L.append("### Gutenberg-Richter parameters (Mc = {:.1f})".format(meta["mc_assumed"]))
    L.append("")
    L.append("| Algorithm | b | SE(b) | N(M>=Mc) | rate(M>=Mc)/yr |")
    L.append("|-----------|-----:|-----:|--------:|---------------:|")
    for k in ("GK", "R", "ZBZ"):
        g = gr[k]
        L.append("| {} | {:.4f} | {:.4f} | {:,} | {:.3f} |".format(
            decl[k]["label"], g["b"], g["se_b"],
            g["n_above_mc"], g["rate_above_mc_per_year"]))
    L.append("")

    L.append("### 2%-in-50-yr PGA across {} sites".format(meta["n_sites"]))
    L.append("")
    L.append("| Algorithm | median (g) | mean (g) | p05 | p25 | p75 | p95 |")
    L.append("|-----------|----------:|--------:|----:|----:|----:|----:|")
    for k in ("GK", "R", "ZBZ"):
        h = haz[k]
        L.append("| {} | {:.4f} | {:.4f} | {:.4f} | {:.4f} | {:.4f} | {:.4f} |".format(
            decl[k]["label"], h["median_g"], h["mean_g"],
            h["p05_g"], h["p25_g"], h["p75_g"], h["p95_g"]))
    L.append("")

    L.append("### Site-wise algorithm spread")
    L.append("")
    L.append("For each site, |max − min| across the three algorithms divided by their mean:")
    L.append("")
    L.append("| Quantile | Relative range |")
    L.append("|----------|---------------:|")
    for q in ("p25_rel_range", "median_rel_range", "p75_rel_range", "p95_rel_range"):
        L.append("| {} | {:.2%} |".format(q, sp[q]))
    L.append("")

    L.append("### Bootstrap (within-GK) vs algorithm spread")
    L.append("")
    L.append("| Source of variability | Relative width (median) | Relative width (p95) |")
    L.append("|-----------------------|----------------------:|--------------------:|")
    L.append("| Bootstrap 95% CI      | {:.2%}                | {:.2%}              |".format(
        bva["bootstrap_rel_95ci_median"], bva["bootstrap_rel_95ci_p95"]))
    L.append("| Algorithm min-max     | {:.2%}                | {:.2%}              |".format(
        bva["algorithm_rel_range_median_at_boot_sites"],
        bva["algorithm_rel_range_p95_at_boot_sites"]))
    L.append("")
    L.append("Ratio (algorithm/median bootstrap) = **{:.2f}×**.".format(
        bva["ratio_alg_to_boot"]))
    L.append("")

    L.append("### Sensitivity: completeness magnitude Mc")
    L.append("")
    L.append("| Mc | GK median (g) | R median (g) | ZBZ median (g) | Alg relative range |")
    L.append("|----|-------------:|-------------:|--------------:|------------------:|")
    for mc_key in sorted(mc_sens):
        row = mc_sens[mc_key]
        def fmt(k):
            if row.get(k) is None:
                return "—"
            return "{:.4f}".format(row[k]["median_pga_g"])
        arr = row.get("alg_rel_range_of_medians", None)
        arr_s = "—" if arr is None else "{:.2%}".format(arr)
        L.append("| {} | {} | {} | {} | {} |".format(
            mc_key, fmt("GK"), fmt("R"), fmt("ZBZ"), arr_s))
    L.append("")

    L.append("### Sensitivity: alternate GMPE (harder attenuation)")
    L.append("")
    L.append("| Algorithm | median (g) | mean (g) |")
    L.append("|-----------|----------:|--------:|")
    for k in ("GK", "R", "ZBZ"):
        h = alt_sens["per_algorithm"][k]
        L.append("| {} | {:.4f} | {:.4f} |".format(decl[k]["label"],
                 h["median_g"], h["mean_g"]))
    L.append("")
    L.append("Algorithm relative range (alt GMPE) = **{:.2%}**.".format(
        alt_sens["alg_rel_range_of_medians"]))
    L.append("")

    L.append("### Sensitivity: restricted era (1995–2023)")
    L.append("")
    if late_sens["alg_rel_range_of_medians"] is not None:
        L.append("| Algorithm | median (g) |")
        L.append("|-----------|----------:|")
        for k in ("GK", "R", "ZBZ"):
            v = late_sens["per_algorithm"].get(k)
            if v is None:
                L.append("| {} | — |".format(decl[k]["label"]))
            else:
                L.append("| {} | {:.4f} |".format(decl[k]["label"], v["median_g"]))
        L.append("")
        L.append("Algorithm relative range (1995–2023) = **{:.2%}**.".format(
            late_sens["alg_rel_range_of_medians"]))
    else:
        L.append("Late-era subset had insufficient data.")
    L.append("")

    with open(path_md, "w") as f:
        f.write("\n".join(L) + "\n")
    print("  Wrote {}".format(path_md))


# ────────────────────────────────────────────────────────────────
# Verification
# ────────────────────────────────────────────────────────────────

def run_verify():
    print("=" * 60)
    print("VERIFICATION MODE")
    print("=" * 60)

    rp = os.path.join(WORKSPACE, RESULTS_JSON)
    mp = os.path.join(WORKSPACE, REPORT_MD)
    passed = 0
    total = 18

    def check(num, desc, cond, detail=""):
        nonlocal passed
        mark = "PASS" if cond else "FAIL"
        print("  [{}/{}] {}: {}{}".format(num, total, mark, desc,
              " ({})".format(detail) if detail else ""))
        if cond:
            passed += 1
        return cond

    if not check(1, "results.json exists", os.path.exists(rp)):
        print("\nVERIFICATION FAILED")
        return False

    with open(rp) as f:
        R = json.load(f)
    has_keys = all(k in R for k in (
        "metadata", "declustering", "gr_parameters", "hazard_summary",
        "site_grid_spread", "bootstrap_vs_algorithm",
        "sensitivity_mc", "sensitivity_alt_gmpe", "sensitivity_late_era"))
    check(2, "Required top-level keys present", has_keys)

    meta = R["metadata"]
    check(3, "n_sites >= 400",
          meta["n_sites"] >= 400,
          "{} sites".format(meta["n_sites"]))
    check(4, "Enough events in catalog (>= 3000)",
          meta["n_events_total"] >= 3000,
          "{:,} events".format(meta["n_events_total"]))
    check(5, "All three algorithms present",
          all(a in R["declustering"] for a in ("GK", "R", "ZBZ")))

    # b-values within physical range
    b_ok = all(0.4 <= R["gr_parameters"][k]["b"] <= 1.6
               for k in ("GK", "R", "ZBZ"))
    check(6, "b-values in [0.4, 1.6] for all algorithms",
          b_ok,
          ", ".join("{}={:.3f}".format(k, R["gr_parameters"][k]["b"])
                    for k in ("GK", "R", "ZBZ")))

    # Algorithm spread is a positive, finite number
    sp = R["site_grid_spread"]
    check(7, "Algorithm median relative range finite and positive",
          sp["median_rel_range"] > 0 and sp["median_rel_range"] < 10,
          "{:.2%}".format(sp["median_rel_range"]))

    # Bootstrap comparison populated
    bva = R["bootstrap_vs_algorithm"]
    check(8, "Bootstrap reps >= {} and >= 20 sites".format(N_BOOTSTRAP - 10),
          meta["n_bootstrap"] >= N_BOOTSTRAP - 10 and
          bva["n_sites_boot"] >= 20,
          "n_boot={}, n_sites_boot={}".format(
              meta["n_bootstrap"], bva["n_sites_boot"]))

    check(9, "Mc sensitivity sweep covers >= 3 Mc values",
          len(R["sensitivity_mc"]) >= 3,
          "{} entries".format(len(R["sensitivity_mc"])))

    check(10, "report.md exists and non-empty",
          os.path.exists(mp) and os.path.getsize(mp) > 200,
          "{} bytes".format(os.path.getsize(mp) if os.path.exists(mp) else 0))

    # NEW: Declustering actually removes events — retention in [0.05, 0.95]
    retention = {k: R["declustering"][k]["fraction_retained"]
                 for k in ("GK", "R", "ZBZ")}
    ret_ok = all(0.05 <= v <= 0.95 for v in retention.values())
    check(11, "Declustering retains 5–95% of events (not a no-op and not empty)",
          ret_ok,
          ", ".join("{}={:.1%}".format(k, retention[k])
                    for k in ("GK", "R", "ZBZ")))

    # NEW: CI width is not pathologically tight — sanity check
    ci_w = bva["bootstrap_rel_95ci_median"]
    check(12, "Bootstrap 95% CI relative width > 1% of estimate (sanity)",
          ci_w > 0.01 and math.isfinite(ci_w),
          "{:.2%}".format(ci_w))

    # NEW: Effect-size plausibility (ratio acts like a Cohen-d-equivalent)
    ratio = bva["ratio_alg_to_boot"]
    check(13, "Effect-size ratio alg/boot in [0, 20] (plausibility ceiling)",
          math.isfinite(ratio) and 0 <= ratio <= 20,
          "ratio={:.2f}".format(ratio))

    # NEW: Sensitivity — finding is robust across Mc
    mc_sens = R["sensitivity_mc"]
    mc_medians_by_alg = {k: [] for k in ("GK", "R", "ZBZ")}
    for mc_key, row in mc_sens.items():
        for k in ("GK", "R", "ZBZ"):
            if row.get(k) is not None:
                mc_medians_by_alg[k].append(row[k]["median_pga_g"])
    # Each algorithm's median PGA across Mc should not swing by >10x
    mc_swing_ok = True
    swings = {}
    for k, vals in mc_medians_by_alg.items():
        if vals and min(vals) > 0:
            swings[k] = max(vals) / min(vals)
            if swings[k] > 10.0:
                mc_swing_ok = False
        else:
            mc_swing_ok = False
            swings[k] = float("inf")
    check(14, "Mc sensitivity: median PGA per algorithm swings <10x across Mc",
          mc_swing_ok,
          ", ".join("{}={:.2f}x".format(k, swings[k]) for k in swings))

    # NEW: Negative/falsification — alternate GMPE does not produce pathological
    # cross-algorithm spread
    alt_rel = R["sensitivity_alt_gmpe"]["alg_rel_range_of_medians"]
    check(15, "Negative control: alt-GMPE alg relative range in [0, 2.0]",
          math.isfinite(alt_rel) and 0.0 <= alt_rel <= 2.0,
          "{:.2%}".format(alt_rel))

    # NEW: Per-algorithm median PGAs are in the physical range (0, 3.0) g
    haz = R["hazard_summary"]
    med_ok = all(0.0 < haz[k]["median_g"] < 3.0 for k in ("GK", "R", "ZBZ"))
    check(16, "All per-algorithm median PGAs in (0, 3.0) g",
          med_ok,
          ", ".join("{}={:.4f}g".format(k, haz[k]["median_g"])
                    for k in ("GK", "R", "ZBZ")))

    # NEW 17: Sensitivity-invariance — algorithm ordering (by median PGA) is
    # stable between the reference Mc=4.0 hazard and the late-era sensitivity.
    # If the ordering flips under a time restriction, the headline finding is
    # not robust. We require the ORDER (by median PGA) to be preserved or at
    # worst a single adjacent swap.
    def _order(d, keys=("GK", "R", "ZBZ")):
        return tuple(sorted(keys, key=lambda k: d.get(k, {}).get("median_g", 0.0)))
    ref_order = _order(haz)
    late_haz = R["sensitivity_late_era"].get("per_algorithm") or {}
    order_ok = True
    if late_haz and all(late_haz.get(k) for k in ("GK", "R", "ZBZ")):
        late_order = _order(late_haz)
        # Count positional agreements; require >= 2/3 (i.e., at most one swap)
        agree = sum(1 for a, b in zip(ref_order, late_order) if a == b)
        order_ok = agree >= 2
    check(17, "Algorithm ordering stable under era restriction (>=2/3 match)",
          order_ok,
          "ref={}, late={}".format(ref_order,
              _order(late_haz) if late_haz and all(late_haz.get(k) for k in
                  ("GK","R","ZBZ")) else "n/a"))

    # NEW 18: Falsification / internal-consistency — the bootstrap 95% CI at
    # its MEDIAN site must bracket the GK median-of-per-site PGA to well within
    # a factor of 3 (a loose plausibility check: a bootstrap that does not
    # overlap the point estimate at all would indicate a coding error).
    ci_mid = bva["bootstrap_rel_95ci_median"]
    # CI width as a fraction of estimate; require < 3 (300%) for sanity.
    falsif_ok = math.isfinite(ci_mid) and ci_mid < 3.0
    check(18, "Falsification: bootstrap CI not implausibly wide (<300% of estimate)",
          falsif_ok,
          "ci_mid_width={:.2%}".format(ci_mid))

    print("")
    print("{}/{} assertions passed.".format(passed, total))
    if passed == total:
        print("ALL CHECKS PASSED")
        print("VERIFICATION PASSED")
        return True
    print("VERIFICATION FAILED")
    return False


# ────────────────────────────────────────────────────────────────
# Main
# ────────────────────────────────────────────────────────────────

LIMITATIONS_TEXT = [
    "1. Simplified declustering implementations — our GK/Reasenberg/ZBZ",
    "   approximate the published algorithms; exact retention fractions may",
    "   differ from reference implementations.",
    "2. Uniform areal source zone — 1-degree cells with truncated GR; no fault",
    "   geometry, characteristic events, or spatially varying b.",
    "3. Single simplified GMPE plus one alternate — does not capture the full",
    "   epistemic GMPE logic tree used in modern NSHMs.",
    "4. Bootstrap is within-algorithm only — does NOT estimate joint",
    "   epistemic uncertainty across algorithm choice (that is the separate",
    "   algorithm-spread metric).",
    "5. M>=3.5 fetch and Mc=4.0 — plausibly CONUS-complete but likely complete",
    "   at lower Mc in California / PNW; a spatially-varying Mc is out of scope.",
    "6. Breaks under: ComCat schema changes, catalog <5 yr, <500 events above Mc,",
    "   or offline execution with empty cache.",
]


def print_limitations():
    print("")
    print("=" * 60)
    print("LIMITATIONS AND ASSUMPTIONS")
    print("=" * 60)
    for line in LIMITATIONS_TEXT:
        print(line)
    print("")
    print("What this analysis does NOT show:")
    print("  - Authoritative site hazard (read as internally-comparable across")
    print("    algorithms, not as a replacement for the 2023 NSHM).")
    print("  - The full GMPE epistemic spread (needs a logic tree, not 2 GMPEs).")
    print("  - Truth-level mainshock/aftershock labels for any single event.")


def _check_disk_space(path, required_bytes):
    """Return (ok, free_bytes). Does not raise on missing path."""
    try:
        usage = shutil.disk_usage(path)
        return usage.free >= required_bytes, usage.free
    except (OSError, FileNotFoundError):
        return True, None


def main():
    random.seed(SEED)
    try:
        os.makedirs(DATA_DIR, exist_ok=True)
    except PermissionError as e:
        print("ERROR: cannot create workspace '{}' (permission denied): {}".format(
            DATA_DIR, e), file=sys.stderr)
        sys.exit(2)
    except OSError as e:
        print("ERROR: cannot create workspace '{}': {}".format(DATA_DIR, e),
              file=sys.stderr)
        sys.exit(2)

    # Disk-space preflight: need ~5 MB for cached chunks + outputs.
    ok, free = _check_disk_space(WORKSPACE, 5 * 1024 * 1024)
    if not ok:
        print("ERROR: insufficient disk space in '{}': {} bytes free, "
              "need ~5 MB.".format(WORKSPACE, free), file=sys.stderr)
        sys.exit(2)

    print("=" * 60)
    print("DECLUSTERING ALGORITHM SENSITIVITY — CONUS HAZARD")
    print("=" * 60)
    print("seed={}, bootstrap={}, m_min_fetch={}, Mc={}".format(
        SEED, N_BOOTSTRAP, MIN_MAG_FETCH, MC_ASSUMED))
    print("")

    print("[1/4] Loading ANSS ComCat CONUS catalog...")
    try:
        events, provenance = load_data()
    except (urllib.error.URLError, urllib.error.HTTPError) as e:
        print("ERROR: failed to fetch catalog (network): {}".format(e),
              file=sys.stderr)
        sys.exit(2)
    except RuntimeError as e:
        print("ERROR: failed to fetch catalog: {}".format(e), file=sys.stderr)
        sys.exit(2)
    except (OSError, IOError) as e:
        print("ERROR: I/O error while loading catalog: {}".format(e),
              file=sys.stderr)
        sys.exit(2)

    if len(events) < 500:
        print("ERROR: catalog has only {} events (need >= 500). Widen the "
              "bounding box, extend the time window, or lower MIN_MAG_FETCH.".format(
                  len(events)), file=sys.stderr)
        sys.exit(2)

    print("")
    print("[2/4] Running declustering + hazard pipeline...")
    try:
        results = run_analysis(events)
    except RuntimeError as e:
        print("ERROR: analysis failed: {}".format(e), file=sys.stderr)
        sys.exit(2)
    results["data_provenance"] = provenance
    results["limitations"] = LIMITATIONS_TEXT

    print("")
    print("[3/4] Writing results & report...")
    try:
        generate_report(results)
    except (OSError, IOError) as e:
        print("ERROR: failed to write output files: {}".format(e),
              file=sys.stderr)
        sys.exit(2)

    print("")
    print("[4/4] Summary:")
    sp = results["site_grid_spread"]
    bva = results["bootstrap_vs_algorithm"]
    print("  Algorithm site-wise relative range (median) = {:.2%}".format(
        sp["median_rel_range"]))
    print("  Bootstrap 95% CI width (median, GK)         = {:.2%}".format(
        bva["bootstrap_rel_95ci_median"]))
    print("  Ratio algorithm / bootstrap                  = {:.2f}×".format(
        bva["ratio_alg_to_boot"]))

    print_limitations()

    print("")
    print("ANALYSIS COMPLETE")


if __name__ == "__main__":
    try:
        if "--verify" in sys.argv:
            ok = run_verify()
            sys.exit(0 if ok else 1)
        main()
    except KeyboardInterrupt:
        print("\nERROR: interrupted by user", file=sys.stderr)
        sys.exit(130)
    except Exception as e:
        print("ERROR: unhandled exception: {}: {}".format(
            type(e).__name__, e), file=sys.stderr)
        sys.exit(2)
SCRIPT_EOF
```

**Expected output:**
- Exit code 0, no stdout.
- File `/tmp/claw4s_auto_declustering-algorithm-sensitivity-of-conus-hazard-estimates/analyze.py` created, ~40 KB.

## Step 3: Run the analysis

Run the script to download the catalog, decluster it three ways, and compute hazard curves.

```bash
cd /tmp/claw4s_auto_declustering-algorithm-sensitivity-of-conus-hazard-estimates && python3 analyze.py
```

**Expected output:**
- stdout progress markers `[1/4]` through `[4/4]`.
- Paginated downloads from USGS ComCat into `data/` (each cached file paired with a `.sha256` sidecar).
- Three declustering passes with fraction retained: `GK ≈ 60–90%`, `R ≈ 40–90%`, `ZBZ ≈ 30–90%`.
- GR MLE `b, SE(b)` for each declustered catalog (b roughly in 0.6–1.4).
- PSHA hazard curves at ~500 CONUS sites for each algorithm.
- Bootstrap 100-rep within-GK comparison at 100 sites.
- Sensitivity: Mc ∈ {3.5, 4.0, 4.5}, alternate GMPE, 1995–2023 restriction.
- A final `LIMITATIONS AND ASSUMPTIONS` block (≥4 caveats).
- Final line `ANALYSIS COMPLETE`.
- Artifacts created: `results.json`, `report.md`, `data/conus_*.csv` with SHA256 sidecars.
- Exit code 0.

## Step 4: Verify the results

Run the script in `--verify` mode to machine-check the produced artifacts.

```bash
cd /tmp/claw4s_auto_declustering-algorithm-sensitivity-of-conus-hazard-estimates && python3 analyze.py --verify
```

**Expected output:**
- 18 verification assertions, each tagged `PASS`.
- Final line `ALL CHECKS PASSED` (alias `VERIFICATION PASSED`).
- Exit code 0.

## Success Criteria

These are measurable, machine-checkable conditions. Items marked `[verify]` are enforced by assertions in `python3 analyze.py --verify`.

1. Script runs to completion and prints `ANALYSIS COMPLETE` (exit code 0).
2. [verify] All 18 verification assertions in `--verify` mode pass, and the script prints `ALL CHECKS PASSED`.
3. [verify] `results.json` contains the top-level keys `metadata`, `declustering`, `gr_parameters`, `hazard_summary`, `site_grid_spread`, `bootstrap_vs_algorithm`, `sensitivity_mc`, `sensitivity_alt_gmpe`, `sensitivity_late_era`, `data_provenance`.
4. [verify] All three declustering algorithms produce plausible `b` in the physical range [0.4, 1.6].
5. [verify] Declustering actually removes events: each algorithm retains between 5% and 95% of input.
6. [verify] Site grid has ≥400 sites and catalog has ≥3,000 events.
7. [verify] Bootstrap ran ≥90 replicates at ≥20 sites; bootstrap 95% CI median width is > 1% of the estimate (sanity: bootstrap is not pathologically tight).
8. [verify] Effect-size plausibility: the headline ratio `algorithm range ÷ bootstrap 95% CI` is finite and lies in [0, 20] (a Cohen-d-equivalent ceiling; ratios >20 would indicate pathological bootstrap collapse).
9. [verify] Mc sensitivity: the median 2%-in-50-yr PGA at Mc=3.5 and Mc=4.5 differ by less than a factor of 10; if the finding (ordering of algorithms by median PGA) depends on Mc by orders of magnitude, that is flagged as non-robust.
10. [verify] Negative/falsification check: when an alternate (harder-attenuation) GMPE is substituted, the algorithm-induced relative range of medians stays in [0, 2.0] — i.e., changing the GMPE does not bleed into an implausible cross-algorithm spread. The three algorithms should still remain ordered consistently relative to the baseline GMPE (a qualitative invariance that a plausibility-preserving analysis must respect).
11. [verify] All per-algorithm median site PGAs are in the physical range (0, 3.0) g.
12. [verify] `data_provenance` is populated with SHA-256 digests for every cached chunk.
13. [verify] Sensitivity-invariance: the algorithm ordering (by median PGA) under the 1995-2023 restricted-era sensitivity matches the reference 1990-2023 ordering in at least 2 of 3 positions (i.e., no order-reversing flip under the era restriction).
14. [verify] Falsification / internal-consistency: the bootstrap 95% CI median width is finite and < 300% of the estimate (a bootstrap that is implausibly wide would indicate a coding error in the resampling loop).
15. `report.md` is generated, ≥200 bytes, and includes sections for declustering, GR, hazard summary, and sensitivity.

## Failure Conditions

Any of the following is treated as a failure; the script prints a clear message to `stderr` and exits with a nonzero code (1 for verification failure, 2 for runtime failure):

1. **USGS ComCat unreachable.** The download helper retries 3× with exponential backoff per chunk. If all three attempts fail, `main()` catches the exception, writes `ERROR: failed to fetch catalog: <detail>` to stderr, and exits with code 2 — no partial `results.json` is written.
2. **Empty / malformed CSV chunk.** The script raises `RuntimeError("Unexpected short response")`, which is caught by `main()` and turned into a clean stderr error plus exit code 2.
3. **Insufficient events above Mc.** If `fit_gr` cannot estimate `b` (fewer than 20 events above Mc), `run_hazard_for_catalog` returns `None` and `run_analysis` raises `RuntimeError` with a remediation hint ("lower Mc or widen the catalog"). Caught and exited with code 2.
4. **Any `--verify` assertion fails.** Exit code 1, with `VERIFICATION FAILED` and the failing assertion number(s) printed to stdout.
5. **Python < 3.8.** Not explicitly checked, but the stdlib features used (f-strings, `math.log2`, type hints in comments) target 3.8+; earlier interpreters will fail to import the script.
6. **Write permission denied on workspace.** `os.makedirs(..., exist_ok=True)` will raise `PermissionError`; caught in `main()` and reported.

## Limitations and Assumptions

The following caveats describe what this analysis does **not** show and under which assumptions the results are trustworthy. The script also echoes a subset of these to stdout in a `LIMITATIONS AND ASSUMPTIONS` block.

1. **Simplified declustering implementations.** Our Gardner-Knopoff uses the 1974 time/space tables directly; our Reasenberg is a simplified link-based variant with a Wells-Coppersmith-like interaction radius, not the exact USGS implementation; our Zaliapin-Ben-Zion uses a sliding `WINDOW=5000` prior-event neighbor search to keep runtime feasible. Results are expected to track the published behavior of these algorithms qualitatively but may differ in exact retention fractions. *What this does NOT show:* the absolute truth of mainshock-vs-aftershock classification for any single event.
2. **Uniform areal source zone, no fault geometry.** We bin declustered events into 1°×1° cells and treat each cell as a uniform Poisson source with truncated GR. This ignores fault segmentation, characteristic-earthquake corrections, and spatially varying `b`. *What this does NOT show:* site-specific hazard consistent with the full 2023 NSHM — our PGA estimates should be read as internally-comparable across algorithms, not as authoritative site hazard.
3. **Single simplified GMPE plus one alternate.** We use a Boore-Joyner-Fumal 1997-style PGA GMPE with a fixed aleatory sigma, plus one harder-attenuation alternate for sensitivity. Epistemic GMPE uncertainty in modern NSHMs involves logic trees over 5–15 GMPEs; the true GMPE-induced hazard spread is expected to be larger than what a two-GMPE sensitivity can reveal. *What this does NOT show:* full ground-motion epistemic uncertainty.
4. **M ≥ 3.5 fetch, Mc = 4.0 fit.** The ANSS ComCat at M ≥ 3.5 over CONUS from 1990-2023 is plausibly complete at Mc = 4.0 nationwide but likely complete at lower Mc in California and the Pacific Northwest. A region-specific Mc would change the b-values and rates; our Mc sweep {3.5, 4.0, 4.5} probes this but a spatially varying Mc is out of scope.
5. **Bootstrap is within-algorithm only.** Our 200 (now 100, for runtime) bootstraps resample events from the GK-declustered catalog to quantify within-algorithm variability. This does not estimate the full joint epistemic uncertainty across algorithm choice — that is precisely what our algorithm-spread metric captures separately.
6. **Breaks under:** (a) USGS ComCat schema change (column renames in CSV); (b) catalog windows shorter than 5 years (GR fits become unstable); (c) bounding boxes with fewer than 500 events above Mc (hazard rates dominated by Poisson noise in single cells); (d) offline execution against an empty cache.