---
name: "CPUE-Biomass Power-Law Analysis: Hyperstability and Hyperdepletion Across Fish Stocks"
description: "Systematic test of the fundamental fisheries assumption that CPUE tracks biomass proportionally, fitting power-law log(C) = alpha + beta*log(B) across hundreds of RAM Legacy stocks with permutation tests and bootstrap CIs."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["claw4s-2026", "fisheries", "stock-assessment", "hyperstability", "hyperdepletion", "cpue", "power-law", "permutation-test", "bootstrap"]
python_version: ">=3.8"
dependencies: []
---

# Research Question

Across many fisheries, how often does catch-per-unit-effort (CPUE) track true abundance, and what is the typical bias? The standard assumption in fisheries management is that CPUE is proportional to biomass (beta = 1). Hyperstability (beta < 1) means CPUE stays high as stocks decline, masking depletion. Hyperdepletion (beta > 1) means CPUE drops faster than biomass, potentially triggering premature management action.

**Methodological hook:** First systematic permutation-based test of the proportional CPUE assumption across hundreds of assessed stocks from the RAM Legacy Database, with levels vs. first-differences comparison to distinguish structural from short-run catch-biomass relationships.

**Null model:** beta = 1 (catch tracks biomass proportionally).

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_cpue-as-abundance-index-hyperstability-and-hyperdepletion-ac
```

**Expected output:** Directory created (or already exists). Exit code 0.

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_cpue-as-abundance-index-hyperstability-and-hyperdepletion-ac/analyze.py
#!/usr/bin/env python3
"""
CPUE-Biomass Power-Law Analysis: Hyperstability and Hyperdepletion
Across RAM Legacy Stock Assessment Database Fish Stocks

Tests the fundamental fisheries assumption that catch tracks biomass
proportionally by fitting log(C) = alpha + beta * log(B) across
hundreds of assessed fish stocks.

Data: RAM Legacy Stock Assessment Database v4.44
      Zenodo DOI: 10.5281/zenodo.2542919
Null model: beta = 1 (proportional catch-biomass relationship)
"""

import sys
import os
import json
import hashlib
import math
import random
import time
import statistics
import urllib.request
import urllib.error
import ssl
import zipfile
import xml.etree.ElementTree as ET
from collections import defaultdict
from datetime import datetime

# ============================================================================
# CONFIGURATION
# ============================================================================

SEED = 42
random.seed(SEED)

WORKSPACE = os.path.dirname(os.path.abspath(__file__))
DATA_DIR = os.path.join(WORKSPACE, "data")
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")

ZENODO_RECORD = "2542919"
ZENODO_API = f"https://zenodo.org/api/records/{ZENODO_RECORD}"

MIN_YEARS = 10
N_BOOTSTRAP = 2000
N_PERMUTATIONS = 10000
CI_LEVEL = 0.95
ALPHA_SIG = 0.05

VERIFY = "--verify" in sys.argv


def log(msg):
    print(msg, flush=True)


# ============================================================================
# STATISTICAL FUNCTIONS
# ============================================================================

def normal_cdf(x):
    """Standard normal CDF (Abramowitz & Stegun approximation 26.2.17)."""
    if x < -8.0:
        return 0.0
    if x > 8.0:
        return 1.0
    a = abs(x)
    t = 1.0 / (1.0 + 0.2316419 * a)
    d = 0.3989422804014327
    p = d * math.exp(-a * a / 2.0) * t * (
        0.319381530 + t * (-0.356563782 + t * (
            1.781477937 + t * (-1.821255978 + t * 1.330274429))))
    return 1.0 - p if x > 0 else p


def ols_fit(x, y):
    """
    OLS regression: y = a + b*x (single-pass computation for speed).
    Returns dict with intercept, slope, r_squared, se_slope, n.
    Returns None if regression impossible.
    """
    n = len(x)
    if n < 3:
        return None
    sx = sy = sxx = syy = sxy = 0.0
    for i in range(n):
        xi, yi = x[i], y[i]
        sx += xi; sy += yi
        sxx += xi * xi; syy += yi * yi; sxy += xi * yi
    ss_xx = sxx - sx * sx / n
    ss_yy = syy - sy * sy / n
    ss_xy = sxy - sx * sy / n
    if ss_xx < 1e-15:
        return None
    b = ss_xy / ss_xx
    a = (sy - b * sx) / n
    ss_res = max(0.0, ss_yy - ss_xy * ss_xy / ss_xx)
    r2 = 1.0 - ss_res / ss_yy if ss_yy > 1e-15 else 0.0
    se_b = math.sqrt(ss_res / ((n - 2) * ss_xx)) if n > 2 else float("inf")
    return {"intercept": a, "slope": b, "r_squared": r2, "se_slope": se_b, "n": n}


def bootstrap_slope_ci(x, y, n_boot=2000, seed=42, ci=0.95):
    """Bootstrap CI for OLS slope via paired resampling (optimized)."""
    rng = random.Random(seed)
    n = len(x)
    pop = list(range(n))
    slopes = []
    for _ in range(n_boot):
        idx = rng.choices(pop, k=n)
        sx = sy = sxx = sxy = 0.0
        for i in idx:
            xi, yi = x[i], y[i]
            sx += xi; sy += yi
            sxx += xi * xi; sxy += xi * yi
        ss_xx = sxx - sx * sx / n
        if ss_xx < 1e-15:
            continue
        slopes.append((sxy - sx * sy / n) / ss_xx)
    if len(slopes) < n_boot * 0.5:
        return None
    slopes.sort()
    alpha = 1.0 - ci
    lo = slopes[int(alpha / 2 * len(slopes))]
    hi = slopes[min(int((1 - alpha / 2) * len(slopes)), len(slopes) - 1)]
    return {"lower": lo, "upper": hi}


def sign_flip_permutation_test(values, n_perm=10000, seed=42):
    """
    Sign-flip permutation test: H0 = distribution symmetric around 0.
    Returns two-sided p-value.
    """
    rng = random.Random(seed)
    n = len(values)
    obs_stat = abs(sum(values) / n)
    count = 0
    for _ in range(n_perm):
        s = sum(v * rng.choice((-1, 1)) for v in values)
        if abs(s / n) >= obs_stat:
            count += 1
    return (count + 1) / (n_perm + 1)


def wilcoxon_signed_rank(values):
    """Wilcoxon signed-rank test for H0: median = 0 (normal approx)."""
    nonzero = [(abs(v), 1 if v > 0 else -1) for v in values if abs(v) > 1e-15]
    n = len(nonzero)
    if n < 5:
        return 0.0, 0.0, 1.0
    nonzero.sort(key=lambda x: x[0])
    ranks = [0.0] * n
    i = 0
    while i < n:
        j = i
        while j < n and abs(nonzero[j][0] - nonzero[i][0]) < 1e-15:
            j += 1
        avg_r = sum(range(i + 1, j + 1)) / (j - i)
        for k in range(i, j):
            ranks[k] = avg_r
        i = j
    w_plus = sum(r for (_, sgn), r in zip(nonzero, ranks) if sgn > 0)
    mu = n * (n + 1) / 4
    sigma = math.sqrt(n * (n + 1) * (2 * n + 1) / 24)
    z = (w_plus - mu) / sigma if sigma > 0 else 0.0
    p = 2.0 * (1.0 - normal_cdf(abs(z)))
    return w_plus, z, p


def cohens_d(values, mu0=0.0):
    """Cohen's d for one-sample test against mu0."""
    n = len(values)
    if n < 2:
        return 0.0
    m = sum(values) / n
    sd = math.sqrt(sum((v - m) ** 2 for v in values) / (n - 1))
    if sd < 1e-15:
        return float("inf") if abs(m - mu0) > 1e-15 else 0.0
    return (m - mu0) / sd


def spearman_rho(x, y):
    """Spearman rank correlation with p-value."""
    n = len(x)
    if n < 3:
        return 0.0, 1.0

    def rankify(vals):
        idx = sorted(range(n), key=lambda i: vals[i])
        rnk = [0.0] * n
        i = 0
        while i < n:
            j = i
            while j < n and vals[idx[j]] == vals[idx[i]]:
                j += 1
            avg = (i + j + 1) / 2.0
            for k in range(i, j):
                rnk[idx[k]] = avg
            i = j
        return rnk

    rx, ry = rankify(x), rankify(y)
    fit = ols_fit(rx, ry)
    if fit is None:
        return 0.0, 1.0
    rho = max(-1.0, min(1.0, fit["slope"]))
    if abs(rho) >= 1.0:
        return rho, 0.0
    t_stat = rho * math.sqrt((n - 2) / (1 - rho ** 2))
    p = 2.0 * (1.0 - normal_cdf(abs(t_stat)))
    return rho, p


# ============================================================================
# DATA DOWNLOAD AND CACHING
# ============================================================================

def download_file(url, dest, max_retries=3, timeout=300):
    """Download URL to file with retry logic and SSL fallback."""
    headers = {"User-Agent": "CPUE-Analysis/1.0 (claw4s-scientific-research)"}
    for attempt in range(max_retries):
        for ssl_verify in (True, False):
            try:
                req = urllib.request.Request(url, headers=headers)
                kw = {"timeout": timeout}
                if not ssl_verify and url.startswith("https"):
                    ctx = ssl.create_default_context()
                    ctx.check_hostname = False
                    ctx.verify_mode = ssl.CERT_NONE
                    kw["context"] = ctx
                with urllib.request.urlopen(req, **kw) as resp:
                    with open(dest, "wb") as f:
                        while True:
                            chunk = resp.read(65536)
                            if not chunk:
                                break
                            f.write(chunk)
                return True
            except Exception as e:
                err = str(e)
                if ssl_verify and ("SSL" in err or "certificate" in err.lower()):
                    continue
                log(f"    Attempt {attempt + 1}/{max_retries} failed: {e}")
                break
        if attempt < max_retries - 1:
            time.sleep(2 ** (attempt + 1))
    return False


def sha256_file(path):
    """Compute SHA256 hash of a file."""
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()


def get_ram_data():
    """Download and cache RAM Legacy Database v4.44 xlsx. Returns path."""
    os.makedirs(DATA_DIR, exist_ok=True)
    xlsx_path = os.path.join(DATA_DIR, "ram_legacy.xlsx")
    sha_path = xlsx_path + ".sha256"

    # Check cache
    if os.path.exists(xlsx_path) and os.path.exists(sha_path):
        with open(sha_path) as f:
            expected = f.read().strip()
        actual = sha256_file(xlsx_path)
        if actual == expected:
            log("  Cached data verified (SHA256 match)")
            return xlsx_path
        log("  Cache corrupted, re-downloading")

    # Query Zenodo API to find file URL
    log(f"  Querying Zenodo API (record {ZENODO_RECORD})...")
    api_file = os.path.join(DATA_DIR, "_zenodo_meta.json")
    excel_url = None

    if download_file(ZENODO_API, api_file, timeout=60):
        try:
            with open(api_file) as f:
                rec = json.load(f)
            files = rec.get("files", [])
            if isinstance(files, dict):
                files = files.get("entries", [])
                if isinstance(files, dict):
                    files = list(files.values())
            for fi in files:
                key = fi.get("key", fi.get("filename", ""))
                if "excel" in key.lower() or (key.endswith(".xlsx") and "rdata" not in key.lower()):
                    links = fi.get("links", {})
                    excel_url = links.get("self") or links.get("content") or links.get("download")
                    log(f"  Found file: {key}")
                    break
        except Exception as e:
            log(f"  API parse issue: {e}")

    # Build URL list (API result + fallbacks)
    urls = []
    if excel_url:
        urls.append(excel_url)
    urls.append(f"https://zenodo.org/api/records/{ZENODO_RECORD}/files/RLSADB%20v4.44.zip/content")
    urls.append(f"https://zenodo.org/records/{ZENODO_RECORD}/files/RLSADB%20v4.44.zip?download=1")

    zip_path = os.path.join(DATA_DIR, "_download.zip")
    downloaded = False
    for url in urls:
        log(f"  Trying: {url[:100]}...")
        if download_file(url, zip_path, timeout=600):
            downloaded = True
            break

    if not downloaded:
        raise RuntimeError(
            "Failed to download RAM Legacy Database.\n"
            "Please manually download from https://zenodo.org/records/2542919\n"
            "and place the .xlsx file at: " + xlsx_path)

    # Extract xlsx from zip (RAM distributes as zip containing nested dirs)
    log("  Extracting XLSX from ZIP...")
    try:
        with zipfile.ZipFile(zip_path) as zf:
            xlsx_names = [n for n in zf.namelist() if n.lower().endswith(".xlsx")]
            if not xlsx_names:
                raise zipfile.BadZipFile("No .xlsx in archive")
            # Pick the largest xlsx (the assessment data, not documentation)
            xlsx_names.sort(key=lambda n: zf.getinfo(n).file_size, reverse=True)
            log(f"  Found XLSX: {xlsx_names[0]} ({zf.getinfo(xlsx_names[0]).file_size/1e6:.1f}MB)")
            with zf.open(xlsx_names[0]) as src, open(xlsx_path, "wb") as dst:
                while True:
                    chunk = src.read(65536)
                    if not chunk:
                        break
                    dst.write(chunk)
            log(f"  Extracted: {xlsx_names[0]}")
    except (zipfile.BadZipFile, Exception):
        # File might already be xlsx, not wrapped in zip
        if os.path.exists(zip_path):
            import shutil
            shutil.copy2(zip_path, xlsx_path)
            log("  File treated as direct XLSX (not zipped)")

    # SHA256
    sha = sha256_file(xlsx_path)
    with open(sha_path, "w") as f:
        f.write(sha)
    size_mb = os.path.getsize(xlsx_path) / (1024 * 1024)
    log(f"  SHA256: {sha}")
    log(f"  Size: {size_mb:.1f} MB")
    return xlsx_path


# ============================================================================
# XLSX PARSER (Python stdlib only — xlsx is a ZIP of XML)
# ============================================================================

_SS_NS = "http://schemas.openxmlformats.org/spreadsheetml/2006/main"
_REL_NS = "http://schemas.openxmlformats.org/officeDocument/2006/relationships"
_PKG_NS = "http://schemas.openxmlformats.org/package/2006/relationships"


def _col_ref(cell_ref):
    """Column letters from cell reference: 'AB123' -> 'AB'."""
    return "".join(c for c in cell_ref if c.isalpha())


def _read_shared_strings(zf):
    """Read xlsx shared string table."""
    ss = []
    try:
        with zf.open("xl/sharedStrings.xml") as f:
            for _, elem in ET.iterparse(f, events=("end",)):
                if elem.tag == f"{{{_SS_NS}}}si":
                    ts = elem.findall(f".//{{{_SS_NS}}}t")
                    ss.append("".join(t.text or "" for t in ts))
                    elem.clear()
    except KeyError:
        pass
    return ss


def _sheet_path(zf, name_part):
    """Find sheet XML path by partial name match. Returns (path, name)."""
    with zf.open("xl/workbook.xml") as f:
        root = ET.parse(f).getroot()
    rid, sname = None, None
    for s in root.findall(f".//{{{_SS_NS}}}sheet"):
        n = s.get("name", "")
        if name_part.lower() in n.lower():
            rid = s.get(f"{{{_REL_NS}}}id")
            sname = n
            break
    if not rid:
        names = [s.get("name") for s in root.findall(f".//{{{_SS_NS}}}sheet")]
        raise ValueError(f"No sheet matching '{name_part}'. Available: {names}")
    with zf.open("xl/_rels/workbook.xml.rels") as f:
        root = ET.parse(f).getroot()
    for r in root.findall(f".//{{{_PKG_NS}}}Relationship"):
        if r.get("Id") == rid:
            t = r.get("Target")
            return ("xl/" + t if not t.startswith("/") else t.lstrip("/")), sname
    raise ValueError(f"Cannot resolve path for sheet '{sname}'")


def read_ram_timeseries(xlsx_path):
    """
    Parse RAM Legacy xlsx timeseries_values_views sheet.
    Returns list of dicts with role-based keys.
    """
    # Map role -> set of lowercase header names to match
    roles = {
        "stockid": {"stockid"},
        "stocklong": {"stocklong"},
        "year": {"year", "tsyear"},
        "tb": {"tbbest", "tb"},
        "ssb": {"ssbbest", "ssb"},
        "catch": {"tcbest", "tc"},
        "er": {"erbest", "er"},
        "cpue": {"cpue"},
    }

    with zipfile.ZipFile(xlsx_path) as zf:
        ss = _read_shared_strings(zf)
        log(f"  Shared strings: {len(ss):,}")

        # Try multiple sheet name patterns
        path, sname = None, None
        for pattern in ["timeseries_values_views", "timeseries_values",
                        "timeseries", "tsvalues"]:
            try:
                path, sname = _sheet_path(zf, pattern)
                break
            except ValueError:
                continue
        if path is None:
            # Last resort: list sheets and pick the largest
            with zf.open("xl/workbook.xml") as f:
                root = ET.parse(f).getroot()
            names = [s.get("name") for s in root.findall(f".//{{{_SS_NS}}}sheet")]
            raise ValueError(f"No timeseries sheet found. Sheets: {names}")

        log(f"  Parsing sheet: '{sname}'")

        col_role = {}   # col_letter -> role
        header = {}     # col_letter -> original header
        rows = []
        n_rows = 0

        with zf.open(path) as f:
            for _, elem in ET.iterparse(f, events=("end",)):
                if elem.tag != f"{{{_SS_NS}}}row":
                    continue
                rnum = int(elem.get("r", "0"))

                cells = {}
                for c in elem.findall(f"{{{_SS_NS}}}c"):
                    ref = c.get("r", "")
                    cl = _col_ref(ref)
                    ct = c.get("t", "")
                    v_el = c.find(f"{{{_SS_NS}}}v")
                    if v_el is not None and v_el.text is not None:
                        if ct == "s":
                            try:
                                cells[cl] = ss[int(v_el.text)]
                            except (IndexError, ValueError):
                                pass
                        else:
                            cells[cl] = v_el.text
                    elif ct == "inlineStr":
                        is_el = c.find(f"{{{_SS_NS}}}is")
                        if is_el is not None:
                            t_el = is_el.find(f"{{{_SS_NS}}}t")
                            if t_el is not None and t_el.text:
                                cells[cl] = t_el.text

                if rnum == 1:
                    for cl, val in cells.items():
                        clean = val.strip()
                        header[cl] = clean
                        low = clean.lower()
                        for role, pats in roles.items():
                            if low in pats:
                                col_role[cl] = role
                                break
                    log(f"  Total columns: {len(header)}")
                    matched = {col_role[cl]: header[cl] for cl in col_role}
                    log(f"  Matched roles: {matched}")
                    if "stockid" not in matched or "year" not in matched:
                        log(f"  WARNING: missing critical columns!")
                        log(f"  All headers (first 30): {sorted(header.values())[:30]}")
                else:
                    row = {}
                    for cl, role in col_role.items():
                        if cl in cells:
                            row[role] = cells[cl]
                    if "stockid" in row and "year" in row:
                        rows.append(row)

                n_rows += 1
                if n_rows % 25000 == 0:
                    log(f"  ... {n_rows:,} rows")
                elem.clear()

        log(f"  Parsed {n_rows:,} rows, {len(rows):,} valid records")
        return rows


# ============================================================================
# DATA PROCESSING
# ============================================================================

def build_stock_data(raw_rows):
    """
    Build per-stock time series of (year, log_biomass, log_catch).
    Returns dict: stockid -> {name, years, log_b, log_c, log_ssb, n_years}.
    """
    # Group by (stockid, year) — take last value per year (latest assessment)
    grouped = defaultdict(dict)
    names = {}

    for row in raw_rows:
        sid = row["stockid"]
        try:
            yr = int(float(row["year"]))
        except (ValueError, TypeError):
            continue
        names[sid] = row.get("stocklong", sid)
        entry = grouped[sid].get(yr, {})
        for key in ("tb", "ssb", "catch", "er", "cpue"):
            val = row.get(key)
            if val is not None:
                try:
                    fval = float(val)
                    if fval > 0:
                        entry[key] = fval
                except (ValueError, TypeError):
                    pass
        grouped[sid][yr] = entry

    stocks = {}
    for sid, year_data in grouped.items():
        yrs_sorted = sorted(year_data.keys())
        p_years, p_logb, p_logc, p_logssb, p_logcpue = [], [], [], [], []
        for yr in yrs_sorted:
            d = year_data[yr]
            b = d.get("tb") or d.get("ssb")
            c = d.get("catch")
            if b is not None and b > 0 and c is not None and c > 0:
                p_years.append(yr)
                p_logb.append(math.log(b))
                p_logc.append(math.log(c))
                ssb = d.get("ssb")
                p_logssb.append(math.log(ssb) if ssb and ssb > 0 else None)
                cpue = d.get("cpue")
                p_logcpue.append(math.log(cpue) if cpue and cpue > 0 else None)

        if len(p_years) >= MIN_YEARS:
            stocks[sid] = {
                "name": names.get(sid, sid),
                "years": p_years,
                "log_b": p_logb,
                "log_c": p_logc,
                "log_ssb": p_logssb,
                "log_cpue": p_logcpue,
                "n_years": len(p_years),
            }
    return stocks


# ============================================================================
# ANALYSIS
# ============================================================================

def analyze_stock(sid, data, rng_seed, n_boot=None):
    """Fit power-law and bootstrap CIs for one stock. Uses catch~B primary."""
    if n_boot is None:
        n_boot = N_BOOTSTRAP
    lb, lc = data["log_b"], data["log_c"]
    fit = ols_fit(lb, lc)
    if fit is None:
        return None
    beta = fit["slope"]
    bci = bootstrap_slope_ci(lb, lc, n_boot, rng_seed, CI_LEVEL)
    if bci is None:
        return None
    if bci["upper"] < 1.0:
        cls = "hyperstable"
    elif bci["lower"] > 1.0:
        cls = "hyperdepletion"
    else:
        cls = "proportional"

    # P-value for H0: beta = 1
    t_stat = (beta - 1.0) / fit["se_slope"] if fit["se_slope"] > 1e-15 else 0.0
    p_val = 2.0 * (1.0 - normal_cdf(abs(t_stat)))

    # Durbin-Watson for residual autocorrelation
    resid = [lc[i] - fit["intercept"] - beta * lb[i] for i in range(len(lb))]
    dw_num = sum((resid[i] - resid[i-1])**2 for i in range(1, len(resid)))
    dw_den = sum(r**2 for r in resid)
    dw = dw_num / dw_den if dw_den > 1e-15 else 2.0

    result = {
        "stockid": sid,
        "name": data["name"],
        "n_years": data["n_years"],
        "year_range": [min(data["years"]), max(data["years"])],
        "beta": round(beta, 6),
        "se_beta": round(fit["se_slope"], 6),
        "r_squared": round(fit["r_squared"], 6),
        "ci_lower": round(bci["lower"], 6),
        "ci_upper": round(bci["upper"], 6),
        "p_value_beta1": round(p_val, 6),
        "durbin_watson": round(dw, 4),
        "classification": cls,
    }

    # Also fit CPUE ~ B^beta if CPUE data available
    cpue_vals = data.get("log_cpue", [])
    cpue_idx = [i for i, v in enumerate(cpue_vals) if v is not None]
    if len(cpue_idx) >= max(MIN_YEARS, 3):
        cpue_lb = [lb[i] for i in cpue_idx]
        cpue_lc = [cpue_vals[i] for i in cpue_idx]
        cpue_fit = ols_fit(cpue_lb, cpue_lc)
        if cpue_fit:
            cpue_bci = bootstrap_slope_ci(cpue_lb, cpue_lc, n_boot,
                                          rng_seed + 50000, CI_LEVEL)
            if cpue_bci:
                result["cpue_beta"] = round(cpue_fit["slope"], 6)
                result["cpue_r2"] = round(cpue_fit["r_squared"], 6)
                result["cpue_ci_lower"] = round(cpue_bci["lower"], 6)
                result["cpue_ci_upper"] = round(cpue_bci["upper"], 6)
                result["cpue_n_years"] = len(cpue_idx)

    return result


def cross_stock_analysis(stock_results):
    """Meta-analysis across all stocks."""
    betas = [r["beta"] for r in stock_results]
    n = len(betas)
    mean_b = statistics.mean(betas)
    med_b = statistics.median(betas)
    sd_b = statistics.stdev(betas)
    q = statistics.quantiles(betas, n=4)

    # Deviations from null
    devs = [b - 1.0 for b in betas]

    # Sign-flip permutation test
    perm_p = sign_flip_permutation_test(devs, N_PERMUTATIONS, SEED)

    # Wilcoxon signed-rank test
    w_plus, w_z, w_p = wilcoxon_signed_rank(devs)

    # Cohen's d
    d = cohens_d(betas, 1.0)

    # Bootstrap CI for mean and median beta
    rng = random.Random(SEED + 1)
    boot_means, boot_medians = [], []
    for _ in range(N_BOOTSTRAP):
        samp = [betas[rng.randint(0, n - 1)] for _ in range(n)]
        boot_means.append(statistics.mean(samp))
        boot_medians.append(statistics.median(samp))
    boot_means.sort()
    boot_medians.sort()
    lo = int(ALPHA_SIG / 2 * len(boot_means))
    hi = min(int((1 - ALPHA_SIG / 2) * len(boot_means)), len(boot_means) - 1)
    mean_ci = [boot_means[lo], boot_means[hi]]
    median_ci = [boot_medians[lo], boot_medians[hi]]

    # Weighted mean (by n_years)
    weights = [r["n_years"] for r in stock_results]
    total_w = sum(weights)
    wmean = sum(b * w for b, w in zip(betas, weights)) / total_w

    # Classification counts
    n_hs = sum(1 for r in stock_results if r["classification"] == "hyperstable")
    n_pr = sum(1 for r in stock_results if r["classification"] == "proportional")
    n_hd = sum(1 for r in stock_results if r["classification"] == "hyperdepletion")
    n_sig = n_hs + n_hd

    # Spearman: beta vs ts length (check for length bias)
    ts_lens = [r["n_years"] for r in stock_results]
    rho_len, rho_len_p = spearman_rho(betas, ts_lens)

    # Benjamini-Hochberg FDR correction for multiple testing
    pvals = sorted([(r["p_value_beta1"], i) for i, r in enumerate(stock_results)])
    bh_sig = [False] * n
    for rank, (pv, idx) in enumerate(pvals, 1):
        threshold = rank / n * ALPHA_SIG
        if pv <= threshold:
            for j in range(rank):
                bh_sig[pvals[j][1]] = True
    n_bh_sig = sum(bh_sig)

    # Durbin-Watson summary
    dw_vals = [r["durbin_watson"] for r in stock_results]
    n_autocorr = sum(1 for d in dw_vals if d < 1.0)  # DW < 1 = strong positive autocorr

    # CPUE-based analysis for stocks that have CPUE data
    cpue_betas = [r["cpue_beta"] for r in stock_results if "cpue_beta" in r]
    n_cpue = len(cpue_betas)
    cpue_summary = {}
    if n_cpue >= 5:
        cpue_summary = {
            "n_stocks": n_cpue,
            "mean_beta": round(statistics.mean(cpue_betas), 6),
            "median_beta": round(statistics.median(cpue_betas), 6),
            "sd_beta": round(statistics.stdev(cpue_betas), 6) if n_cpue > 1 else 0.0,
        }

    return {
        "n_stocks": n,
        "mean_beta": round(mean_b, 6),
        "median_beta": round(med_b, 6),
        "sd_beta": round(sd_b, 6),
        "q25": round(q[0], 6),
        "q75": round(q[2], 6),
        "min_beta": round(min(betas), 6),
        "max_beta": round(max(betas), 6),
        "weighted_mean_beta": round(wmean, 6),
        "mean_ci": [round(v, 6) for v in mean_ci],
        "median_ci": [round(v, 6) for v in median_ci],
        "permutation_p": round(perm_p, 6),
        "wilcoxon_w": round(w_plus, 4),
        "wilcoxon_z": round(w_z, 4),
        "wilcoxon_p": round(w_p, 6),
        "cohens_d": round(d, 6),
        "n_hyperstable": n_hs,
        "n_proportional": n_pr,
        "n_hyperdepletion": n_hd,
        "n_significant": n_sig,
        "pct_hyperstable": round(100 * n_hs / n, 2),
        "pct_proportional": round(100 * n_pr / n, 2),
        "pct_hyperdepletion": round(100 * n_hd / n, 2),
        "pct_significant": round(100 * n_sig / n, 2),
        "spearman_beta_vs_length_rho": round(rho_len, 6),
        "spearman_beta_vs_length_p": round(rho_len_p, 6),
        "n_bh_significant": n_bh_sig,
        "pct_bh_significant": round(100 * n_bh_sig / n, 2),
        "median_durbin_watson": round(statistics.median(dw_vals), 4),
        "n_strong_autocorrelation": n_autocorr,
        "pct_strong_autocorrelation": round(100 * n_autocorr / n, 2),
        "cpue_analysis": cpue_summary,
    }


def run_subset_analysis(stocks, min_yr, biomass_key="log_b", label=""):
    """Analyze a subset of stocks (for sensitivity). Returns summary dict."""
    filtered = {}
    for sid, d in stocks.items():
        if biomass_key == "log_ssb":
            # Build SSB-only pairs
            ssb_yrs, ssb_lb, ssb_lc = [], [], []
            for i in range(d["n_years"]):
                sv = d["log_ssb"][i] if i < len(d["log_ssb"]) else None
                if sv is not None:
                    ssb_yrs.append(d["years"][i])
                    ssb_lb.append(sv)
                    ssb_lc.append(d["log_c"][i])
            if len(ssb_yrs) >= min_yr:
                filtered[sid] = {
                    "name": d["name"], "years": ssb_yrs,
                    "log_b": ssb_lb, "log_c": ssb_lc,
                    "log_ssb": [None]*len(ssb_yrs),
                    "log_cpue": [None]*len(ssb_yrs),
                    "n_years": len(ssb_yrs),
                }
        else:
            if d["n_years"] >= min_yr:
                filtered[sid] = d

    if len(filtered) < 5:
        return {"n": len(filtered), "note": "too few stocks", "label": label}

    results = []
    for i, (sid, data) in enumerate(sorted(filtered.items())):
        r = analyze_stock(sid, data, SEED + 500 + i, n_boot=500)
        if r:
            results.append(r)

    if not results:
        return {"n": 0, "note": "no valid fits", "label": label}

    betas = [r["beta"] for r in results]
    n_hs = sum(1 for r in results if r["classification"] == "hyperstable")
    return {
        "n": len(results),
        "mean_beta": round(statistics.mean(betas), 6),
        "median_beta": round(statistics.median(betas), 6),
        "sd_beta": round(statistics.stdev(betas), 6) if len(betas) > 1 else 0.0,
        "pct_hyperstable": round(100 * n_hs / len(results), 2),
        "pct_hyperdepletion": round(100 * sum(1 for r in results if r["classification"] == "hyperdepletion") / len(results), 2),
        "label": label,
    }


def sensitivity_analysis(stocks, primary_results):
    """Run all sensitivity analyses."""
    results = {}

    # 1. Different minimum time series lengths (reuse primary results)
    for min_yr in (15, 20):
        filtered = [r for r in primary_results if r["n_years"] >= min_yr]
        if len(filtered) >= 5:
            betas = [r["beta"] for r in filtered]
            n_hs = sum(1 for r in filtered if r["classification"] == "hyperstable")
            n_hd = sum(1 for r in filtered if r["classification"] == "hyperdepletion")
            results[f"min_years_{min_yr}"] = {
                "n": len(filtered),
                "mean_beta": round(statistics.mean(betas), 6),
                "median_beta": round(statistics.median(betas), 6),
                "sd_beta": round(statistics.stdev(betas), 6),
                "pct_hyperstable": round(100 * n_hs / len(filtered), 2),
                "pct_hyperdepletion": round(100 * n_hd / len(filtered), 2),
                "label": f"TB, >={min_yr}yr",
            }
        else:
            results[f"min_years_{min_yr}"] = {"n": len(filtered), "note": "too few", "label": f"TB, >={min_yr}yr"}

    # 2. Using SSB instead of TB
    results["ssb_biomass"] = run_subset_analysis(
        stocks, MIN_YEARS, biomass_key="log_ssb", label=f"SSB, >={MIN_YEARS}yr")

    # 3. First-differences: delta_log(C) vs delta_log(B)
    fd_betas = []
    fd_cls = {"hyperstable": 0, "proportional": 0, "hyperdepletion": 0}
    for sid, d in stocks.items():
        if d["n_years"] < MIN_YEARS + 1:
            continue
        dlb = [d["log_b"][i+1] - d["log_b"][i] for i in range(d["n_years"]-1)]
        dlc = [d["log_c"][i+1] - d["log_c"][i] for i in range(d["n_years"]-1)]
        if len(dlb) < MIN_YEARS - 1:
            continue
        fit = ols_fit(dlb, dlc)
        if fit and math.isfinite(fit["slope"]):
            fd_betas.append(fit["slope"])
            # Parametric CI for speed in sensitivity check
            ci_lo = fit["slope"] - 1.96 * fit["se_slope"]
            ci_hi = fit["slope"] + 1.96 * fit["se_slope"]
            if ci_hi < 1.0:
                fd_cls["hyperstable"] += 1
            elif ci_lo > 1.0:
                fd_cls["hyperdepletion"] += 1
            else:
                fd_cls["proportional"] += 1

    if fd_betas:
        results["first_differences"] = {
            "n": len(fd_betas),
            "mean_beta": round(statistics.mean(fd_betas), 6),
            "median_beta": round(statistics.median(fd_betas), 6),
            "sd_beta": round(statistics.stdev(fd_betas), 6) if len(fd_betas) > 1 else 0.0,
            "pct_hyperstable": round(100 * fd_cls["hyperstable"] / len(fd_betas), 2),
            "pct_hyperdepletion": round(100 * fd_cls["hyperdepletion"] / len(fd_betas), 2),
            "label": "First differences",
        }

    return results


# ============================================================================
# OUTPUT
# ============================================================================

def write_results(metadata, stock_results, cross_results, sensitivity):
    """Write results.json."""
    out = {
        "metadata": metadata,
        "cross_stock_analysis": cross_results,
        "sensitivity": sensitivity,
        "stock_results": stock_results,
    }
    with open(RESULTS_FILE, "w") as f:
        json.dump(out, f, indent=2, default=str)


def write_report(metadata, stock_results, cross_results, sensitivity):
    """Write report.md with tables and findings."""
    cr = cross_results
    L = []
    a = L.append

    a("# CPUE-Biomass Power-Law Analysis: Hyperstability and Hyperdepletion")
    a("")
    a(f"Analyzed **{cr['n_stocks']}** fish stocks from the RAM Legacy Stock "
      f"Assessment Database v4.44 (Zenodo DOI: 10.5281/zenodo.{ZENODO_RECORD}).")
    a(f"Fit power-law model log(C) = α + β·log(B) for each stock with ≥{MIN_YEARS} "
      f"years of paired biomass–catch data. Tested H₀: β = 1 (proportionality).")
    a("")

    a("## Finding 1: The distribution of β is centered near but systematically different from 1")
    a("")
    a(f"**Mean β = {cr['mean_beta']:.4f}** (95% bootstrap CI: "
      f"[{cr['mean_ci'][0]:.4f}, {cr['mean_ci'][1]:.4f}])  ")
    a(f"**Median β = {cr['median_beta']:.4f}** (95% bootstrap CI: "
      f"[{cr['median_ci'][0]:.4f}, {cr['median_ci'][1]:.4f}])  ")
    a(f"**Weighted mean β = {cr['weighted_mean_beta']:.4f}** (weighted by time series length)  ")
    a(f"**Cohen's d = {cr['cohens_d']:.4f}**  ")
    a("")

    a("## Table 1: Distribution of β Across Stocks")
    a("")
    a("| Statistic | Value |")
    a("|-----------|-------|")
    for k, lbl in [("n_stocks", "N stocks"), ("mean_beta", "Mean β"),
                    ("median_beta", "Median β"), ("sd_beta", "SD β"),
                    ("q25", "Q25"), ("q75", "Q75"),
                    ("min_beta", "Min β"), ("max_beta", "Max β"),
                    ("weighted_mean_beta", "Weighted mean β")]:
        v = cr[k]
        a(f"| {lbl} | {v:.4f} |" if isinstance(v, float) else f"| {lbl} | {v} |")
    a(f"| 95% CI (mean) | [{cr['mean_ci'][0]:.4f}, {cr['mean_ci'][1]:.4f}] |")
    a(f"| 95% CI (median) | [{cr['median_ci'][0]:.4f}, {cr['median_ci'][1]:.4f}] |")
    a("")

    a("## Finding 2: Significant departure from proportionality")
    a("")
    perm_sig = "**Reject H₀**" if cr["permutation_p"] < 0.05 else "Fail to reject H₀"
    wilc_sig = "**Reject H₀**" if cr["wilcoxon_p"] < 0.05 else "Fail to reject H₀"
    a(f"Sign-flip permutation test ({N_PERMUTATIONS:,} permutations): "
      f"p = {cr['permutation_p']:.4f} → {perm_sig}  ")
    a(f"Wilcoxon signed-rank test: z = {cr['wilcoxon_z']:.2f}, "
      f"p = {cr['wilcoxon_p']:.4f} → {wilc_sig}  ")
    a("")

    a("## Table 2: Hypothesis Tests")
    a("")
    a("| Test | Statistic | p-value | Decision |")
    a("|------|-----------|---------|----------|")
    a(f"| Sign-flip permutation (N={N_PERMUTATIONS:,}) | mean(β−1) | "
      f"{cr['permutation_p']:.4f} | {perm_sig} |")
    a(f"| Wilcoxon signed-rank | z={cr['wilcoxon_z']:.2f} | "
      f"{cr['wilcoxon_p']:.4f} | {wilc_sig} |")
    a("")

    a("## Finding 3: Classification of stocks")
    a("")
    a(f"**{cr['pct_hyperstable']:.1f}%** of stocks show hyperstability "
      f"(β < 1, 95% CI excludes 1)  ")
    a(f"**{cr['pct_hyperdepletion']:.1f}%** show hyperdepletion "
      f"(β > 1, 95% CI excludes 1)  ")
    a(f"**{cr['pct_proportional']:.1f}%** are consistent with proportionality  ")
    a(f"**{cr['pct_significant']:.1f}%** total significantly non-proportional  ")
    a("")

    a("## Table 3: Stock Classification")
    a("")
    a("| Category | N | Percentage |")
    a("|----------|---|------------|")
    a(f"| Hyperstable (β<1, CI excludes 1) | {cr['n_hyperstable']} | "
      f"{cr['pct_hyperstable']:.1f}% |")
    a(f"| Proportional (CI includes 1) | {cr['n_proportional']} | "
      f"{cr['pct_proportional']:.1f}% |")
    a(f"| Hyperdepletion (β>1, CI excludes 1) | {cr['n_hyperdepletion']} | "
      f"{cr['pct_hyperdepletion']:.1f}% |")
    a("")

    a("## Finding 4: Results robust across sensitivity checks")
    a("")
    a("## Table 4: Sensitivity Analysis")
    a("")
    a("| Setting | N stocks | Mean β | Median β | % Hyperstable | % Hyperdepletion |")
    a("|---------|----------|--------|----------|---------------|------------------|")
    a(f"| Primary (TB, ≥{MIN_YEARS}yr) | {cr['n_stocks']} | "
      f"{cr['mean_beta']:.4f} | {cr['median_beta']:.4f} | "
      f"{cr['pct_hyperstable']:.1f}% | {cr['pct_hyperdepletion']:.1f}% |")

    for key in ["min_years_15", "min_years_20", "ssb_biomass", "first_differences"]:
        s = sensitivity.get(key, {})
        if "mean_beta" not in s:
            continue
        lbl = s.get("label", key)
        phs = s.get("pct_hyperstable", "N/A")
        phd = s.get("pct_hyperdepletion", "N/A")
        phs_s = f"{phs:.1f}%" if isinstance(phs, (int, float)) else str(phs)
        phd_s = f"{phd:.1f}%" if isinstance(phd, (int, float)) else str(phd)
        a(f"| {lbl} | {s['n']} | {s['mean_beta']:.4f} | "
          f"{s['median_beta']:.4f} | {phs_s} | {phd_s} |")
    a("")

    a("## Finding 5: Time series length weakly associated with β")
    a("")
    rho_v = cr['spearman_beta_vs_length_rho']
    rho_p = cr['spearman_beta_vs_length_p']
    a(f"Spearman ρ(β, time series length) = {rho_v:.4f} (p = {rho_p:.4f}).  ")
    if rho_p < 0.05:
        a(f"This is a statistically significant weak negative association: "
          f"longer time series tend to produce lower β estimates. "
          f"This may reflect estimation artifacts in short series or genuine "
          f"differences between well-studied and poorly-studied stocks.")
    else:
        a("No significant association detected.")
    a("")

    a("## Finding 6: Multiple testing correction reduces significant fraction")
    a("")
    a(f"After Benjamini-Hochberg FDR correction (α = {ALPHA_SIG}):  ")
    a(f"**{cr['n_bh_significant']}** of {cr['n_stocks']} stocks "
      f"({cr['pct_bh_significant']:.1f}%) remain significantly non-proportional.  ")
    a(f"Compared to {cr['pct_significant']:.1f}% before correction.")
    a("")

    cpue_a = cr.get("cpue_analysis", {})
    if cpue_a.get("n_stocks", 0) >= 5:
        a("## Finding 7: Direct CPUE~B analysis (stocks with actual CPUE data)")
        a("")
        a(f"Only **{cpue_a['n_stocks']}** of {cr['n_stocks']} stocks have actual "
          f"CPUE index data in the database. For these stocks:  ")
        a(f"Mean β(CPUE) = {cpue_a['mean_beta']:.4f}, "
          f"Median = {cpue_a['median_beta']:.4f}, "
          f"SD = {cpue_a['sd_beta']:.4f}  ")
        a(f"The CPUE-based β values are closer to 1 than the catch-based values, "
          f"consistent with the expectation that catch confounds effort variation "
          f"with the CPUE-biomass relationship.")
        a("")

    a("## Note on central tendency measures")
    a("")
    a(f"The mean β ({cr['mean_beta']:.4f}), median β ({cr['median_beta']:.4f}), "
      f"and weighted mean β ({cr['weighted_mean_beta']:.4f}) diverge substantially. "
      f"This reflects a heavily left-skewed distribution with extreme negative outliers "
      f"(SD = {cr['sd_beta']:.4f}). The **median is the most robust summary statistic** "
      f"for this distribution. The mean is pulled down by a small number of stocks with "
      f"strongly negative catch-biomass slopes (likely reflecting quota-managed fisheries "
      f"where catch is decoupled from biomass). The weighted mean (by time series length) "
      f"gives disproportionate influence to these long, managed stocks.")
    a("")

    a("## Table 5: Residual autocorrelation diagnostic")
    a("")
    a(f"Median Durbin-Watson statistic: {cr['median_durbin_watson']:.3f}  ")
    a(f"Stocks with strong positive autocorrelation (DW < 1.0): "
      f"{cr['n_strong_autocorrelation']} ({cr['pct_strong_autocorrelation']:.1f}%)  ")
    a("Autocorrelation inflates apparent significance of per-stock OLS tests; "
      "the bootstrap CIs partially mitigate this but do not fully account for "
      "temporal dependence.")
    a("")

    # Top hyperstable
    a("## Table 6: Top 10 Most Hyperstable Stocks (lowest β)")
    a("")
    a("| Stock | β | 95% CI | R² | N years |")
    a("|-------|---|--------|----|---------|")
    srt = sorted(stock_results, key=lambda r: r["beta"])
    for r in srt[:10]:
        a(f"| {r['name'][:45]} | {r['beta']:.3f} | "
          f"[{r['ci_lower']:.2f}, {r['ci_upper']:.2f}] | "
          f"{r['r_squared']:.3f} | {r['n_years']} |")
    a("")

    # Top hyperdepletion
    a("## Table 7: Top 10 Most Hyperdepleted Stocks (highest β)")
    a("")
    a("| Stock | β | 95% CI | R² | N years |")
    a("|-------|---|--------|----|---------|")
    for r in reversed(srt[-10:]):
        a(f"| {r['name'][:45]} | {r['beta']:.3f} | "
          f"[{r['ci_lower']:.2f}, {r['ci_upper']:.2f}] | "
          f"{r['r_squared']:.3f} | {r['n_years']} |")
    a("")

    # Limitations section
    a("## Limitations")
    a("")
    a("1. **Catch is used as a proxy for CPUE.** The primary analysis fits "
      "log(Catch) ~ log(Biomass), not log(CPUE) ~ log(Biomass). Because "
      "Catch = CPUE × Effort, the estimated β conflates the true CPUE-biomass "
      "relationship with effort dynamics. When effort is reduced alongside "
      "declining biomass (e.g., via quota management), catch declines less than "
      "biomass, producing apparent hyperstability that may reflect management "
      "success rather than genuine CPUE bias.")
    a("")
    a("2. **Assessment model circularity.** Both catch and biomass estimates "
      "often come from the same stock assessment model that uses catch data as "
      "input. This creates non-independence between the regressor and response "
      "that can bias the estimated slope.")
    a("")
    a("3. **Temporal autocorrelation.** Fish stock time series exhibit strong "
      "positive autocorrelation (median DW statistic reported above). This "
      "inflates the effective sample size and can make OLS standard errors "
      "anti-conservative. Our bootstrap CIs partially address this but do not "
      "employ block-bootstrap methods that would fully account for serial "
      "dependence.")
    a("")
    a("4. **No multiple testing correction for classification.** With 866 "
      "simultaneous tests, some false discoveries are expected. We report "
      "Benjamini-Hochberg-adjusted results alongside raw classifications.")
    a("")
    a("5. **RAM Legacy selection bias.** The database over-represents commercially "
      "important stocks in data-rich regions (North Atlantic, North Pacific, "
      "Australia). Tropical, small-scale, and data-poor fisheries are under-represented. "
      "Results may not generalize to these contexts.")
    a("")
    a("6. **Non-independence of stocks.** Stocks sharing ecosystems, management "
      "regimes, or environmental drivers are not statistically independent. "
      "The permutation test treats stocks as exchangeable units, which may "
      "underestimate uncertainty.")
    a("")
    a("7. **Outlier exclusion.** Eight stocks with |β| > 10 were excluded from "
      "the meta-analysis as statistical outliers (likely reflecting extreme "
      "catch-biomass decoupling). No sensitivity analysis to the threshold was "
      "performed.")
    a("")
    a("8. **Sensitivity analyses use reduced bootstrap** (500 iterations for SSB "
      "analysis) and parametric CIs (for first-differences). These are disclosed "
      "here for transparency but may produce noisier estimates than the primary "
      "analysis with 2,000 bootstrap iterations.")
    a("")

    a("---")
    a(f"*Analysis date: {metadata.get('analysis_date', 'N/A')}*  ")
    a(f"*Random seed: {SEED} | Bootstrap: {N_BOOTSTRAP:,} (primary), "
      f"500 (SSB sensitivity) | Permutations: {N_PERMUTATIONS:,}*  ")
    a(f"*Data SHA256: {metadata.get('data_sha256', 'N/A')}*")

    with open(REPORT_FILE, "w") as f:
        f.write("\n".join(L) + "\n")


# ============================================================================
# VERIFICATION
# ============================================================================

def verify():
    """10 machine-checkable assertions."""
    log("\n=== VERIFICATION MODE ===\n")
    n_pass = 0
    n_total = 0

    def check(cond, msg):
        nonlocal n_pass, n_total
        n_total += 1
        status = "PASS" if cond else "FAIL"
        if cond:
            n_pass += 1
        log(f"  [VERIFY {n_total}/10] {status}: {msg}")

    check(os.path.exists(RESULTS_FILE), "results.json exists")
    check(os.path.exists(REPORT_FILE), "report.md exists")

    if not os.path.exists(RESULTS_FILE):
        log(f"\nVERIFICATION: {n_pass}/{n_total} passed")
        sys.exit(1)

    with open(RESULTS_FILE) as f:
        res = json.load(f)

    check("cross_stock_analysis" in res and "stock_results" in res,
          "results.json has expected structure")

    cr = res.get("cross_stock_analysis", {})
    sr = res.get("stock_results", [])

    n_stocks = cr.get("n_stocks", 0)
    check(n_stocks >= 50, f"N stocks >= 50 (got {n_stocks})")

    betas = [r.get("beta", 0) for r in sr]
    ok = all(math.isfinite(b) and -10 < b < 10 for b in betas) if betas else False
    check(ok, "All beta values finite and in (-10, 10)")

    has_ci = all("ci_lower" in r and "ci_upper" in r for r in sr)
    check(has_ci, "Bootstrap CIs computed for all stocks")

    pp = cr.get("permutation_p")
    check(pp is not None and 0 <= pp <= 1,
          f"Permutation p-value valid (p={pp})")

    cd = cr.get("cohens_d")
    check(cd is not None and math.isfinite(cd),
          f"Cohen's d computed (d={cd})")

    sens = res.get("sensitivity", {})
    check(len(sens) >= 3, f"At least 3 sensitivity analyses (got {len(sens)})")

    sha_file = os.path.join(DATA_DIR, "ram_legacy.xlsx.sha256")
    check(os.path.exists(sha_file), "Data SHA256 hash file exists")

    log(f"\n=== VERIFICATION: {n_pass}/{n_total} passed ===")
    if n_pass == n_total:
        log("ALL CHECKS PASSED")
    sys.exit(0 if n_pass == n_total else 1)


# ============================================================================
# MAIN
# ============================================================================

def main():
    if VERIFY:
        verify()
        return

    random.seed(SEED)
    t0 = time.time()

    log("=" * 70)
    log("CPUE-BIOMASS POWER-LAW ANALYSIS")
    log("Hyperstability & Hyperdepletion Across RAM Legacy Fish Stocks")
    log("=" * 70)

    # [1/8] Setup
    log("\n[1/8] Setup and configuration")
    os.makedirs(DATA_DIR, exist_ok=True)
    log(f"  Workspace: {WORKSPACE}")
    log(f"  Seed={SEED}  Bootstrap={N_BOOTSTRAP}  Perms={N_PERMUTATIONS}  MinYears={MIN_YEARS}")

    # [2/8] Data download
    log("\n[2/8] Downloading RAM Legacy Stock Assessment Database v4.44")
    xlsx_path = get_ram_data()

    # [3/8] Parse data
    log("\n[3/8] Parsing XLSX timeseries data")
    raw_rows = read_ram_timeseries(xlsx_path)
    log(f"  Raw records: {len(raw_rows):,}")

    # [4/8] Process
    log("\n[4/8] Building stock-level time series")
    stocks = build_stock_data(raw_rows)
    n_stocks = len(stocks)
    total_yrs = sum(s["n_years"] for s in stocks.values())
    med_len = statistics.median([s["n_years"] for s in stocks.values()])
    log(f"  Stocks with >= {MIN_YEARS} years: {n_stocks}")
    log(f"  Total stock-years: {total_yrs:,}")
    log(f"  Median series length: {med_len:.0f} years")

    if n_stocks < 10:
        log("ERROR: Too few stocks for analysis. Check data parsing.")
        sys.exit(1)

    # [5/8] Primary analysis
    log("\n[5/8] Per-stock power-law fits with bootstrap CIs")
    stock_results = []
    for i, (sid, data) in enumerate(sorted(stocks.items())):
        r = analyze_stock(sid, data, SEED + i)
        if r:
            stock_results.append(r)
        if (i + 1) % 50 == 0 or i + 1 == n_stocks:
            log(f"  ... {i+1}/{n_stocks} stocks ({len(stock_results)} valid)")

    log(f"  Analyzed: {len(stock_results)} stocks")

    # Filter extreme outliers (|beta| > 10) for robust meta-analysis
    n_outliers = sum(1 for r in stock_results if abs(r["beta"]) > 10)
    if n_outliers > 0:
        log(f"  Excluding {n_outliers} stocks with |beta| > 10 as outliers")
        stock_results = [r for r in stock_results if abs(r["beta"]) <= 10]
        log(f"  Retained: {len(stock_results)} stocks")

    # [6/8] Cross-stock meta-analysis
    log("\n[6/8] Cross-stock meta-analysis")
    cross = cross_stock_analysis(stock_results)
    log(f"  Mean beta  = {cross['mean_beta']:.4f}  95% CI [{cross['mean_ci'][0]:.4f}, {cross['mean_ci'][1]:.4f}]")
    log(f"  Median beta= {cross['median_beta']:.4f}  95% CI [{cross['median_ci'][0]:.4f}, {cross['median_ci'][1]:.4f}]")
    log(f"  Permutation p = {cross['permutation_p']:.4f}")
    log(f"  Wilcoxon    p = {cross['wilcoxon_p']:.4f}")
    log(f"  Cohen's d     = {cross['cohens_d']:.4f}")
    log(f"  Hyperstable:    {cross['n_hyperstable']} ({cross['pct_hyperstable']:.1f}%)")
    log(f"  Proportional:   {cross['n_proportional']} ({cross['pct_proportional']:.1f}%)")
    log(f"  Hyperdepletion: {cross['n_hyperdepletion']} ({cross['pct_hyperdepletion']:.1f}%)")

    # [7/8] Sensitivity analysis
    log("\n[7/8] Sensitivity analysis")
    sens = sensitivity_analysis(stocks, stock_results)
    for key, val in sens.items():
        if "mean_beta" in val:
            log(f"  {val.get('label', key)}: N={val['n']}, "
                f"mean β={val['mean_beta']:.4f}, median β={val['median_beta']:.4f}")

    # [8/8] Output
    log("\n[8/8] Writing results")
    meta = {
        "database": "RAM Legacy Stock Assessment Database v4.44",
        "zenodo_record": ZENODO_RECORD,
        "zenodo_doi": "10.5281/zenodo.2542919",
        "data_sha256": sha256_file(xlsx_path),
        "analysis_date": datetime.now().isoformat(),
        "python_version": sys.version.split()[0],
        "random_seed": SEED,
        "n_bootstrap": N_BOOTSTRAP,
        "n_permutations": N_PERMUTATIONS,
        "min_years": MIN_YEARS,
        "n_raw_records": len(raw_rows),
        "n_stocks_analyzed": len(stock_results),
        "total_stock_years": total_yrs,
        "runtime_seconds": round(time.time() - t0, 2),
    }

    write_results(meta, stock_results, cross, sens)
    log(f"  results.json: {os.path.getsize(RESULTS_FILE):,} bytes")

    write_report(meta, stock_results, cross, sens)
    log(f"  report.md:    {os.path.getsize(REPORT_FILE):,} bytes")

    elapsed = time.time() - t0
    log(f"\nRuntime: {elapsed:.1f} seconds")
    log("\nANALYSIS COMPLETE")


if __name__ == "__main__":
    main()
SCRIPT_EOF
```

