code/reproduce_theory.py

#!/usr/bin/env python3
"""Exact finite-sample numerical audits for Claims 1--5.

This script reproduces the two-component Gaussian-mixture model from
"The Catastrophic Failure of the k-Means Algorithm in High Dimensions,
and How Hartigan's Algorithm Avoids It" without depending on the authors'
implementation.

The key reduction is exact.  Conditional on fixed ground-truth labels, each
feature column of X is Gaussian with covariance

    Sigma = tau**2 * Z Z.T + sigma**2 * I.

Every Lloyd/Hartigan fixed-point comparison depends on X only through the
Gram matrix G = X X.T, and therefore G is Wishart(d, Sigma).  Bartlett draws
make exhaustive partition checks inexpensive even when d is very large.

Only the Python standard library and NumPy are required.
"""

from __future__ import annotations

import argparse
import csv
import html
import json
import math
import os
import platform
import sys
import time
from dataclasses import asdict, dataclass
from pathlib import Path
from typing import Any, Iterable, Sequence

import numpy as np


PAPER_URL = "https://arxiv.org/abs/2602.09936"
CODE_URL = (
    "https://github.com/Lederman-Group/Catastrophic_Failure_KMeans/"
    "tree/e460f6198c42685e75efa7877bb66a55cab4d460"
)
ZENODO_URL = "https://zenodo.org/records/20435365"
FLOAT_EPS = np.finfo(np.float64).eps
TIE_ULPS = 128.0


@dataclass(frozen=True)
class Profile:
    name: str
    high_dimensions: tuple[int, ...]
    low_dimensions: tuple[int, ...]
    high_trials: int
    low_trials: int
    partition_batch_size: int
    claim5_dimensions: tuple[int, ...]
    claim5_ratios: tuple[float, ...]
    claim5_trials: int
    wishart_batch_size: int


PROFILES = {
    "smoke": Profile(
        name="smoke",
        high_dimensions=(64, 1024, 4096, 4245, 4631, 5913, 8192, 32768, 255164, 349345),
        low_dimensions=(64, 1024, 8192, 349345),
        high_trials=240,
        low_trials=160,
        partition_batch_size=80,
        claim5_dimensions=(8, 128, 1024, 8192, 100000),
        claim5_ratios=(0.5, 0.9, 1.0, 1.1, 1.5, 2.0),
        claim5_trials=1200,
        wishart_batch_size=1200,
    ),
    "full": Profile(
        name="full",
        high_dimensions=(
            64,
            256,
            1024,
            4096,
            4245,
            4631,
            5913,
            8192,
            32768,
            131072,
            226812,
            255164,
            349345,
            524288,
        ),
        low_dimensions=(64, 1024, 8192, 131072, 524288),
        high_trials=3000,
        low_trials=1200,
        partition_batch_size=150,
        claim5_dimensions=(8, 32, 128, 512, 1024, 2048, 8192, 32768, 100000),
        claim5_ratios=(0.1, 0.5, 0.9, 1.0, 1.1, 1.5, 2.0),
        claim5_trials=10000,
        wishart_batch_size=2500,
    ),
}


def _native(value: Any) -> Any:
    """Convert NumPy values recursively into strict JSON-compatible values."""

    if isinstance(value, dict):
        return {str(key): _native(val) for key, val in value.items()}
    if isinstance(value, (list, tuple)):
        return [_native(val) for val in value]
    if isinstance(value, np.ndarray):
        return _native(value.tolist())
    if isinstance(value, (np.integer,)):
        return int(value)
    if isinstance(value, (np.floating, float)):
        val = float(value)
        return val if math.isfinite(val) else None
    if isinstance(value, (np.bool_, bool)):
        return bool(value)
    return value


def wilson_interval(successes: int, trials: int, alpha: float = 0.05) -> tuple[float, float]:
    """Two-sided Wilson score interval using the standard-normal quantile.

    Python's standard library has no inverse normal CDF, but NormalDist does.
    Importing locally keeps compatibility with older minimal Python installs.
    """

    if trials <= 0:
        raise ValueError("trials must be positive")
    if not 0 <= successes <= trials:
        raise ValueError("successes must lie in [0, trials]")
    if not 0.0 < alpha < 1.0:
        raise ValueError("alpha must lie in (0, 1)")
    from statistics import NormalDist

    z = NormalDist().inv_cdf(1.0 - alpha / 2.0)
    phat = successes / trials
    z2 = z * z
    denominator = 1.0 + z2 / trials
    center = (phat + z2 / (2.0 * trials)) / denominator
    half_width = (
        z
        * math.sqrt(phat * (1.0 - phat) / trials + z2 / (4.0 * trials * trials))
        / denominator
    )
    return max(0.0, center - half_width), min(1.0, center + half_width)


def enumerate_bipartitions(n: int) -> np.ndarray:
    """Enumerate each nonempty unlabeled bipartition exactly once.

    Point zero is fixed in cluster 0, removing complement symmetry.  Rows are
    binary membership indicators for cluster 1.
    """

    if n < 2:
        raise ValueError("n must be at least 2")
    values = np.arange(1, 1 << (n - 1), dtype=np.uint64)
    bit_positions = np.arange(n - 1, dtype=np.uint64)
    remaining = ((values[:, None] >> bit_positions[None, :]) & 1).astype(np.float64)
    return np.concatenate((np.zeros((remaining.shape[0], 1)), remaining), axis=1)


def q_balanced_mask(partitions: np.ndarray, q: float) -> np.ndarray:
    """Definition 2.5, including strict bounds and the |C_k| > 2 clause."""

    if q <= 1.0:
        raise ValueError("Corollary 3.8 requires q > 1")
    n = partitions.shape[1]
    sizes = partitions.sum(axis=1)
    lower = n / 2.0 - q * math.sqrt(n / 4.0)
    upper = n / 2.0 + q * math.sqrt(n / 4.0)
    other_sizes = n - sizes
    return (
        (sizes > 2.0)
        & (other_sizes > 2.0)
        & (sizes > lower)
        & (sizes < upper)
        & (other_sizes > lower)
        & (other_sizes < upper)
    )


def model_covariance(true_labels: np.ndarray, tau2: float, sigma2: float) -> np.ndarray:
    """Covariance across observations for one independent feature column."""

    labels = np.asarray(true_labels)
    if labels.ndim != 1:
        raise ValueError("true_labels must be one-dimensional")
    if tau2 <= 0.0 or sigma2 <= 0.0:
        raise ValueError("tau2 and sigma2 must be positive")
    shared_center = labels[:, None] == labels[None, :]
    return tau2 * shared_center.astype(np.float64) + sigma2 * np.eye(labels.size)


