#!/usr/bin/env python3
"""
Reproducible analysis: does help-desk case resolution time differ by service level?

Pre-registered question
------------------------
In a help-desk process event log, does the case-level time-to-resolution differ
across service_level tiers?

Design
------
- Unit of analysis is the CASE (the event log has multiple event rows per case),
  so the log is aggregated to one value per case before any test. This guards
  against pseudoreplication.
- Primary outcome: resolution time = hours from a case's first event to its first
  "Resolve ticket" activity.
- Grouping variable: modal service_level per case, restricted to tiers with
  adequate sample size (Value 1, Value 2, Value 3). Tier "Value 4" (n=3 cases) is
  excluded a priori.
- Primary test: Kruskal-Wallis H test (nonparametric; resolution time is heavily
  right-skewed) across the three tiers, with epsilon-squared effect size.
- Follow-up: pairwise Mann-Whitney U tests (3 comparisons) with Holm correction
  for multiple comparisons, and rank-biserial correlation as the pairwise effect
  size.

Robustness / sensitivity checks (pre-registered)
------------------------------------------------
- R1 (alternative outcome): Kruskal-Wallis on full throughput-to-Closed time, to
  show whether any tier difference is specific to handling time or driven by
  administrative closure timing.
- R2 (attribute ambiguity): re-run the primary test after dropping the cases whose
  service_level is not constant within the case.
- R3 (attribute definition): re-run the primary test using the first-observed
  service_level per case instead of the modal value.

Outputs (written under the working directory)
----------------------------------------------
- figures/fig-1-resolution-by-service-level.png
- figures/fig-2-robustness-sensitivity.png
- tables/tbl-1-resolution-summary-by-tier.csv
- tables/tbl-2-robustness-sensitivity.csv
- tables/datapackage.json   (Frictionless)
- results.json              (every statistic)

Re-running this script reproduces every number, figure, and table.
"""

import os
import sys
import json
import hashlib
import random
from itertools import combinations

import numpy as np
import pandas as pd
from scipy import stats

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib as mpl

# --------------------------------------------------------------------------- #
# 0. Configuration and reproducibility
# --------------------------------------------------------------------------- #
SEED = 20240517
random.seed(SEED)
np.random.seed(SEED)

DATA_URL = "https://ndownloader.figshare.com/files/23993303"
DATA_MD5 = "a91d74377c55be5cbdfda05dae69db30"
# md5-gated local fallback (used only when the canonical download is unreachable)
LOCAL_FALLBACK = "/Users/markhahnel/datasetpapers-handfeed/italian-helpdesk/finale.csv"
CACHE_NAME = "finale.csv"

TIERS = ["Value 1", "Value 2", "Value 3"]   # service_level tiers analysed
EXCLUDED_TIER = "Value 4"                    # n=3 cases, excluded a priori
RESOLVE_ACTIVITY = "Resolve ticket"
CLOSE_ACTIVITY = "Closed"
TS_FORMAT = "%Y/%m/%d %H:%M:%S.%f"

FIG_DIR = "figures"
TBL_DIR = "tables"

# Colourblind-safe palette (Okabe-Ito subset), one hue per tier
TIER_COLORS = {
    "Value 1": "#0072B2",   # blue
    "Value 2": "#E69F00",   # orange
    "Value 3": "#009E73",   # bluish green
}


# --------------------------------------------------------------------------- #
# 1. Data acquisition: canonical download + md5-gated local fallback
# --------------------------------------------------------------------------- #
def _md5(path, chunk=1 << 20):
    h = hashlib.md5()
    with open(path, "rb") as fh:
        for block in iter(lambda: fh.read(chunk), b""):
            h.update(block)
    return h.hexdigest()


