"""
Deeper analysis of the curated Barnes & Noble gift/gift-adjacent ocean-import sample.

Important interpretation rules
------------------------------
1. This dataset is a curated public-record sample, not a complete census of Barnes & Noble imports.
2. `source_country` follows the public trade/customs source and is not automatically proof of manufacturing origin.
3. Weight-based shares describe this sample only. They should not be generalized to B&N's entire supply chain.
4. Date clustering may reflect both actual shipment timing and the way public sources expose records.
"""

from pathlib import Path
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

ROOT = Path(__file__).resolve().parents[1]
DATA = ROOT / "data" / "bn_gift_import_sample.csv"
DATA_DIR = ROOT / "data"
FIGURES = ROOT / "figures"
FIGURES.mkdir(exist_ok=True)


def hhi_from_weights(weights: pd.Series) -> float:
    """HHI on 0-1 scale using weight shares."""
    weights = weights.dropna().astype(float)
    if weights.sum() == 0:
        return float("nan")
    shares = weights / weights.sum()
    return float((shares**2).sum())


def effective_number(hhi: float) -> float:
    """Inverse HHI; intuitive 'effective number' of equally sized groups."""
    return float(1 / hhi) if hhi and not np.isnan(hhi) else float("nan")


def gini(values: pd.Series) -> float:
    """Gini coefficient for non-negative shipment weights."""
    x = np.sort(values.dropna().astype(float).to_numpy())
    x = x[x >= 0]
    if len(x) == 0 or x.sum() == 0:
        return float("nan")
    n = len(x)
    ranks = np.arange(1, n + 1)
    return float((2 * np.sum(ranks * x) / (n * x.sum())) - ((n + 1) / n))


def product_count(description: str) -> int:
    """Approximate number of product types from semicolon-delimited manifest text."""
    if pd.isna(description) or not str(description).strip():
        return 0
    return str(description).count(";") + 1


def save_svg(name: str):
    plt.tight_layout()
    plt.savefig(FIGURES / f"{name}.svg", format="svg", bbox_inches="tight")
    plt.savefig(FIGURES / f"{name}.png", dpi=180, bbox_inches="tight")
    plt.close()


df = pd.read_csv(DATA, parse_dates=["arrival_date"])
df["product_count_est"] = df["description"].map(product_count)
df["weight_tons"] = df["weight_kg"] / 1000
known = df.dropna(subset=["weight_kg"]).copy()
known = known[known["weight_kg"] >= 0].copy()

total_weight = float(known["weight_kg"].sum())

# -----------------------------
# Country-level concentration
# -----------------------------
country = (
    known.groupby("source_country", as_index=False)
    .agg(shipments=("weight_kg", "size"), weight_kg=("weight_kg", "sum"))
    .sort_values("weight_kg", ascending=False)
)
country["weight_share"] = country["weight_kg"] / total_weight
country_hhi = hhi_from_weights(country["weight_kg"])
country["cumulative_weight_share"] = country["weight_share"].cumsum()
country.to_csv(DATA_DIR / "country_metrics.csv", index=False)

# -----------------------------
# Supplier-level concentration
# -----------------------------
supplier = (
    known.groupby("supplier_group", as_index=False)
    .agg(
        shipments=("weight_kg", "size"),
        weight_kg=("weight_kg", "sum"),
        categories=("category", "nunique"),
        approx_product_types=("product_count_est", "sum"),
    )
    .sort_values("weight_kg", ascending=False)
)
supplier["weight_share"] = supplier["weight_kg"] / total_weight
supplier["cumulative_weight_share"] = supplier["weight_share"].cumsum()
supplier_hhi = hhi_from_weights(supplier["weight_kg"])
supplier.to_csv(DATA_DIR / "supplier_metrics.csv", index=False)

# -----------------------------
# Category composition
# -----------------------------
category = (
    known.groupby("category", as_index=False)
    .agg(shipments=("weight_kg", "size"), weight_kg=("weight_kg", "sum"))
    .sort_values("weight_kg", ascending=False)
)
category["weight_share"] = category["weight_kg"] / total_weight
category.to_csv(DATA_DIR / "category_metrics.csv", index=False)

