"""
04-03  Descriptive Statistics Dashboard
Chapters 03-01, 03-02, 04-01, 04-02

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

pd.set_option("display.width", 130)
pd.set_option("display.max_columns", 30)

df = pd.read_csv("data/employees.csv")
sal = df["salary"]
n = len(sal)


def describe_one(x):
    x = pd.Series(x)
    q1, q3 = x.quantile([0.25, 0.75])
    return pd.Series({
        "n":       x.size,
        "mean":    x.mean(),
        "median":  x.median(),
        "sd":      x.std(ddof=1),
        "var":     x.var(ddof=1),
        "cv_pct":  100 * x.std(ddof=1) / x.mean(),
        "min":     x.min(),
        "Q1":      q1,
        "Q3":      q3,
        "max":     x.max(),
        "IQR":     q3 - q1,
        "skew_SK": 3 * (x.mean() - x.median()) / x.std(ddof=1),
    })


# ── 1. OVERALL SUMMARY ──────────────────────────────────────────────────────
print("\n===== 1. SALARY — OVERALL =====")
print(describe_one(sal).round(3).to_string())
print(f"mode(s): {list(sal.mode())}")
print(f"trimmed mean (10% each end): {stats.trim_mean(sal, 0.1):.2f}")

print("\nPopulation vs sample spread (for contrast):")
print(f"  sample sd     {sal.std(ddof=1):.2f}   (divides by n-1)")
print(f"  population sd {sal.std(ddof=0):.2f}   (divides by N)")

# ── 2. BY DEPARTMENT ────────────────────────────────────────────────────────
print("\n===== 2. SALARY BY DEPARTMENT =====")
by_dept = df.groupby("department")["salary"].apply(describe_one).unstack()
print(by_dept.round(2).to_string())

cv = by_dept["cv_pct"]
print("\nMost / least variable department (by CV, not by sd):")
print(f"  most   {cv.idxmax()} ({cv.max():.1f}%)")
print(f"  least  {cv.idxmin()} ({cv.min():.1f}%)")

# ── 3. THE WEIGHTED-MEAN CHECK ──────────────────────────────────────────────
print("\n===== 3. WEIGHTED MEAN CHECK =====")
means, counts = by_dept["mean"], by_dept["n"]
print(f"overall mean           {sal.mean():.2f}")
print(f"weighted mean of means {np.average(means, weights=counts):.2f}   <- must match")
print(f"PLAIN average of means {means.mean():.2f}   <- only equal when every n is equal")

# ── 4. OUTLIERS, BOTH RULES ─────────────────────────────────────────────────
print("\n===== 4. OUTLIERS =====")

q1, q3 = sal.quantile([0.25, 0.75])
iqr = q3 - q1
lower, upper = q1 - 1.5 * iqr, q3 + 1.5 * iqr
ext_lo, ext_hi = q1 - 3 * iqr, q3 + 3 * iqr

df["z"] = stats.zscore(sal, ddof=1)
df["iqr_flag"] = (sal < lower) | (sal > upper)
df["z_flag"] = df["z"].abs() > 3

print(f"IQR fences:     {lower:,.0f}  to  {upper:,.0f}")
print(f"Extreme fences: {ext_lo:,.0f}  to  {ext_hi:,.0f}")
print(f"IQR-rule outliers: {df['iqr_flag'].sum()}")
print(f"z-rule  outliers: {df['z_flag'].sum()}")

flagged = df.loc[df["iqr_flag"] | df["z_flag"],
                 ["employee_id", "department", "salary", "z", "iqr_flag", "z_flag"]]
print(flagged.sort_values("salary", ascending=False).to_string(index=False)
      if len(flagged) else "(none)")

if df["iqr_flag"].sum() != df["z_flag"].sum():
    print("\nThe two rules DISAGREE. The extreme value inflates s, shrinking its own")
    print("z-score -- the masking effect. The resistant IQR rule is the one to trust.")

print("\nEffect of the extreme value on each statistic:")
comp = pd.DataFrame({
    "with_outlier":    describe_one(sal),
    "without_outlier": describe_one(sal[~df["iqr_flag"]]),
}).T[["mean", "median", "sd", "IQR"]]
print(comp.round(2).to_string())

# ── 5. MOST UNUSUAL EMPLOYEES ───────────────────────────────────────────────
print("\n===== 5. FIVE MOST UNUSUAL SALARIES (by |z|) =====")
print(df.reindex(df["z"].abs().sort_values(ascending=False).index)
        .head(5)[["employee_id", "department", "salary", "z"]]
        .round(3).to_string(index=False))

# ── 6. DASHBOARD ────────────────────────────────────────────────────────────
fig, axes = plt.subplots(2, 2, figsize=(12, 9))

axes[0, 0].hist(sal, bins=15, color="#8A5FBF", edgecolor="white")
axes[0, 0].axvline(sal.mean(), color="#0FA3A3", lw=2, label="mean")
axes[0, 0].axvline(sal.median(), color="#0B7A7A", lw=2, ls="--", label="median")
axes[0, 0].set(title="Salary distribution", xlabel="Salary")
axes[0, 0].legend()

groups = [g["salary"].values for _, g in df.groupby("department")]
labels = list(df.groupby("department").groups.keys())
axes[0, 1].boxplot(groups, tick_labels=labels)
axes[0, 1].set(title="Salary by department", ylabel="Salary")
axes[0, 1].tick_params(axis="x", rotation=30)

means_sorted = by_dept["mean"].sort_values()
axes[1, 0].bar(means_sorted.index, means_sorted.values, color="#5B2A86")
axes[1, 0].set(title="Mean salary by department", ylabel="Salary")
axes[1, 0].tick_params(axis="x", rotation=30)

axes[1, 1].scatter(df["years_service"], sal, color="#5B2A86")
m, b = np.polyfit(df["years_service"], sal, 1)
xs = np.linspace(df["years_service"].min(), df["years_service"].max(), 50)
axes[1, 1].plot(xs, m * xs + b, color="#0FA3A3", lw=2)
axes[1, 1].set(title="Salary vs years of service",
               xlabel="Years of service", ylabel="Salary")

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