def acquire_data():
    """Return a path to an md5-verified copy of finale.csv.

    Strategy:
      1. If a cached copy already exists and matches DATA_MD5, use it.
      2. Otherwise attempt the canonical figshare download; verify md5.
      3. If the download fails (offline / blocked), fall back to the local copy,
         which is *also* md5-gated: a mismatch is a hard failure.
    """
    if os.path.exists(CACHE_NAME) and _md5(CACHE_NAME) == DATA_MD5:
        print(f"[data] using verified cache: {CACHE_NAME}")
        return CACHE_NAME

    # 2. canonical download
    try:
        import urllib.request
        print(f"[data] downloading canonical source: {DATA_URL}")
        urllib.request.urlretrieve(DATA_URL, CACHE_NAME)
        got = _md5(CACHE_NAME)
        if got != DATA_MD5:
            raise ValueError(f"downloaded md5 {got} != expected {DATA_MD5}")
        print("[data] download verified against expected md5")
        return CACHE_NAME
    except Exception as exc:  # noqa: BLE001
        print(f"[data] canonical download unavailable ({exc}); trying local fallback")

    # 3. md5-gated local fallback
    if os.path.exists(LOCAL_FALLBACK):
        got = _md5(LOCAL_FALLBACK)
        if got != DATA_MD5:
            raise SystemExit(
                f"[data] FATAL: local fallback md5 {got} != expected {DATA_MD5}"
            )
        print(f"[data] using md5-verified local fallback: {LOCAL_FALLBACK}")
        return LOCAL_FALLBACK

    raise SystemExit(
        "[data] FATAL: no verified copy of finale.csv available "
        "(canonical download failed and no local fallback found)."
    )


# --------------------------------------------------------------------------- #
# 2. Aggregate event log -> one row per case
# --------------------------------------------------------------------------- #
def build_case_table(df):
    df = df.copy()
    df["ts"] = pd.to_datetime(df["Complete Timestamp"], format=TS_FORMAT)

    def per_case(sub):
        sub = sub.sort_values("ts")
        start = sub["ts"].min()
        end = sub["ts"].max()
        resolves = sub.loc[sub["Activity"] == RESOLVE_ACTIVITY, "ts"]
        first_resolve = resolves.min() if len(resolves) else pd.NaT
        sl_mode = sub["service_level"].mode().iloc[0]
        sl_first = sub["service_level"].iloc[0]
        return pd.Series({
            "n_events": len(sub),
            "throughput_h": (end - start).total_seconds() / 3600.0,
            "resolve_h": ((first_resolve - start).total_seconds() / 3600.0
                          if pd.notna(first_resolve) else np.nan),
            "service_level": sl_mode,
            "service_level_first": sl_first,
            "sl_ambiguous": sub["service_level"].nunique() > 1,
        })

    case = df.groupby("Case ID", sort=False).apply(per_case, include_groups=False)
    return case


# --------------------------------------------------------------------------- #
# 3. Statistics helpers
# --------------------------------------------------------------------------- #
def epsilon_squared(H, n, k):
    """Epsilon-squared effect size for Kruskal-Wallis."""
    return (H - k + 1) / (n - k)


def holm_correction(pvals):
    """Holm-Bonferroni corrected p-values, preserving input order."""
    m = len(pvals)
    order = np.argsort(pvals)
    adj = [None] * m
    running = 0.0
    for rank, idx in enumerate(order):
        val = (m - rank) * pvals[idx]
        running = max(running, val)
        adj[idx] = min(running, 1.0)
    return adj


def rank_biserial(x, y, U):
    """Rank-biserial correlation from a Mann-Whitney U (x vs y)."""
    return 1.0 - 2.0 * U / (len(x) * len(y))


def kruskal_block(frame, outcome, group_col, tiers):
    sub = frame.dropna(subset=[outcome])
    sub = sub[sub[group_col].isin(tiers)]
    groups = [sub.loc[sub[group_col] == t, outcome].values for t in tiers]
    groups = [g for g in groups if len(g) > 0]
    H, p = stats.kruskal(*groups)
    k = len(groups)
    n = sum(len(g) for g in groups)
    eps2 = epsilon_squared(H, n, k)
    return {
        "H": float(H), "df": int(k - 1), "p": float(p),
        "epsilon_squared": float(eps2), "n": int(n), "k": int(k),
        "group_n": {t: int((sub[group_col] == t).sum()) for t in tiers},
        "group_median_h": {t: float(np.median(sub.loc[sub[group_col] == t, outcome]))
                           for t in tiers if (sub[group_col] == t).any()},
    }


