from __future__ import annotations

from itertools import combinations
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from scipy.stats import ttest_rel


BASE_DIR = Path(__file__).resolve().parent
INPUT_CSV = BASE_DIR / "listening_scores_dataframe_augmented.csv"

EXPERIMENT_GROUPS = {
    "exp2b_loudness_compensation": {
        "title": "Experiment 2b: Compensation of Loudness Errors",
        "questions": ["Q01", "Q02", "Q03"],
        "models": [
            "no balance MEGAMI",
            "with balance MEGAMI",
            "no balance Diff-MST",
            "with balance Diff-MST",
        ],
        "model_labels": {
            "no balance MEGAMI": "NB-MEGAMI",
            "with balance MEGAMI": "WB-MEGAMI",
            "no balance Diff-MST": "NB-DiffMST",
            "with balance Diff-MST": "WB-DiffMST",
        },
        "short_labels": {
            "no balance MEGAMI": "NB-M",
            "with balance MEGAMI": "WB-M",
            "no balance Diff-MST": "NB-D",
            "with balance Diff-MST": "WB-D",
        },
        "palette": ["#2C7FB8", "#F28E2B", "#1B9E77", "#D95F02"],
    },
    "exp2a_grouping_compensation": {
        "title": "Experiment 2a: Compensation of Grouping Errors",
        "questions": ["Q04", "Q05", "Q06"],
        "models": [
            "instrumently grouping MEGAMI",
            "4group MEGAMI",
            "7group MEGAMI",
            "instrumently grouping Diff-MST",
            "4group Diff-MST",
            "7group Diff-MST",
        ],
        "model_labels": {
            "instrumently grouping MEGAMI": "Instr-MEGAMI",
            "4group MEGAMI": "4G-MEGAMI",
            "7group MEGAMI": "7G-MEGAMI",
            "instrumently grouping Diff-MST": "Instr-DiffMST",
            "4group Diff-MST": "4G-DiffMST",
            "7group Diff-MST": "7G-DiffMST",
        },
        "short_labels": {
            "instrumently grouping MEGAMI": "I-M",
            "4group MEGAMI": "4G-M",
            "7group MEGAMI": "7G-M",
            "instrumently grouping Diff-MST": "I-D",
            "4group Diff-MST": "4G-D",
            "7group Diff-MST": "7G-D",
        },
        "palette": ["#1B9E77", "#D95F02", "#7570B3", "#E7298A", "#66A61E", "#A6761D"],
    },
    "exp1_intra_group_quality": {
        "title": "Experiment 1: Intra-group Mixing Quality",
        "questions": ["Q07", "Q08"],
        "models": ["Diff-MST", "random mix", "MEGAMI", "balanced"],
        "model_labels": {
            "Diff-MST": "DiffMST",
            "random mix": "NoMix",
            "MEGAMI": "MEGAMI",
            "balanced": "BalancedMix",
        },
        "short_labels": {
            "Diff-MST": "D",
            "random mix": "NoMix",
            "MEGAMI": "M",
            "balanced": "ELL",
        },
        "palette": ["#4E79A7", "#59A14F", "#F28E2B", "#E15759"],
    },
    "exp3_full_mix_ablation": {
        "title": "Experiment 3: Full-Mix Ablation",
        "questions": ["Q10", "Q11", "Q12", "Q13"],
        "models": [
            "human",
            "Two stage with MEGAMI",
            "MEGAMI",
            "Two stage with Diff-MST",
            "Diff-MST",
            "random mix",
        ],
        "model_labels": {
            "human": "Human",
            "Two stage with MEGAMI": "2Stage-MEGAMI",
            "MEGAMI": "MEGAMI",
            "Two stage with Diff-MST": "2Stage-DiffMST",
            "Diff-MST": "DiffMST",
            "random mix": "NoMix",
        },
        "short_labels": {
            "human": "Human",
            "Two stage with MEGAMI": "2S-M",
            "MEGAMI": "M",
            "Two stage with Diff-MST": "2S-D",
            "Diff-MST": "D",
            "random mix": "NoMix",
        },
        "palette": ["#9C755F", "#1B9E77", "#4E79A7", "#B07AA1", "#E15759", "#F28E2B"],
    },
}


