import argparse
import csv
import hashlib
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np

SEED = 20260925
M = 1000
REPLICATES = 1000
ACCEPTANCE_RATE = 0.20
REVIEWERS = 3
OMEGAS = (0.0, 0.25, 0.50, 0.75, 1.0)
ALPHA = 0.80
SIGMA_AI = 0.40
KAPPA = 0.10
DELTA = 0.10
SEVERITY_SD = 0.20
SPECIFIC_NOISE_SD = 0.40
MU = 0.0
ALLOCATION_OMEGA = 0.50
EXTRA_FRACTIONS = (0.25, 0.50)
MECHANISM_SIGMAS = (0.0, 0.20, 0.40, 0.80, 1.20)
ABILITY_PROFILES = {
    "low": (0.20, 0.35, 0.50),
    "middle": (0.35, 0.55, 0.75),
    "high": (0.50, 0.70, 0.90),
}
SEVERITY_SDS = (0.0, 0.10, 0.20, 0.50)
SPECIFIC_NOISE_SDS = (0.10, 0.40, 0.80, 1.00)
STRATEGIES = ("uniform_extra", "boundary_min_one", "boundary_reallocating")
ROOT = Path(__file__).resolve().parents[1]
RESULTS_DIR = ROOT / "results"
FIGURES_DIR = ROOT / "figures"


def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument("--seed", type=int, default=SEED)
    return parser.parse_args()


def rng_for(seed_sequence):
    return np.random.default_rng(seed_sequence)