# -----------------------------
# Shipment-level diagnostics
# -----------------------------
shipment_metrics = df[
    [
        "arrival_date",
        "bol",
        "supplier_group",
        "source_country",
        "category",
        "mixed_container",
        "weight_kg",
        "packages",
        "product_count_est",
        "description",
    ]
].copy()
shipment_metrics["weight_share_of_known_sample"] = shipment_metrics["weight_kg"] / total_weight
shipment_metrics["kg_per_package"] = shipment_metrics["weight_kg"] / shipment_metrics["packages"]
shipment_metrics.to_csv(DATA_DIR / "shipment_metrics.csv", index=False)

# -----------------------------
# Robustness checks
# Remove largest k records by weight and recompute headline concentration metrics.
# -----------------------------
rows = []
ordered = known.sort_values("weight_kg", ascending=False).reset_index(drop=True)
for k in [0, 1, 3, 5]:
    sub = ordered.iloc[k:].copy()
    subtotal = float(sub["weight_kg"].sum())
    country_w = sub.groupby("source_country")["weight_kg"].sum()
    supplier_w = sub.groupby("supplier_group")["weight_kg"].sum()
    rows.append(
        {
            "largest_shipments_removed": k,
            "remaining_shipments": len(sub),
            "remaining_weight_kg": subtotal,
            "china_weight_share": float(country_w.get("China", 0) / subtotal),
            "form_pac_weight_share": float(supplier_w.get("FORM + PAC", 0) / subtotal),
            "country_hhi": hhi_from_weights(country_w),
            "supplier_hhi": hhi_from_weights(supplier_w),
            "mixed_weight_share": float(
                sub.loc[sub["mixed_container"], "weight_kg"].sum() / subtotal
            ),
        }
    )
robustness = pd.DataFrame(rows)
robustness.to_csv(DATA_DIR / "robustness_checks.csv", index=False)

# -----------------------------
# Descriptive comparisons: mixed vs single-SKU records
# -----------------------------
mixed_summary = (
    known.groupby("mixed_container")
    .agg(
        shipments=("weight_kg", "size"),
        total_weight_kg=("weight_kg", "sum"),
        mean_weight_kg=("weight_kg", "mean"),
        median_weight_kg=("weight_kg", "median"),
        median_product_types=("product_count_est", "median"),
        mean_product_types=("product_count_est", "mean"),
    )
    .reset_index()
)
mixed_summary["weight_share"] = mixed_summary["total_weight_kg"] / total_weight
mixed_summary.to_csv(DATA_DIR / "consolidation_metrics.csv", index=False)

single_known = known.loc[~known["mixed_container"]].sort_values("weight_kg", ascending=False)
single_without_largest = single_known.iloc[1:] if len(single_known) > 1 else single_known

# -----------------------------
# Time clustering (descriptive only)
# -----------------------------
monthly = (
    known.assign(month=known["arrival_date"].dt.to_period("M").astype(str))
    .groupby("month", as_index=False)
    .agg(shipments=("weight_kg", "size"), weight_kg=("weight_kg", "sum"))
)
monthly["weight_share"] = monthly["weight_kg"] / total_weight
monthly.to_csv(DATA_DIR / "monthly_metrics.csv", index=False)

jul_sep_2026 = (df["arrival_date"] >= "2026-07-01") & (df["arrival_date"] <= "2026-09-30")
jul_sep_2026_known = (known["arrival_date"] >= "2026-07-01") & (
    known["arrival_date"] <= "2026-09-30"
)