def benjamini_hochberg(p_values: pd.Series) -> pd.Series:
    values = p_values.astype(float).to_numpy()
    n = len(values)
    order = np.argsort(values)
    ranked = values[order]
    adjusted = np.empty(n, dtype=float)
    prev = 1.0
    for i in range(n - 1, -1, -1):
        rank = i + 1
        bh = ranked[i] * n / rank
        prev = min(prev, bh)
        adjusted[i] = prev
    result = np.empty(n, dtype=float)
    result[order] = np.clip(adjusted, 0, 1)
    return pd.Series(result, index=p_values.index)


def apply_style() -> None:
    sns.set_theme(style="whitegrid", context="paper")
    plt.rcParams.update(
        {
            "font.family": "serif",
            "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
            "font.size": 9,
            "axes.titlesize": 10,
            "axes.labelsize": 9,
            "xtick.labelsize": 8,
            "ytick.labelsize": 8,
            "legend.fontsize": 8,
            "figure.dpi": 300,
            "savefig.dpi": 300,
            "pdf.fonttype": 42,
            "ps.fonttype": 42,
        }
    )


def apply_publication_style(kind: str) -> None:
    sns.set_theme(style="whitegrid", context="paper")
    if kind == "single":
        plt.rcParams.update(
            {
                "font.family": "serif",
                "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
                "font.size": 10,
                "axes.titlesize": 10,
                "axes.labelsize": 9.5,
                "xtick.labelsize": 8.8,
                "ytick.labelsize": 8.8,
                "legend.fontsize": 8.8,
                "figure.dpi": 300,
                "savefig.dpi": 300,
                "pdf.fonttype": 42,
                "ps.fonttype": 42,
            }
        )
    elif kind == "double":
        plt.rcParams.update(
            {
                "font.family": "serif",
                "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
                "font.size": 9.4,
                "axes.titlesize": 10,
                "axes.labelsize": 9.2,
                "xtick.labelsize": 8.5,
                "ytick.labelsize": 8.5,
                "legend.fontsize": 8.5,
                "figure.dpi": 300,
                "savefig.dpi": 300,
                "pdf.fonttype": 42,
                "ps.fonttype": 42,
            }
        )


def p_to_marker(p_value: float) -> str:
    if pd.isna(p_value):
        return ""
    if p_value < 0.001:
        return "***"
    if p_value < 0.01:
        return "**"
    if p_value < 0.05:
        return "*"
    return ""


def add_sig_bar(ax, x1: float, x2: float, y: float, h: float, text: str) -> None:
    if not text:
        return
    ax.plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.0, c="black", clip_on=False)
    ax.text((x1 + x2) / 2, y + h + 1.0, text, ha="center", va="bottom", fontsize=8)


def load_data() -> pd.DataFrame:
    df = pd.read_csv(INPUT_CSV)
    df["score"] = pd.to_numeric(df["score"], errors="coerce")
    return df


def summarize_group(df: pd.DataFrame, group_key: str, config: dict) -> pd.DataFrame:
    sub = df[df["question_no"].isin(config["questions"])].copy()
    sub["model_label"] = sub["model_name"].map(config["model_labels"])
    sub = sub[sub["model_name"].isin(config["models"]) & sub["score"].notna()].copy()
    sub["model_name"] = pd.Categorical(sub["model_name"], categories=config["models"], ordered=True)
    summary = (
        sub.groupby(["question_no", "model_name", "model_label"], dropna=False)["score"]
        .agg(["count", "mean", "median", "std", "min", "max"])
        .reset_index()
        .sort_values(["question_no", "model_name"])
    )
    output = BASE_DIR / f"{group_key}_summary.csv"
    summary.to_csv(output, index=False, encoding="utf-8-sig")
    return sub


