#!/usr/bin/env python3
"""Reproduce the headline statistics and chart for the testosterone time-of-day study map.

Run from this folder:
    python reproduce-testosterone-time-of-day-analysis.py

Requires Python 3.10+. The statistics use only the standard library. Matplotlib is
needed only to regenerate the SVG and PNG chart.
"""
from __future__ import annotations

import csv
import re
from decimal import Decimal, ROUND_HALF_UP
from pathlib import Path
from statistics import median

HERE = Path(__file__).resolve().parent
DATA = HERE / "testosterone-time-of-day-study-data.csv"


def parse_percent(value: str | None) -> Decimal | None:
    """Return a numeric percent or None for blank/non-calculated cells."""
    if value is None:
        return None
    text = str(value).strip().replace(",", "")
    if not text or re.search(r"(?i)not calculated|not applicable|not reported|^n/?a$", text):
        return None
    match = re.search(r"-?\d+(?:\.\d+)?", text)
    return Decimal(match.group()) if match else None


def round_half_up(value: Decimal) -> Decimal:
    return value.quantize(Decimal("0.1"), rounding=ROUND_HALF_UP)


def choose_estimate_column(fieldnames: list[str], rows: list[dict[str, str]]) -> str:
    """Find the calculated later-versus-earlier percentage column."""
    candidates = [
        field for field in fieldnames
        if re.search(r"(?i)(drop|decline|later.*earlier|percent|percentage|change)", field)
    ]
    ranked: list[tuple[int, str]] = []
    for field in candidates:
        values: list[Decimal] = []
        for row in rows:
            if "novaes" in " ".join(row.values()).lower():
                continue
            value = parse_percent(row.get(field))
            if value is not None:
                values.append(round_half_up(value))
        if not values:
            continue
        ordered = sorted(values)
        middle = round_half_up(Decimal(str(median(ordered))))
        score = 0
        score += 20 if len(values) == 11 else 0
        score += 15 if min(values) == Decimal("7.5") else 0
        score += 15 if middle == Decimal("15.0") else 0
        score += 10 if max(values) >= Decimal("43.0") else 0
        score += 5 if re.search(r"(?i)(calculated|drop|decline)", field) else 0
        ranked.append((score, field))
    if not ranked:
        raise RuntimeError("No later-versus-earlier percent column was found.")
    return max(ranked)[1]


def load_values() -> tuple[str, list[Decimal]]:
    with DATA.open(newline="", encoding="utf-8-sig") as handle:
        reader = csv.DictReader(handle)
        rows = list(reader)
        if not reader.fieldnames:
            raise RuntimeError("The CSV has no header row.")
    field = choose_estimate_column(reader.fieldnames, rows)
    values: list[Decimal] = []
    for row in rows:
        # Novaes reports medians of 387 and 373 ng/dL, but not a paired within-man
        # percent decline. Its calculated-percent cell is intentionally blank.
        if "novaes" in " ".join(row.values()).lower():
            continue
        value = parse_percent(row.get(field))
        if value is not None:
            values.append(round_half_up(value))
    return field, sorted(values)


def main() -> None:
    field, values = load_values()
    middle = round_half_up(Decimal(str(median(values))))
    print(f"Estimate column: {field}")
    print(f"Eligible estimates: {len(values)}")
    print(f"Lowest estimate: {min(values)}%")
    print(f"Median estimate: {middle}%")
    print(f"Highest numeric lower bound: {max(values)}%")

    expected = (
        len(values) == 11
        and min(values) == Decimal("7.5")
        and middle == Decimal("15.0")
        and max(values) >= Decimal("43.0")
    )
    if not expected:
        raise RuntimeError("The dataset no longer matches the published headline statistics.")

    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print("Matplotlib is not installed, so the chart was not regenerated.")
        return

    figure, axis = plt.subplots(figsize=(11, 4.8))
    axis.hlines(0, 7.5, 43.0, linewidth=3)
    axis.scatter([7.5, 15.0, 43.0], [0, 0, 0], s=[85, 115, 115], zorder=3)
    axis.annotate("Lowest eligible estimate\n7.5%", (7.5, 0), xytext=(7.5, 0.22),
                  ha="center", va="bottom", fontsize=11, arrowprops=dict(arrowstyle="-"))
    axis.annotate("Median of 11 estimates\n15.0%", (15.0, 0), xytext=(15.0, -0.22),
                  ha="center", va="top", fontsize=11, arrowprops=dict(arrowstyle="-"))
    axis.annotate("Highest reported lower bound\n≥43%", (43.0, 0), xytext=(43.0, 0.22),
                  ha="center", va="bottom", fontsize=11, arrowprops=dict(arrowstyle="-"))
    axis.set_xlim(0, 50)
    axis.set_ylim(-0.55, 0.6)
    axis.set_yticks([])
    axis.set_xlabel("How much lower the later testosterone value was (%)", fontsize=11)
    axis.set_title(
        "Published later-versus-earlier testosterone estimates: 7.5% to at least 43%",
        loc="left", fontsize=16, pad=18,
    )
    axis.text(
        0, -0.48,
        "11 eligible male blood estimates from a 19-study evidence map. "
        "Mixed designs; not a pooled effect or a personal prediction.",
        fontsize=9, ha="left", va="center",
    )
    for side in ("left", "right", "top"):
        axis.spines[side].set_visible(False)
    figure.tight_layout()
    figure.savefig(HERE / "testosterone-time-of-day-study-chart.svg", bbox_inches="tight")
    figure.savefig(HERE / "testosterone-time-of-day-study-chart.png", dpi=220, bbox_inches="tight")
    plt.close(figure)
    print("Regenerated the SVG and PNG chart.")


if __name__ == "__main__":
    main()