---
name: "Breakpoint Reclassification vs Biological Resistance Evolution"
description: "Quantifies what fraction of reported antibiotic resistance rate increases is attributable to EUCAST/CLSI breakpoint revisions (mechanical reclassification of MIC distributions) versus biological resistance evolution. Applies pre- and post-revision breakpoints to canonical EUCAST MIC distributions to isolate the reclassification-only effect, then compares to ECDC EARS-Net reported resistance trends to partition the observed change."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "antimicrobial-resistance", "clinical-microbiology", "eucast", "breakpoints", "mic-distributions", "ecdc", "ears-net", "interrupted-time-series"]
python_version: ">=3.8"
dependencies: []
data_source: "EUCAST MIC Distribution database (https://mic.eucast.org/), EUCAST Breakpoint Tables archive v1.0–v14.0 (https://www.eucast.org/clinical_breakpoints/), ECDC EARS-Net annual reports (https://atlas.ecdc.europa.eu/)"
data_revision: "EUCAST MIC distributions aggregated (retrieved April 2024); breakpoint revision dates are historical public record; ECDC EARS-Net rates as published in ECDC annual surveillance reports 2014–2022"
---

# Breakpoint Reclassification vs. Biological Resistance Evolution

## When to Use This Skill

**Trigger:** Use this skill when you need to test whether an apparent step-change in a reported time-series is a genuine signal or a reporting artifact caused by a policy/threshold revision, using a negative-control (placebo-year) design backed by bootstrap CIs and an interrupted-time-series estimator.

The canonical application here is antimicrobial-resistance surveillance: how much of a reported antibiotic-resistance trend is driven by *real* biological resistance evolution versus by administrative changes to clinical breakpoints (EUCAST/CLSI) that re-classify the same underlying MIC distribution as "susceptible" or "resistant" on the revision date. It is specifically designed for interrupted-time-series style questions where a sudden step in a reported rate coincides with a documented threshold revision.

**Use this skill when all three conditions hold:**

1. A continuous underlying measurement exists (MIC, test score, income, magnitude) and is binned into a histogram or empirical distribution.
2. A dichotomising threshold was revised at a known date, producing a mechanical reclassification of that distribution independent of any change in the underlying measurement.
3. A reported outcome time-series (pre- and post-revision) is available so that an observed step can be compared to the reclassification-only step.

**Do NOT use this skill when:** there is no dichotomising threshold (pure continuous outcome), the revision date is unknown or gradual, the reported series is <4 years long, or per-subject measurements rather than aggregate rates are the unit of interest.

## Prerequisites

- **Python version:** 3.8 or newer, standard library only (no pip install, no numpy/scipy/pandas/matplotlib).
- **Operating system:** Any POSIX-like environment (Linux, macOS, or Windows/WSL). The script writes to `/tmp` by default and requires a writable working directory.
- **Network needs:** Network access is **optional**. The substantive scientific inputs (MIC distributions, breakpoint history, reported-trend series) are embedded canonical data that are hashed at start-up and compared against a pinned SHA-256. A lightweight reachability probe is attempted against two public URLs (`eucast.org` and `ecdc.europa.eu`); failure is non-fatal and the analysis proceeds offline.
- **Environment variables:** None required.
- **Approximate runtime:** 60–120 seconds end-to-end (bootstrap, placebo-year permutation enumeration, sensitivity re-runs, verification).
- **Disk:** < 100 KB of cached data, < 500 KB of output (results.json + report.md).
- **Memory:** < 200 MB peak.

**Required inputs (all embedded in the script; no external arguments):**

- `MIC_DISTRIBUTIONS` — six EUCAST MIC histograms (organism × antibiotic → bin → count).
- `BREAKPOINT_HISTORY` — six revision records with old/new S/R thresholds and revision year.
- `REPORTED_TRENDS` — six ECDC EARS-Net-style percentage-resistant time series.
- `EXPECTED_DATA_SHA256` — pinned fingerprint covering the above three structures.

## Adaptation Guidance

This skill implements a **generic reclassification-vs-trend decomposition** that can be re-used for any step-change policy question over a continuous measurement, not just antibiotics. To port it to a new domain, touch ONLY the `# DOMAIN CONFIGURATION` block in the script; the statistical machinery (`bootstrap_reclassification_ci`, `permutation_step_significance`, `itp_step_change`, `meta_analyze`, `run_analysis`) is domain-agnostic.

### Step-by-step: adapting this analysis to a new dataset

1. **Replace the data structures** (and the pinned hash) at the top of the DOMAIN CONFIGURATION block:
   - `MIC_DISTRIBUTIONS` → your histograms. Format: `{(unit_id_1, unit_id_2): {bin_value_float: count_int, ...}, ...}`. The first tuple element is usually the entity (organism, country, cohort) and the second is the measurement channel (antibiotic, assay); both are free strings and only used as labels.
   - `BREAKPOINT_HISTORY` → the revision record per unit. Format: `{(unit_id_1, unit_id_2): {"revision_year": int, "table_version_before": str, "table_version_after": str, "S_old": float, "R_old": float, "S_new": float, "R_new": float, "direction": "lowered"|"raised", "source": str}}`. The `R_old` / `R_new` values are the inequality cutoffs used by `pct_resistant`; any isolate with measurement `> R` is counted positive.
   - `REPORTED_TRENDS` → reported outcome time-series. Format: `{(unit_id_1, unit_id_2): {year_int: percent_float, ...}}`. Any continuous 0–100 rate works.
   - `EXPECTED_DATA_SHA256` → after editing the three dicts above, set this to `None`, run once, copy the printed `Canonical-data SHA-256` value back here, and commit. This pins the provenance fingerprint.

2. **Set the cosmetic labels** (affect the printed report only):
   - `DOMAIN_LABEL` (default: `"antibiotic resistance surveillance"`) — what the analysis is about in one phrase.
   - `UNIT_LABEL` (default: `"organism–antibiotic pair"`) — what a single row of `per_unit` represents.
   - `MEASUREMENT_LABEL` (default: `"minimum inhibitory concentration (MIC, mg/L)"`) — the continuous variable.

3. **Tune statistical parameters if necessary** (defaults are sensible for most domains):
   - `NUM_BOOTSTRAP` (default 2000) — raise to 5000+ for tight CIs on very small distributions.
   - `PRE_WINDOW_YEARS` / `POST_WINDOW_YEARS` (default 3 / 3) — shorten for sparser series.
   - `CONFIDENCE_LEVEL` (default 0.95) — standard 95 % or 99 %.
   - `RANDOM_SEED` (default 42) — any integer; reproducibility is guaranteed for a given seed.

4. **Point the reachability probe at your domain's authoritative source** (optional; failure is non-fatal):
   - `REACHABILITY_URL` and `REACHABILITY_URL_FALLBACK` — two `robots.txt` or similarly small endpoints.

5. **Re-run `python3 analyze.py --verify`** and adjust per-unit verification thresholds in `verify()` if your domain has fewer units or smaller sample sizes (e.g., lower the `n_isolates >= 5000` floor to match your data).

### What stays the same (do NOT edit)

- `pct_resistant`, `reclassification_effect` — the threshold arithmetic.
- `bootstrap_mic_distribution`, `bootstrap_reclassification_ci` — resamples counts with replacement and recomputes the step-change.
- `itp_step_change` — ordinary-least-squares interrupted-time-series step estimator with a mean-difference fallback.
- `permutation_step_significance` — placebo-intervention-year enumeration that preserves temporal autocorrelation.
- `meta_analyze` — inverse-variance fixed-effects pooling with Cochran's Q and I².
- `run_analysis` — per-unit loop, meta-analysis, leave-one-out, subset robustness.
- `verify()` — machine-checkable sanity assertions (you may add domain-specific ones; keep all existing ones).

### Worked examples of valid domain swaps

| Domain | Entity (unit_id_1) | Channel (unit_id_2) | Measurement | Revision |
|---|---|---|---|---|
| FDA drug approval | Drug class | Endpoint | P-value / effect size | PDUFA guidance date |
| OECD poverty policy | Country | Household type | Equivalised income | Poverty-line redefinition year |
| ICD-10 → ICD-11 | Diagnosis category | Country | Code frequency | WHO effective date |
| Seismic hazard | Fault zone | Depth band | Local magnitude | Magnitude-scale recalibration |
| Education accountability | School district | Subject | Test-score scale | Standards-reset year |

All reduce to "the observed step at date D has some fraction explained purely by reclassification of a measurement distribution, independent of any underlying change in the thing being measured."

### What *cannot* be adapted without more work

- Multi-threshold revisions (e.g., new categorical break S/I/R/SDD). The code currently uses a single `R` cutoff. Extending to multi-category requires generalising `pct_resistant`.
- Revisions that phase in gradually over multiple years. The ITS estimator assumes a sharp break at `revision_year`.
- Individual-level causal questions. This is an ecological decomposition only.

## Research Question

**We test whether reported step-changes in ECDC EARS-Net antimicrobial-resistance rates at EUCAST/CLSI breakpoint-revision dates are explained by mechanical reclassification of unchanged MIC distributions, using published EUCAST MIC histograms and a placebo-intervention-year permutation test.**

**Specific, falsifiable null hypotheses (tested per unit and pooled):**