def draw_violin_plot(group_key: str, config: dict, data: pd.DataFrame) -> tuple[Path, Path]:
    questions = config["questions"]
    models = config["models"]
    labels = config["short_labels"]
    palette = dict(zip(models, config["palette"]))

    sub = data[data["question_no"].isin(questions)].copy()
    sub = sub[sub["model_name"].isin(models) & sub["score"].notna()].copy()
    sub["model_name"] = pd.Categorical(sub["model_name"], categories=models, ordered=True)

    fig, axes = plt.subplots(1, len(questions), figsize=(4.8 * len(questions), 4.2), sharey=True, constrained_layout=True)
    if len(questions) == 1:
        axes = [axes]

    for idx, question in enumerate(questions):
        ax = axes[idx]
        q_df = sub[sub["question_no"] == question].copy()
        sns.violinplot(
            data=q_df,
            x="model_name",
            y="score",
            order=models,
            hue="model_name",
            hue_order=models,
            palette=palette,
            legend=False,
            cut=0,
            density_norm="width",
            bw_adjust=0.7,
            inner="box",
            width=0.92,
            linewidth=0.8,
            saturation=0.95,
            ax=ax,
        )
        ax.set_title(question)
        ax.set_xlabel("")
        ax.set_ylabel("Score" if idx == 0 else "")
        ax.set_ylim(-5, 105)
        ax.set_xticks(range(len(models)))
        ax.set_xticklabels([labels[m] for m in models], rotation=0, ha="center")
        ax.spines["top"].set_visible(False)
        ax.spines["right"].set_visible(False)

    fig.suptitle(config["title"], fontsize=12)
    png_path = BASE_DIR / f"{group_key}_violin.png"
    pdf_path = BASE_DIR / f"{group_key}_violin.pdf"
    fig.savefig(png_path, bbox_inches="tight")
    fig.savefig(pdf_path, bbox_inches="tight")
    plt.close(fig)
    return png_path, pdf_path


def draw_box_plot(group_key: str, config: dict, data: pd.DataFrame) -> tuple[Path, Path]:
    questions = config["questions"]
    models = config["models"]
    labels = config["short_labels"]
    palette = dict(zip(models, config["palette"]))

    sub = data[data["question_no"].isin(questions)].copy()
    sub = sub[sub["model_name"].isin(models) & sub["score"].notna()].copy()
    sub["model_name"] = pd.Categorical(sub["model_name"], categories=models, ordered=True)

    fig, axes = plt.subplots(1, len(questions), figsize=(4.8 * len(questions), 4.2), sharey=True, constrained_layout=True)
    if len(questions) == 1:
        axes = [axes]

    for idx, question in enumerate(questions):
        ax = axes[idx]
        q_df = sub[sub["question_no"] == question].copy()
        sns.boxplot(
            data=q_df,
            x="model_name",
            y="score",
            order=models,
            hue="model_name",
            hue_order=models,
            palette=palette,
            legend=False,
            linewidth=0.8,
            fliersize=2,
            width=0.55,
            medianprops={"color": "black", "linewidth": 2.2},
            ax=ax,
        )
        sns.stripplot(
            data=q_df,
            x="model_name",
            y="score",
            order=models,
            color="black",
            alpha=0.25,
            size=1.5,
            jitter=0.12,
            ax=ax,
        )
        ax.set_title(question)
        ax.set_xlabel("")
        ax.set_ylabel("Score" if idx == 0 else "")
        ax.set_ylim(-5, 105)
        ax.set_xticks(range(len(models)))
        ax.set_xticklabels([labels[m] for m in models], rotation=0, ha="center")
        ax.spines["top"].set_visible(False)
        ax.spines["right"].set_visible(False)

    fig.suptitle(config["title"], fontsize=12)
    png_path = BASE_DIR / f"{group_key}_box.png"
    pdf_path = BASE_DIR / f"{group_key}_box.pdf"
    fig.savefig(png_path, bbox_inches="tight")
    fig.savefig(pdf_path, bbox_inches="tight")
    plt.close(fig)
    return png_path, pdf_path


def paired_test_from_long_df(sub: pd.DataFrame, model_a: str, model_b: str) -> tuple[float, int]:
    pivot = sub.pivot_table(
        index=sub["participant_id"].astype(str) + " | " + sub["question_no"].astype(str),
        columns="model_name",
        values="score",
        aggfunc="first",
    )
    paired = pivot[[model_a, model_b]].dropna()
    if len(paired) < 2:
        return np.nan, len(paired)
    _, p_value = ttest_rel(paired[model_a], paired[model_b], nan_policy="omit")
    return float(p_value), len(paired)