# -----------------------------
# Compact machine-readable summary
# -----------------------------
summary = {
    "sample_shipments": int(len(df)),
    "known_weight_shipments": int(len(known)),
    "known_weight_kg": total_weight,
    "known_weight_lb": float(total_weight * 2.2046226218),
    "shipment_weight_mean_kg": float(known["weight_kg"].mean()),
    "shipment_weight_median_kg": float(known["weight_kg"].median()),
    "shipment_weight_q1_kg": float(known["weight_kg"].quantile(0.25)),
    "shipment_weight_q3_kg": float(known["weight_kg"].quantile(0.75)),
    "largest_shipment_kg": float(known["weight_kg"].max()),
    "largest_shipment_weight_share": float(known["weight_kg"].max() / total_weight),
    "top_5_shipments_weight_share": float(
        known.nlargest(5, "weight_kg")["weight_kg"].sum() / total_weight
    ),
    "shipment_weight_gini": gini(known["weight_kg"]),
    "china_weight_share": float(country.set_index("source_country").loc["China", "weight_share"]),
    "country_hhi": country_hhi,
    "effective_source_countries": effective_number(country_hhi),
    "form_pac_weight_share": float(
        supplier.set_index("supplier_group").loc["FORM + PAC", "weight_share"]
    ),
    "top_3_suppliers_weight_share": float(supplier.head(3)["weight_share"].sum()),
    "supplier_hhi": supplier_hhi,
    "effective_suppliers": effective_number(supplier_hhi),
    "mixed_container_row_share": float(df["mixed_container"].mean()),
    "mixed_container_weight_share": float(
        known.loc[known["mixed_container"], "weight_kg"].sum() / total_weight
    ),
    "mixed_median_weight_kg": float(
        known.loc[known["mixed_container"], "weight_kg"].median()
    ),
    "single_median_weight_kg": float(
        known.loc[~known["mixed_container"], "weight_kg"].median()
    ),
    "single_mean_weight_excluding_largest_kg": float(
        single_without_largest["weight_kg"].mean()
    ),
    "mixed_mean_weight_kg": float(
        known.loc[known["mixed_container"], "weight_kg"].mean()
    ),
    "mixed_median_product_types": float(
        df.loc[df["mixed_container"], "product_count_est"].median()
    ),
    "single_median_product_types": float(
        df.loc[~df["mixed_container"], "product_count_est"].median()
    ),
    "max_approx_product_types_one_record": int(df["product_count_est"].max()),
    "jul_sep_2026_row_share": float(jul_sep_2026.mean()),
    "jul_sep_2026_weight_share": float(
        known.loc[jul_sep_2026_known, "weight_kg"].sum() / total_weight
    ),
    "us_port_known_rows": int(df["us_port"].notna().sum()),
    "foreign_port_known_rows": int(df["foreign_port"].notna().sum()),
    "hs_code_known_rows": int(df["hs_code"].notna().sum()),
}

(DATA_DIR / "analysis_summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")

# -----------------------------
# Figures
# -----------------------------
# 1. Source-country weight
plt.figure(figsize=(8.5, 5.2))
plot_country = country.head(6)
plt.bar(plot_country["source_country"], plot_country["weight_kg"] / 1000)
plt.ylabel("Known shipment weight (metric tons)")
plt.title("Recorded source-country concentration in the sample")
plt.xticks(rotation=25, ha="right")
save_svg("weight_by_source_country")

# 2. Supplier Pareto chart
plt.figure(figsize=(9.5, 5.7))
pareto = supplier.head(10).copy()
x = np.arange(len(pareto))
plt.bar(x, pareto["weight_kg"] / 1000)
plt.ylabel("Known weight (metric tons)")
plt.xticks(x, pareto["supplier_group"], rotation=40, ha="right")
plt.title("Supplier concentration: top suppliers by known sample weight")
save_svg("supplier_concentration")

# 3. Shipment weight distribution (sorted)
plt.figure(figsize=(9.0, 5.3))
ranked = known.sort_values("weight_kg", ascending=False).reset_index(drop=True)
plt.bar(np.arange(1, len(ranked) + 1), ranked["weight_kg"] / 1000)
plt.axhline(known["weight_kg"].median() / 1000, linestyle="--", label="Median")
plt.xlabel("Shipment rank (largest to smallest)")
plt.ylabel("Shipment weight (metric tons)")
plt.title("Shipment-size distribution reveals a strong upper tail")
plt.legend()
save_svg("shipment_weight_distribution")