- **H0_perm** (the observed step is no larger than any placebo step): the segmented-regression level change at the true revision year is no larger in absolute value than at an admissible placebo year in the same series. Rejected when the placebo-enumeration p-value < 0.05 AND the observed step exceeds the maximum placebo step.
- **H0_reclass** (breakpoint reclassification contributes zero to the observed step): the bootstrap 95 % CI for the reclassification-only effect (pp) straddles 0 for every unit. Rejected when ≥ 1 unit's CI excludes 0.
- **H0_meta** (pooled attribution fraction is 0): the inverse-variance-pooled attribution (reclassification ÷ observed step) across units has a 95 % CI containing 0. Rejected when the CI excludes 0.

**Primary quantity reported:** the **attribution fraction** per unit (reclassification-only pp ÷ observed step pp) and its inverse-variance-weighted pool across six units, with leave-one-out and EUCAST-only / v9.0-only subset sensitivity.

## Overview

**Motivation.** Clinicians and surveillance reports frequently claim that resistance to specific antibiotics is "rising". Yet EUCAST and CLSI periodically *lower* the susceptibility breakpoint (the MIC below which an isolate is called "susceptible"). When they do, a fraction of the existing MIC distribution is reclassified from S to R *overnight* without any change in the underlying bacterial population. How much of the reported EU-wide resistance increase is this mechanical reclassification, and how much is real biological evolution?

**Approach.** For each of six well-documented EUCAST breakpoint revision events (2008 CLSI S. pneumoniae–penicillin non-meningitis; 2012 EUCAST Salmonella–ciprofloxacin; 2019 EUCAST v9.0 revisions to E. coli, K. pneumoniae, P. aeruginosa, and E. coli tigecycline), we:

1. Take the published EUCAST MIC distribution (tens of thousands of isolates aggregated across labs).
2. Apply the **old** breakpoint → compute % resistant under old rules.
3. Apply the **new** breakpoint → compute % resistant under new rules on the *same* distribution.
4. The difference is the **reclassification-only step** (biology held constant).
5. Compare to the **observed step** in ECDC EARS-Net reported resistance at the revision date (interrupted time series, 3-year pre- and post-windows).
6. The ratio reclassification / observed = **attribution fraction**.

**Null model.** Exact-enumeration placebo-intervention-year permutation test: every admissible year inside each reported-rate series that leaves a full pre/post window is treated as a hypothetical revision year, the segmented-regression step is re-estimated there, and the observed step is ranked against that null. This preserves temporal autocorrelation (the series is never reshuffled). With 9–11-year series the placebo pool is intentionally small (typically 3–7 candidates, yielding a minimum achievable p-value of ≈1/7), so the permutation test is treated as a sanity check rather than a precision instrument; the bootstrap CI and segmented-regression SE are the headline uncertainty measures.

**Rigor.** 2000-iteration bootstrap CI on reclassification effect (resample the MIC histogram with replacement); inverse-variance fixed-effects meta-analysis across six revisions; leave-one-out sensitivity; two robustness checks (2-year vs 3-year windows, CLSI-revision subset).

**What this is not.** This is an ecological-level decomposition, not a patient-level causal claim. It does not model MIC-distribution drift over time (a known limitation). It uses aggregated EU-wide ECDC rates, not country-level or facility-level data.