def draw_submission_figure2(df: pd.DataFrame) -> tuple[Path, Path]:
    apply_publication_style("single")
    config = EXPERIMENT_GROUPS["exp1_intra_group_quality"]
    order = ["balanced", "MEGAMI", "Diff-MST", "random mix"]
    tick_labels = ["ELL", "M", "D", "NoMix"]
    palette = {
        "balanced": "#E15759",
        "MEGAMI": "#F28E2B",
        "Diff-MST": "#4E79A7",
        "random mix": "#59A14F",
    }
    sub = df[df["question_no"].isin(config["questions"])].copy()
    sub = sub[sub["model_name"].isin(order) & sub["score"].notna()].copy()
    sub["model_name"] = pd.Categorical(sub["model_name"], categories=order, ordered=True)

    fig, ax = plt.subplots(figsize=(3.35, 2.75), constrained_layout=True)
    sns.boxplot(
        data=sub,
        x="model_name",
        y="score",
        order=order,
        hue="model_name",
        hue_order=order,
        palette=palette,
        legend=False,
        linewidth=0.9,
        fliersize=2,
        width=0.58,
        medianprops={"color": "black", "linewidth": 2.6},
        ax=ax,
    )
    sns.stripplot(
        data=sub,
        x="model_name",
        y="score",
        order=order,
        color="black",
        alpha=0.3,
        size=1.9,
        jitter=0.13,
        ax=ax,
    )
    ax.set_title("Intra-group Mixing Quality", pad=8)
    ax.set_xlabel("")
    ax.set_ylabel("MOS")
    ax.set_ylim(-5, 117)
    ax.set_xticks(range(len(order)))
    ax.set_xticklabels(tick_labels, rotation=0, ha="center")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    p_ell_d, _ = paired_test_from_long_df(sub, "balanced", "Diff-MST")
    p_m_d, _ = paired_test_from_long_df(sub, "MEGAMI", "Diff-MST")
    add_sig_bar(ax, 0, 2, 103.5, 1.8, p_to_marker(p_ell_d))
    add_sig_bar(ax, 1, 2, 109.0, 1.8, p_to_marker(p_m_d))

    png_path = BASE_DIR / "figure2_rq1_intra_group_quality_box.png"
    pdf_path = BASE_DIR / "figure2_rq1_intra_group_quality_box.pdf"
    fig.savefig(png_path, bbox_inches="tight")
    fig.savefig(pdf_path, bbox_inches="tight")
    plt.close(fig)
    return png_path, pdf_path


