"""
04-02  Frequency Distribution Builder
Chapters 02-01, 02-02

Run from the project folder:   python python/analysis.py
"""
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

x = pd.read_csv("data/service_times.csv")["minutes"].to_numpy()
n = x.size

# ── 1. CLASS DESIGN ─────────────────────────────────────────────────────────
print("\n===== 1. CLASS DESIGN =====")

rng = x.max() - x.min()
k_sturges = int(np.ceil(1 + 3.322 * np.log10(n)))
k_sqrt = int(np.ceil(np.sqrt(n)))
w_raw = int(np.ceil(rng / k_sturges))
w = int(np.ceil(w_raw / 10) * 10)          # round up to a clean multiple of 10
start = int(np.floor(x.min() / w) * w)

print(f"n                 {n}")
print(f"min / max         {x.min():.1f} / {x.max():.1f}")
print(f"range R           {rng:.1f}")
print(f"k (Sturges)       {k_sturges}      k (sqrt n) = {k_sqrt}")
print(f"raw width         {w_raw}  ->  chosen width w = {w}")
print(f"first lower limit {start}")

n_classes = int(np.ceil((x.max() - start) / w)) + 1
lower = start + w * np.arange(n_classes)
upper = lower + w - 1
edges = np.append(lower - 0.5, upper[-1] + 0.5)

# ── 2. THE FREQUENCY TABLE ──────────────────────────────────────────────────
print("\n===== 2. FREQUENCY DISTRIBUTION =====")

f, _ = np.histogram(x, bins=edges)

dist = pd.DataFrame({
    "class":     [f"{lo}-{up}" for lo, up in zip(lower, upper)],
    "lower_bnd": edges[:-1],
    "upper_bnd": edges[1:],
    "midpoint":  (lower + upper) / 2,
    "f":         f,
    "rel_f":     (f / n).round(4),
    "cum_f":     f.cumsum(),
    "cum_pct":   (100 * f.cumsum() / n).round(1),
})
print(dist.to_string(index=False))

print(f"\nCHECK  sum(f) = {f.sum()}  (must equal n = {n})")
print(f"CHECK  sum(rel_f) = {dist['rel_f'].sum():.4f}  (must equal 1)")

# ── 3. GROUPED vs EXACT STATISTICS ──────────────────────────────────────────
print("\n===== 3. GROUPED vs EXACT =====")

mids = dist["midpoint"].to_numpy()
xbar_g = np.average(mids, weights=f)
s_g = np.sqrt((f * (mids - xbar_g) ** 2).sum() / (n - 1))
xbar, s = x.mean(), x.std(ddof=1)

print(f"grouped mean {xbar_g:7.3f}    exact mean {xbar:7.3f}    "
      f"error {100*abs(xbar_g-xbar)/xbar:.2f}%")
print(f"grouped sd   {s_g:7.3f}    exact sd   {s:7.3f}    "
      f"error {100*abs(s_g-s)/s:.2f}%")

# ── 4. SHAPE ────────────────────────────────────────────────────────────────
print("\n===== 4. SHAPE =====")

med = np.median(x)
q1, q3 = np.percentile(x, [25, 75])
print(f"mean   {xbar:.2f}")
print(f"median {med:.2f}")
print(f"Pearson SK = 3(mean - median)/s = {3*(xbar-med)/s:.3f}")
print(f"five-number summary: {np.percentile(x, [0,25,50,75,100]).round(2)}")
print(f"IQR: {q3-q1:.2f}")

if xbar > med * 1.02:
    shape = "RIGHT-SKEWED"
elif xbar < med * 0.98:
    shape = "LEFT-SKEWED"
else:
    shape = "roughly SYMMETRIC"
print(f"\nShape: {shape}")
print("=> report the "
      + ("MEAN and STANDARD DEVIATION" if shape == "roughly SYMMETRIC"
         else "MEDIAN and IQR"))

# ── 5. THE FOUR PLOTS ───────────────────────────────────────────────────────
fig, axes = plt.subplots(2, 2, figsize=(11, 8.5))

# histogram — bars touch by construction
axes[0, 0].hist(x, bins=edges, color="#8A5FBF", edgecolor="white")
axes[0, 0].set(title="Histogram", xlabel="Service time (minutes)", ylabel="Frequency")

# frequency polygon — MIDPOINTS, closed on the axis
poly_x = np.concatenate(([mids[0] - w], mids, [mids[-1] + w]))
poly_y = np.concatenate(([0], f, [0]))
axes[0, 1].plot(poly_x, poly_y, "o-", color="#5B2A86", lw=2)
axes[0, 1].axhline(0, color="grey", lw=0.8)
axes[0, 1].set(title="Frequency polygon", xlabel="Class midpoint", ylabel="Frequency")

# ogive — upper BOUNDARIES, starting at 0
axes[1, 0].plot(edges, np.concatenate(([0], f.cumsum())), "o-", color="#0FA3A3", lw=2)
axes[1, 0].axhline(n / 2, ls="--", color="grey")
axes[1, 0].axvline(med, ls="--", color="grey")
axes[1, 0].set(title="Ogive (cumulative frequency)",
               xlabel="Upper class boundary", ylabel="Cumulative frequency")

# boxplot
axes[1, 1].boxplot(x, vert=False, showmeans=True)
axes[1, 1].set(title="Boxplot", xlabel="Service time (minutes)", yticks=[])

plt.tight_layout()
plt.savefig("frequency_plots.png", dpi=110)
print("\nWrote frequency_plots.png")