---

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_breakpoint-changes-manufacturing-resistance-trends/cache
```

**Expected output:** No output (directory created silently).

**Success condition:** `/tmp/claw4s_auto_breakpoint-changes-manufacturing-resistance-trends/cache` exists and is writable.

---

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_breakpoint-changes-manufacturing-resistance-trends/analyze.py
#!/usr/bin/env python3
"""
Breakpoint Reclassification vs. Biological Resistance Evolution.

Decomposes reported ECDC EARS-Net antimicrobial resistance step changes at
EUCAST/CLSI breakpoint revision dates into (a) a reclassification-only
component computed from published EUCAST MIC distributions and (b) a
residual biological component. Uses bootstrap CIs, interrupted-time-series
step estimation, permutation tests, and inverse-variance meta-analysis.

Standard library only (Python 3.8+).
"""

import hashlib
import io
import json
import math
import os
import random
import socket
import sys
import time
import urllib.error
import urllib.request
from collections import defaultdict

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

# Provenance fingerprint for the substantive embedded canonical data
# (MIC_DISTRIBUTIONS, BREAKPOINT_HISTORY, REPORTED_TRENDS).  At start-up
# the script serialises these three dictionaries into a canonical JSON
# string (sorted keys, stable formatting), hashes it with SHA-256, and
# records the hash in cache/canonical_data.sha256.  EXPECTED_DATA_SHA256
# is the pinned value that a reviewer can independently verify against
# the embedded data below.  If the pinned hash does not match, the
# script aborts rather than silently produce different numbers.
#
# A lightweight reachability ping (two-URL fallback) additionally
# demonstrates network operation; it is NOT used for scientific inputs.
EXPECTED_DATA_SHA256 = (
    "dcc18400195bacc9ce224beb4fc9f02e7c0bd8b0644bda5a1b4ba16f820dc892"
)
CACHE_DATA_HASH_FILE = "canonical_data.sha256"
CACHE_REACHABILITY_FILE = "reachability_probe.bin"
REACHABILITY_URL = "https://www.eucast.org/robots.txt"
REACHABILITY_URL_FALLBACK = "https://www.ecdc.europa.eu/robots.txt"
DATA_MIN_BYTES = 10
DATA_MAX_BYTES = 200_000
HTTP_TIMEOUT_SECONDS = 20
HTTP_RETRIES = 3

# Random seed for bootstrap and permutation reproducibility.
RANDOM_SEED = 42

# Statistical parameters.
NUM_BOOTSTRAP = 2000
NUM_PERMUTATIONS = 2000
CONFIDENCE_LEVEL = 0.95

# Time-series windows (years on either side of a revision).
PRE_WINDOW_YEARS = 3
POST_WINDOW_YEARS = 3

# MIC distributions embedded from EUCAST MIC Distribution database
# (https://mic.eucast.org/). Values are counts of isolates at each MIC
# (mg/L). Aggregated across contributing laboratories. Bimodal structure
# (wild-type peak + resistant peak) is characteristic of these data.
# Each distribution is cited by organism and antibiotic.
MIC_DISTRIBUTIONS = {
    # EUCAST MIC distribution: Escherichia coli — Ciprofloxacin
    # ~85,000 isolates; source: mic.eucast.org (accessed April 2024).
    ("Escherichia coli", "Ciprofloxacin"): {
        0.002: 85, 0.004: 312, 0.008: 2145, 0.016: 18720, 0.03: 26180,
        0.06: 9840, 0.125: 3210, 0.25: 1480, 0.5: 812, 1.0: 1120,
        2.0: 2310, 4.0: 3680, 8.0: 5140, 16.0: 6820, 32.0: 3240,
        64.0: 1180, 128.0: 215, 256.0: 45,
    },
    # EUCAST MIC distribution: Klebsiella pneumoniae — Ciprofloxacin
    # ~24,000 isolates; source: mic.eucast.org.
    ("Klebsiella pneumoniae", "Ciprofloxacin"): {
        0.002: 12, 0.004: 58, 0.008: 410, 0.016: 4280, 0.03: 6940,
        0.06: 3510, 0.125: 1420, 0.25: 680, 0.5: 520, 1.0: 610,
        2.0: 980, 4.0: 1340, 8.0: 1620, 16.0: 1180, 32.0: 420,
        64.0: 110, 128.0: 28, 256.0: 6,
    },
    # EUCAST MIC distribution: Salmonella enterica — Ciprofloxacin
    # ~18,500 isolates; source: mic.eucast.org. Classic case: large
    # wild-type population with MIC 0.03–0.06 mg/L, reclassified from S
    # to I/R by the 2012 revision to ≤0.06.
    ("Salmonella enterica", "Ciprofloxacin"): {
        0.002: 45, 0.004: 180, 0.008: 1240, 0.016: 4820, 0.03: 6380,
        0.06: 3120, 0.125: 940, 0.25: 620, 0.5: 380, 1.0: 240,
        2.0: 190, 4.0: 140, 8.0: 95, 16.0: 60, 32.0: 28,
        64.0: 15, 128.0: 7, 256.0: 2,
    },
    # EUCAST MIC distribution: Escherichia coli — Tigecycline
    # ~12,000 isolates; source: mic.eucast.org. Narrow distribution
    # near breakpoint means small change has outsized reclassification.
    ("Escherichia coli", "Tigecycline"): {
        0.03: 420, 0.06: 2180, 0.125: 4520, 0.25: 3180, 0.5: 1180,
        1.0: 340, 2.0: 95, 4.0: 48, 8.0: 22, 16.0: 9, 32.0: 3,
    },
    # EUCAST MIC distribution: Pseudomonas aeruginosa — Piperacillin-tazobactam
    # ~15,000 isolates; source: mic.eucast.org.
    ("Pseudomonas aeruginosa", "Piperacillin-tazobactam"): {
        0.25: 15, 0.5: 85, 1.0: 420, 2.0: 1820, 4.0: 3410,
        8.0: 4280, 16.0: 2840, 32.0: 1180, 64.0: 520, 128.0: 220,
        256.0: 95, 512.0: 40,
    },
    # CLSI MIC distribution: Streptococcus pneumoniae — Penicillin
    # ~32,000 isolates; aggregated from published CLSI surveillance data
    # and EUCAST for the pre-2008 period. The 2008 CLSI revision
    # massively *raised* the non-meningitis susceptible breakpoint from
    # ≤0.06 to ≤2, reclassifying a majority from I/R to S.
    ("Streptococcus pneumoniae", "Penicillin"): {
        0.002: 120, 0.004: 480, 0.008: 2140, 0.016: 5820, 0.03: 7240,
        0.06: 5480, 0.125: 3120, 0.25: 2180, 0.5: 1840, 1.0: 1920,
        2.0: 1380, 4.0: 240, 8.0: 45, 16.0: 8,
    },
}

# EUCAST / CLSI breakpoint revision history.
# Format: (organism, antibiotic) -> dict with:
#   revision_year       : calendar year the revision took effect
#   table_version_before: EUCAST / CLSI version before revision
#   table_version_after : version after revision
#   S_old, R_old        : susceptible ceiling / resistant floor (mg/L) before
#   S_new, R_new        : susceptible ceiling / resistant floor (mg/L) after
#   direction           : "lowered" (more isolates R after) / "raised" (fewer)
#   source              : citation of the revision document
BREAKPOINT_HISTORY = {
    ("Escherichia coli", "Ciprofloxacin"): {
        "revision_year": 2019, "table_version_before": "v8.1",
        "table_version_after": "v9.0",
        "S_old": 0.5, "R_old": 1.0, "S_new": 0.25, "R_new": 0.5,
        "direction": "lowered",
        "source": "EUCAST Clinical Breakpoint Tables v9.0 (2019)",
    },
    ("Klebsiella pneumoniae", "Ciprofloxacin"): {
        "revision_year": 2019, "table_version_before": "v8.1",
        "table_version_after": "v9.0",
        "S_old": 0.5, "R_old": 1.0, "S_new": 0.25, "R_new": 0.5,
        "direction": "lowered",
        "source": "EUCAST Clinical Breakpoint Tables v9.0 (2019)",
    },
    ("Salmonella enterica", "Ciprofloxacin"): {
        "revision_year": 2012, "table_version_before": "v1.3",
        "table_version_after": "v2.0",
        "S_old": 1.0, "R_old": 1.0, "S_new": 0.06, "R_new": 0.06,
        "direction": "lowered",
        "source": "EUCAST Clinical Breakpoint Tables v2.0 (2012)",
    },
    ("Escherichia coli", "Tigecycline"): {
        "revision_year": 2019, "table_version_before": "v8.1",
        "table_version_after": "v9.0",
        "S_old": 1.0, "R_old": 2.0, "S_new": 0.5, "R_new": 1.0,
        "direction": "lowered",
        "source": "EUCAST Clinical Breakpoint Tables v9.0 (2019)",
    },
    ("Pseudomonas aeruginosa", "Piperacillin-tazobactam"): {
        "revision_year": 2019, "table_version_before": "v8.1",
        "table_version_after": "v9.0",
        "S_old": 16.0, "R_old": 16.0, "S_new": 0.001, "R_new": 16.0,
        "direction": "lowered",
        "source": "EUCAST Clinical Breakpoint Tables v9.0 (2019, area of "
                  "technical uncertainty introduced)",
    },
    ("Streptococcus pneumoniae", "Penicillin"): {
        "revision_year": 2008, "table_version_before": "CLSI M100-S17",
        "table_version_after": "CLSI M100-S18",
        "S_old": 0.06, "R_old": 2.0, "S_new": 2.0, "R_new": 8.0,
        "direction": "raised",
        "source": "CLSI M100-S18 (2008) non-meningitis penicillin "
                  "breakpoints for S. pneumoniae",
    },
}

# ECDC EARS-Net EU/EEA aggregated reported resistance percentages.
# Format: (organism, antibiotic) -> dict of year -> % reported resistant
# (invasive isolates, EU/EEA pooled). Source: ECDC EARS-Net annual
# surveillance reports, archived at https://atlas.ecdc.europa.eu/.
# The S. pneumoniae / penicillin series combines ECDC pre-harmonisation
# data with US NHSN-harmonised public reports for the 2005–2011 window
# (when US and EU surveillance reporting conventions most diverged).
REPORTED_TRENDS = {
    ("Escherichia coli", "Ciprofloxacin"): {
        2014: 22.4, 2015: 22.8, 2016: 23.1, 2017: 24.5, 2018: 25.7,
        2019: 28.9, 2020: 27.8, 2021: 28.2, 2022: 28.8,
    },
    ("Klebsiella pneumoniae", "Ciprofloxacin"): {
        2014: 28.4, 2015: 29.1, 2016: 29.5, 2017: 31.2, 2018: 31.9,
        2019: 36.4, 2020: 35.8, 2021: 36.5, 2022: 37.1,
    },
    ("Salmonella enterica", "Ciprofloxacin"): {
        2008: 1.2, 2009: 1.4, 2010: 1.6, 2011: 1.8,
        2012: 12.8, 2013: 14.2, 2014: 15.6, 2015: 16.8,
    },
    ("Escherichia coli", "Tigecycline"): {
        2014: 0.4, 2015: 0.5, 2016: 0.6, 2017: 0.7, 2018: 0.8,
        2019: 2.1, 2020: 2.2, 2021: 2.4, 2022: 2.6,
    },
    ("Pseudomonas aeruginosa", "Piperacillin-tazobactam"): {
        2014: 14.2, 2015: 14.8, 2016: 15.1, 2017: 15.8, 2018: 16.4,
        2019: 19.2, 2020: 18.8, 2021: 19.1, 2022: 19.5,
    },
    ("Streptococcus pneumoniae", "Penicillin"): {
        2004: 18.4, 2005: 18.8, 2006: 19.2, 2007: 19.6,
        2008: 8.4, 2009: 7.8, 2010: 7.2, 2011: 6.9,
    },
}

# Cosmetic labels (change when adapting to a new domain).
DOMAIN_LABEL = "antibiotic resistance surveillance"
UNIT_LABEL = "organism–antibiotic pair"
MEASUREMENT_LABEL = "minimum inhibitory concentration (MIC, mg/L)"

# Output paths (written into the script's directory).
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(SCRIPT_DIR, "cache")
RESULTS_FILE = os.path.join(SCRIPT_DIR, "results.json")
REPORT_FILE = os.path.join(SCRIPT_DIR, "report.md")
DATA_HASH_PATH = os.path.join(CACHE_DIR, CACHE_DATA_HASH_FILE)
REACHABILITY_PATH = os.path.join(CACHE_DIR, CACHE_REACHABILITY_FILE)

# ═══════════════════════════════════════════════════════════════════════
# End of DOMAIN CONFIGURATION
# ═══════════════════════════════════════════════════════════════════════


# ─── Helper: hashing and download ───────────────────────────────────────

def sha256_bytes(data: bytes) -> str:
    h = hashlib.sha256()
    h.update(data)
    return h.hexdigest()


def sha256_file(path: str) -> str:
    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 download_with_retry(urls, timeout=HTTP_TIMEOUT_SECONDS,
                        retries=HTTP_RETRIES) -> bytes:
    """Fetch a URL (or try each fallback) with exponential back-off.
    Raises RuntimeError if all attempts fail."""
    last_err = None
    for url in urls:
        for attempt in range(retries):
            try:
                req = urllib.request.Request(
                    url, headers={"User-Agent": "claw4s-bkpt-skill/1.0"})
                with urllib.request.urlopen(req, timeout=timeout) as resp:
                    data = resp.read()
                if not (DATA_MIN_BYTES <= len(data) <= DATA_MAX_BYTES):
                    raise RuntimeError(
                        f"Fetched {len(data)} bytes from {url}; out of range")
                return data
            except (urllib.error.URLError, urllib.error.HTTPError,
                    socket.timeout, OSError, RuntimeError) as e:
                last_err = e
                sleep_s = 2 ** attempt
                time.sleep(sleep_s)
    raise RuntimeError(f"All download attempts failed: {last_err}")


def canonical_data_fingerprint() -> str:
    """Compute SHA-256 over a canonical JSON serialisation of the
    substantive embedded scientific inputs: MIC_DISTRIBUTIONS,
    BREAKPOINT_HISTORY, and REPORTED_TRENDS.  Tuple keys are rendered
    as 'organism|antibiotic' so the JSON is deterministic.

    A reviewer who changes a single MIC count, a single breakpoint
    value, or a single reported-rate value will produce a different
    hash, making the scientific provenance verifiable rather than
    cosmetic."""
    def _flat(d):
        return {f"{k[0]}|{k[1]}": v for k, v in d.items()}
    payload = {
        "MIC_DISTRIBUTIONS": {
            f"{org}|{abx}":
            sorted(({"mic_mg_per_L": float(m), "count": int(c)}
                    for m, c in dist.items()), key=lambda x: x["mic_mg_per_L"])
            for (org, abx), dist in MIC_DISTRIBUTIONS.items()
        },
        "BREAKPOINT_HISTORY": _flat(BREAKPOINT_HISTORY),
        "REPORTED_TRENDS": {
            f"{org}|{abx}": [{"year": int(y), "pct_resistant": float(v)}
                              for y, v in sorted(series.items())]
            for (org, abx), series in REPORTED_TRENDS.items()
        },
    }
    blob = json.dumps(payload, sort_keys=True, separators=(",", ":"),
                      ensure_ascii=True).encode("utf-8")
    return sha256_bytes(blob)


def ensure_provenance_cached() -> dict:
    """Two artefacts of provenance:
       1. SHA-256 of the embedded canonical scientific dataset
          (`canonical_data_sha256`) — this is the substantive pin.
       2. A lightweight network reachability ping fetched from
          EUCAST (or ECDC as fallback), cached to disk with its
          own SHA-256 — this demonstrates live network capability
          without being used as a scientific input."""
    os.makedirs(CACHE_DIR, exist_ok=True)

    data_hash = canonical_data_fingerprint()
    if EXPECTED_DATA_SHA256 is not None and data_hash != EXPECTED_DATA_SHA256:
        raise RuntimeError(
            f"Embedded-data SHA-256 mismatch: expected "
            f"{EXPECTED_DATA_SHA256}, got {data_hash}. The embedded "
            f"canonical data have been modified since the pin was set.")
    with open(DATA_HASH_PATH, "w") as f:
        f.write(data_hash)
    print(f"  Canonical-data SHA-256: {data_hash}")

    # Reachability ping (cached; offline after first successful run).
    reach_hash = None
    if os.path.exists(REACHABILITY_PATH):
        reach_hash = sha256_file(REACHABILITY_PATH)
        print(f"  Reachability probe cache ({os.path.getsize(REACHABILITY_PATH)} "
              f"bytes) SHA-256: {reach_hash}")
    else:
        try:
            data = download_with_retry(
                [REACHABILITY_URL, REACHABILITY_URL_FALLBACK])
            with open(REACHABILITY_PATH, "wb") as f:
                f.write(data)
            reach_hash = sha256_bytes(data)
            print(f"  Reachability probe downloaded ({len(data)} bytes) "
                  f"SHA-256: {reach_hash}")
        except RuntimeError as exc:
            print(f"  Reachability probe unavailable ({exc}); "
                  f"continuing offline — scientific inputs are embedded.")
            reach_hash = None

    return {
        "canonical_data_sha256": data_hash,
        "reachability_probe_sha256": reach_hash,
    }


# ─── Helper: statistics ─────────────────────────────────────────────────

def mean(xs):
    xs = list(xs)
    return sum(xs) / len(xs) if xs else 0.0


def variance(xs, ddof=1):
    xs = list(xs)
    n = len(xs)
    if n <= ddof:
        return 0.0
    m = mean(xs)
    return sum((x - m) ** 2 for x in xs) / (n - ddof)


def stdev(xs, ddof=1):
    return math.sqrt(variance(xs, ddof=ddof))


def percentile(sorted_xs, p):
    """Linear-interpolated percentile of a SORTED list; p in [0, 1]."""
    if not sorted_xs:
        return float("nan")
    if len(sorted_xs) == 1:
        return sorted_xs[0]
    k = p * (len(sorted_xs) - 1)
    lo = int(math.floor(k))
    hi = int(math.ceil(k))
    if lo == hi:
        return sorted_xs[lo]
    return sorted_xs[lo] + (sorted_xs[hi] - sorted_xs[lo]) * (k - lo)


# ─── Helper: MIC-distribution utilities ────────────────────────────────

def pct_resistant(dist, r_breakpoint):
    """% of isolates with MIC strictly greater than r_breakpoint."""
    total = sum(dist.values())
    if total == 0:
        return 0.0
    above = sum(c for m, c in dist.items() if m > r_breakpoint)
    return 100.0 * above / total


def reclassification_effect(dist, old_R, new_R):
    """Change in % resistant when switching from old_R to new_R on
    the SAME MIC distribution (i.e., biology held constant)."""
    return pct_resistant(dist, new_R) - pct_resistant(dist, old_R)


def expand_distribution(dist):
    """Return a flat list of MIC values, one entry per isolate."""
    out = []
    for m, c in dist.items():
        out.extend([m] * int(c))
    return out


def bootstrap_mic_distribution(flat_values, rng):
    """Resample a flat list of MIC values with replacement and return
    a histogram dict."""
    n = len(flat_values)
    hist = defaultdict(int)
    for _ in range(n):
        v = flat_values[rng.randrange(n)]
        hist[v] += 1
    return dict(hist)


def bootstrap_reclassification_ci(dist, old_R, new_R,
                                   n_boot=NUM_BOOTSTRAP,
                                   seed=RANDOM_SEED,
                                   ci=CONFIDENCE_LEVEL):
    """Bootstrap CI (percentile method) for reclassification effect."""
    rng = random.Random(seed)
    flat = expand_distribution(dist)
    effects = []
    for _ in range(n_boot):
        boot_dist = bootstrap_mic_distribution(flat, rng)
        effects.append(reclassification_effect(boot_dist, old_R, new_R))
    effects.sort()
    alpha = (1.0 - ci) / 2.0
    lo = percentile(effects, alpha)
    hi = percentile(effects, 1.0 - alpha)
    return {
        "mean": mean(effects),
        "sd": stdev(effects),
        "lo": lo,
        "hi": hi,
        "n_boot": n_boot,
    }


# ─── Helper: interrupted-time-series step estimation ───────────────────

def _solve_linear_system(A, b):
    """Solve A x = b by Gauss–Jordan elimination. A is an n×n list-of-lists,
    b is length n. Returns x as a length-n list.  Raises ValueError if
    singular."""
    n = len(A)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for i in range(n):
        # Partial pivot.
        max_row = max(range(i, n), key=lambda r: abs(M[r][i]))
        if abs(M[max_row][i]) < 1e-14:
            raise ValueError("singular matrix")
        M[i], M[max_row] = M[max_row], M[i]
        piv = M[i][i]
        M[i] = [v / piv for v in M[i]]
        for r in range(n):
            if r == i:
                continue
            factor = M[r][i]
            M[r] = [M[r][c] - factor * M[i][c] for c in range(n + 1)]
    return [M[i][n] for i in range(n)]


def itp_step_change(series, revision_year,
                    pre_window=PRE_WINDOW_YEARS,
                    post_window=POST_WINDOW_YEARS):
    """Segmented-regression interrupted-time-series estimator:

        y_t = β0 + β1·(t − t_rev) + β2·post_t + β3·(t − t_rev)·post_t + ε

    where post_t = 1 for t ≥ t_rev, 0 otherwise. β2 is the level
    change (the "step") at the revision year. Uses only bounded years
    [revision_year − pre_window, revision_year + post_window]. Returns
    the step β2 and its standard error from the sandwich (X'X)^{-1}·σ²
    residual-variance estimator.

    Falls back to a pre/post mean-difference estimator when the design
    matrix is rank-deficient (too few pre or post points to identify
    β1 or β3)."""
    yrs = sorted(series.keys())
    keep = [y for y in yrs
            if revision_year - pre_window <= y <= revision_year + post_window]
    if len(keep) < 4:
        return {"step": float("nan"), "slope_pre": float("nan"),
                "slope_change": float("nan"), "pre_mean": float("nan"),
                "post_mean": float("nan"), "n_pre": 0, "n_post": 0,
                "se_step": float("nan"), "method": "insufficient_data"}
    pre_years = [y for y in keep if y < revision_year]
    post_years = [y for y in keep if y >= revision_year]
    pre_vals = [series[y] for y in pre_years]
    post_vals = [series[y] for y in post_years]

    # Build segmented-regression design.
    t_centered = [y - revision_year for y in keep]
    post_flags = [1 if y >= revision_year else 0 for y in keep]
    y = [series[t] for t in keep]

    def _design(use_slope_change):
        if use_slope_change:
            X = [[1.0, t_centered[i], float(post_flags[i]),
                  t_centered[i] * post_flags[i]]
                 for i in range(len(keep))]
        else:
            X = [[1.0, t_centered[i], float(post_flags[i])]
                 for i in range(len(keep))]
        return X

    def _fit(X, y):
        # OLS via normal equations: (X'X) β = X'y.
        p = len(X[0])
        XtX = [[sum(X[k][i] * X[k][j] for k in range(len(X)))
                for j in range(p)] for i in range(p)]
        Xty = [sum(X[k][i] * y[k] for k in range(len(X))) for i in range(p)]
        beta = _solve_linear_system(XtX, Xty)
        resid = [y[k] - sum(beta[j] * X[k][j] for j in range(p))
                 for k in range(len(X))]
        dof = max(1, len(X) - p)
        sigma2 = sum(r * r for r in resid) / dof
        # Var(β) = σ² (X'X)^{-1}; we need only diag for SEs.
        inv_diag = []
        for i in range(p):
            e_i = [1.0 if k == i else 0.0 for k in range(p)]
            inv_col = _solve_linear_system(XtX, e_i)
            inv_diag.append(inv_col[i])
        se = [math.sqrt(sigma2 * d) if d >= 0 else float("nan")
              for d in inv_diag]
        return beta, se, sigma2

    # Try the full four-parameter segmented regression.
    try:
        if len(pre_vals) >= 2 and len(post_vals) >= 2:
            beta, se, _ = _fit(_design(True), y)
            step, slope_pre, slope_change = beta[2], beta[1], beta[3]
            se_step = se[2]
            method = "segmented_regression_full"
        elif len(pre_vals) >= 1 and len(post_vals) >= 1 and len(keep) >= 3:
            beta, se, _ = _fit(_design(False), y)
            step, slope_pre, slope_change = beta[2], beta[1], float("nan")
            se_step = se[2]
            method = "segmented_regression_no_slope_change"
        else:
            raise ValueError("rank_deficient")
    except (ValueError, ZeroDivisionError):
        if len(pre_vals) == 0 or len(post_vals) == 0:
            return {"step": float("nan"), "slope_pre": float("nan"),
                    "slope_change": float("nan"),
                    "pre_mean": mean(pre_vals) if pre_vals else float("nan"),
                    "post_mean": mean(post_vals) if post_vals else float("nan"),
                    "n_pre": len(pre_vals), "n_post": len(post_vals),
                    "se_step": float("nan"), "method": "failed"}
        step = mean(post_vals) - mean(pre_vals)
        var_pre = variance(pre_vals) if len(pre_vals) > 1 else 0.0
        var_post = variance(post_vals) if len(post_vals) > 1 else 0.0
        se_step = math.sqrt((var_pre / max(len(pre_vals), 1)) +
                            (var_post / max(len(post_vals), 1)))
        slope_pre = float("nan"); slope_change = float("nan")
        method = "mean_difference_fallback"

    return {
        "step": step, "slope_pre": slope_pre, "slope_change": slope_change,
        "pre_mean": mean(pre_vals), "post_mean": mean(post_vals),
        "n_pre": len(pre_vals), "n_post": len(post_vals),
        "se_step": se_step, "method": method,
    }


# ─── Helper: permutation test for step-size at revision year ──────────

def permutation_step_significance(series, revision_year,
                                  pre_window=PRE_WINDOW_YEARS,
                                  post_window=POST_WINDOW_YEARS,
                                  n_perm=NUM_PERMUTATIONS,
                                  seed=RANDOM_SEED):
    """Null: the revision year is uninformative — any year is a priori
    equally likely to host the 'revision'.  **Placebo-intervention
    permutation**: preserve the series values in their observed
    temporal order, but treat each admissible non-revision year as a
    candidate revision year.  Compute the ITS step at each such
    placebo year and compare |observed| to the empirical distribution
    of placebo |steps|.  This preserves all autocorrelation.

    Because the series is short (typically 9 years), we enumerate
    admissible placebo years rather than randomly sample; n_perm
    then caps at the number of admissible placebos.  p-value uses
    +1/+1 smoothing."""
    yrs = sorted(series.keys())
    if len(yrs) < 4:
        return {"p_value": float("nan"), "n_perm": 0,
                "observed_step": float("nan"),
                "null_mean_abs_step": float("nan"),
                "null": "placebo_intervention_years"}
    observed = itp_step_change(series, revision_year,
                               pre_window, post_window)["step"]
    if observed is None or math.isnan(observed):
        return {"p_value": float("nan"), "n_perm": 0,
                "observed_step": float("nan"),
                "null_mean_abs_step": float("nan"),
                "null": "placebo_intervention_years"}
    # Enumerate admissible placebo years: those with at least one
    # pre-year and one post-year within the pre/post windows, and not
    # the true revision year.
    candidates = []
    for y in yrs:
        has_pre = any((y - pre_window) <= yy < y for yy in yrs)
        has_post = any(y <= yy <= (y + post_window) for yy in yrs)
        if has_pre and has_post and y != revision_year:
            candidates.append(y)
    if not candidates:
        return {"p_value": float("nan"), "n_perm": 0,
                "observed_step": observed,
                "null_mean_abs_step": float("nan"),
                "null": "placebo_intervention_years"}
    # If we have fewer candidates than n_perm, enumerate each
    # exactly; otherwise sample with replacement.  (In practice,
    # for 9-year series we always enumerate.)
    rng = random.Random(seed)
    null_steps = []
    if len(candidates) <= n_perm:
        perm_years = candidates
    else:
        perm_years = [candidates[rng.randrange(len(candidates))]
                      for _ in range(n_perm)]
    for y in perm_years:
        step = itp_step_change(series, y, pre_window, post_window)["step"]
        if step is not None and not math.isnan(step):
            null_steps.append(step)
    if not null_steps:
        return {"p_value": float("nan"), "n_perm": 0,
                "observed_step": observed,
                "null_mean_abs_step": float("nan"),
                "null": "placebo_intervention_years"}
    obs_abs = abs(observed)
    tail = sum(1 for a in null_steps if abs(a) >= obs_abs)
    p = (tail + 1) / (len(null_steps) + 1)
    return {
        "p_value": p,
        "n_perm": len(null_steps),
        "observed_step": observed,
        "null_mean_abs_step": mean(abs(a) for a in null_steps),
        "null": "placebo_intervention_years",
    }


# ─── Helper: inverse-variance fixed-effects meta-analysis ──────────────

def meta_analyze(point_estimates, ses):
    """Inverse-variance fixed-effects pooling with Cochran's Q + I²."""
    paired = [(p, s) for p, s in zip(point_estimates, ses)
              if s is not None and s > 0 and not math.isnan(p)
              and not math.isnan(s)]
    if not paired:
        return {"pooled": float("nan"), "se": float("nan"),
                "ci_lo": float("nan"), "ci_hi": float("nan"), "k": 0,
                "Q": float("nan"), "Q_df": 0, "I2": float("nan")}
    weights = [1.0 / (s ** 2) for _, s in paired]
    w_sum = sum(weights)
    pooled = sum(w * p for w, (p, _) in zip(weights, paired)) / w_sum
    pooled_se = math.sqrt(1.0 / w_sum)
    # Cochran's Q for heterogeneity; I² = max(0, (Q - df) / Q).
    Q = sum(w * (p - pooled) ** 2 for w, (p, _) in zip(weights, paired))
    df = len(paired) - 1
    I2 = max(0.0, (Q - df) / Q) if Q > 0 else 0.0
    z = 1.959963984540054
    return {
        "pooled": pooled, "se": pooled_se,
        "ci_lo": pooled - z * pooled_se,
        "ci_hi": pooled + z * pooled_se,
        "k": len(paired),
        "Q": Q, "Q_df": df, "I2": I2,
    }


# ─── Load data ─────────────────────────────────────────────────────────

def load_data():
    """Compute canonical-data provenance hash, optionally fetch network
    reachability probe, and package the embedded canonical data into
    analysis structures."""
    provenance = ensure_provenance_cached()

    units = []
    for key, bkpt in BREAKPOINT_HISTORY.items():
        if key not in MIC_DISTRIBUTIONS:
            continue
        if key not in REPORTED_TRENDS:
            continue
        units.append({
            "organism": key[0],
            "antibiotic": key[1],
            "mic_distribution": MIC_DISTRIBUTIONS[key],
            "breakpoint": bkpt,
            "reported_trend": REPORTED_TRENDS[key],
        })
    return {
        "canonical_data_sha256": provenance["canonical_data_sha256"],
        "reachability_probe_sha256": provenance["reachability_probe_sha256"],
        "units": units,
        "num_bootstrap": NUM_BOOTSTRAP,
        "num_permutations": NUM_PERMUTATIONS,
        "pre_window_years": PRE_WINDOW_YEARS,
        "post_window_years": POST_WINDOW_YEARS,
        "confidence_level": CONFIDENCE_LEVEL,
        "random_seed": RANDOM_SEED,
    }


# ─── Analyse ───────────────────────────────────────────────────────────

def analyse_unit(unit, seed_offset=0):
    """Full analysis for a single organism–antibiotic unit."""
    dist = unit["mic_distribution"]
    bkpt = unit["breakpoint"]
    series = unit["reported_trend"]

    # % resistant under each breakpoint regime.
    pct_r_old = pct_resistant(dist, bkpt["R_old"])
    pct_r_new = pct_resistant(dist, bkpt["R_new"])
    reclass = pct_r_new - pct_r_old

    # Bootstrap CI on reclassification effect.
    boot = bootstrap_reclassification_ci(
        dist, bkpt["R_old"], bkpt["R_new"],
        n_boot=NUM_BOOTSTRAP,
        seed=RANDOM_SEED + seed_offset,
    )

    # Observed step in the reported time series at the revision year.
    step = itp_step_change(series, bkpt["revision_year"])

    # Attribution = reclassification / observed step.
    observed = step["step"]
    if observed and not math.isnan(observed) and observed != 0:
        attribution = reclass / observed
    else:
        attribution = float("nan")

    # Uncertainty in attribution (delta method using boot SE on reclass
    # and ITS SE on step).
    if (not math.isnan(observed) and observed != 0
            and not math.isnan(step["se_step"])
            and not math.isnan(boot["sd"])):
        # Var(R/O) ~= (1/O)^2 Var(R) + (R/O^2)^2 Var(O)
        var_attr = ((boot["sd"] / observed) ** 2 +
                    (reclass * step["se_step"] / (observed ** 2)) ** 2)
        se_attr = math.sqrt(var_attr)
        ci_attr = [attribution - 1.96 * se_attr,
                   attribution + 1.96 * se_attr]
    else:
        se_attr = float("nan")
        ci_attr = [float("nan"), float("nan")]

    # Permutation test: is the observed step at the revision year
    # larger than expected under a null shuffle of the same time series?
    perm = permutation_step_significance(series, bkpt["revision_year"])

    # Sensitivity: 2-year windows.
    step_2y = itp_step_change(series, bkpt["revision_year"],
                              pre_window=2, post_window=2)
    attr_2y = (reclass / step_2y["step"]
               if step_2y["step"] and not math.isnan(step_2y["step"])
               and step_2y["step"] != 0 else float("nan"))

    return {
        "organism": unit["organism"],
        "antibiotic": unit["antibiotic"],
        "revision_year": bkpt["revision_year"],
        "version_before": bkpt["table_version_before"],
        "version_after": bkpt["table_version_after"],
        "R_old_mg_per_L": bkpt["R_old"],
        "R_new_mg_per_L": bkpt["R_new"],
        "direction": bkpt["direction"],
        "source": bkpt["source"],
        "n_isolates": sum(dist.values()),
        "pct_resistant_under_old_R": pct_r_old,
        "pct_resistant_under_new_R": pct_r_new,
        "reclassification_effect_pp": reclass,
        "reclassification_bootstrap": boot,
        "observed_step_pp": observed,
        "observed_step_se": step["se_step"],
        "observed_pre_mean": step["pre_mean"],
        "observed_post_mean": step["post_mean"],
        "n_pre_years": step["n_pre"],
        "n_post_years": step["n_post"],
        "its_method": step["method"],
        "its_slope_pre": step["slope_pre"],
        "its_slope_change": step["slope_change"],
        "attribution_fraction": attribution,
        "attribution_se": se_attr,
        "attribution_ci95": ci_attr,
        "permutation_test": perm,
        "sensitivity_2yr_window": {
            "observed_step_pp": step_2y["step"],
            "attribution_fraction": attr_2y,
            "its_method": step_2y["method"],
        },
    }


def run_analysis(data):
    """Per-unit analysis plus meta-analysis across units."""
    per_unit = []
    for i, unit in enumerate(data["units"]):
        per_unit.append(analyse_unit(unit, seed_offset=i))

    # Meta-analyse attribution across units (inverse-variance weights).
    attrs = [u["attribution_fraction"] for u in per_unit]
    ses = [u["attribution_se"] for u in per_unit]
    meta_all = meta_analyze(attrs, ses)

    # Leave-one-out sensitivity for meta-analysis.
    loo = []
    for skip in range(len(per_unit)):
        a = [u["attribution_fraction"]
             for i, u in enumerate(per_unit) if i != skip]
        s = [u["attribution_se"]
             for i, u in enumerate(per_unit) if i != skip]
        m = meta_analyze(a, s)
        loo.append({
            "skipped": f"{per_unit[skip]['organism']} / "
                       f"{per_unit[skip]['antibiotic']}",
            "pooled": m["pooled"], "ci_lo": m["ci_lo"], "ci_hi": m["ci_hi"],
        })

    # Subset sensitivity: EUCAST-only (exclude CLSI S. pneumoniae).
    eucast_mask = [
        u for u in per_unit
        if u["version_before"].startswith("v")
    ]
    meta_eucast = meta_analyze(
        [u["attribution_fraction"] for u in eucast_mask],
        [u["attribution_se"] for u in eucast_mask],
    )

    # Subset sensitivity: 2019 EUCAST v9.0 cohort only.
    v9_mask = [u for u in per_unit if u["version_after"] == "v9.0"]
    meta_v9 = meta_analyze(
        [u["attribution_fraction"] for u in v9_mask],
        [u["attribution_se"] for u in v9_mask],
    )

    # Global permutation p-value: combine per-unit p-values via Fisher.
    ps = [u["permutation_test"]["p_value"] for u in per_unit
          if not math.isnan(u["permutation_test"]["p_value"])]
    if ps:
        fisher_stat = -2.0 * sum(math.log(max(p, 1e-12)) for p in ps)
        # Chi-square survival with 2k df via series expansion would be
        # complex; report the statistic and df, defer p to table.
        fisher_df = 2 * len(ps)
    else:
        fisher_stat, fisher_df = float("nan"), 0

    # Mean and median attribution across units (descriptive).
    mean_attr = mean(a for a in attrs if not math.isnan(a))
    sorted_attrs = sorted(a for a in attrs if not math.isnan(a))
    median_attr = percentile(sorted_attrs, 0.5) if sorted_attrs else float("nan")

    return {
        "per_unit": per_unit,
        "meta_all": meta_all,
        "meta_eucast_only": meta_eucast,
        "meta_2019_v9_only": meta_v9,
        "leave_one_out": loo,
        "fisher_combined": {
            "statistic": fisher_stat, "df": fisher_df, "k": len(ps),
        },
        "descriptive": {
            "mean_attribution": mean_attr,
            "median_attribution": median_attr,
            "n_units": len(per_unit),
        },
        "config": {
            "num_bootstrap": data["num_bootstrap"],
            "num_permutations": data["num_permutations"],
            "pre_window_years": data["pre_window_years"],
            "post_window_years": data["post_window_years"],
            "confidence_level": data["confidence_level"],
            "random_seed": data["random_seed"],
        },
        "canonical_data_sha256": data["canonical_data_sha256"],
        "reachability_probe_sha256": data["reachability_probe_sha256"],
    }


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

def generate_report(results):
    """Write structured JSON and human-readable Markdown."""
    with open(RESULTS_FILE, "w") as f:
        json.dump(results, f, indent=2, default=str)

    lines = []
    lines.append("# Breakpoint Reclassification vs. Biological Resistance "
                 "Evolution — Results")
    lines.append("")
    cfg = results["config"]
    lines.append(f"- Bootstrap iterations: {cfg['num_bootstrap']}")
    lines.append(f"- Permutation iterations: {cfg['num_permutations']}")
    lines.append(f"- Pre/post window (years): "
                 f"{cfg['pre_window_years']}/{cfg['post_window_years']}")
    lines.append(f"- Random seed: {cfg['random_seed']}")
    lines.append(f"- Units analysed: {results['descriptive']['n_units']}")
    lines.append(f"- Canonical data SHA-256: "
                 f"{results['canonical_data_sha256']}")
    rh = results['reachability_probe_sha256']
    lines.append(f"- Reachability probe SHA-256: "
                 f"{rh if rh else '(offline — probe skipped)'}")
    lines.append("")

    lines.append("## Per-unit results")
    lines.append("")
    lines.append("| Organism / Antibiotic | Year | Old R | New R | "
                 "Reclass (pp) | Observed step (pp) | Attribution (frac) |")
    lines.append("|---|---|---|---|---|---|---|")
    for u in results["per_unit"]:
        lines.append(
            f"| {u['organism']} / {u['antibiotic']} | "
            f"{u['revision_year']} | "
            f"{u['R_old_mg_per_L']} | {u['R_new_mg_per_L']} | "
            f"{u['reclassification_effect_pp']:+.2f} "
            f"(95% CI {u['reclassification_bootstrap']['lo']:+.2f}, "
            f"{u['reclassification_bootstrap']['hi']:+.2f}) | "
            f"{u['observed_step_pp']:+.2f} | "
            f"{u['attribution_fraction']:+.3f} |")
    lines.append("")

    lines.append("## Meta-analysis")
    lines.append("")
    m = results["meta_all"]
    lines.append(f"- **All units (k={m['k']}):** pooled attribution = "
                 f"{m['pooled']:.3f} (95% CI {m['ci_lo']:.3f}, "
                 f"{m['ci_hi']:.3f})")
    m = results["meta_eucast_only"]
    lines.append(f"- **EUCAST only (k={m['k']}):** "
                 f"{m['pooled']:.3f} (95% CI {m['ci_lo']:.3f}, "
                 f"{m['ci_hi']:.3f})")
    m = results["meta_2019_v9_only"]
    lines.append(f"- **2019 v9.0 cohort (k={m['k']}):** "
                 f"{m['pooled']:.3f} (95% CI {m['ci_lo']:.3f}, "
                 f"{m['ci_hi']:.3f})")
    lines.append("")

    lines.append("## Leave-one-out sensitivity")
    lines.append("")
    lines.append("| Unit removed | Pooled | 95% CI |")
    lines.append("|---|---|---|")
    for l in results["leave_one_out"]:
        lines.append(f"| {l['skipped']} | {l['pooled']:.3f} | "
                     f"[{l['ci_lo']:.3f}, {l['ci_hi']:.3f}] |")
    lines.append("")

    lines.append("## Permutation-test summary")
    lines.append("")
    lines.append("Null: the revision year is uninformative (label-shuffle "
                 "of the reported time series); test statistic = |observed "
                 "step at revision year|.")
    lines.append("")
    lines.append("| Unit | |Observed step| | p-value | n perms |")
    lines.append("|---|---|---|---|")
    for u in results["per_unit"]:
        p = u["permutation_test"]
        lines.append(f"| {u['organism']} / {u['antibiotic']} | "
                     f"{abs(p['observed_step']):.2f} | "
                     f"{p['p_value']:.4f} | {p['n_perm']} |")
    lines.append("")

    lines.append("## Heterogeneity (Cochran's Q, I²)")
    lines.append("")
    m = results["meta_all"]
    lines.append(f"- **All units:** Q = {m['Q']:.2f} (df = {m['Q_df']}), "
                 f"I² = {100.0 * m['I2']:.1f}%")
    m = results["meta_eucast_only"]
    lines.append(f"- **EUCAST only:** Q = {m['Q']:.2f} (df = {m['Q_df']}), "
                 f"I² = {100.0 * m['I2']:.1f}%")
    m = results["meta_2019_v9_only"]
    lines.append(f"- **2019 v9.0 cohort:** Q = {m['Q']:.2f} "
                 f"(df = {m['Q_df']}), I² = {100.0 * m['I2']:.1f}%")
    lines.append("")

    lines.append("## Outlier flag")
    lines.append("")
    loo = results["leave_one_out"]
    if loo:
        overall = results["meta_all"]["pooled"]
        deltas = [(abs(l["pooled"] - overall), l) for l in loo]
        deltas.sort(reverse=True)
        max_delta, worst = deltas[0]
        lines.append(f"- Removing **{worst['skipped']}** shifts the pooled "
                     f"estimate by {max_delta:.3f} (from {overall:.3f} to "
                     f"{worst['pooled']:.3f}). Treat any meta-analytic "
                     f"interpretation with awareness of this leverage.")
    lines.append("")

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


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

def verify():
    """Machine-checkable sanity assertions on the produced results."""
    if not os.path.exists(RESULTS_FILE):
        print("FAIL: results.json missing")
        sys.exit(1)
    with open(RESULTS_FILE) as f:
        r = json.load(f)

    checks = []

    # 1. Expected number of units analysed.
    n = r["descriptive"]["n_units"]
    checks.append(("num_units ≥ 6", n >= 6))

    # 2. Each unit has a non-trivial MIC sample size.
    small = [u for u in r["per_unit"] if u["n_isolates"] < 5000]
    checks.append(("every unit has ≥ 5000 isolates", len(small) == 0))

    # 3. Bootstrap iterations as configured.
    checks.append(("bootstrap iters = 2000",
                   r["config"]["num_bootstrap"] == 2000))

    # 4. Permutation sampling budget as configured (note: placebo-year
    # enumeration yields far fewer actual draws; see check #15 for the
    # actual floor).
    checks.append(("permutation budget = 2000 (config only)",
                   r["config"]["num_permutations"] == 2000))

    # 5. All per-unit bootstrap CIs have lo ≤ mean ≤ hi.
    bad_ci = [u for u in r["per_unit"]
              if not (u["reclassification_bootstrap"]["lo"]
                      <= u["reclassification_bootstrap"]["mean"]
                      <= u["reclassification_bootstrap"]["hi"])]
    checks.append(("all bootstrap CIs contain their point estimate",
                   len(bad_ci) == 0))

    # 6. Reclassification effects are finite real numbers.
    bad_re = [u for u in r["per_unit"]
              if (u["reclassification_effect_pp"] is None or
                  (isinstance(u["reclassification_effect_pp"], float) and
                   math.isnan(u["reclassification_effect_pp"])))]
    checks.append(("all reclassification effects are finite",
                   len(bad_re) == 0))

    # 7. Meta-analysis pooled attribution is in plausible range [-5, 5].
    pooled = r["meta_all"]["pooled"]
    checks.append(("meta_all pooled attribution in [-5, 5]",
                   -5.0 <= pooled <= 5.0))

    # 8. Meta-analysis CI contains the point estimate.
    mlo, mhi = r["meta_all"]["ci_lo"], r["meta_all"]["ci_hi"]
    checks.append(("meta_all CI contains the point estimate",
                   mlo <= pooled <= mhi))

    # 9. Canonical-data SHA-256 is valid hex and matches the re-computed hash.
    ph = r["canonical_data_sha256"]
    valid_ph = (isinstance(ph, str) and len(ph) == 64
                and all(c in "0123456789abcdef" for c in ph))
    # Recompute to confirm determinism.
    recomputed = canonical_data_fingerprint()
    checks.append(("canonical data SHA-256 is valid hex",
                   valid_ph and ph == recomputed))

    # 10. At least one unit has reclassification effect |x| > 1 pp.
    big = [u for u in r["per_unit"]
           if abs(u["reclassification_effect_pp"]) > 1.0]
    checks.append(("≥ 1 unit has |reclassification| > 1 pp", len(big) >= 1))

    # 11. Leave-one-out has as many rows as units.
    checks.append(("leave-one-out has n_units rows",
                   len(r["leave_one_out"]) == n))

    # 12. Seed is 42.
    checks.append(("random seed = 42", r["config"]["random_seed"] == 42))

    # 13. Every bootstrap CI has non-trivial width.
    trivial = [u for u in r["per_unit"]
               if u["reclassification_effect_pp"] != 0
               and (u["reclassification_bootstrap"]["hi"]
                    - u["reclassification_bootstrap"]["lo"]) <= 0]
    checks.append(("bootstrap CI width > 0 where reclass != 0",
                   len(trivial) == 0))

    # 14. Meta-analysis exposes heterogeneity statistics.
    has_Q = ("Q" in r["meta_all"] and "I2" in r["meta_all"]
             and r["meta_all"]["Q_df"] == r["meta_all"]["k"] - 1)
    checks.append(("meta_all reports Cochran's Q and I²", has_Q))

    # 15. Permutation test has at least 3 admissible placebo years per unit
    # (enumerated; short time-series permit fewer than NUM_PERMUTATIONS).
    bad_perm = [u for u in r["per_unit"]
                if u["permutation_test"]["n_perm"] < 3]
    checks.append(("permutation test ≥ 3 placebo years per unit",
                   len(bad_perm) == 0))

    # 16. All ITS fits are non-NaN and use segmented regression or a
    # documented fallback.
    valid_methods = {"segmented_regression_full",
                     "segmented_regression_no_slope_change",
                     "mean_difference_fallback"}
    bad_its = [u for u in r["per_unit"]
               if (u.get("its_method") not in valid_methods)
               or math.isnan(u["observed_step_pp"])]
    checks.append(("every unit has a non-NaN ITS step with a "
                   "documented method", len(bad_its) == 0))

    # 17. Every unit's permutation test has a non-NaN p-value.
    nan_p = [u for u in r["per_unit"]
             if math.isnan(u["permutation_test"]["p_value"])]
    checks.append(("every permutation p-value is non-NaN", len(nan_p) == 0))

    # 18. Null type is placebo_intervention_years (preserves order and AC).
    wrong_null = [u for u in r["per_unit"]
                  if u["permutation_test"].get("null")
                  != "placebo_intervention_years"]
    checks.append(("permutation null is placebo_intervention_years",
                   len(wrong_null) == 0))

    passed = sum(1 for _, ok in checks if ok)
    for name, ok in checks:
        print(f"  {'PASS' if ok else 'FAIL'}  {name}")
    if passed == len(checks):
        print(f"\nVERIFICATION PASSED ({passed}/{len(checks)})")
        sys.exit(0)
    else:
        print(f"\nVERIFICATION FAILED ({passed}/{len(checks)})")
        sys.exit(1)


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

def main():
    if len(sys.argv) > 1 and sys.argv[1] == "--verify":
        verify()
        return

    print("[1/5] Fetching provenance fingerprint (download + cache)...")
    random.seed(RANDOM_SEED)
    data = load_data()

    print(f"[2/5] Loaded {len(data['units'])} organism–antibiotic units "
          f"from embedded canonical data.")
    for u in data["units"]:
        print(f"    - {u['organism']} / {u['antibiotic']} "
              f"(revision {u['breakpoint']['revision_year']}, "
              f"{sum(u['mic_distribution'].values())} isolates)")

    print(f"[3/5] Running {NUM_BOOTSTRAP}-iter bootstrap + "
          f"placebo-year permutation test per unit...")
    results = run_analysis(data)

    print("[4/5] Writing results.json and report.md...")
    generate_report(results)
    print(f"    - {RESULTS_FILE}")
    print(f"    - {REPORT_FILE}")

    print("[5/5] Summary:")
    print(f"    - Units analysed:             "
          f"{results['descriptive']['n_units']}")
    m = results["meta_all"]
    print(f"    - Pooled attribution (all):   "
          f"{m['pooled']:.3f} (95% CI {m['ci_lo']:.3f}, {m['ci_hi']:.3f})")
    m = results["meta_eucast_only"]
    print(f"    - Pooled attribution (EUCAST): "
          f"{m['pooled']:.3f} (95% CI {m['ci_lo']:.3f}, {m['ci_hi']:.3f})")
    m = results["meta_2019_v9_only"]
    print(f"    - Pooled attribution (v9.0):  "
          f"{m['pooled']:.3f} (95% CI {m['ci_lo']:.3f}, {m['ci_hi']:.3f})")

    print("\nANALYSIS COMPLETE")


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

**Expected output:** No stdout. The file `/tmp/claw4s_auto_breakpoint-changes-manufacturing-resistance-trends/analyze.py` now exists.

**Success condition:** `analyze.py` is a valid Python file (syntactically loadable by `python3 -c "import ast; ast.parse(open('/tmp/claw4s_auto_breakpoint-changes-manufacturing-resistance-trends/analyze.py').read())"` — absolute path; does not depend on current directory).

---

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_breakpoint-changes-manufacturing-resistance-trends && python3 analyze.py
```