def balanced_groups(rng, count):
    groups = np.empty(count, dtype=np.int8)
    groups[: count // 2] = -1
    groups[(count + 1) // 2 :] = 1
    rng.shuffle(groups)
    return groups


def threshold_and_truth(quality):
    accepted_count = int(round(M * ACCEPTANCE_RATE))
    threshold = float(np.partition(quality, M - accepted_count)[M - accepted_count])
    accepted = quality >= threshold
    if int(accepted.sum()) != accepted_count:
        order = np.argsort(quality, kind="mergesort")
        accepted = np.zeros(M, dtype=bool)
        accepted[order[-accepted_count:]] = True
    return threshold, accepted


def base_review_data(rng, ability_values, severity_sd, specific_noise_sd, ai_sigma, error_unit=None):
    quality = rng.standard_normal(M)
    groups = balanced_groups(rng, M)
    abilities = rng.permutation(np.asarray(ability_values, dtype=float))
    severity = rng.normal(0.0, severity_sd, REVIEWERS)
    noise = rng.normal(0.0, specific_noise_sd, (M, REVIEWERS))
    base_human = quality[:, None] * abilities[None, :] + severity[None, :] + noise
    if error_unit is None:
        error_unit = rng.standard_normal(M)
    ai = ALPHA * quality + KAPPA * groups + ai_sigma * error_unit
    return quality, groups, abilities, severity, base_human, ai


def summarize(values):
    array = np.asarray(values, dtype=float)
    return {
        "mean": float(np.mean(array)),
        "median": float(np.median(array)),
        "q025": float(np.quantile(array, 0.025)),
        "q975": float(np.quantile(array, 0.975)),
    }


def metric_columns(prefix, values):
    summary = summarize(values)
    return {
        f"{prefix}_mean": summary["mean"],
        f"{prefix}_median": summary["median"],
        f"{prefix}_q025": summary["q025"],
        f"{prefix}_q975": summary["q975"],
    }


def review_metrics(assisted, final_score, accepted, groups, quality, theta):
    consensus = np.mean(np.all((assisted >= theta) == (assisted[:, :1] >= theta), axis=1))
    plus = groups == 1
    return {
        "mae": np.mean(np.abs(final_score - quality)),
        "decision_error": np.mean((final_score >= theta) != accepted),
        "wsd": np.mean(np.std(assisted, axis=1, ddof=1)),
        "decision_consensus": consensus,
        "group_mae_plus": np.mean(np.abs(final_score[plus] - quality[plus])),
        "group_mae_minus": np.mean(np.abs(final_score[~plus] - quality[~plus])),
        "decision_acceptance_rate": np.mean(final_score >= theta),
    }


def group_bias_for_full_review(base_human, quality, error_unit, omega):
    human_plus = (1.0 - omega) * (np.mean(base_human, axis=1) + DELTA)
    human_minus = (1.0 - omega) * (np.mean(base_human, axis=1) - DELTA)
    ai_plus = ALPHA * quality + KAPPA + SIGMA_AI * error_unit
    ai_minus = ALPHA * quality - KAPPA + SIGMA_AI * error_unit
    return np.mean((human_plus + omega * ai_plus) - (human_minus + omega * ai_minus))


def run_omega(replicates, seed_sequence):
    rng = rng_for(seed_sequence)
    metric_names = (
        "mae",
        "decision_error",
        "wsd",
        "decision_consensus",
        "group_bias",
        "group_mae_plus",
        "group_mae_minus",
        "decision_acceptance_rate",
        "theta",
    )
    values = {name: np.empty((len(OMEGAS), replicates), dtype=float) for name in metric_names}
    base_abilities = (0.20, 0.50, 0.80)
    for replicate in range(replicates):
        quality, groups, _, _, base_human, ai = base_review_data(
            rng, base_abilities, SEVERITY_SD, SPECIFIC_NOISE_SD, SIGMA_AI
        )
        theta, accepted = threshold_and_truth(quality)
        human = base_human + DELTA * groups[:, None]
        for index, omega in enumerate(OMEGAS):
            assisted = (1.0 - omega) * human + omega * ai[:, None]
            final_score = np.mean(assisted, axis=1)
            metrics = review_metrics(assisted, final_score, accepted, groups, quality, theta)
            metrics["group_bias"] = 2.0 * ((1.0 - omega) * DELTA + omega * KAPPA)
            metrics["theta"] = theta
            for name in metric_names:
                values[name][index, replicate] = metrics[name]
    rows = []
    for index, omega in enumerate(OMEGAS):
        row = {"omega": float(omega)}
        for name in metric_names:
            row.update(metric_columns(name, values[name][index]))
        row["human_coverage_mean"] = 1.0
        row["avg_reviews_per_manuscript_mean"] = float(REVIEWERS)
        rows.append(row)
    return rows


def run_parameter_sensitivity(replicates, seed_sequence):
    rng = rng_for(seed_sequence)
    cell_names = [
        (profile_name, severity_sd, noise_sd)
        for profile_name in ABILITY_PROFILES
        for severity_sd in SEVERITY_SDS
        for noise_sd in SPECIFIC_NOISE_SDS
    ]
    metric_names = (
        "mae",
        "decision_error",
        "wsd",
        "decision_consensus",
        "group_bias",
        "group_mae_plus",
        "group_mae_minus",
    )
    values = {name: np.empty((len(cell_names), replicates), dtype=float) for name in metric_names}
    cell_index = {cell: index for index, cell in enumerate(cell_names)}
    for replicate in range(replicates):
        quality = rng.standard_normal(M)
        groups = balanced_groups(rng, M)
        theta, accepted = threshold_and_truth(quality)
        severity_unit = rng.standard_normal(REVIEWERS)
        noise_unit = rng.standard_normal((M, REVIEWERS))
        ai_error = SIGMA_AI * rng.standard_normal(M)
        ai = ALPHA * quality + KAPPA * groups + ai_error
        for profile_name, ability_values in ABILITY_PROFILES.items():
            abilities = rng.permutation(np.asarray(ability_values, dtype=float))
            for severity_sd in SEVERITY_SDS:
                severity = severity_sd * severity_unit
                for noise_sd in SPECIFIC_NOISE_SDS:
                    base_human = (
                        quality[:, None] * abilities[None, :]
                        + severity[None, :]
                        + noise_sd * noise_unit
                    )
                    human = base_human + DELTA * groups[:, None]
                    assisted = (1.0 - ALLOCATION_OMEGA) * human + ALLOCATION_OMEGA * ai[:, None]
                    final_score = np.mean(assisted, axis=1)
                    metrics = review_metrics(assisted, final_score, accepted, groups, quality, theta)
                    metrics["group_bias"] = 2.0 * (
                        (1.0 - ALLOCATION_OMEGA) * DELTA + ALLOCATION_OMEGA * KAPPA
                    )
                    index = cell_index[(profile_name, severity_sd, noise_sd)]
                    for name in metric_names:
                        values[name][index, replicate] = metrics[name]
    rows = []
    for index, (profile_name, severity_sd, noise_sd) in enumerate(cell_names):
        profile = ABILITY_PROFILES[profile_name]
        row = {
            "ability_profile": profile_name,
            "ability_min": min(profile),
            "ability_mean": float(np.mean(profile)),
            "ability_max": max(profile),
            "severity_sd": severity_sd,
            "specific_noise_sd": noise_sd,
            "omega": ALLOCATION_OMEGA,
        }
        for name in metric_names:
            row.update(metric_columns(name, values[name][index]))
        row["human_coverage_mean"] = 1.0
        row["avg_reviews_per_manuscript_mean"] = float(REVIEWERS)
        rows.append(row)
    return rows


def run_mechanism(replicates, seed_sequence):
    rng = rng_for(seed_sequence)
    metric_names = ("mae", "decision_error", "wsd", "decision_consensus", "group_bias")
    values = {
        name: np.empty((len(MECHANISM_SIGMAS), len(OMEGAS), replicates), dtype=float)
        for name in metric_names
    }
    base_abilities = (0.20, 0.50, 0.80)
    for replicate in range(replicates):
        quality = rng.standard_normal(M)
        groups = balanced_groups(rng, M)
        abilities = rng.permutation(np.asarray(base_abilities, dtype=float))
        severity = rng.normal(0.0, SEVERITY_SD, REVIEWERS)
        noise = rng.normal(0.0, SPECIFIC_NOISE_SD, (M, REVIEWERS))
        base_human = quality[:, None] * abilities[None, :] + severity[None, :] + noise
        human = base_human + DELTA * groups[:, None]
        theta, accepted = threshold_and_truth(quality)
        error_unit = rng.standard_normal(M)
        for sigma_index, sigma_eta in enumerate(MECHANISM_SIGMAS):
            ai = ALPHA * quality + KAPPA * groups + sigma_eta * error_unit
            for omega_index, omega in enumerate(OMEGAS):
                assisted = (1.0 - omega) * human + omega * ai[:, None]
                final_score = np.mean(assisted, axis=1)
                metrics = review_metrics(assisted, final_score, accepted, groups, quality, theta)
                metrics["group_bias"] = 2.0 * ((1.0 - omega) * DELTA + omega * KAPPA)
                for name in metric_names:
                    values[name][sigma_index, omega_index, replicate] = metrics[name]
    rows = []
    for sigma_index, sigma_eta in enumerate(MECHANISM_SIGMAS):
        sigma_rows = []
        for omega_index, omega in enumerate(OMEGAS):
            row = {"sigma_ai_common_error": sigma_eta, "omega": float(omega)}
            for name in metric_names:
                row.update(metric_columns(name, values[name][sigma_index, omega_index]))
            sigma_rows.append(row)
        zero = sigma_rows[0]
        full = sigma_rows[-1]
        min_mae_omega = float(min(sigma_rows, key=lambda item: item["mae_mean"])["omega"])
        wsd_drop = float(zero["wsd_mean"] - full["wsd_mean"])
        mae_rise_to_full = float(full["mae_mean"] - zero["mae_mean"])
        pseudo_consensus = wsd_drop > 0.0 and full["mae_mean"] > zero["mae_mean"]
        for row in sigma_rows:
            row["wsd_drop_omega0_to_1_mean"] = wsd_drop
            row["mae_change_omega0_to_1_mean"] = mae_rise_to_full
            row["min_mae_omega"] = min_mae_omega
            row["pseudo_consensus_pattern"] = int(pseudo_consensus)
        rows.extend(sigma_rows)
    return rows


def evenly_spaced_ids(count, target):
    ids = np.floor((np.arange(target, dtype=float) + 0.5) * count / target).astype(int)
    return np.unique(ids)


def strategy_counts(strategy, distance, extra_count, total_manuscripts):
    counts = np.ones(total_manuscripts, dtype=np.int16)
    if extra_count == 0:
        return counts
    if strategy == "uniform_extra":
        counts[evenly_spaced_ids(total_manuscripts, extra_count)] = 2
    elif strategy == "boundary_min_one":
        counts[np.argsort(distance, kind="mergesort")[:extra_count]] = 2
    elif strategy == "boundary_reallocating":
        ordered = np.argsort(distance, kind="mergesort")
        zero_count = min(extra_count, (total_manuscripts - extra_count) // 2)
        counts[ordered[: extra_count + zero_count]] = 2
        if zero_count:
            counts[ordered[-zero_count:]] = 0
    else:
        raise ValueError(strategy)
    if int(counts.sum()) != total_manuscripts + extra_count:
        raise RuntimeError(
            f"allocation budget mismatch: strategy={strategy}, extra={extra_count}, "
            f"actual={int(counts.sum())}, target={total_manuscripts + extra_count}"
        )
    return counts


def evaluate_allocation(counts, assisted, ai, base_human, quality, groups, theta, accepted):
    total = int(counts.sum())
    one = counts == 1
    two = counts == 2
    zero = counts == 0
    final_score = ai.copy()
    final_score[one] = assisted[one, 0]
    final_score[two] = np.mean(assisted[two, :2], axis=1)
    wsd = np.mean(np.std(assisted[two, :2], axis=1, ddof=1)) if np.any(two) else np.nan
    consensus_all = np.mean(
        np.concatenate(
            (
                np.ones(int(one.sum()), dtype=bool),
                (assisted[two, 0] >= theta) == (assisted[two, 1] >= theta),
            )
        )
    )
    consensus_multi = np.mean(
        (assisted[two, 0] >= theta) == (assisted[two, 1] >= theta)
    ) if np.any(two) else np.nan
    plus = groups == 1
    minus = ~plus
    ai_error = ai - KAPPA * groups
    human_plus = (1.0 - ALLOCATION_OMEGA) * (base_human + DELTA)
    human_plus = human_plus + ALLOCATION_OMEGA * (ai_error + KAPPA)[:, None]
    human_minus = (1.0 - ALLOCATION_OMEGA) * (base_human - DELTA)
    human_minus = human_minus + ALLOCATION_OMEGA * (ai_error - KAPPA)[:, None]
    final_plus = ai_error + KAPPA
    final_minus = ai_error - KAPPA
    final_plus[one] = human_plus[one, 0]
    final_plus[two] = np.mean(human_plus[two, :2], axis=1)
    final_minus[one] = human_minus[one, 0]
    final_minus[two] = np.mean(human_minus[two, :2], axis=1)
    group_bias = np.mean(final_plus - final_minus)
    avg_plus = float(np.mean(counts[plus]))
    avg_minus = float(np.mean(counts[minus]))
    return {
        "mae": np.mean(np.abs(final_score - quality)),
        "decision_error": np.mean((final_score >= theta) != accepted),
        "multi_review_wsd": wsd,
        "decision_consensus_reviewed": consensus_all,
        "decision_consensus_multi_review": consensus_multi,
        "group_bias": group_bias,
        "group_mae_plus": np.mean(np.abs(final_score[plus] - quality[plus])),
        "group_mae_minus": np.mean(np.abs(final_score[minus] - quality[minus])),
        "human_coverage": np.mean(~zero),
        "multi_review_coverage": np.mean(two),
        "ai_only_coverage": np.mean(zero),
        "avg_reviews": total / M,
        "total_reviews": total,
        "group_avg_reviews_plus": avg_plus,
        "group_avg_reviews_minus": avg_minus,
        "group_resource_gap": abs(avg_plus - avg_minus) / (avg_plus + avg_minus),
    }


def run_allocation(replicates, seed_sequence):
    rng = rng_for(seed_sequence)
    records = []
    base_abilities = (0.20, 0.50, 0.80)
    for extra_fraction in EXTRA_FRACTIONS:
        extra_count = int(round(M * extra_fraction))
        for replicate in range(replicates):
            quality = rng.standard_normal(M)
            groups = balanced_groups(rng, M)
            abilities = rng.permutation(np.asarray(base_abilities, dtype=float))
            severity = rng.normal(0.0, SEVERITY_SD, REVIEWERS)
            noise = rng.normal(0.0, SPECIFIC_NOISE_SD, (M, REVIEWERS))
            base_human = quality[:, None] * abilities[None, :] + severity[None, :] + noise
            error_unit = rng.standard_normal(M)
            ai = ALPHA * quality + KAPPA * groups + SIGMA_AI * error_unit
            human = base_human + DELTA * groups[:, None]
            assisted = (1.0 - ALLOCATION_OMEGA) * human + ALLOCATION_OMEGA * ai[:, None]
            theta, accepted = threshold_and_truth(quality)
            distance = np.abs(ai - theta)
            uniform_counts = strategy_counts("uniform_extra", distance, extra_count, M)
            uniform_metrics = evaluate_allocation(
                uniform_counts, assisted, ai, base_human, quality, groups, theta, accepted
            )
            for strategy in STRATEGIES:
                counts = strategy_counts(strategy, distance, extra_count, M)
                metrics = evaluate_allocation(
                    counts, assisted, ai, base_human, quality, groups, theta, accepted
                )
                metrics["coverage_delta_vs_uniform"] = metrics["human_coverage"] - uniform_metrics["human_coverage"]
                metrics["total_review_delta_vs_uniform"] = metrics["total_reviews"] - uniform_metrics["total_reviews"]
                metrics["avg_review_delta_vs_uniform"] = metrics["avg_reviews"] - uniform_metrics["avg_reviews"]
                metrics["group_resource_gap_delta_vs_uniform"] = (
                    metrics["group_resource_gap"] - uniform_metrics["group_resource_gap"]
                )
                metrics["mae_delta_vs_uniform"] = metrics["mae"] - uniform_metrics["mae"]
                metrics["decision_error_delta_vs_uniform"] = (
                    metrics["decision_error"] - uniform_metrics["decision_error"]
                )
                records.append(
                    {
                        "extra_fraction": extra_fraction,
                        "strategy": strategy,
                        "omega": ALLOCATION_OMEGA,
                        "target_total_reviews": M + extra_count,
                        **metrics,
                    }
                )
    grouped = {}
    for record in records:
        key = (record["extra_fraction"], record["strategy"])
        grouped.setdefault(key, []).append(record)
    rows = []
    metric_names = tuple(key for key in records[0] if key not in {"extra_fraction", "strategy", "omega", "target_total_reviews"})
    for (extra_fraction, strategy), items in grouped.items():
        row = {
            "extra_fraction": extra_fraction,
            "strategy": strategy,
            "omega": ALLOCATION_OMEGA,
            "target_total_reviews": M + int(round(M * extra_fraction)),
        }
        for name in metric_names:
            row.update(metric_columns(name, [item[name] for item in items]))
        rows.append(row)
    return rows


def parameter_rows(seed):
    rows = [
        {"parameter": "seed", "value": seed, "description": "fixed_default_seed"},
        {"parameter": "manuscripts_M", "value": M, "description": "manuscripts_per_replicate"},
        {"parameter": "monte_carlo_replicates", "value": REPLICATES, "description": "paired_replicates_per_cell"},
        {"parameter": "acceptance_rate", "value": ACCEPTANCE_RATE, "description": "exact_top_fraction_by_quality"},
        {"parameter": "threshold_rule", "value": "top_20_percent_quality", "description": "synthetic_oracle_threshold"},
        {"parameter": "reviewers", "value": REVIEWERS, "description": "main_sensitivity_mechanism"},
        {"parameter": "base_abilities", "value": "0.20|0.50|0.80", "description": "permuted_within_replicate"},
        {"parameter": "severity_sd", "value": SEVERITY_SD, "description": "main_reviewer_severity"},
        {"parameter": "specific_noise_sd", "value": SPECIFIC_NOISE_SD, "description": "main_reviewer_manuscript_noise"},
        {"parameter": "human_mu", "value": MU, "description": "global_score_offset"},
        {"parameter": "human_group_delta", "value": DELTA, "description": "synthetic_group_effect"},
        {"parameter": "ai_quality_alpha", "value": ALPHA, "description": "ai_quality_loading"},
        {"parameter": "ai_group_kappa", "value": KAPPA, "description": "synthetic_ai_group_effect"},
        {"parameter": "ai_common_error_sd", "value": SIGMA_AI, "description": "main_shared_ai_error"},
        {"parameter": "omegas", "value": "|".join(str(value) for value in OMEGAS), "description": "human_ai_blend_grid"},
        {"parameter": "allocation_omega", "value": ALLOCATION_OMEGA, "description": "allocation_sensitivity_blend"},
        {"parameter": "allocation_extra_fractions", "value": "|".join(str(value) for value in EXTRA_FRACTIONS), "description": "same_budget_extra_review_grid"},
        {"parameter": "allocation_strategies", "value": "|".join(STRATEGIES), "description": "paired_same_budget_strategies"},
        {"parameter": "mechanism_common_error_sds", "value": "|".join(str(value) for value in MECHANISM_SIGMAS), "description": "shared_ai_error_mechanism"},
        {"parameter": "ability_profiles", "value": "low|middle|high", "description": "sensitivity_profiles"},
        {"parameter": "severity_sensitivity_sds", "value": "|".join(str(value) for value in SEVERITY_SDS), "description": "sensitivity_grid"},
        {"parameter": "specific_noise_sensitivity_sds", "value": "|".join(str(value) for value in SPECIFIC_NOISE_SDS), "description": "sensitivity_grid"},
        {"parameter": "group_labels", "value": "-1|+1", "description": "abstract_random_groups_independent_of_quality"},
        {"parameter": "gittins_index_computed", "value": False, "description": "intentionally_not_computed"},
    ]
    return rows


def write_csv(path, rows):
    if not rows:
        raise RuntimeError(f"no rows for {path}")
    fieldnames = list(rows[0])
    with path.open("w", encoding="utf-8", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=fieldnames, lineterminator="\n")
        writer.writeheader()
        writer.writerows(rows)


def finite_check(rows):
    for row in rows:
        for value in row.values():
            if isinstance(value, (bool, np.bool_)):
                continue
            if isinstance(value, (int, float, np.integer, np.floating)):
                if not np.isfinite(float(value)):
                    raise RuntimeError(f"nonfinite_output: {key}={value}")


def plot_omega(rows, path):
    figure, axes = plt.subplots(2, 2, figsize=(11.5, 8.0), constrained_layout=True)
    panels = (
        (axes[0, 0], "mae", "MAE", False),
        (axes[0, 1], "wsd", "Within-manuscript WSD", False),
        (axes[1, 0], "decision_consensus", "Decision consensus", True),
        (axes[1, 1], "group_bias", "Counterfactual group bias", False),
    )
    for axis, metric, label, percent in panels:
        x = np.asarray([row["omega"] for row in rows])
        mean = np.asarray([row[f"{metric}_mean"] for row in rows])
        low = np.asarray([row[f"{metric}_q025"] for row in rows])
        high = np.asarray([row[f"{metric}_q975"] for row in rows])
        if percent:
            mean = 100.0 * mean
            low = 100.0 * low
            high = 100.0 * high
        errors = np.vstack((np.maximum(mean - low, 0.0), np.maximum(high - mean, 0.0)))
        axis.errorbar(x, mean, yerr=errors, marker="o", capsize=3)
        axis.set_xlabel("omega")
        ylabel = f"{label} (%)" if percent else label
        axis.set_ylabel(ylabel)
        axis.set_xticks(OMEGAS)
        axis.grid(alpha=0.25)
    figure.suptitle("Synthetic peer-review mechanism simulation")
    figure.savefig(path, dpi=180, metadata={"Software": "peer_review.py"})
    plt.close(figure)


def plot_mechanism(rows, path):
    figure, axes = plt.subplots(1, 2, figsize=(11.0, 4.3), constrained_layout=True)
    colors = {0.0: "#3366cc", 0.4: "#777777", 1.2: "#cc3333"}
    for sigma in (0.0, 0.4, 1.2):
        selected = sorted(
            (row for row in rows if row["sigma_ai_common_error"] == sigma),
            key=lambda row: row["omega"],
        )
        x = [row["omega"] for row in selected]
        mae = [row["mae_mean"] for row in selected]
        wsd = [row["wsd_mean"] for row in selected]
        axes[0].plot(x, mae, marker="o", color=colors[sigma], label=f"shared error SD={sigma:g}")
        axes[1].plot(x, wsd, marker="o", color=colors[sigma], label=f"shared error SD={sigma:g}")
    axes[0].set_ylabel("MAE")
    axes[1].set_ylabel("Within-manuscript WSD")
    for axis in axes:
        axis.set_xlabel("omega")
        axis.set_xticks(OMEGAS)
        axis.grid(alpha=0.25)
        axis.legend(frameon=False)
    axes[0].set_title("Accuracy can worsen with AI reliance")
    axes[1].set_title("Shared AI signal compresses dispersion")
    figure.savefig(path, dpi=180, metadata={"Software": "peer_review.py"})
    plt.close(figure)


def plot_allocation(rows, path):
    figure, axes = plt.subplots(1, 2, figsize=(11.0, 4.3), constrained_layout=True)
    markers = {"uniform_extra": "o", "boundary_min_one": "s", "boundary_reallocating": "^"}
    labels = {
        "uniform_extra": "Uniform extras",
        "boundary_min_one": "Boundary, min 1",
        "boundary_reallocating": "Boundary, reallocating",
    }
    for extra_fraction in EXTRA_FRACTIONS:
        for strategy in STRATEGIES:
            row = next(
                item
                for item in rows
                if item["extra_fraction"] == extra_fraction and item["strategy"] == strategy
            )
            label = f"{labels[strategy]} (+{extra_fraction:.0%})"
            axes[0].scatter(
                row["human_coverage_mean"],
                row["mae_mean"],
                marker=markers[strategy],
                s=65,
                label=label,
            )
            axes[1].scatter(
                row["human_coverage_mean"],
                row["decision_error_mean"],
                marker=markers[strategy],
                s=65,
                label=label,
            )
    axes[0].set_title("Same-review budget: MAE")
    axes[1].set_title("Same-review budget: decision error")
    for axis in axes:
        axis.set_xlabel("Human-review coverage")
        axis.set_xlim(0.40, 1.03)
        axis.grid(alpha=0.25)
        axis.legend(frameon=False, fontsize=8)
    axes[0].set_ylabel("MAE")
    axes[1].set_ylabel("20% threshold decision error")
    figure.savefig(path, dpi=180, metadata={"Software": "peer_review.py"})
    plt.close(figure)


def validate_outputs(csv_paths, png_paths):
    validation = {
        "finite_numeric_csv": True,
        "csv_rows": {},
        "png_shapes": {},
        "sha256": {},
    }
    for path in csv_paths:
        with path.open("r", encoding="utf-8", newline="") as handle:
            rows = list(csv.DictReader(handle))
        if not rows:
            raise RuntimeError(f"empty CSV: {path}")
        for row in rows:
            for value in row.values():
                if value is None or value == "":
                    continue
                try:
                    numeric = float(value)
                except ValueError:
                    continue
                if not np.isfinite(numeric):
                    raise RuntimeError(f"nonfinite CSV value: {path}")
        validation["csv_rows"][path.name] = len(rows)
    for path in png_paths:
        if not path.is_file() or path.stat().st_size == 0:
            raise RuntimeError(f"missing PNG: {path}")
        image = mpimg.imread(path)
        if image.ndim not in (2, 3) or image.shape[0] < 100 or image.shape[1] < 100:
            raise RuntimeError(f"invalid PNG: {path}")
        if not np.all(np.isfinite(image)):
            raise RuntimeError(f"nonfinite PNG: {path}")
        validation["png_shapes"][path.name] = list(image.shape)
    for path in csv_paths + png_paths:
        validation["sha256"][path.name] = hashlib.sha256(path.read_bytes()).hexdigest()
    return validation


def print_summary(seed, omega_rows, mechanism_rows, allocation_rows, validation):
    omega_zero = omega_rows[0]
    omega_full = omega_rows[-1]
    high_error_zero = next(
        row
        for row in mechanism_rows
        if row["sigma_ai_common_error"] == 1.2 and row["omega"] == 0.0
    )
    high_error_full = next(
        row
        for row in mechanism_rows
        if row["sigma_ai_common_error"] == 1.2 and row["omega"] == 1.0
    )
    print(f"status=pass")
    print(f"seed={seed}")
    print(f"M={M}")
    print(f"replicates={REPLICATES}")
    print(f"omega_zero_mae={omega_zero['mae_mean']:.6f}")
    print(f"omega_full_mae={omega_full['mae_mean']:.6f}")
    print(f"omega_zero_wsd={omega_zero['wsd_mean']:.6f}")
    print(f"omega_full_wsd={omega_full['wsd_mean']:.6f}")
    print(f"omega_zero_consensus={omega_zero['decision_consensus_mean']:.6f}")
    print(f"omega_full_consensus={omega_full['decision_consensus_mean']:.6f}")
    print(f"omega_zero_decision_error={omega_zero['decision_error_mean']:.6f}")
    print(f"omega_full_decision_error={omega_full['decision_error_mean']:.6f}")
    print(f"omega_full_group_bias={omega_full['group_bias_mean']:.6f}")
    print(
        "high_shared_error_wsd_drop="
        f"{high_error_zero['wsd_mean'] - high_error_full['wsd_mean']:.6f}"
    )
    print(
        "high_shared_error_mae_change="
        f"{high_error_full['mae_mean'] - high_error_zero['mae_mean']:.6f}"
    )
    for fraction in EXTRA_FRACTIONS:
        for strategy in STRATEGIES:
            row = next(
                item
                for item in allocation_rows
                if item["extra_fraction"] == fraction and item["strategy"] == strategy
            )
            print(
                f"allocation_{fraction:.2f}_{strategy}="
                f"coverage:{row['human_coverage_mean']:.6f},"
                f"reviews:{row['avg_reviews_mean']:.6f},"
                f"review_delta:{row['avg_review_delta_vs_uniform_mean']:.6f},"
                f"mae:{row['mae_mean']:.6f},"
                f"de:{row['decision_error_mean']:.6f},"
                f"group_resource_gap:{row['group_resource_gap_mean']:.6f}"
            )
    print(f"csv_rows={validation['csv_rows']}")
    print(f"png_shapes={validation['png_shapes']}")
    print(f"sha256={validation['sha256']}")


def main():
    args = parse_args()
    master = np.random.SeedSequence(args.seed)
    omega_seed, sensitivity_seed, mechanism_seed, allocation_seed = master.spawn(4)
    parameter_data = parameter_rows(args.seed)
    omega_rows = run_omega(REPLICATES, omega_seed)
    sensitivity_rows = run_parameter_sensitivity(REPLICATES, sensitivity_seed)
    mechanism_rows = run_mechanism(REPLICATES, mechanism_seed)
    allocation_rows = run_allocation(REPLICATES, allocation_seed)
    for rows in (parameter_data, omega_rows, sensitivity_rows, mechanism_rows, allocation_rows):
        finite_check(rows)
    parameter_path = RESULTS_DIR / "peer_review_parameters.csv"
    omega_path = RESULTS_DIR / "peer_review_omega_summary.csv"
    sensitivity_path = RESULTS_DIR / "peer_review_parameter_sensitivity.csv"
    mechanism_path = RESULTS_DIR / "peer_review_common_error_mechanism.csv"
    allocation_path = RESULTS_DIR / "peer_review_allocation_sensitivity.csv"
    write_csv(parameter_path, parameter_data)
    write_csv(omega_path, omega_rows)
    write_csv(sensitivity_path, sensitivity_rows)
    write_csv(mechanism_path, mechanism_rows)
    write_csv(allocation_path, allocation_rows)
    plot_omega(omega_rows, FIGURES_DIR / "peer_review_summary.png")
    plot_mechanism(mechanism_rows, FIGURES_DIR / "peer_review_common_error_mechanism.png")
    plot_allocation(allocation_rows, FIGURES_DIR / "peer_review_allocation_sensitivity.png")
    csv_paths = [parameter_path, omega_path, sensitivity_path, mechanism_path, allocation_path]
    png_paths = [
        FIGURES_DIR / "peer_review_summary.png",
        FIGURES_DIR / "peer_review_common_error_mechanism.png",
        FIGURES_DIR / "peer_review_allocation_sensitivity.png",
    ]
    validation = validate_outputs(csv_paths, png_paths)
    print_summary(args.seed, omega_rows, mechanism_rows, allocation_rows, validation)


if __name__ == "__main__":
    main()