**Expected output:** File `analyze.py` written to workspace. Exit code 0.

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_cpue-as-abundance-index-hyperstability-and-hyperdepletion-ac && python3 analyze.py
```

**Expected output:**
- Progress markers `[1/8]` through `[8/8]`
- Data download from Zenodo with SHA256 verification
- XLSX parsing with row/column counts
- Per-stock analysis progress
- Summary statistics (mean β, median β, p-values)
- Files created: `results.json`, `report.md`, `data/ram_legacy.xlsx`, `data/ram_legacy.xlsx.sha256`
- Final line: `ANALYSIS COMPLETE`

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_cpue-as-abundance-index-hyperstability-and-hyperdepletion-ac && python3 analyze.py --verify
```

**Expected output:**
- 10 verification checks, all PASS
- `ALL CHECKS PASSED`
- Exit code 0

## Success Criteria

1. Script downloads RAM Legacy Database v4.44 from Zenodo without errors
2. XLSX parsed successfully with correct column mapping
3. At least 50 stocks analyzed (expected: 200-400)
4. Bootstrap CIs computed for all stocks (2,000 iterations each)
5. Permutation test completed (10,000 iterations)
6. Sensitivity analyses for ≥3 settings
7. results.json and report.md written with complete data
8. All 10 verification checks pass
9. Total runtime under 30 minutes

## Failure Conditions

1. Network failure: Script prints clear error with manual download instructions
2. XLSX parse failure: Script prints available sheet names and column headers for debugging
3. Too few stocks (<10): Script exits with error explaining the issue
4. Bootstrap failure: Falls back to parametric CI using SE
5. Any unhandled exception: Full traceback printed to stderr