**Expected output:**

```
[1/5] Fetching provenance fingerprint (download + cache)...
  Canonical-data SHA-256: dcc18400195bacc9ce224beb4fc9f02e7c0bd8b0644bda5a1b4ba16f820dc892
  Reachability probe downloaded (... bytes) SHA-256: <64-hex string>
[2/5] Loaded 6 organism–antibiotic units from embedded canonical data.
    - Escherichia coli / Ciprofloxacin (revision 2019, 86534 isolates)
    - Klebsiella pneumoniae / Ciprofloxacin (revision 2019, 24124 isolates)
    - Salmonella enterica / Ciprofloxacin (revision 2012, 18502 isolates)
    - Escherichia coli / Tigecycline (revision 2019, 11997 isolates)
    - Pseudomonas aeruginosa / Piperacillin-tazobactam (revision 2019, 14925 isolates)
    - Streptococcus pneumoniae / Penicillin (revision 2008, 32013 isolates)
[3/5] Running 2000-iter bootstrap + placebo-year permutation test per unit...
[4/5] Writing results.json and report.md...
    - /tmp/.../results.json
    - /tmp/.../report.md
[5/5] Summary:
    - Units analysed:             6
    - Pooled attribution (all):   0.113 (95% CI 0.104, 0.122)
    - Pooled attribution (EUCAST): 1.047 (95% CI 1.001, 1.093)
    - Pooled attribution (v9.0):  0.698 (95% CI 0.555, 0.842)

ANALYSIS COMPLETE
```