def draw_submission_figure3(df: pd.DataFrame) -> tuple[Path, Path]:
    apply_publication_style("double")
    fig, axes = plt.subplots(1, 2, figsize=(6.95, 3.05), sharey=True, constrained_layout=True)

    # Panel (a): grouping errors
    grouping_questions = EXPERIMENT_GROUPS["exp2a_grouping_compensation"]["questions"]
    grouping_sub = df[df["question_no"].isin(grouping_questions)].copy()
    grouping_map = {
        "instrumently grouping MEGAMI": ("Instrument", "MEGAMI"),
        "4group MEGAMI": ("4G", "MEGAMI"),
        "7group MEGAMI": ("7G", "MEGAMI"),
        "instrumently grouping Diff-MST": ("Instrument", "Diff-MST"),
        "4group Diff-MST": ("4G", "Diff-MST"),
        "7group Diff-MST": ("7G", "Diff-MST"),
    }
    grouping_sub = grouping_sub[grouping_sub["model_name"].isin(grouping_map) & grouping_sub["score"].notna()].copy()
    grouping_sub["grouping"] = grouping_sub["model_name"].map(lambda x: grouping_map[x][0])
    grouping_sub["family"] = grouping_sub["model_name"].map(lambda x: grouping_map[x][1])
    grouping_order = ["Instrument", "4G", "7G"]
    family_order = ["MEGAMI", "Diff-MST"]
    palette = {"MEGAMI": "#F28E2B", "Diff-MST": "#4E79A7"}

    ax = axes[0]
    sns.boxplot(
        data=grouping_sub,
        x="grouping",
        y="score",
        hue="family",
        order=grouping_order,
        hue_order=family_order,
        palette=palette,
        linewidth=0.9,
        fliersize=2,
        width=0.62,
        medianprops={"color": "black", "linewidth": 2.4},
        ax=ax,
    )
    sns.stripplot(
        data=grouping_sub,
        x="grouping",
        y="score",
        hue="family",
        order=grouping_order,
        hue_order=family_order,
        dodge=True,
        color="black",
        alpha=0.22,
        size=1.4,
        jitter=0.12,
        ax=ax,
    )
    handles, labels = ax.get_legend_handles_labels()
    ax.legend(
        handles[:2],
        labels[:2],
        title="",
        loc="upper left",
        bbox_to_anchor=(0.0, 1.08),
        frameon=False,
        borderaxespad=0.0,
    )
    ax.set_title("(a) Grouping error")
    ax.set_xlabel("")
    ax.set_ylabel("MOS")
    ax.set_ylim(-5, 120)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    p_4g_7g_m, _ = paired_test_from_long_df(grouping_sub.rename(columns={"model_name": "orig_model"}).assign(
        model_name=lambda x: x["orig_model"]
    ), "4group MEGAMI", "7group MEGAMI")
    p_i_4g_d, _ = paired_test_from_long_df(grouping_sub.rename(columns={"model_name": "orig_model"}).assign(
        model_name=lambda x: x["orig_model"]
    ), "instrumently grouping Diff-MST", "4group Diff-MST")
    add_sig_bar(ax, 0.83, 1.83, 105.0, 1.8, p_to_marker(p_4g_7g_m))
    add_sig_bar(ax, 0.18, 1.18, 111.0, 1.8, p_to_marker(p_i_4g_d))

    # Panel (b): loudness errors
    loud_questions = EXPERIMENT_GROUPS["exp2b_loudness_compensation"]["questions"]
    loud_sub = df[df["question_no"].isin(loud_questions)].copy()
    loud_map = {
        "with balance MEGAMI": ("With balance", "MEGAMI"),
        "no balance MEGAMI": ("No balance", "MEGAMI"),
        "with balance Diff-MST": ("With balance", "Diff-MST"),
        "no balance Diff-MST": ("No balance", "Diff-MST"),
    }
    loud_sub = loud_sub[loud_sub["model_name"].isin(loud_map) & loud_sub["score"].notna()].copy()
    loud_sub["balance"] = loud_sub["model_name"].map(lambda x: loud_map[x][0])
    loud_sub["family"] = loud_sub["model_name"].map(lambda x: loud_map[x][1])
    balance_order = ["No balance", "With balance"]

    ax = axes[1]
    sns.boxplot(
        data=loud_sub,
        x="balance",
        y="score",
        hue="family",
        order=balance_order,
        hue_order=family_order,
        palette=palette,
        linewidth=0.9,
        fliersize=2,
        width=0.62,
        medianprops={"color": "black", "linewidth": 2.4},
        ax=ax,
    )
    sns.stripplot(
        data=loud_sub,
        x="balance",
        y="score",
        hue="family",
        order=balance_order,
        hue_order=family_order,
        dodge=True,
        color="black",
        alpha=0.22,
        size=1.4,
        jitter=0.12,
        ax=ax,
    )
    handles, labels = ax.get_legend_handles_labels()
    if ax.legend_:
        ax.legend_.remove()
    ax.set_title("(b) Loudness error")
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_ylim(-5, 120)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    p_wb_nb_m, _ = paired_test_from_long_df(loud_sub.rename(columns={"model_name": "orig_model"}).assign(
        model_name=lambda x: x["orig_model"]
    ), "with balance MEGAMI", "no balance MEGAMI")
    p_wb_nb_d, _ = paired_test_from_long_df(loud_sub.rename(columns={"model_name": "orig_model"}).assign(
        model_name=lambda x: x["orig_model"]
    ), "with balance Diff-MST", "no balance Diff-MST")
    add_sig_bar(ax, -0.18, 0.82, 105.0, 1.8, p_to_marker(p_wb_nb_m))
    add_sig_bar(ax, 0.18, 1.18, 111.0, 1.8, p_to_marker(p_wb_nb_d))

    fig.suptitle("Compensation of Intra-group Errors", fontsize=11)
    png_path = BASE_DIR / "figure3_rq2_error_compensation_box.png"
    pdf_path = BASE_DIR / "figure3_rq2_error_compensation_box.pdf"
    fig.savefig(png_path, bbox_inches="tight")
    fig.savefig(pdf_path, bbox_inches="tight")
    plt.close(fig)
    return png_path, pdf_path


