import numpy as np, pandas as pd, pingouin as pg
from scipy import stats
d = pd.read_csv("gtt.csv")

def fmt(p):
    return "< 0.001" if p < 0.001 else f"= {p:.3f}"

def auc(g):                      # trapezoids over the actual time points (handles uneven spacing)
    g = g.sort_values("time")
    return np.trapezoid(g.glucose, g.time)

def iauc(g):                     # net incremental AUC: glucose minus the mouse's own 0-min value
    g = g.sort_values("time")
    return np.trapezoid(g.glucose - g.glucose.iloc[0], g.time)

per = d.groupby(["group", "mouse"]).apply(lambda g: pd.Series({"auc": auc(g), "iauc": iauc(g), "fasting": g.glucose.iloc[0]}), include_groups=False).reset_index()
for col, label in [("fasting", "Fasting glucose (mmol/L)"), ("auc", "Total AUC (mmol/L x min)"), ("iauc", "Incremental AUC (mmol/L x min)")]:
    h, c = per.loc[per.group == "HFD", col], per.loc[per.group == "chow", col]
    r = stats.ttest_ind(h, c, equal_var=False)
    ci = r.confidence_interval()
    print(f"{label}: chow {c.mean():.1f}, HFD {h.mean():.1f}, difference {h.mean() - c.mean():.1f} "
          f"(95% CI {ci.low:.1f} to {ci.high:.1f}), Welch t = {r.statistic:.2f}, df = {r.df:.1f}, p {fmt(r.pvalue)}")

# Time course: mixed ANOVA (group between, time within) with the Greenhouse-Geisser correction
aov = pg.mixed_anova(data=d, dv="glucose", within="time", between="group", subject="mouse", correction=True)
row = aov[aov.Source == "Interaction"].iloc[0]
grp = aov[aov.Source == "group"].iloc[0]
print(f"Group (curve height): F({grp.DF1:.0f}, {grp.DF2:.0f}) = {grp.F:.2f}, p {fmt(grp['p_unc'])}")
print(f"Group x time: F({row.DF1:.0f}, {row.DF2:.0f}) = {row.F:.2f}, Greenhouse-Geisser epsilon = {row.eps:.2f}, "
      f"corrected p = {row['p_GG_corr']:.3f} (uncorrected {row['p_unc']:.3f})")

# What a t-test at every time point would show (shown for comparison; not recommended)
for t, g in d.groupby("time"):
    h, c = g.loc[g.group == "HFD", "glucose"], g.loc[g.group == "chow", "glucose"]
    print(f"{t:>3} min: HFD - chow = {h.mean() - c.mean():.2f} mmol/L, Welch p {fmt(stats.ttest_ind(h, c, equal_var=False).pvalue)}")