Note on permutation semantics: the permutation null is a **placebo-intervention-year enumeration**: each candidate year inside the observed series that admits a full pre/post window is treated as a hypothetical revision year, the ITS step is re-estimated there, and the observed step is ranked against that null. The number of admissible placebo years is typically 3–7 (≥ 3 enforced). The `NUM_PERMUTATIONS` config value is an upper sampling budget, **not** the number of placebo years actually evaluated; `results.json[unit].permutation_test.n_perm` records the actual count.

**Success condition:** stdout contains `ANALYSIS COMPLETE`. Files `results.json` and `report.md` are written in the workspace directory. `cache/canonical_data.sha256` is written and contains the pinned hash above; `cache/reachability_probe.bin` is written when the reachability probe succeeds (non-fatal otherwise).

**Failure modes:**
- **Canonical-hash mismatch:** the script aborts if the SHA-256 of the embedded `MIC_DISTRIBUTIONS + BREAKPOINT_HISTORY + REPORTED_TRENDS` does not equal `EXPECTED_DATA_SHA256`. This means the embedded data were modified; restore the original values or update the pinned hash.
- `URLError`/`HTTPError` on the reachability probe: both primary (`https://www.eucast.org/robots.txt`) and fallback (`https://www.ecdc.europa.eu/robots.txt`) failed. The script continues offline and records an `unavailable` tag in `results.json`. No user action required for the scientific output.
- `FileNotFoundError` on write: the workspace path was not created in Step 1.
- `KeyError` inside `analyse_unit`: a `BREAKPOINT_HISTORY` key was not present in both `MIC_DISTRIBUTIONS` and `REPORTED_TRENDS`; the `load_data` filter should prevent this.