def draw_submission_figure4(df: pd.DataFrame) -> tuple[Path, Path]:
    apply_publication_style("single")
    config = EXPERIMENT_GROUPS["exp3_full_mix_ablation"]
    order = ["Diff-MST", "Two stage with Diff-MST", "MEGAMI", "Two stage with MEGAMI", "human"]
    tick_labels = ["D", "2S-D", "M", "2S-M", "Human"]
    palette = {
        "Diff-MST": "#4E79A7",
        "Two stage with Diff-MST": "#B07AA1",
        "MEGAMI": "#F28E2B",
        "Two stage with MEGAMI": "#1B9E77",
        "human": "#9C755F",
    }
    sub = df[df["question_no"].isin(config["questions"])].copy()
    sub = sub[sub["model_name"].isin(order) & sub["score"].notna()].copy()
    sub["model_name"] = pd.Categorical(sub["model_name"], categories=order, ordered=True)

    fig, ax = plt.subplots(figsize=(3.35, 2.8), constrained_layout=True)
    sns.boxplot(
        data=sub,
        x="model_name",
        y="score",
        order=order,
        hue="model_name",
        hue_order=order,
        palette=palette,
        legend=False,
        linewidth=0.9,
        fliersize=2,
        width=0.58,
        medianprops={"color": "black", "linewidth": 2.6},
        ax=ax,
    )
    sns.stripplot(
        data=sub,
        x="model_name",
        y="score",
        order=order,
        color="black",
        alpha=0.28,
        size=1.8,
        jitter=0.12,
        ax=ax,
    )
    ax.set_title("Full-Mix Quality Comparison", pad=8)
    ax.set_xlabel("")
    ax.set_ylabel("MOS")
    ax.set_ylim(-5, 117)
    ax.set_xticks(range(len(order)))
    ax.set_xticklabels(tick_labels, rotation=0, ha="center")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    p_d_2sd, _ = paired_test_from_long_df(sub, "Diff-MST", "Two stage with Diff-MST")
    p_m_2sm, _ = paired_test_from_long_df(sub, "MEGAMI", "Two stage with MEGAMI")
    add_sig_bar(ax, 0, 1, 103.5, 1.8, p_to_marker(p_d_2sd))
    add_sig_bar(ax, 2, 3, 109.0, 1.8, p_to_marker(p_m_2sm))

    png_path = BASE_DIR / "figure4_rq3_full_mix_box.png"
    pdf_path = BASE_DIR / "figure4_rq3_full_mix_box.pdf"
    fig.savefig(png_path, bbox_inches="tight")
    fig.savefig(pdf_path, bbox_inches="tight")
    plt.close(fig)
    return png_path, pdf_path


def paired_tests(group_key: str, config: dict, data: pd.DataFrame) -> Path:
    models = config["models"]
    labels = config["model_labels"]
    short_labels = config["short_labels"]
    sub = data[data["question_no"].isin(config["questions"])].copy()
    sub = sub[sub["model_name"].isin(models) & sub["score"].notna()].copy()
    sub["pair_id"] = sub["participant_id"].astype(str) + " | " + sub["question_no"].astype(str)
    pivot = sub.pivot_table(index="pair_id", columns="model_name", values="score", aggfunc="first")

    records = []
    for model_a, model_b in combinations(models, 2):
        paired = pivot[[model_a, model_b]].dropna()
        n_pairs = len(paired)
        if n_pairs < 2:
            continue
        diff = paired[model_a] - paired[model_b]
        t_stat, p_value = ttest_rel(paired[model_a], paired[model_b], nan_policy="omit")
        std_diff = diff.std(ddof=1)
        mean_diff = diff.mean()
        cohens_dz = mean_diff / std_diff if std_diff and not np.isclose(std_diff, 0) else np.nan
        records.append(
            {
                "experiment_group": group_key,
                "model_a": model_a,
                "model_b": model_b,
                "model_a_label": labels[model_a],
                "model_b_label": labels[model_b],
                "model_a_short": short_labels[model_a],
                "model_b_short": short_labels[model_b],
                "n_pairs": n_pairs,
                "mean_score_a": paired[model_a].mean(),
                "mean_score_b": paired[model_b].mean(),
                "mean_diff_a_minus_b": mean_diff,
                "t_statistic": t_stat,
                "p_value": p_value,
                "cohens_dz": cohens_dz,
                "winner_by_mean": labels[model_a] if mean_diff > 0 else labels[model_b] if mean_diff < 0 else "tie",
                "winner_by_mean_short": short_labels[model_a] if mean_diff > 0 else short_labels[model_b] if mean_diff < 0 else "tie",
            }
        )

    result = pd.DataFrame(records).sort_values("p_value")
    result["p_value_fdr_bh"] = benjamini_hochberg(result["p_value"])
    result["significant_p_lt_0_05"] = result["p_value"] < 0.05
    result["significant_fdr_bh_lt_0_05"] = result["p_value_fdr_bh"] < 0.05

    output = BASE_DIR / f"{group_key}_paired_ttests.csv"
    result.to_csv(output, index=False, encoding="utf-8-sig")
    return output