def sample_wishart_normalized(
    rng: np.random.Generator,
    covariance: np.ndarray,
    degrees_of_freedom: int,
    draws: int,
) -> np.ndarray:
    """Draw Wishart matrices with Bartlett's decomposition and divide by d."""

    covariance = np.asarray(covariance, dtype=np.float64)
    p = covariance.shape[0]
    if covariance.shape != (p, p):
        raise ValueError("covariance must be square")
    if degrees_of_freedom < p:
        raise ValueError("Bartlett decomposition here requires d >= matrix dimension")
    if draws <= 0:
        raise ValueError("draws must be positive")

    chol = np.linalg.cholesky(covariance)
    bartlett = np.zeros((draws, p, p), dtype=np.float64)
    row, col = np.tril_indices(p, k=-1)
    bartlett[:, row, col] = rng.standard_normal((draws, row.size))
    for idx in range(p):
        bartlett[:, idx, idx] = np.sqrt(
            rng.chisquare(degrees_of_freedom - idx, size=draws)
        )
    transformed = np.einsum("ij,tjk->tik", chol, bartlett, optimize=True)
    gram = np.einsum("tik,tjk->tij", transformed, transformed, optimize=True)
    gram /= float(degrees_of_freedom)
    return gram


def gram_partition_distances(
    grams: np.ndarray, partitions: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
    """Squared distances to both empirical centroids from Gram matrices.

    Returned arrays have shape (draw, partition, point).  Since input Gram
    matrices may be divided by d, distances carry the same harmless scaling.
    """

    grams = np.asarray(grams, dtype=np.float64)
    partitions = np.asarray(partitions, dtype=np.float64)
    if grams.ndim != 3 or grams.shape[1] != grams.shape[2]:
        raise ValueError("grams must have shape (draw, n, n)")
    if partitions.ndim != 2 or partitions.shape[1] != grams.shape[1]:
        raise ValueError("partition width must match Gram dimension")

    n = grams.shape[1]
    sizes1 = partitions.sum(axis=1)
    sizes0 = n - sizes1
    if np.any(sizes1 <= 0.0) or np.any(sizes0 <= 0.0):
        raise ValueError("all partitions must contain two nonempty clusters")
    complements = 1.0 - partitions

    gb1 = np.einsum("tij,mj->tmi", grams, partitions, optimize=True)
    gb0 = np.einsum("tij,mj->tmi", grams, complements, optimize=True)
    quadratic1 = np.einsum("mi,tmi->tm", partitions, gb1, optimize=True)
    quadratic0 = np.einsum("mi,tmi->tm", complements, gb0, optimize=True)
    diagonal = np.diagonal(grams, axis1=1, axis2=2)[:, None, :]

    distance1 = (
        diagonal
        - 2.0 * gb1 / sizes1[None, :, None]
        + quadratic1[:, :, None] / (sizes1[None, :, None] ** 2)
    )
    distance0 = (
        diagonal
        - 2.0 * gb0 / sizes0[None, :, None]
        + quadratic0[:, :, None] / (sizes0[None, :, None] ** 2)
    )
    return distance0, distance1


def _near_tie_mask(margin: np.ndarray, current: np.ndarray, alternative: np.ndarray) -> np.ndarray:
    scale = np.maximum(np.abs(current) + np.abs(alternative), np.finfo(np.float64).tiny)
    return np.abs(margin) <= TIE_ULPS * FLOAT_EPS * scale


def fixed_point_tests(
    grams: np.ndarray, partitions: np.ndarray
) -> dict[str, np.ndarray | int | float]:
    """Apply exact Lloyd and Hartigan fixed-point definitions to every row."""

    distance0, distance1 = gram_partition_distances(grams, partitions)
    labels = partitions[None, :, :].astype(bool)
    current = np.where(labels, distance1, distance0)
    alternative = np.where(labels, distance0, distance1)

    # A move is accepted only for a strict improvement.  Thus a fixed point
    # uses current <= alternative (ties occur with probability zero in-model).
    lloyd_margin = alternative - current
    lloyd_fixed = np.all(lloyd_margin >= 0.0, axis=2)
    lloyd_near_ties = _near_tie_mask(lloyd_margin, current, alternative)

    n = partitions.shape[1]
    sizes1 = partitions.sum(axis=1)
    sizes0 = n - sizes1
    current_sizes = np.where(partitions.astype(bool), sizes1[:, None], sizes0[:, None])
    alternative_sizes = n - current_sizes
    movable = current_sizes > 1.0
    current_factor = np.zeros_like(current_sizes, dtype=np.float64)
    current_factor[movable] = current_sizes[movable] / (current_sizes[movable] - 1.0)
    alternative_factor = alternative_sizes / (alternative_sizes + 1.0)
    hartigan_current = current * current_factor[None, :, :]
    hartigan_alternative = alternative * alternative_factor[None, :, :]
    hartigan_margin = hartigan_alternative - hartigan_current
    hartigan_fixed = np.all((~movable[None, :, :]) | (hartigan_margin >= 0.0), axis=2)
    hartigan_near_ties = movable[None, :, :] & _near_tie_mask(
        hartigan_margin, hartigan_current, hartigan_alternative
    )

    finite_lloyd_margin = np.abs(lloyd_margin).ravel()
    finite_hart_margin = np.abs(hartigan_margin[:, movable]).ravel()
    return {
        "lloyd_fixed": lloyd_fixed,
        "hartigan_fixed": hartigan_fixed,
        "lloyd_near_ties": int(lloyd_near_ties.sum()),
        "hartigan_near_ties": int(hartigan_near_ties.sum()),
        "lloyd_comparisons": int(lloyd_margin.size),
        "hartigan_comparisons": int(grams.shape[0] * movable.sum()),
        "minimum_abs_lloyd_margin": float(finite_lloyd_margin.min()),
        "minimum_abs_hartigan_margin": float(finite_hart_margin.min()),
        "negative_distance_count": int(((distance0 < 0.0) | (distance1 < 0.0)).sum()),
    }


def corollary_38_parameters(n: int, q: float, beta: float) -> tuple[float, float]:
    """Return sigma^2 and rho_q from Theorem 3.6 / Corollary 3.8."""

    if n < 1 or q <= 1.0 or beta <= 1.0:
        raise ValueError("requires n >= 1, q > 1, and beta > 1")
    root_n = math.sqrt(n)
    sigma = beta * (root_n * q + n - 2.0) / (
        math.sqrt(2.0) * math.sqrt(root_n * q + n)
    )
    sigma2 = sigma * sigma
    numerator = (
        sigma2
        * (root_n * q + n - 2.0)
        * (root_n * q + n)
        * (root_n * q + n + 2.0)
        * (root_n * (sigma2 + 2.0) * (root_n + q) - 4.0)
    )
    denominator = (
        n * sigma2 * (root_n + q) ** 2 + (root_n * q + n - 2.0) ** 2
    ) ** 2
    rho_q = numerator / denominator
    if not 0.0 < rho_q < 1.0:
        raise ArithmeticError(f"rho_q must lie in (0, 1), got {rho_q}")
    return sigma2, rho_q


def rho_h_uniform(n: int, tau2: float, sigma2: float, r_star: float) -> float:
    term = 4.0 * tau2 * r_star * r_star / n / (3.0 * tau2 + sigma2)
    rho = 1.0 - term * term
    if not 0.0 < rho < 1.0:
        raise ArithmeticError(f"rho_h must lie in (0, 1), got {rho}")
    return rho


def probability_bound(log_prefactor: float, rho: float, dimension: int) -> tuple[float, float]:
    """Return capped bound and its uncapped natural logarithm."""

    log_bound = log_prefactor + (dimension / 4.0) * math.log(rho)
    if log_bound >= 0.0:
        return 1.0, log_bound
    return math.exp(log_bound), log_bound


def _correct_partition_mask(partitions: np.ndarray, true_labels: np.ndarray) -> np.ndarray:
    labels = np.asarray(true_labels, dtype=np.float64)
    return np.all(partitions == labels[None, :], axis=1) | np.all(
        partitions == (1.0 - labels)[None, :], axis=1
    )


def run_claims_1_to_4(
    rng: np.random.Generator,
    profile: Profile,
    *,
    n: int = 8,
    tau2: float = 1.0,
    q: float = 1.5,
    beta: float = 1.5,
) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    """Exhaustively audit fixed points and both corollary events."""

    if n % 2:
        raise ValueError("the default audit expects balanced, even n")
    true_labels = np.repeat(np.arange(2), n // 2)
    partitions = enumerate_bipartitions(n)
    balanced = q_balanced_mask(partitions, q)
    correct = _correct_partition_mask(partitions, true_labels)
    incorrect = ~correct
    high_sigma2, rho_q = corollary_38_parameters(n, q, beta)
    low_sigma2 = 1.0
    r_star = 0.5

    settings = [
        ("high_noise", high_sigma2, profile.high_dimensions, profile.high_trials, True),
        ("low_noise_control", low_sigma2, profile.low_dimensions, profile.low_trials, False),
    ]
    rows: list[dict[str, Any]] = []
    aggregate_diagnostics = {
        "lloyd_near_ties": 0,
        "hartigan_near_ties": 0,
        "lloyd_comparisons": 0,
        "hartigan_comparisons": 0,
        "negative_distance_count": 0,
        "hartigan_not_lloyd_violations": 0,
        "minimum_abs_lloyd_margin": math.inf,
        "minimum_abs_hartigan_margin": math.inf,
    }

    for regime, sigma2, dimensions, trial_count, cor38_applies in settings:
        covariance = model_covariance(true_labels, tau2, sigma2)
        rho_h = rho_h_uniform(n, tau2, sigma2, r_star)
        for dimension in dimensions:
            all_balanced_successes = 0
            any_incorrect_hartigan_successes = 0
            correct_hartigan_successes = 0
            balanced_fixed_fraction_sum = 0.0
            incorrect_hartigan_fraction_sum = 0.0
            all_lloyd_fixed_fraction_sum = 0.0
            trials_done = 0
            while trials_done < trial_count:
                batch = min(profile.partition_batch_size, trial_count - trials_done)
                grams = sample_wishart_normalized(rng, covariance, dimension, batch)
                tested = fixed_point_tests(grams, partitions)
                lloyd_fixed = np.asarray(tested["lloyd_fixed"], dtype=bool)
                hartigan_fixed = np.asarray(tested["hartigan_fixed"], dtype=bool)

                balanced_fixed = lloyd_fixed[:, balanced]
                incorrect_hartigan = hartigan_fixed[:, incorrect]
                all_balanced_successes += int(np.all(balanced_fixed, axis=1).sum())
                any_incorrect_hartigan_successes += int(
                    np.any(incorrect_hartigan, axis=1).sum()
                )
                correct_hartigan_successes += int(np.all(hartigan_fixed[:, correct], axis=1).sum())
                balanced_fixed_fraction_sum += float(balanced_fixed.mean(axis=1).sum())
                incorrect_hartigan_fraction_sum += float(
                    incorrect_hartigan.mean(axis=1).sum()
                )
                all_lloyd_fixed_fraction_sum += float(lloyd_fixed.mean(axis=1).sum())
                subset_violations = hartigan_fixed & ~lloyd_fixed
                aggregate_diagnostics["hartigan_not_lloyd_violations"] += int(
                    subset_violations.sum()
                )
                for key in (
                    "lloyd_near_ties",
                    "hartigan_near_ties",
                    "lloyd_comparisons",
                    "hartigan_comparisons",
                    "negative_distance_count",
                ):
                    aggregate_diagnostics[key] += int(tested[key])
                aggregate_diagnostics["minimum_abs_lloyd_margin"] = min(
                    aggregate_diagnostics["minimum_abs_lloyd_margin"],
                    float(tested["minimum_abs_lloyd_margin"]),
                )
                aggregate_diagnostics["minimum_abs_hartigan_margin"] = min(
                    aggregate_diagnostics["minimum_abs_hartigan_margin"],
                    float(tested["minimum_abs_hartigan_margin"]),
                )
                trials_done += batch

            exists_not_lloyd_fixed = trial_count - all_balanced_successes
            lloyd_event_ci = wilson_interval(exists_not_lloyd_fixed, trial_count)
            hartigan_event_ci = wilson_interval(
                any_incorrect_hartigan_successes, trial_count
            )
            cor38_bound, cor38_log = probability_bound(
                n * math.log(2.0) + math.log(n), rho_q, dimension
            )
            cor312_bound, cor312_log = probability_bound(
                n * math.log(2.0), rho_h, dimension
            )
            row = {
                "claim_group": "claims_1_to_4",
                "regime": regime,
                "n": n,
                "dimension": dimension,
                "tau2": tau2,
                "sigma2": sigma2,
                "q": q,
                "beta": beta,
                "trials": trial_count,
                "total_unlabeled_partitions": int(partitions.shape[0]),
                "q_balanced_partitions": int(balanced.sum()),
                "incorrect_partitions": int(incorrect.sum()),
                "prob_all_q_balanced_lloyd_fixed": all_balanced_successes / trial_count,
                "prob_all_q_balanced_lloyd_fixed_ci_low": wilson_interval(
                    all_balanced_successes, trial_count
                )[0],
                "prob_all_q_balanced_lloyd_fixed_ci_high": wilson_interval(
                    all_balanced_successes, trial_count
                )[1],
                "mean_fraction_q_balanced_lloyd_fixed": balanced_fixed_fraction_sum
                / trial_count,
                "mean_fraction_all_lloyd_fixed": all_lloyd_fixed_fraction_sum / trial_count,
                "prob_exists_q_balanced_not_lloyd_fixed": exists_not_lloyd_fixed
                / trial_count,
                "prob_exists_q_balanced_not_lloyd_fixed_ci_low": lloyd_event_ci[0],
                "prob_exists_q_balanced_not_lloyd_fixed_ci_high": lloyd_event_ci[1],
                "corollary_38_applies": cor38_applies,
                "rho_q": rho_q if cor38_applies else None,
                "corollary_38_bound": cor38_bound if cor38_applies else None,
                "corollary_38_log_uncapped_bound": cor38_log if cor38_applies else None,
                "prob_any_incorrect_hartigan_fixed": any_incorrect_hartigan_successes
                / trial_count,
                "prob_any_incorrect_hartigan_fixed_ci_low": hartigan_event_ci[0],
                "prob_any_incorrect_hartigan_fixed_ci_high": hartigan_event_ci[1],
                "mean_fraction_incorrect_hartigan_fixed": incorrect_hartigan_fraction_sum
                / trial_count,
                "prob_correct_partition_hartigan_fixed": correct_hartigan_successes
                / trial_count,
                "rho_h": rho_h,
                "corollary_312_bound": cor312_bound,
                "corollary_312_log_uncapped_bound": cor312_log,
            }
            rows.append(row)

    if aggregate_diagnostics["hartigan_not_lloyd_violations"]:
        raise AssertionError("Hartigan fixed points must be a subset of Lloyd fixed points")
    return rows, {
        "true_labels": true_labels.tolist(),
        "high_sigma2": high_sigma2,
        "low_sigma2": low_sigma2,
        "rho_q": rho_q,
        "partition_count": int(partitions.shape[0]),
        "q_balanced_partition_count": int(balanced.sum()),
        "diagnostics": aggregate_diagnostics,
    }


def theorem_34_sigma2_threshold(c: int, c_bar: int, tau2: float) -> float:
    """Square of Eq. (13): sigma > sqrt(2*c_bar)*tau*(c-1)/sqrt(c(c+c_bar))."""

    return 2.0 * c_bar * tau2 * (c - 1.0) ** 2 / (c * (c + c_bar))


def theorem_34_rho(sigma2: float, tau2: float, c: int, c_bar: int) -> float:
    numerator = (
        4.0
        * sigma2
        * (c - 1.0)
        * c**2
        * c_bar
        * (c_bar + 1.0)
        * (c * (sigma2 + 2.0 * tau2) - 2.0 * tau2)
    )
    denominator = (
        -c * (sigma2 + 4.0 * tau2) * c_bar
        + c**2 * (sigma2 + 2.0 * (sigma2 + tau2) * c_bar)
        + 2.0 * tau2 * c_bar
    ) ** 2
    rho = numerator / denominator
    # At the (excluded) equality sigma^2=sigma_0^2, rho is exactly one.
    # Permit that control row while retaining a guard against real drift.
    if -1e-14 <= rho <= 1.0 + 1e-14:
        return min(1.0, max(0.0, rho))
    raise ArithmeticError(f"Theorem 3.4 rho outside [0,1]: {rho}")


def theorem_39_rho(
    sigma2: float,
    tau2: float,
    c: int,
    c_bar: int,
    purity_current: float,
    purity_alternative: float,
) -> float:
    current_term = c / (c - 1.0) * (1.0 - purity_current) ** 2
    alternative_term = c_bar / (c_bar + 1.0) * (1.0 - purity_alternative) ** 2
    ratio = tau2 * (current_term - alternative_term) / (
        tau2 * (current_term + alternative_term) + sigma2
    )
    rho = 1.0 - ratio * ratio
    if not 0.0 <= rho < 1.0:
        raise ArithmeticError(f"Theorem 3.9 rho outside [0,1): {rho}")
    return rho


def _claim5_designs() -> list[dict[str, Any]]:
    n = 40
    true_labels = np.repeat(np.arange(2), n // 2)
    official = np.concatenate(
        (
            np.zeros(5, dtype=int),
            np.ones(15, dtype=int),
            np.ones(5, dtype=int),
            np.zeros(15, dtype=int),
        )
    )
    single_error = true_labels.copy()
    single_error[0] = 1
    return [
        {
            "name": "paper_balanced_purity_025_075",
            "true_labels": true_labels,
            "current_labels": official,
            "point": 0,
        },
        {
            "name": "single_misclassified_point",
            "true_labels": true_labels,
            "current_labels": single_error,
            "point": 0,
        },
    ]


def _single_point_covariance(
    true_labels: np.ndarray,
    current_labels: np.ndarray,
    point: int,
    tau2: float,
    sigma2: float,
) -> tuple[np.ndarray, dict[str, Any]]:
    current_label = int(current_labels[point])
    alternative_label = 1 - current_label
    current_members = current_labels == current_label
    alternative_members = current_labels == alternative_label
    c = int(current_members.sum())
    c_bar = int(alternative_members.sum())
    if c <= 1 or c_bar <= 0:
        raise ValueError("Claim 5 designs require a movable point and nonempty alternative")
    coefficients = np.zeros((2, true_labels.size), dtype=np.float64)
    coefficients[:, point] = 1.0
    coefficients[0, current_members] -= 1.0 / c
    coefficients[1, alternative_members] -= 1.0 / c_bar
    covariance = model_covariance(true_labels, tau2, sigma2)
    reduced = coefficients @ covariance @ coefficients.T
    true_class = true_labels[point]
    purity_current = float(np.mean(true_labels[current_members] == true_class))
    purity_alternative = float(np.mean(true_labels[alternative_members] == true_class))
    return reduced, {
        "c": c,
        "c_bar": c_bar,
        "purity_current": purity_current,
        "purity_alternative": purity_alternative,
    }


def run_claim_5(
    rng: np.random.Generator,
    profile: Profile,
    *,
    tau2: float = 1.0,
) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    rows: list[dict[str, Any]] = []
    diagnostics = {
        "lloyd_near_ties": 0,
        "hartigan_near_ties": 0,
        "comparisons": 0,
    }
    design_summaries: list[dict[str, Any]] = []
    for design in _claim5_designs():
        # Geometry-derived cluster counts/purities do not depend on sigma2.
        _, geometry = _single_point_covariance(
            design["true_labels"],
            design["current_labels"],
            design["point"],
            tau2,
            1.0,
        )
        c = geometry["c"]
        c_bar = geometry["c_bar"]
        threshold_sigma2 = theorem_34_sigma2_threshold(c, c_bar, tau2)
        design_summaries.append(
            {
                "design": design["name"],
                **geometry,
                "sigma2_threshold": threshold_sigma2,
            }
        )
        for ratio in profile.claim5_ratios:
            sigma2 = ratio * threshold_sigma2
            reduced_covariance, checked_geometry = _single_point_covariance(
                design["true_labels"],
                design["current_labels"],
                design["point"],
                tau2,
                sigma2,
            )
            if checked_geometry != geometry:
                raise AssertionError("cluster geometry unexpectedly changed")
            rho_lloyd = theorem_34_rho(sigma2, tau2, c, c_bar)
            rho_hartigan = theorem_39_rho(
                sigma2,
                tau2,
                c,
                c_bar,
                geometry["purity_current"],
                geometry["purity_alternative"],
            )
            for dimension in profile.claim5_dimensions:
                lloyd_stays = 0
                hartigan_stays = 0
                done = 0
                while done < profile.claim5_trials:
                    batch = min(profile.wishart_batch_size, profile.claim5_trials - done)
                    wishart2 = sample_wishart_normalized(
                        rng, reduced_covariance, dimension, batch
                    )
                    distance_current = wishart2[:, 0, 0]
                    distance_alternative = wishart2[:, 1, 1]
                    lloyd_margin = distance_alternative - distance_current
                    hartigan_current = c / (c - 1.0) * distance_current
                    hartigan_alternative = c_bar / (c_bar + 1.0) * distance_alternative
                    hartigan_margin = hartigan_alternative - hartigan_current
                    # Current assignment is wrong in both designs.  Staying means
                    # no strict improvement under the relevant criterion.
                    lloyd_stays += int((lloyd_margin >= 0.0).sum())
                    hartigan_stays += int((hartigan_margin >= 0.0).sum())
                    diagnostics["lloyd_near_ties"] += int(
                        _near_tie_mask(
                            lloyd_margin, distance_current, distance_alternative
                        ).sum()
                    )
                    diagnostics["hartigan_near_ties"] += int(
                        _near_tie_mask(
                            hartigan_margin, hartigan_current, hartigan_alternative
                        ).sum()
                    )
                    diagnostics["comparisons"] += 2 * batch
                    done += batch

                lloyd_ci = wilson_interval(lloyd_stays, profile.claim5_trials)
                hartigan_ci = wilson_interval(hartigan_stays, profile.claim5_trials)
                log_lloyd_tail = (dimension / 4.0) * math.log(rho_lloyd)
                lloyd_lower_bound = -math.expm1(log_lloyd_tail)
                hartigan_upper_bound = math.exp((dimension / 4.0) * math.log(rho_hartigan))
                rows.append(
                    {
                        "claim_group": "claim_5",
                        "design": design["name"],
                        "n": int(design["true_labels"].size),
                        "dimension": dimension,
                        "tau2": tau2,
                        "sigma2": sigma2,
                        "sigma2_threshold": threshold_sigma2,
                        "sigma2_threshold_ratio": ratio,
                        "theorem_34_strict_condition": ratio > 1.0,
                        "c": c,
                        "c_bar": c_bar,
                        "purity_current": geometry["purity_current"],
                        "purity_alternative": geometry["purity_alternative"],
                        "trials": profile.claim5_trials,
                        "prob_lloyd_stays_wrong": lloyd_stays / profile.claim5_trials,
                        "prob_lloyd_stays_wrong_ci_low": lloyd_ci[0],
                        "prob_lloyd_stays_wrong_ci_high": lloyd_ci[1],
                        "rho_lloyd": rho_lloyd,
                        "theorem_34_lloyd_stay_lower_bound": lloyd_lower_bound
                        if ratio > 1.0
                        else None,
                        "prob_hartigan_stays_wrong": hartigan_stays
                        / profile.claim5_trials,
                        "prob_hartigan_stays_wrong_ci_low": hartigan_ci[0],
                        "prob_hartigan_stays_wrong_ci_high": hartigan_ci[1],
                        "rho_hartigan": rho_hartigan,
                        "theorem_39_hartigan_stay_upper_bound": hartigan_upper_bound,
                    }
                )
    return rows, {"designs": design_summaries, "diagnostics": diagnostics}


def _direct_fixed_tests(data: np.ndarray, partitions: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    lloyd_results = []
    hartigan_results = []
    for partition in partitions.astype(int):
        centers = np.stack([data[partition == k].mean(axis=0) for k in (0, 1)])
        squared = ((data[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
        current = squared[np.arange(data.shape[0]), partition]
        alternative = squared[np.arange(data.shape[0]), 1 - partition]
        lloyd_results.append(bool(np.all(current <= alternative)))
        sizes = np.bincount(partition, minlength=2)
        current_size = sizes[partition]
        alternative_size = sizes[1 - partition]
        movable = current_size > 1
        current_weighted = np.zeros_like(current)
        current_weighted[movable] = (
            current[movable] * current_size[movable] / (current_size[movable] - 1.0)
        )
        alternative_weighted = (
            alternative * alternative_size / (alternative_size + 1.0)
        )
        hartigan_results.append(
            bool(np.all((~movable) | (current_weighted <= alternative_weighted)))
        )
    return np.asarray(lloyd_results), np.asarray(hartigan_results)


def run_self_checks(seed: int = 29589) -> dict[str, Any]:
    """Algebraic and numerical checks independent of Monte Carlo outcomes."""

    rng = np.random.default_rng(seed + 991)
    n = 8
    d = 64
    tau2 = 1.0
    sigma2 = 9.64053853187796
    true_labels = np.repeat(np.arange(2), n // 2)
    means = rng.normal(scale=math.sqrt(tau2), size=(2, d))
    data = means[true_labels] + rng.normal(scale=math.sqrt(sigma2), size=(n, d))
    partitions = enumerate_bipartitions(n)
    gram = (data @ data.T / d)[None, :, :]
    distance0, distance1 = gram_partition_distances(gram, partitions)

    maximum_distance_error = 0.0
    for pidx, partition in enumerate(partitions.astype(int)):
        for cluster in (0, 1):
            center = data[partition == cluster].mean(axis=0)
            direct = ((data - center) ** 2).sum(axis=1) / d
            formula = distance0[0, pidx] if cluster == 0 else distance1[0, pidx]
            maximum_distance_error = max(
                maximum_distance_error, float(np.max(np.abs(direct - formula)))
            )
    if maximum_distance_error > 2e-12:
        raise AssertionError(f"Gram/direct centroid mismatch: {maximum_distance_error}")

    gram_fixed = fixed_point_tests(gram, partitions)
    direct_lloyd, direct_hartigan = _direct_fixed_tests(data, partitions)
    if not np.array_equal(np.asarray(gram_fixed["lloyd_fixed"])[0], direct_lloyd):
        raise AssertionError("direct and Gram Lloyd fixed-point tests disagree")
    if not np.array_equal(np.asarray(gram_fixed["hartigan_fixed"])[0], direct_hartigan):
        raise AssertionError("direct and Gram Hartigan fixed-point tests disagree")

    # Check the Hartigan weighted-distance identity against direct WCSS delta.
    partition = partitions[37].astype(int)
    point = next(i for i in range(n) if np.sum(partition == partition[i]) > 1)
    current_label = partition[point]
    alternative_label = 1 - current_label
    sizes = np.bincount(partition, minlength=2)
    centers = np.stack([data[partition == k].mean(axis=0) for k in (0, 1)])
    before = float(((data - centers[partition]) ** 2).sum())
    current_distance = float(((data[point] - centers[current_label]) ** 2).sum())
    alternative_distance = float(((data[point] - centers[alternative_label]) ** 2).sum())
    expected_delta = (
        sizes[alternative_label] / (sizes[alternative_label] + 1.0) * alternative_distance
        - sizes[current_label] / (sizes[current_label] - 1.0) * current_distance
    )
    moved = partition.copy()
    moved[point] = alternative_label
    moved_centers = np.stack([data[moved == k].mean(axis=0) for k in (0, 1)])
    after = float(((data - moved_centers[moved]) ** 2).sum())
    wcss_delta_error = abs((after - before) - expected_delta)
    if wcss_delta_error > 2e-10:
        raise AssertionError(f"Hartigan WCSS identity mismatch: {wcss_delta_error}")

    # Regression constants calculated from the paper's displayed equations.
    calculated_sigma2, rho_q = corollary_38_parameters(8, 1.5, 1.5)
    if not math.isclose(calculated_sigma2, 9.64053853187796, rel_tol=2e-14):
        raise AssertionError("Corollary 3.8 sigma regression failed")
    if not math.isclose(rho_q, 0.9928404727106878, rel_tol=2e-14):
        raise AssertionError("Corollary 3.8 rho regression failed")
    rho_h = rho_h_uniform(8, 1.0, calculated_sigma2, 0.5)
    if not math.isclose(rho_h, 0.9999022112550369, rel_tol=2e-14):
        raise AssertionError("Corollary 3.12 rho regression failed")
    balanced_count = int(q_balanced_mask(partitions, 1.5).sum())
    if partitions.shape != (127, 8) or balanced_count != 91:
        raise AssertionError(
            f"partition enumeration regression failed: {partitions.shape}, {balanced_count}"
        )

    lo, hi = wilson_interval(0, 100)
    if lo != 0.0 or not (0.036 < hi < 0.038):
        raise AssertionError("Wilson interval regression failed")

    return {
        "status": "passed",
        "direct_x_dimension": d,
        "partitions_checked": int(partitions.shape[0]),
        "maximum_distance_error": maximum_distance_error,
        "wcss_delta_error": wcss_delta_error,
        "rho_q_regression": rho_q,
        "rho_h_regression": rho_h,
        "q_balanced_partitions": balanced_count,
    }


def write_csv(path: Path, rows: Sequence[dict[str, Any]]) -> None:
    if not rows:
        raise ValueError("cannot write an empty CSV")
    fieldnames: list[str] = []
    seen = set()
    for row in rows:
        for key in row:
            if key not in seen:
                fieldnames.append(key)
                seen.add(key)
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.DictWriter(handle, fieldnames=fieldnames)
        writer.writeheader()
        writer.writerows(rows)


def _svg_polyline(
    points: Sequence[tuple[float, float]], color: str, *, dashed: bool = False
) -> str:
    encoded = " ".join(f"{x:.2f},{y:.2f}" for x, y in points)
    dash = ' stroke-dasharray="7 5"' if dashed else ""
    circles = "".join(
        f'<circle cx="{x:.2f}" cy="{y:.2f}" r="3.0" fill="{color}"/>'
        for x, y in points
    )
    return (
        f'<polyline points="{encoded}" fill="none" stroke="{color}" '
        f'stroke-width="2.3"{dash}/>{circles}'
    )


def _panel_svg(
    *,
    x: float,
    y: float,
    width: float,
    height: float,
    title: str,
    series: Sequence[dict[str, Any]],
    x_label: str = "dimension d (log scale)",
    y_label: str = "probability",
) -> str:
    left, right, top, bottom = 62.0, 18.0, 42.0, 52.0
    plot_x = x + left
    plot_y = y + top
    plot_w = width - left - right
    plot_h = height - top - bottom
    all_x = [point[0] for item in series for point in item["values"]]
    if not all_x:
        return ""
    log_min = math.log10(min(all_x))
    log_max = math.log10(max(all_x))
    if log_max == log_min:
        log_max += 1.0

    def sx(value: float) -> float:
        return plot_x + (math.log10(value) - log_min) / (log_max - log_min) * plot_w

    def sy(value: float) -> float:
        return plot_y + (1.0 - min(1.0, max(0.0, value))) * plot_h

    parts = [
        f'<g><rect x="{x}" y="{y}" width="{width}" height="{height}" rx="8" '
        'fill="#ffffff" stroke="#d8dee9"/>',
        f'<text x="{x + width / 2}" y="{y + 25}" text-anchor="middle" '
        f'font-size="15" font-weight="700">{html.escape(title)}</text>',
    ]
    for tick in (0.0, 0.25, 0.5, 0.75, 1.0):
        py = sy(tick)
        parts.append(
            f'<line x1="{plot_x}" y1="{py}" x2="{plot_x + plot_w}" y2="{py}" '
            'stroke="#e5e9f0"/>'
        )
        parts.append(
            f'<text x="{plot_x - 8}" y="{py + 4}" text-anchor="end" '
            f'font-size="11">{tick:g}</text>'
        )
    tick_powers = range(math.floor(log_min), math.ceil(log_max) + 1)
    for power in tick_powers:
        value = 10.0**power
        if min(all_x) <= value <= max(all_x):
            px = sx(value)
            parts.append(
                f'<line x1="{px}" y1="{plot_y}" x2="{px}" y2="{plot_y + plot_h}" '
                'stroke="#edf0f5"/>'
            )
            parts.append(
                f'<text x="{px}" y="{plot_y + plot_h + 18}" text-anchor="middle" '
                f'font-size="11">10^{power}</text>'
            )
    for item in series:
        plotted = [(sx(px), sy(py)) for px, py in item["values"]]
        parts.append(_svg_polyline(plotted, item["color"], dashed=item.get("dashed", False)))
    legend_x = plot_x + 5
    legend_y = plot_y + 15
    for index, item in enumerate(series):
        ly = legend_y + index * 17
        dash = ' stroke-dasharray="6 4"' if item.get("dashed", False) else ""
        parts.append(
            f'<line x1="{legend_x}" y1="{ly}" x2="{legend_x + 22}" y2="{ly}" '
            f'stroke="{item["color"]}" stroke-width="2.3"{dash}/>'
        )
        parts.append(
            f'<text x="{legend_x + 28}" y="{ly + 4}" font-size="10">'
            f'{html.escape(item["label"])}</text>'
        )
    parts.extend(
        [
            f'<text x="{plot_x + plot_w / 2}" y="{y + height - 12}" '
            f'text-anchor="middle" font-size="12">{html.escape(x_label)}</text>',
            f'<text x="{x + 15}" y="{plot_y + plot_h / 2}" text-anchor="middle" '
            f'font-size="12" transform="rotate(-90 {x + 15} {plot_y + plot_h / 2})">'
            f'{html.escape(y_label)}</text>',
            "</g>",
        ]
    )
    return "".join(parts)


def build_claims_1_4_svg(rows: Sequence[dict[str, Any]]) -> str:
    high = sorted(
        (row for row in rows if row["regime"] == "high_noise"),
        key=lambda row: row["dimension"],
    )
    low = sorted(
        (row for row in rows if row["regime"] == "low_noise_control"),
        key=lambda row: row["dimension"],
    )
    panel1 = [
        {
            "label": "high noise: empirical all fixed",
            "color": "#2563eb",
            "values": [(row["dimension"], row["prob_all_q_balanced_lloyd_fixed"]) for row in high],
        },
        {
            "label": "Cor. 3.8 guaranteed lower bound",
            "color": "#1d4ed8",
            "dashed": True,
            "values": [(row["dimension"], 1.0 - row["corollary_38_bound"]) for row in high],
        },
        {
            "label": "low-noise control",
            "color": "#94a3b8",
            "values": [(row["dimension"], row["prob_all_q_balanced_lloyd_fixed"]) for row in low],
        },
    ]
    panel2 = [
        {
            "label": "high noise: empirical any incorrect",
            "color": "#ea580c",
            "values": [(row["dimension"], row["prob_any_incorrect_hartigan_fixed"]) for row in high],
        },
        {
            "label": "Cor. 3.12 upper bound",
            "color": "#c2410c",
            "dashed": True,
            "values": [(row["dimension"], row["corollary_312_bound"]) for row in high],
        },
        {
            "label": "low-noise control",
            "color": "#94a3b8",
            "values": [(row["dimension"], row["prob_any_incorrect_hartigan_fixed"]) for row in low],
        },
    ]
    body = _panel_svg(
        x=18, y=55, width=555, height=350, title="Claims 1 & 3 — Lloyd fixed-point proliferation", series=panel1
    ) + _panel_svg(
        x=587, y=55, width=555, height=350, title="Claims 2 & 4 — incorrect Hartigan fixed points", series=panel2
    )
    return (
        '<svg xmlns="http://www.w3.org/2000/svg" width="1160" height="425" '
        'viewBox="0 0 1160 425" role="img" aria-label="Claims 1 to 4 numerical audit">'
        '<rect width="1160" height="425" fill="#f8fafc"/>'
        '<text x="580" y="30" text-anchor="middle" font-family="system-ui,sans-serif" '
        'font-size="20" font-weight="700">Exact exhaustive Wishart/Gram audit</text>'
        f'<g font-family="system-ui,sans-serif" fill="#172033">{body}</g></svg>'
    )


def build_claim_5_svg(rows: Sequence[dict[str, Any]]) -> str:
    design = "paper_balanced_purity_025_075"
    selected_ratios = sorted(
        set(row["sigma2_threshold_ratio"] for row in rows if row["design"] == design)
    )
    if len(selected_ratios) > 4:
        selected_ratios = [selected_ratios[0], 1.0, selected_ratios[-2], selected_ratios[-1]]
    palette = ["#64748b", "#7c3aed", "#2563eb", "#dc2626"]
    lloyd_series = []
    hartigan_series = []
    for color, ratio in zip(palette, selected_ratios):
        subset = sorted(
            (
                row
                for row in rows
                if row["design"] == design
                and row["sigma2_threshold_ratio"] == ratio
            ),
            key=lambda row: row["dimension"],
        )
        lloyd_series.append(
            {
                "label": f"sigma²/sigma0²={ratio:g}",
                "color": color,
                "values": [(row["dimension"], row["prob_lloyd_stays_wrong"]) for row in subset],
            }
        )
        hartigan_series.append(
            {
                "label": f"sigma²/sigma0²={ratio:g}",
                "color": color,
                "values": [(row["dimension"], row["prob_hartigan_stays_wrong"]) for row in subset],
            }
        )
    body = _panel_svg(
        x=18,
        y=55,
        width=555,
        height=350,
        title="Claim 5 — Lloyd keeps the wrong assignment",
        series=lloyd_series,
    ) + _panel_svg(
        x=587,
        y=55,
        width=555,
        height=350,
        title="Control — Hartigan keeps the wrong assignment",
        series=hartigan_series,
    )
    return (
        '<svg xmlns="http://www.w3.org/2000/svg" width="1160" height="425" '
        'viewBox="0 0 1160 425" role="img" aria-label="Claim 5 numerical audit">'
        '<rect width="1160" height="425" fill="#f8fafc"/>'
        '<text x="580" y="30" text-anchor="middle" font-family="system-ui,sans-serif" '
        'font-size="20" font-weight="700">Single-sample persistence across dimension</text>'
        f'<g font-family="system-ui,sans-serif" fill="#172033">{body}</g></svg>'
    )


def _html_table(rows: Sequence[dict[str, Any]], columns: Sequence[str], limit: int | None = None) -> str:
    shown = rows if limit is None else rows[:limit]
    header = "".join(f"<th>{html.escape(column)}</th>" for column in columns)
    body = []
    for row in shown:
        cells = []
        for column in columns:
            value = row.get(column)
            if isinstance(value, float):
                rendered = f"{value:.6g}"
            elif value is None:
                rendered = "—"
            else:
                rendered = str(value)
            cells.append(f"<td>{html.escape(rendered)}</td>")
        body.append("<tr>" + "".join(cells) + "</tr>")
    return f"<table><thead><tr>{header}</tr></thead><tbody>{''.join(body)}</tbody></table>"


def build_html_report(
    summary: dict[str, Any],
    claims_1_4_rows: Sequence[dict[str, Any]],
    claim5_rows: Sequence[dict[str, Any]],
    claims_svg: str,
    claim5_svg: str,
) -> str:
    high = [row for row in claims_1_4_rows if row["regime"] == "high_noise"]
    claim5_paper = [
        row for row in claim5_rows if row["design"] == "paper_balanced_purity_025_075"
    ]
    columns_14 = (
        "regime",
        "dimension",
        "trials",
        "prob_all_q_balanced_lloyd_fixed",
        "corollary_38_bound",
        "prob_any_incorrect_hartigan_fixed",
        "corollary_312_bound",
    )
    columns_5 = (
        "design",
        "sigma2_threshold_ratio",
        "dimension",
        "trials",
        "prob_lloyd_stays_wrong",
        "theorem_34_lloyd_stay_lower_bound",
        "prob_hartigan_stays_wrong",
        "theorem_39_hartigan_stay_upper_bound",
    )
    diagnostics = summary["claims_1_to_4"]["diagnostics"]
    return f"""<!doctype html>
<html lang="en"><head><meta charset="utf-8"><meta name="viewport" content="width=device-width,initial-scale=1">
<title>Claims 1–5 exact theory audit</title>
<style>
body{{font:15px/1.5 system-ui,sans-serif;color:#172033;background:#f1f5f9;margin:0}}
main{{max-width:1200px;margin:auto;padding:28px}} section{{background:white;border-radius:12px;padding:22px;margin:18px 0;box-shadow:0 1px 3px #0001}}
h1,h2{{line-height:1.2}} code{{background:#eef2f7;padding:2px 5px;border-radius:4px}} table{{border-collapse:collapse;width:100%;font-size:12px}}
th,td{{border:1px solid #d8dee9;padding:6px;text-align:right}} th:first-child,td:first-child{{text-align:left}} th{{background:#eef2f7;position:sticky;top:0}}
.scroll{{overflow:auto;max-height:480px}} .note{{border-left:4px solid #2563eb;padding:8px 12px;background:#eff6ff}} a{{color:#1d4ed8}}
svg{{max-width:100%;height:auto}}
</style></head><body><main>
<h1>Exact Wishart/Gram numerical audit: Claims 1–5</h1>
<p>Profile <code>{html.escape(summary['profile'])}</code>, seed <code>{summary['seed']}</code>, NumPy {html.escape(summary['environment']['numpy'])}, runtime {summary['runtime_seconds']:.2f}s.</p>
<p><a href="{PAPER_URL}">Paper</a> · <a href="{CODE_URL}">author code (audited commit)</a> · <a href="{ZENODO_URL}">author archive</a></p>
<section><h2>Method and scope</h2>
<p>For fixed labels, each feature column is Gaussian with covariance <code>Σ=τ²ZZᵀ+σ²I</code>, so <code>G=XXᵀ ~ Wishart(d,Σ)</code>. Since both fixed-point criteria use only distances to empirical centroids, exhaustive Gram-space tests are distributionally exact, not a proxy. Complement-symmetric partitions are enumerated once.</p>
<p class="note">A numerical audit is not a replacement for a proof. Corollary bounds are probabilities over independent datasets. The exact event is evaluated once per dataset; candidate partitions are not treated as independent trials.</p>
<p>Built-in self-check: <strong>{html.escape(summary['self_check']['status'])}</strong>; {summary['self_check']['partitions_checked']} direct-X partitions checked, maximum Gram/direct distance error {summary['self_check']['maximum_distance_error']:.3g}.</p>
</section>
<section><h2>Claims 1–4</h2>{claims_svg}
<p>High-noise parameters: n=8, τ²=1, q=1.5, β=1.5, σ²={summary['claims_1_to_4']['high_sigma2']:.8g}. The q-balanced filter uses strict Definition 2.5 inequalities and cluster sizes &gt;2.</p>
<p>Near-tie diagnostics: {diagnostics['lloyd_near_ties']} Lloyd and {diagnostics['hartigan_near_ties']} Hartigan comparisons within {TIE_ULPS:g} scaled ulps; {diagnostics['negative_distance_count']} negative computed distances. Hartigan-fixed/not-Lloyd-fixed violations: {diagnostics['hartigan_not_lloyd_violations']}.</p>
<div class="scroll">{_html_table(claims_1_4_rows, columns_14)}</div></section>
<section><h2>Claim 5</h2>{claim5_svg}
<p>The paper-matched design uses n=40, τ²=1, c=c̄=20 and purities 0.25/0.75, giving σ₀²=18.05. A second design contains exactly one misclassified point. The Theorem 3.4 lower bound is shown only when the strict condition σ²&gt;σ₀² holds.</p>
<div class="scroll">{_html_table(claim5_paper, columns_5)}</div></section>
<section><h2>Artifacts</h2><p>Machine-readable values are in <code>claims_1_4_results.csv</code>, <code>claim_5_results.csv</code>, and <code>summary.json</code>. Standalone figures are <code>claims_1_4.svg</code> and <code>claim_5.svg</code>.</p></section>
</main></body></html>"""


def _summary_findings(rows_14: Sequence[dict[str, Any]], rows_5: Sequence[dict[str, Any]]) -> dict[str, Any]:
    high = sorted(
        (row for row in rows_14 if row["regime"] == "high_noise"),
        key=lambda row: row["dimension"],
    )
    high_last = high[-1]
    claim5_above = [
        row
        for row in rows_5
        if row["design"] == "paper_balanced_purity_025_075"
        and row["sigma2_threshold_ratio"] > 1.0
    ]
    claim5_last = max(claim5_above, key=lambda row: row["dimension"])
    return {
        "claim_1_highest_dimension_prob_all_balanced_fixed": high_last[
            "prob_all_q_balanced_lloyd_fixed"
        ],
        "claim_2_highest_dimension_prob_any_incorrect_hartigan_fixed": high_last[
            "prob_any_incorrect_hartigan_fixed"
        ],
        "claim_3_highest_dimension_empirical_bad_event": high_last[
            "prob_exists_q_balanced_not_lloyd_fixed"
        ],
        "claim_3_highest_dimension_bound": high_last["corollary_38_bound"],
        "claim_4_highest_dimension_empirical_bad_event": high_last[
            "prob_any_incorrect_hartigan_fixed"
        ],
        "claim_4_highest_dimension_bound": high_last["corollary_312_bound"],
        "claim_5_high_dimension_lloyd_stays_wrong": claim5_last[
            "prob_lloyd_stays_wrong"
        ],
        "claim_5_high_dimension_hartigan_stays_wrong": claim5_last[
            "prob_hartigan_stays_wrong"
        ],
    }


def parse_args(argv: Sequence[str] | None = None) -> argparse.Namespace:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--profile", choices=sorted(PROFILES), default="smoke")
    parser.add_argument("--output-dir", type=Path, default=Path("outputs/theory"))
    parser.add_argument("--seed", type=int, default=29589)
    parser.add_argument(
        "--self-check-only",
        action="store_true",
        help="run algebra/regression checks without creating result files",
    )
    return parser.parse_args(argv)


def main(argv: Sequence[str] | None = None) -> int:
    args = parse_args(argv)
    started = time.perf_counter()
    self_check = run_self_checks(args.seed)
    if args.self_check_only:
        print(json.dumps(_native(self_check), indent=2, allow_nan=False))
        return 0

    profile = PROFILES[args.profile]
    # Independent deterministic streams keep one claim's profile changes from
    # silently changing another claim's samples.
    seed_sequence = np.random.SeedSequence(args.seed)
    stream_14, stream_5 = seed_sequence.spawn(2)
    rows_14, metadata_14 = run_claims_1_to_4(
        np.random.default_rng(stream_14), profile
    )
    rows_5, metadata_5 = run_claim_5(np.random.default_rng(stream_5), profile)
    runtime = time.perf_counter() - started

    args.output_dir.mkdir(parents=True, exist_ok=True)
    claims_svg = build_claims_1_4_svg(rows_14)
    claim5_svg = build_claim_5_svg(rows_5)
    summary = {
        "schema_version": 1,
        "paper": PAPER_URL,
        "author_code": CODE_URL,
        "author_archive": ZENODO_URL,
        "profile": profile.name,
        "profile_parameters": asdict(profile),
        "seed": args.seed,
        "runtime_seconds": runtime,
        "environment": {
            "python": platform.python_version(),
            "numpy": np.__version__,
            "platform": platform.platform(),
            "pid": os.getpid(),
        },
        "self_check": self_check,
        "claims_1_to_4": metadata_14,
        "claim_5": metadata_5,
        "findings": _summary_findings(rows_14, rows_5),
        "definitions": {
            "lloyd_fixed": "every point has current squared distance <= alternative squared distance",
            "hartigan_fixed": "every allowed source-nonsingleton move has weighted current distance <= weighted alternative distance",
            "corollary_38_event": "exists q-approximately balanced partition that is not Lloyd-fixed",
            "corollary_312_event": "exists nonempty incorrect partition that is Hartigan-fixed",
            "ties": "strict improvements move; exact ties stay and near-ties are diagnosed without altering classification",
        },
    }
    report = build_html_report(summary, rows_14, rows_5, claims_svg, claim5_svg)

    write_csv(args.output_dir / "claims_1_4_results.csv", rows_14)
    write_csv(args.output_dir / "claim_5_results.csv", rows_5)
    (args.output_dir / "claims_1_4.svg").write_text(claims_svg, encoding="utf-8")
    (args.output_dir / "claim_5.svg").write_text(claim5_svg, encoding="utf-8")
    (args.output_dir / "theory_reproduction.html").write_text(report, encoding="utf-8")
    (args.output_dir / "summary.json").write_text(
        json.dumps(_native(summary), indent=2, sort_keys=True, allow_nan=False) + "\n",
        encoding="utf-8",
    )

    print(f"Self-check: {self_check['status']}")
    print(f"Profile: {profile.name}; seed: {args.seed}; runtime: {runtime:.2f}s")
    print(f"Results: {args.output_dir.resolve()}")
    print(json.dumps(_native(summary["findings"]), indent=2, allow_nan=False))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())