"""Reproducible analysis for the migraine LLM study.

Model mapping used in the dataset:
- ChatGPT columns: GPT-4o (OpenAI)
- DeepSeek columns: DeepSeek-V3 (DeepSeek)
Responses were collected through the public web interfaces in May 2026.

Input: Additional_file_1_Migraine_LLM_dataset.xlsx (sheet: Dataset)
Python: 3.11+
Packages: pandas, numpy, scipy, scikit-learn, statsmodels, openpyxl
"""
from __future__ import annotations

import math
import re
from pathlib import Path

import numpy as np
import pandas as pd
from scipy.stats import norm, wilcoxon
from sklearn.metrics import cohen_kappa_score
from statsmodels.stats.multitest import multipletests

INPUT = Path("Additional_file_1_Migraine_LLM_dataset.xlsx")
OUTPUT = Path("analysis_results.csv")

VOWELS = set("aeıioöuüAEIİOÖUÜ")


def word_tokens(text: str) -> list[str]:
    return re.findall(
        r"[A-Za-zÇĞİÖŞÜçğıöşüÂÎÛâîû]+(?:['’-][A-Za-zÇĞİÖŞÜçğıöşüÂÎÛâîû]+)*",
        str(text),
    )


def sentence_count(text: str) -> int:
    return max(1, len(re.findall(r"[.!?]+(?:[\"'”’)\]]*)", str(text))))


def atesman_score(text: str) -> float:
    words = word_tokens(text)
    n_words = len(words)
    if n_words == 0:
        return np.nan
    n_syllables = sum(sum(char in VOWELS for char in word) for word in words)
    n_sentences = sentence_count(text)
    raw = 198.825 - 40.175 * (n_syllables / n_words) - 2.610 * (n_words / n_sentences)
    return float(np.clip(raw, 0, 100))


def paired_test(deepseek: pd.Series, chatgpt: pd.Series) -> dict[str, float]:
    x = deepseek.astype(float).to_numpy()
    y = chatgpt.astype(float).to_numpy()
    differences = x - y
    nonzero = differences[differences != 0]
    result = wilcoxon(x, y, zero_method="wilcox", alternative="two-sided", method="auto")
    direction = np.sign(np.mean(differences))
    z_value = float(direction * norm.isf(result.pvalue / 2)) if result.pvalue > 0 else np.inf
    effect_r = abs(z_value) / math.sqrt(len(nonzero)) if len(nonzero) else 0.0
    return {
        "W": float(result.statistic),
        "p": float(result.pvalue),
        "z": z_value,
        "effect_r": effect_r,
        "nonzero_pairs": int(len(nonzero)),
    }


def descriptive(series: pd.Series) -> dict[str, float]:
    values = series.astype(float)
    return {
        "n": int(values.count()),
        "mean": float(values.mean()),
        "sd": float(values.std(ddof=1)),
        "median": float(values.median()),
        "q1": float(values.quantile(0.25)),
        "q3": float(values.quantile(0.75)),
        "min": float(values.min()),
        "max": float(values.max()),
    }


def main() -> None:
    data = pd.read_excel(INPUT, sheet_name="Dataset")

    # Recalculate objective readability from each response text.
    data["ChatGPT_Atesman_recalculated"] = data["ChatGPT_response"].map(atesman_score)
    data["DeepSeek_Atesman_recalculated"] = data["DeepSeek_response"].map(atesman_score)

    outcomes = ["Accuracy", "Completeness", "Empathy", "Safety", "Atesman"]
    results: list[dict[str, object]] = []

    # Descriptive results.
    for model in ("ChatGPT", "DeepSeek"):
        for outcome in outcomes:
            col = (
                f"{model}_Atesman_recalculated"
                if outcome == "Atesman"
                else f"{model}_mean_{outcome}"
            )
            results.append({
                "analysis": "overall descriptive",
                "model_or_comparison": f"{model} - {outcome}",
                **descriptive(data[col]),
            })

    # Primary paired model comparisons, with Holm correction across five outcomes.
    primary_tests = []
    for outcome in outcomes:
        ds_col = "DeepSeek_Atesman_recalculated" if outcome == "Atesman" else f"DeepSeek_mean_{outcome}"
        cg_col = "ChatGPT_Atesman_recalculated" if outcome == "Atesman" else f"ChatGPT_mean_{outcome}"
        test = paired_test(data[ds_col], data[cg_col])
        primary_tests.append((outcome, test))
    adjusted = multipletests([test["p"] for _, test in primary_tests], method="holm")[1]
    for (outcome, test), adjusted_p in zip(primary_tests, adjusted):
        results.append({
            "analysis": "overall paired comparison",
            "model_or_comparison": f"DeepSeek vs ChatGPT - {outcome}",
            **test,
            "holm_p": float(adjusted_p),
        })

    # Prespecified red-flag subgroup.
    red = data[data["Red_flag"].astype(bool)].copy()
    red_tests = []
    for outcome in outcomes:
        ds_col = "DeepSeek_Atesman_recalculated" if outcome == "Atesman" else f"DeepSeek_mean_{outcome}"
        cg_col = "ChatGPT_Atesman_recalculated" if outcome == "Atesman" else f"ChatGPT_mean_{outcome}"
        test = paired_test(red[ds_col], red[cg_col])
        red_tests.append((outcome, test))
    adjusted_red = multipletests([test["p"] for _, test in red_tests], method="holm")[1]
    for (outcome, test), adjusted_p in zip(red_tests, adjusted_red):
        results.append({
            "analysis": "red-flag paired comparison",
            "model_or_comparison": f"DeepSeek vs ChatGPT - {outcome}",
            **test,
            "holm_p": float(adjusted_p),
        })

    # Exact agreement and quadratic-weighted kappa for the four clinician-rated outcomes.
    for model in ("ChatGPT", "DeepSeek"):
        for outcome in ("Accuracy", "Completeness", "Empathy", "Safety"):
            first = data[f"{model}_Burak_{outcome}"].astype(int)
            second = data[f"{model}_Hasan_{outcome}"].astype(int)
            results.append({
                "analysis": "inter-rater agreement",
                "model_or_comparison": f"{model} - {outcome}",
                "exact_agreement": float((first == second).mean()),
                "quadratic_weighted_kappa": float(cohen_kappa_score(first, second, weights="quadratic")),
            })

    pd.DataFrame(results).to_csv(OUTPUT, index=False)
    print(f"Saved {OUTPUT.resolve()}")


if __name__ == "__main__":
    main()