# --------------------------------------------------------------------------- #
# 4. Figures
# --------------------------------------------------------------------------- #
def make_fig1(prim, results, path):
    fig, ax = plt.subplots(figsize=(6.0, 4.2))
    data = [prim.loc[prim.service_level == t, "resolve_h"].values for t in TIERS]
    positions = range(1, len(TIERS) + 1)

    # violin for shape + jittered points behind + median tick
    vp = ax.violinplot(data, positions=positions, widths=0.8,
                       showmeans=False, showmedians=False, showextrema=False)
    for body, t in zip(vp["bodies"], TIERS):
        body.set_facecolor(TIER_COLORS[t])
        body.set_alpha(0.25)
        body.set_edgecolor(TIER_COLORS[t])

    rng = np.random.default_rng(SEED)
    for pos, t, d in zip(positions, TIERS, data):
        jit = rng.uniform(-0.12, 0.12, size=len(d))
        ax.scatter(pos + jit, d, s=6, color=TIER_COLORS[t], alpha=0.30,
                   edgecolors="none", zorder=2)
        med = np.median(d)
        ax.hlines(med, pos - 0.30, pos + 0.30, color="black", lw=2.2, zorder=4)
        # median value labelled centred just above the tick, on a white box,
        # so labels never collide even when medians are close in value
        ax.text(pos, med * 1.35, f"median\n{med:.0f} h",
                va="bottom", ha="center", fontsize=6.8, zorder=5,
                bbox=dict(boxstyle="round,pad=0.15", fc="white",
                          ec="0.6", lw=0.5, alpha=0.9))

    ax.set_yscale("log")
    ax.set_ylim(0.05, 2000)
    ax.set_xticks(list(positions))
    ax.set_xticklabels([f"{t}\n(n={results['primary']['group_n'][t]})" for t in TIERS])
    ax.set_xlabel("Service level")
    ax.set_ylabel("Time to first resolution (hours, log scale)")
    kw = results["primary"]
    ax.set_title(
        "Resolution time barely differs across service tiers\n"
        f"Kruskal-Wallis H={kw['H']:.2f}, p={kw['p']:.3f}, "
        f"$\\epsilon^2$={kw['epsilon_squared']:.4f} (negligible), N={kw['n']}",
        fontsize=9, loc="left")
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    fig.tight_layout()
    fig.savefig(path, dpi=200)
    plt.close(fig)


def make_fig2(results, path):
    """Compare the KW effect (epsilon-squared) and p across analysis variants."""
    variants = [
        ("Primary\n(resolution time)", results["primary"]),
        ("R1 alt outcome\n(throughput to Closed)", results["robustness"]["R1_throughput"]),
        ("R2 drop ambiguous\ncases", results["robustness"]["R2_drop_ambiguous"]),
        ("R3 first-observed\nservice level", results["robustness"]["R3_first_observed"]),
    ]
    labels = [v[0] for v in variants]
    eps = [v[1]["epsilon_squared"] for v in variants]
    ps = [v[1]["p"] for v in variants]
    ns = [v[1]["n"] for v in variants]
    colors = ["#0072B2", "#D55E00", "#56B4E9", "#56B4E9"]

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8.4, 4.2))

    ypos = range(len(labels))
    ax1.hlines(list(ypos), 0, eps, color=colors, lw=2, zorder=1)
    ax1.scatter(eps, list(ypos), color=colors, s=60, zorder=2)
    for y, e, n in zip(ypos, eps, ns):
        ax1.text(e + max(eps) * 0.03, y, f"{e:.4f} (n={n})", va="center", fontsize=7)
    ax1.set_yticks(list(ypos))
    ax1.set_yticklabels(labels, fontsize=7)
    ax1.invert_yaxis()
    ax1.set_xlabel("Effect size (epsilon-squared)")
    ax1.set_xlim(0, max(eps) * 1.35)
    ax1.set_title("Effect size stays small in every variant", fontsize=9, loc="left")
    ax1.axvspan(0, 0.01, color="grey", alpha=0.12)
    ax1.text(0.005, len(labels) - 0.5, "negligible", fontsize=6.5,
             ha="center", color="grey")

    ax2.scatter(ps, list(ypos), color=colors, s=60, zorder=2)
    for y, p in zip(ypos, ps):
        ax2.text(p, y + 0.18, f"p={p:.1e}" if p < 1e-3 else f"p={p:.3f}",
                 va="bottom", ha="center", fontsize=7)
    ax2.axvline(0.05, color="black", ls="--", lw=1)
    ax2.text(0.058, len(labels) - 0.5, "p = 0.05", fontsize=7, ha="left",
             va="center", rotation=90)
    ax2.set_xscale("log")
    ax2.set_xlim(1e-14, 0.3)
    ax2.set_yticks(list(ypos))
    ax2.set_yticklabels([])
    ax2.invert_yaxis()
    ax2.set_xlabel("p-value (log scale)")
    ax2.set_title("Significance is driven by the closure-time outcome", fontsize=9, loc="left")

    for ax in (ax1, ax2):
        for s in ("top", "right"):
            ax.spines[s].set_visible(False)
    fig.suptitle("Robustness and sensitivity of the service-level comparison",
                 fontsize=10, x=0.02, ha="left")
    fig.tight_layout(rect=(0, 0, 1, 0.96))
    fig.savefig(path, dpi=200)
    plt.close(fig)