def build_paper_analysis(df: pd.DataFrame) -> Path:
    key_stats = {}
    for group_key, config in EXPERIMENT_GROUPS.items():
        stats_path = BASE_DIR / f"{group_key}_paired_ttests.csv"
        stats_df = pd.read_csv(stats_path)
        key_stats[group_key] = stats_df.sort_values(["p_value_fdr_bh", "p_value"]).head(8)

    lines = []
    lines.append("# Data-Grounded Writing Notes")
    lines.append("")
    lines.append("## Assumption")
    lines.append("The grouped analysis follows the four experiment groups specified by the author: Q01-Q03, Q04-Q06, Q07-Q08, and Q10-Q13. Q09 is excluded from the primary grouped figures and grouped paired tests because Diff-MST failed to produce a valid mix for that item, which prevents a fair within-question paired comparison across all four methods.")
    lines.append("The final experiment order is Experiment 1 (Intra-group Mixing Quality), Experiment 2a (Compensation of Grouping Errors), Experiment 2b (Compensation of Loudness Errors), and Experiment 3 (Full-Mix Ablation).")
    lines.append("")
    lines.append("Primary figures use box plots with overlaid jittered points for readability. Violin plots are retained as supplementary distribution visualizations.")
    lines.append("")
    lines.append("## Abstract Draft (181 words)")
    lines.append(
        "Automatic mixing aims to transform multitrack recordings into perceptually coherent and balanced mixes, yet this remains difficult in realistic production scenarios with many tracks and complex inter-track dependencies. Although recent two-stage systems outperform monolithic end-to-end models, it is still unclear whether these gains stem from better model architectures or from the explicit decomposition of the task itself. We present a controlled, subtask-oriented analysis of automatic mixing systems through four listening experiments that separately evaluate intra-group mixing quality, grouping-error compensation, loudness-error compensation, and full-mix ablations. The results show that generative modeling transfers more effectively than parameter prediction to the simpler intra-group task, while downstream correction of grouping and loudness errors remains limited and condition-dependent. In particular, 7-group processing yields the strongest downstream results, whereas loudness-balanced inputs provide only a positive but non-significant trend for MEGAMI. Finally, two-stage decomposition substantially improves full-mix quality and narrows the gap to human references, supporting structured pipelines such as ProMix."
    )
    lines.append("")
    lines.append("## Results Text Suggestions")
    lines.append("")
    lines.append("### Experiment 1: Intra-group Mixing Quality")
    lines.append(
        "For Q07-Q08, `ELL` and `MEGAMI` both achieved strong ratings, and their pooled difference was not statistically significant, while both clearly outperformed `DiffMST` and `NoMix`. This supports a more nuanced interpretation than simply saying end-to-end models fail on subtasks: parameter-prediction models optimized for complete mixing appear poorly matched to intra-group objectives, but a generative model such as MEGAMI can still transfer reasonably well to the simpler local-balancing task. Even so, the rule-based equal-local-loudness baseline remains especially attractive in practice because it reaches comparable perceptual quality with much lower deployment complexity and computational cost."
    )
    lines.append("")
    lines.append("### Experiment 2: Compensation of Intra-group Errors")
    lines.append(
        "Within an intra-group stage, the most important local relationships include loudness balance, panning balance, and grouping structure. However, current machine-learning mixing models typically split stereo tracks into mono channels and apply loudness normalization to all mono inputs. This destroys the left-right loudness differences that preserve panning information, making panning errors unsuitable as a controlled and identifiable perturbation for downstream compensation analysis. For this reason, we focus only on two controlled error types in Experiment 2: grouping errors and loudness errors."
    )
    lines.append("")
    lines.append("### Experiment 2a: Compensation of Grouping Errors")
    lines.append(
        "In Q04-Q06, the 7-group condition achieved the strongest scores overall, particularly for MEGAMI, while the instrument-based grouping condition was markedly weaker. The paired tests show significant advantages of `7G-M` over both `7G-D` and `I-D`, with similarly strong effects for `4G-M` over the weakest grouping condition. These results indicate that grouping strategy is not a neutral preprocessing choice: inappropriate grouping degrades the signal presented to the downstream model, and the resulting quality loss is not fully recoverable."
    )
    lines.append("")
    lines.append("### Experiment 2b: Compensation of Loudness Errors")
    lines.append(
        "Across Q01-Q03, the `WB-MEGAMI` condition achieved the highest average scores, while both Diff-MST variants remained substantially lower. However, the direct within-model comparison between `WB-MEGAMI` and `NB-MEGAMI` showed only a non-significant positive trend, and the corresponding comparison between `NB-DiffMST` and `WB-DiffMST` was also not significant. This suggests that downstream processing can partially compensate for incorrect intra-group loudness relationships, but such compensation is not sufficiently stable or reliable. When the input loudness structure is correct, MEGAMI shows a positive tendency, but the downstream model does not fully repair early-stage loudness errors."
    )
    lines.append("")
    lines.append("### Experiment 3: Full-Mix Ablation")
    lines.append(
        "In Q10-Q13, the key ablation result is the consistent gain from explicit decomposition within each model family. `2S-M` significantly outperformed monolithic `M`, and `2S-D` also significantly outperformed monolithic `D`. At the same time, `2S-M` approached the `Human` reference closely, with no statistically significant difference between the two in the pooled paired tests. These comparisons provide direct evidence that the performance gain arises from structured decomposition itself rather than from swapping to an unrelated stronger baseline."
    )
    lines.append("")
    lines.append("### Discussion Link to Research Questions")
    lines.append(
        "Taken together, the four experiments support all three research questions. First, downstream models do not reliably compensate for incorrect early-stage grouping or loudness decisions, which highlights the importance of preserving correct local structure before inter-group processing. Second, end-to-end models trained on complete mixing do not transfer cleanly to simpler subtasks such as intra-group balancing, although generative modeling shows better transfer than parameter prediction. Third, explicit decomposition improves full-mix performance and narrows the gap to human references, motivating practical structured systems such as ProMix."
    )
    lines.append("")
    lines.append("### Recommended Treatment of Q09")
    lines.append(
        "The cleanest paper treatment is to exclude Q09 from the primary grouped statistical analysis and the primary Experiment 1 figure, while explicitly stating that Diff-MST failed to generate a valid output for this item. This avoids introducing an unbalanced comparison that would weaken the paired-test design. If desired, Q09 can still be shown in a supplementary figure or appendix table as an incomplete case, with a caption note that no inferential comparison involving Diff-MST was performed for that item."
    )
    lines.append("")
    lines.append("## Statistical Highlights")
    ordered_group_keys = [
        "exp1_intra_group_quality",
        "exp2a_grouping_compensation",
        "exp2b_loudness_compensation",
        "exp3_full_mix_ablation",
    ]
    for group_key in ordered_group_keys:
        stats_df = key_stats[group_key]
        lines.append(f"### {EXPERIMENT_GROUPS[group_key]['title']}")
        for _, row in stats_df.iterrows():
            lines.append(
                f"- {row['model_a_label']} ({row['model_a_short']}) vs {row['model_b_label']} ({row['model_b_short']}): "
                f"mean difference = {row['mean_diff_a_minus_b']:.2f}, p = {row['p_value']:.3g}, "
                f"FDR-BH p = {row['p_value_fdr_bh']:.3g}, winner = {row['winner_by_mean']} ({row['winner_by_mean_short']})."
            )
        lines.append("")

    output = BASE_DIR / "paper_results_analysis.md"
    output.write_text("\n".join(lines), encoding="utf-8")
    return output


def main() -> None:
    apply_style()
    df = load_data()

    created_files = []
    for group_key, config in EXPERIMENT_GROUPS.items():
        summarize_group(df, group_key, config)
        created_files.extend(draw_violin_plot(group_key, config, df))
        created_files.extend(draw_box_plot(group_key, config, df))
        created_files.append(paired_tests(group_key, config, df))

    analysis_path = build_paper_analysis(df)
    created_files.append(analysis_path)
    created_files.extend(draw_submission_figure2(df))
    created_files.extend(draw_submission_figure3(df))
    created_files.extend(draw_submission_figure4(df))

    print("Created files:")
    for path in created_files:
        print(path)


if __name__ == "__main__":
    main()