# 4. Mixed vs single-SKU median weight
plt.figure(figsize=(7.4, 5.2))
labels = ["Mixed-SKU", "Single-SKU"]
medians = [
    known.loc[known["mixed_container"], "weight_kg"].median() / 1000,
    known.loc[~known["mixed_container"], "weight_kg"].median() / 1000,
]
plt.bar(labels, medians)
plt.ylabel("Median known shipment weight (metric tons)")
plt.title("Mixed-SKU records are typically larger in this sample")
save_svg("mixed_vs_single_median")

# 5. Robustness: concentration after removing largest records
plt.figure(figsize=(8.4, 5.2))
plt.plot(
    robustness["largest_shipments_removed"],
    robustness["china_weight_share"] * 100,
    marker="o",
    label="China-linked weight share",
)
plt.plot(
    robustness["largest_shipments_removed"],
    robustness["form_pac_weight_share"] * 100,
    marker="o",
    label="FORM + PAC weight share",
)
plt.xlabel("Largest shipments removed")
plt.ylabel("Share of remaining known weight (%)")
plt.title("Robustness check: are concentration results driven by outliers?")
plt.xticks(robustness["largest_shipments_removed"])
plt.legend()
save_svg("concentration_robustness")

# 6. Category composition
plt.figure(figsize=(9.4, 6.0))
cat = category.head(10).sort_values("weight_kg")
plt.barh(cat["category"], cat["weight_kg"] / 1000)
plt.xlabel("Known shipment weight (metric tons)")
plt.title("Gift-related sample composition by record-level category")
save_svg("category_weight")

# 7. Timeline
plt.figure(figsize=(9.2, 5.2))
timeline = known.sort_values("arrival_date")
plt.scatter(timeline["arrival_date"], timeline["weight_kg"] / 1000, s=45)
plt.ylabel("Shipment weight (metric tons)")
plt.title("Observed gift-related shipments over time")
plt.xticks(rotation=30, ha="right")
save_svg("shipment_timeline")

# Console output designed to be understandable when run locally.
print("\n=== B&N GIFT-SUPPLY SAMPLE: DEEPER ANALYSIS ===")
print(f"Rows: {len(df)} ({len(known)} with known weight)")
print(f"Known sample weight: {total_weight:,.0f} kg ({total_weight/1000:,.1f} metric tons)")
print(f"Median shipment: {known['weight_kg'].median():,.0f} kg")
print(f"Largest shipment share: {summary['largest_shipment_weight_share']:.1%}")
print(f"Top 5 shipment share: {summary['top_5_shipments_weight_share']:.1%}")
print(f"Shipment-weight Gini: {summary['shipment_weight_gini']:.3f}")
print(f"China-linked weight share: {summary['china_weight_share']:.1%}")
print(f"Country HHI: {country_hhi:.3f} (effective countries: {effective_number(country_hhi):.2f})")
print(f"FORM + PAC weight share: {summary['form_pac_weight_share']:.1%}")
print(f"Top 3 suppliers: {summary['top_3_suppliers_weight_share']:.1%} of known weight")
print(f"Supplier HHI: {supplier_hhi:.3f} (effective suppliers: {effective_number(supplier_hhi):.2f})")
print(f"Mixed-SKU records: {summary['mixed_container_row_share']:.1%} of rows")
print(f"Mixed-SKU weight: {summary['mixed_container_weight_share']:.1%} of known weight")
print(
    "Median shipment weight — mixed vs single: "
    f"{summary['mixed_median_weight_kg']:,.0f} vs {summary['single_median_weight_kg']:,.0f} kg"
)
print(
    "Robustness (remove largest shipment): China share = "
    f"{robustness.loc[robustness['largest_shipments_removed']==1, 'china_weight_share'].iloc[0]:.1%}; "
    "FORM + PAC share = "
    f"{robustness.loc[robustness['largest_shipments_removed']==1, 'form_pac_weight_share'].iloc[0]:.1%}"
)
print("\nDerived tables written to data/*.csv; figures written to figures/*.svg and *.png")