# --------------------------------------------------------------------------- #
# 5. Main
# --------------------------------------------------------------------------- #
def main():
    os.makedirs(FIG_DIR, exist_ok=True)
    os.makedirs(TBL_DIR, exist_ok=True)

    data_path = acquire_data()
    df = pd.read_csv(data_path)

    n_events = len(df)
    n_cases = df["Case ID"].nunique()
    time_min = pd.to_datetime(df["Complete Timestamp"], format=TS_FORMAT).min()
    time_max = pd.to_datetime(df["Complete Timestamp"], format=TS_FORMAT).max()

    case = build_case_table(df)

    # ---- primary analysis frame ----
    prim = case.dropna(subset=["resolve_h"]).copy()
    prim = prim[prim.service_level.isin(TIERS)]

    results = {
        "meta": {
            "seed": SEED,
            "data_md5": DATA_MD5,
            "n_events": int(n_events),
            "n_cases": int(n_cases),
            "time_span": [str(time_min), str(time_max)],
            "outcome_primary": "hours from first case event to first 'Resolve ticket'",
            "grouping": "modal service_level per case",
            "tiers_analysed": TIERS,
            "tier_excluded_a_priori": {EXCLUDED_TIER: int((case.service_level == EXCLUDED_TIER).sum())},
            "n_cases_no_resolve_event": int(case["resolve_h"].isna().sum()),
            "n_cases_ambiguous_service_level": int(case["sl_ambiguous"].sum()),
        }
    }

    # ---- primary: Kruskal-Wallis ----
    results["primary"] = kruskal_block(prim, "resolve_h", "service_level", TIERS)

    # ---- pairwise Mann-Whitney U + Holm ----
    pairs = list(combinations(TIERS, 2))
    praw, blocks = [], []
    for a, b in pairs:
        xa = prim.loc[prim.service_level == a, "resolve_h"].values
        xb = prim.loc[prim.service_level == b, "resolve_h"].values
        U, p = stats.mannwhitneyu(xa, xb, alternative="two-sided")
        praw.append(p)
        blocks.append({"group_a": a, "group_b": b, "U": float(U), "p_raw": float(p),
                       "rank_biserial": float(rank_biserial(xa, xb, U)),
                       "n_a": int(len(xa)), "n_b": int(len(xb))})
    padj = holm_correction(praw)
    for blk, pa in zip(blocks, padj):
        blk["p_holm"] = float(pa)
    results["pairwise"] = blocks

    # ---- robustness / sensitivity ----
    R1 = kruskal_block(case, "throughput_h", "service_level", TIERS)
    prim_noamb = prim[~prim["sl_ambiguous"]]
    R2 = kruskal_block(prim_noamb, "resolve_h", "service_level", TIERS)
    R3 = kruskal_block(case.dropna(subset=["resolve_h"]),
                       "resolve_h", "service_level_first", TIERS)
    results["robustness"] = {
        "R1_throughput": R1,
        "R2_drop_ambiguous": R2,
        "R3_first_observed": R3,
    }

    # ---- tables ----
    # tbl-1: per-tier summary of resolution time
    rows = []
    for t in TIERS:
        d = prim.loc[prim.service_level == t, "resolve_h"].values
        rows.append({
            "service_level": t,
            "n_cases": len(d),
            "median_h": np.median(d),
            "mean_h": np.mean(d),
            "q25_h": np.percentile(d, 25),
            "q75_h": np.percentile(d, 75),
            "min_h": np.min(d),
            "max_h": np.max(d),
        })
    tbl1 = pd.DataFrame(rows)
    tbl1.to_csv(os.path.join(TBL_DIR, "tbl-1-resolution-summary-by-tier.csv"),
                index=False)

    # tbl-2: robustness/sensitivity results
    tbl2 = pd.DataFrame([
        {"variant": "primary_resolution_time", "outcome": "time to first Resolve ticket",
         "H": results["primary"]["H"], "df": results["primary"]["df"],
         "p": results["primary"]["p"], "epsilon_squared": results["primary"]["epsilon_squared"],
         "n": results["primary"]["n"]},
        {"variant": "R1_throughput_to_closed", "outcome": "time to Closed",
         "H": R1["H"], "df": R1["df"], "p": R1["p"],
         "epsilon_squared": R1["epsilon_squared"], "n": R1["n"]},
        {"variant": "R2_drop_ambiguous_cases", "outcome": "time to first Resolve ticket",
         "H": R2["H"], "df": R2["df"], "p": R2["p"],
         "epsilon_squared": R2["epsilon_squared"], "n": R2["n"]},
        {"variant": "R3_first_observed_service_level", "outcome": "time to first Resolve ticket",
         "H": R3["H"], "df": R3["df"], "p": R3["p"],
         "epsilon_squared": R3["epsilon_squared"], "n": R3["n"]},
    ])
    tbl2.to_csv(os.path.join(TBL_DIR, "tbl-2-robustness-sensitivity.csv"), index=False)

    # ---- Frictionless datapackage ----
    datapackage = build_datapackage(tbl1, tbl2)
    with open(os.path.join(TBL_DIR, "datapackage.json"), "w") as fh:
        json.dump(datapackage, fh, indent=2)

    # ---- figures ----
    make_fig1(prim, results, os.path.join(FIG_DIR, "fig-1-resolution-by-service-level.png"))
    make_fig2(results, os.path.join(FIG_DIR, "fig-2-robustness-sensitivity.png"))

    # ---- results.json ----
    with open("results.json", "w") as fh:
        json.dump(results, fh, indent=2)

    print("[done] wrote results.json, figures/, tables/")
    print(json.dumps({
        "primary_p": results["primary"]["p"],
        "primary_eps2": results["primary"]["epsilon_squared"],
        "pairwise_holm": {f"{b['group_a']} vs {b['group_b']}": b["p_holm"] for b in blocks},
    }, indent=2))
    return results