---

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_breakpoint-changes-manufacturing-resistance-trends && python3 analyze.py --verify
```

**Expected output (tail):**

```
  PASS  num_units ≥ 6
  PASS  every unit has ≥ 5000 isolates
  PASS  bootstrap iters = 2000
  PASS  permutation budget = 2000 (config only)
  PASS  all bootstrap CIs contain their point estimate
  PASS  all reclassification effects are finite
  PASS  meta_all pooled attribution in [-5, 5]
  PASS  meta_all CI contains the point estimate
  PASS  provenance SHA256 is valid hex
  PASS  ≥ 1 unit has |reclassification| > 1 pp
  PASS  leave-one-out has n_units rows
  PASS  random seed = 42
  PASS  bootstrap CI width > 0 where reclass != 0
  PASS  meta_all reports Cochran's Q and I²
  PASS  permutation test ≥ 3 placebo years per unit
  PASS  every unit has a non-NaN ITS step with a documented method
  PASS  every permutation p-value is non-NaN
  PASS  permutation null is placebo_intervention_years

VERIFICATION PASSED (18/18)
```

**Success condition:** Exit code 0 and the string `VERIFICATION PASSED` in stdout.

**Failure condition:** Any `FAIL` line in stdout, or exit code 1.

---

## Success Criteria (overall skill)

**Structural (must hold for the run to be considered complete):**

1. Steps 1–4 complete without intervention and without any Python exception.
2. Step 3 stdout ends with the literal string `ANALYSIS COMPLETE`.
3. Step 4 stdout ends with the literal string `VERIFICATION PASSED (18/18)` and exits with code 0.
4. `results.json` is well-formed JSON that parses via `json.load()` without error.
5. `report.md` contains a per-unit table, a meta-analysis summary, a leave-one-out sensitivity table, a Cochran's Q / I² heterogeneity block, and an explicit outlier flag.
6. Cached files in `cache/` survive a second run (`cache/canonical_data.sha256` identical across runs; reachability probe re-used from cache if network is unavailable).

**Quantitative (bounded ranges for plausibility):**

7. Exactly 6 organism–antibiotic units are analysed.
8. Every unit has ≥ 5 000 isolates in its embedded MIC distribution.
9. Every unit's bootstrap 95 % CI half-width on the reclassification effect is < 5 pp (i.e., CIs are informative, not vacuous).
10. Every unit's permutation test enumerates ≥ 3 admissible placebo years.
11. Meta-analytic pooled attribution (all units) lies in the plausible interval [-5, 5].
12. Meta-analytic pooled CI contains the point estimate.
13. Canonical-data SHA-256 matches the pinned `EXPECTED_DATA_SHA256` (otherwise the script aborts before analysis).

**Scientific (substantive interpretation tests):**

14. At least one unit has |reclassification effect| > 1 pp (demonstrates the reclassification mechanism is material, not negligible).
15. Attribution fractions are heterogeneous across units (Cochran's Q >> df in the full-cohort meta-analysis) — the null of "uniform attribution" is rejected, confirming the science is a decomposition, not a single parameter.

All 18 assertions in `verify()` check subsets of criteria 1–15 mechanically.

## Failure Conditions

A run is considered **failed** if any of the following occur:

1. **Exception during analysis** — any unhandled Python exception in Steps 2 or 3. Inspect stderr for the traceback and the specific line number.
2. **Verification fails** — `python3 analyze.py --verify` prints fewer than 18 `PASS` lines or exits with code 1. Inspect which named check failed and consult the matching section of the code (e.g., `bootstrap CI width > 0` points to `bootstrap_reclassification_ci`).
3. **Missing outputs** — after a successful run, `results.json` or `report.md` is not present in the workspace directory. This indicates a write-permission or path error.
4. **Canonical-hash mismatch** — the script aborts with `RuntimeError: Embedded-data SHA-256 mismatch` when the embedded `MIC_DISTRIBUTIONS + BREAKPOINT_HISTORY + REPORTED_TRENDS` do not hash to `EXPECTED_DATA_SHA256`. This means the scientific inputs were modified since the pin was set. Either revert the modification or deliberately update the pin and document why.
5. **Degenerate ITS fit** — any unit reports `its_method = "failed"` or `"insufficient_data"`, or `observed_step_pp` is `NaN`. This indicates a reported-trend series too short for a pre/post step estimate.
6. **Singleton permutation support** — any unit has `permutation_test.n_perm < 3`. This fails verification check #15 and indicates the series does not admit enough placebo years to form even a minimal null distribution.
7. **Implausible pool** — meta-analytic pooled attribution outside [-5, 5]. This indicates that at least one unit has a very small observed step dividing a moderate reclassification effect (attribution fraction explodes); consult the leave-one-out table and exclude that unit if appropriate.

## Limitations and Assumptions

The interpretation of the output is bounded by the following:

1. **Ecological, not patient-level.** Results are EU/EEA-aggregated ECDC EARS-Net rates. The analysis does NOT make any patient-level causal claim about treatment outcomes.
2. **Static MIC distribution assumption.** The reclassification-only component holds biology *exactly* constant (one MIC histogram applied at two breakpoints). If the true MIC distribution is drifting over the pre/post window, the reclassification term is biased upward or downward depending on the direction of drift. This is a known and documented simplification.
3. **Sharp-break ITS.** `itp_step_change` assumes the policy revision takes effect at a single calendar year. Revisions that phase in over several years (as some EUCAST changes do) are modelled approximately; the level change at the nominal revision year absorbs transition dynamics.
4. **Placebo-year permutation has low resolution.** Reported-trend series are 9–11 years; the admissible placebo-year pool is typically 3–7 candidates, so the minimum achievable p-value is ≈ 1/7 ≈ 0.14. The permutation test is a sanity check; the bootstrap CI and segmented-regression SE are the headline uncertainty measures.
5. **Inverse-variance pooling assumes SEs are comparable.** The delta-method SE on the attribution fraction can be inflated when the observed step is small, leaking weight across units. Subset and leave-one-out checks are included to flag this.
6. **The embedded data are curated.** MIC distributions, breakpoint history, and reported rates are manually transcribed from public EUCAST / CLSI / ECDC sources and fingerprinted. A reviewer should verify a sample of the counts against the cited source before citing any numeric result beyond the attribution fraction itself.
7. **No multiple-testing correction across units.** Per-unit p-values are reported as-is; the Fisher combined statistic is reported without a final p because chi-square survival with 2k df is left to the consumer.
8. **What the results do NOT show:** (a) whether any individual patient's treatment changed; (b) whether resistance to *un-revised* antibiotics rose or fell; (c) laboratory-level variation in MIC determination; (d) the country-level heterogeneity inside the EU/EEA aggregate.

## Failure Modes Cheat-Sheet (for quick debugging)

| Symptom | Likely cause | Fix |
|---|---|---|
| `SHA-256 mismatch` abort | Embedded data edited | Revert or re-pin `EXPECTED_DATA_SHA256` |
| `VERIFICATION FAILED` with `n_perm < 3` | Trend series too short | Lengthen series or reduce `PRE_/POST_WINDOW_YEARS` |
| `its_method = "failed"` | Fewer than 2 pre- or post-points | Check `REPORTED_TRENDS` keys around the revision year |
| Bootstrap CI width = 0 | Distribution has a single MIC bin | Inspect histogram; add finer bins if available |
| Network probe error | Firewall / offline | Non-fatal; analysis still completes |
| `ANALYSIS COMPLETE` missing from stdout | Early exception | Check stderr traceback |
- Bootstrap iteration count below 2000 (indicates truncated run). The placebo-year permutation enumerates admissible revision years in the observed series; fewer than 3 placebo years for any unit is a failure (check #15).