def build_datapackage(tbl1, tbl2):
    def fields_for(df, descriptions):
        out = []
        for col in df.columns:
            dtype = df[col].dtype
            if pd.api.types.is_integer_dtype(dtype):
                ftype = "integer"
            elif pd.api.types.is_float_dtype(dtype):
                ftype = "number"
            else:
                ftype = "string"
            out.append({"name": col, "type": ftype,
                        "description": descriptions.get(col, "")})
        return out

    desc1 = {
        "service_level": "Service-level tier (modal value per case).",
        "n_cases": "Number of cases in the tier.",
        "median_h": "Median time to first resolution, hours.",
        "mean_h": "Mean time to first resolution, hours.",
        "q25_h": "25th percentile of resolution time, hours.",
        "q75_h": "75th percentile of resolution time, hours.",
        "min_h": "Minimum resolution time, hours.",
        "max_h": "Maximum resolution time, hours.",
    }
    desc2 = {
        "variant": "Analysis variant (primary or robustness/sensitivity check).",
        "outcome": "Outcome measure used in the variant.",
        "H": "Kruskal-Wallis H statistic.",
        "df": "Degrees of freedom.",
        "p": "Uncorrected p-value.",
        "epsilon_squared": "Epsilon-squared effect size.",
        "n": "Number of cases in the variant.",
    }
    return {
        "profile": "tabular-data-package",
        "name": "helpdesk-resolution-by-service-level",
        "title": "Help-desk case resolution time by service level",
        "licenses": [{"name": "4TU-General-Terms-of-Use",
                      "title": "4TU.ResearchData General Terms of Use"}],
        "resources": [
            {
                "name": "tbl-1-resolution-summary-by-tier",
                "path": "tbl-1-resolution-summary-by-tier.csv",
                "format": "csv", "mediatype": "text/csv", "encoding": "utf-8",
                "profile": "tabular-data-resource",
                "schema": {"fields": fields_for(tbl1, desc1)},
            },
            {
                "name": "tbl-2-robustness-sensitivity",
                "path": "tbl-2-robustness-sensitivity.csv",
                "format": "csv", "mediatype": "text/csv", "encoding": "utf-8",
                "profile": "tabular-data-resource",
                "schema": {"fields": fields_for(tbl2, desc2)},
            },
        ],
    }


if __name__ == "__main__":
    main()