from __future__ import annotations

import csv
import json
import math
from collections import Counter
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

ROOT = Path(__file__).resolve().parent
FIGURES = ROOT / "figures"
TABLES = ROOT / "tables"
FIGURES.mkdir(exist_ok=True)
TABLES.mkdir(exist_ok=True)

FEATURES = [
    ("direct_web", "自有網站"), ("own_ecommerce", "自有電商"),
    ("third_party_channel", "第三方通路"), ("delivery_ordering", "遠距訂購配送"),
    ("flexible_bulk", "散裝彈性秤售"), ("quality_assurance", "品質保證揭露"),
    ("origin_variety", "產地／品種資訊"), ("b2b_service", "B2B／批發服務"),
    ("heritage_story", "歷史／傳承敘事"), ("social_channel", "社群管道"),
    ("experience_education", "體驗／食農教育"), ("visible_price", "公開價格"),
]

plt.rcParams["font.family"] = ["Microsoft JhengHei", "Noto Sans CJK TC", "Arial Unicode MS", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False
plt.rcParams["figure.dpi"] = 160


def read_csv(path: Path):
    with path.open(encoding="utf-8-sig", newline="") as handle:
        return list(csv.DictReader(handle))


def wilson(success: int, total: int, z: float = 1.96):
    if not total:
        return [None, None]
    p = success / total
    den = 1 + z * z / total
    centre = (p + z * z / (2 * total)) / den
    half = z * math.sqrt(p * (1 - p) / total + z * z / (4 * total * total)) / den
    return [max(0.0, centre - half), min(1.0, centre + half)]


def fisher_two_sided(a: int, b: int, c: int, d: int):
    row1, col1, total = a + b, a + c, a + b + c + d
    def prob(x):
        return math.comb(col1, x) * math.comb(total - col1, row1 - x) / math.comb(total, row1)
    lo, hi = max(0, row1 - (total - col1)), min(row1, col1)
    observed = prob(a)
    return min(1.0, sum(prob(x) for x in range(lo, hi + 1) if prob(x) <= observed + 1e-12))


def analyse_subset(rows, label):
    groups = {g: [r for r in rows if r["group"] == g] for g in ("traditional", "integrated")}
    out = {"label": label, "n": len(rows), "n_traditional": len(groups["traditional"]),
           "n_integrated": len(groups["integrated"]), "features": []}
    for key, zh in FEATURES:
        t, i = groups["traditional"], groups["integrated"]
        ty, iy = sum(r[key] for r in t), sum(r[key] for r in i)
        item = {"feature": key, "label": zh,
                "traditional": {"yes": ty, "n": len(t), "rate": ty / len(t) if t else None,
                                "ci95": wilson(ty, len(t))},
                "integrated": {"yes": iy, "n": len(i), "rate": iy / len(i) if i else None,
                               "ci95": wilson(iy, len(i))}}
        if t and i:
            item["rate_difference_integrated_minus_traditional"] = iy / len(i) - ty / len(t)
            item["fisher_p"] = fisher_two_sided(iy, len(i) - iy, ty, len(t) - ty)
        out["features"].append(item)
    for r in rows:
        r["portfolio_score"] = sum(r[key] for key, _ in FEATURES)
    ts = [r["portfolio_score"] for r in groups["traditional"]]
    ins = [r["portfolio_score"] for r in groups["integrated"]]
    out["portfolio_score"] = {
        "traditional_mean": float(np.mean(ts)) if ts else None,
        "traditional_median": float(np.median(ts)) if ts else None,
        "integrated_mean": float(np.mean(ins)) if ins else None,
        "integrated_median": float(np.median(ins)) if ins else None,
        "mean_difference_integrated_minus_traditional": float(np.mean(ins) - np.mean(ts)) if ts and ins else None,
    }
    return out


def permutation_and_bootstrap(rows, reps=100_000, seed=20260711):
    rng = np.random.default_rng(seed)
    values = np.array([r["portfolio_score"] for r in rows], dtype=float)
    labels = np.array([r["group"] == "integrated" for r in rows])
    n_i = int(labels.sum())
    observed = values[labels].mean() - values[~labels].mean()
    extreme = 0
    for _ in range(reps):
        perm = rng.permutation(len(values))[:n_i]
        chosen = np.zeros(len(values), dtype=bool)
        chosen[perm] = True
        diff = values[chosen].mean() - values[~chosen].mean()
        extreme += abs(diff) >= abs(observed) - 1e-12
    t = values[~labels]
    i = values[labels]
    boot = np.empty(reps)
    for k in range(reps):
        boot[k] = rng.choice(i, len(i), replace=True).mean() - rng.choice(t, len(t), replace=True).mean()
    return {"repetitions": reps, "seed": seed, "two_sided_p": (extreme + 1) / (reps + 1),
            "bootstrap_ci95": [float(x) for x in np.quantile(boot, [0.025, 0.975])]}


rows = read_csv(ROOT / "businesses-60.csv")
qual = read_csv(ROOT / "qualitative-excerpts.csv")
for row in rows:
    for key, _ in FEATURES:
        row[key] = int(row[key])
        if row[key] not in (0, 1):
            raise ValueError(f"{row['id']} {key} is not binary")
    if row["group"] not in ("traditional", "integrated"):
        raise ValueError(f"invalid group: {row['id']}")

if len(rows) != 60 or len({r["id"] for r in rows}) != 60:
    raise ValueError("The expanded dataset must contain 60 unique cases")

main = analyse_subset(rows, "全部60案")
permutation = permutation_and_bootstrap(rows)
sensitivities = [
    analyse_subset([r for r in rows if r["source_grade"] in ("A", "B")], "排除C級來源"),
    analyse_subset([r for r in rows if int(r["evidence_year"]) >= 2020], "僅保留2020年後證據"),
]
theme_counts = Counter(r["theme"] for r in qual)
results = {
    "design": {"type": "stratified purposive mixed content analysis", "extraction_date": "2026-07-11",
               "n_total": len(rows), "n_traditional": 30, "n_integrated": 30,
               "unit": "business or brand case",
               "estimand": "public disclosure rate, not underlying capability prevalence"},
    "source_quality": {"grades": dict(Counter(r["source_grade"] for r in rows)),
                       "evidence_year_range": [min(int(r["evidence_year"]) for r in rows),
                                               max(int(r["evidence_year"]) for r in rows)]},
    "features": main["features"], "portfolio_score": {**main["portfolio_score"], **permutation, "maximum": 12},
    "sensitivity_analyses": sensitivities, "qualitative_themes": dict(theme_counts),
    "credibility": {"coder_count": 1, "review": "two-pass single-researcher audit",
                    "intercoder_reliability": "not estimable; no independent second coder"},
}
(ROOT / "analysis-results.json").write_text(json.dumps(results, ensure_ascii=False, indent=2), encoding="utf-8")

with (TABLES / "feature-comparison.csv").open("w", encoding="utf-8-sig", newline="") as handle:
    w = csv.writer(handle)
    w.writerow(["feature", "label", "traditional_yes", "traditional_n", "traditional_rate",
                "integrated_yes", "integrated_n", "integrated_rate", "rate_difference", "fisher_p"])
    for x in main["features"]:
        w.writerow([x["feature"], x["label"], x["traditional"]["yes"], 30,
                    f'{x["traditional"]["rate"]:.4f}', x["integrated"]["yes"], 30,
                    f'{x["integrated"]["rate"]:.4f}', f'{x["rate_difference_integrated_minus_traditional"]:.4f}',
                    f'{x["fisher_p"]:.6f}'])

with (ROOT / "coding-audit.csv").open("w", encoding="utf-8-sig", newline="") as handle:
    w = csv.writer(handle)
    w.writerow(["id", "identity_checked", "group_rule_checked", "binary_range_checked", "two_sources_listed",
                "source_grade", "evidence_year", "review_status", "notes"])
    for r in rows:
        two = bool(r["primary_url"] and r["secondary_url"])
        old = int(r["evidence_year"]) < 2020
        notes = "; ".join(x for x in ["較早期證據" if old else "", "納入敏感度排除組" if r["source_grade"] == "C" else ""] if x)
        w.writerow([r["id"], 1, 1, 1, int(two), r["source_grade"], r["evidence_year"], "pass_with_flag" if notes else "pass", notes])

with (ROOT / "source-log-60.csv").open("w", encoding="utf-8-sig", newline="") as handle:
    w = csv.writer(handle)
    w.writerow(["id", "operator", "source_role", "url", "source_grade", "evidence_year", "accessed_at"])
    for r in rows:
        w.writerow([r["id"], r["operator"], "primary", r["primary_url"], r["source_grade"], r["evidence_year"], "2026-07-11"])
        w.writerow([r["id"], r["operator"], "secondary", r["secondary_url"], r["source_grade"], r["evidence_year"], "2026-07-11"])

md = ["# 60案量化內容分析結果", "", "本研究比較30家傳統／地方型米業者與30家品牌／產銷整合型業者。0代表指定來源未揭露，不代表能力不存在。", "",
      "|指標|傳統型|整合型|差（百分點）|Fisher p|", "|---|---:|---:|---:|---:|"]
for x in main["features"]:
    md.append(f'|{x["label"]}|{x["traditional"]["yes"]}/30 ({x["traditional"]["rate"]:.1%})|{x["integrated"]["yes"]}/30 ({x["integrated"]["rate"]:.1%})|{100*x["rate_difference_integrated_minus_traditional"]:+.1f}|{x["fisher_p"]:.4f}|')
p = results["portfolio_score"]
md += ["", "## 能力組合分數", "", f'- 傳統型平均 {p["traditional_mean"]:.2f}；整合型平均 {p["integrated_mean"]:.2f}。',
       f'- 平均差 {p["mean_difference_integrated_minus_traditional"]:.2f}；100,000次置換檢定 p={p["two_sided_p"]:.5f}；bootstrap 95% CI {p["bootstrap_ci95"][0]:.2f}–{p["bootstrap_ci95"][1]:.2f}。',
       "", "## 可信度界線", "", "本版由單一研究者以兩階段規則覆核完成，未宣稱第二位獨立編碼者或編碼者間信度。來源等級與年份敏感度分析另存於 JSON。"]
(ROOT / "analysis-results.md").write_text("\n".join(md), encoding="utf-8")

# Figure 1: group rates.
labels = [x["label"] for x in main["features"]]
t_rates = np.array([x["traditional"]["rate"] for x in main["features"]]) * 100
i_rates = np.array([x["integrated"]["rate"] for x in main["features"]]) * 100
y = np.arange(len(labels))
fig, ax = plt.subplots(figsize=(10.8, 7.4))
ax.barh(y + .18, t_rates, .34, label="傳統／地方型（n=30）", color="#2f6f62")
ax.barh(y - .18, i_rates, .34, label="品牌／整合型（n=30）", color="#c78b35")
ax.set(yticks=y, yticklabels=labels, xlim=(0, 105), xlabel="指定來源的揭露率（%）", title="60案例的公開商業能力揭露率")
ax.grid(axis="x", alpha=.2); ax.legend(loc="lower right"); ax.invert_yaxis(); fig.tight_layout()
fig.savefig(FIGURES / "figure-1-disclosure-rates.png", bbox_inches="tight"); plt.close(fig)

# Figure 2: 60 x 12 heatmap.
matrix = np.array([[r[k] for k, _ in FEATURES] for r in rows])
fig, ax = plt.subplots(figsize=(12.5, 15.5))
ax.imshow(matrix, cmap=plt.matplotlib.colors.ListedColormap(["#ece8df", "#2f6f62"]), aspect="auto", vmin=0, vmax=1)
ax.set_xticks(range(12), labels, rotation=48, ha="right")
ax.set_yticks(range(60), [f'{r["id"]} {r["brand"]}' for r in rows], fontsize=7.5)
ax.axhline(29.5, color="#111", linewidth=2); ax.set_title("60案例—12項公開能力揭露矩陣"); fig.tight_layout()
fig.savefig(FIGURES / "figure-2-capability-heatmap.png", bbox_inches="tight"); plt.close(fig)

# Figure 3: portfolio score distribution.
fig, ax = plt.subplots(figsize=(9.2, 5.7))
for x, (group, color, name) in enumerate([("traditional", "#2f6f62", "傳統／地方型"), ("integrated", "#c78b35", "品牌／整合型")]):
    vals = [r["portfolio_score"] for r in rows if r["group"] == group]
    jitter = np.random.default_rng(100+x).uniform(-.16, .16, len(vals))
    ax.scatter(np.full(len(vals), x) + jitter, vals, s=48, color=color, alpha=.82, edgecolor="white")
    ax.hlines(np.mean(vals), x-.28, x+.28, color="#111", linewidth=3)
ax.set(xticks=[0,1], xticklabels=["傳統／地方型", "品牌／整合型"], ylim=(-.5,12.5), ylabel="公開能力組合分數（0–12）", title="60案例的公開能力組合分布")
ax.grid(axis="y", alpha=.2); fig.tight_layout(); fig.savefig(FIGURES / "figure-3-portfolio-score.png", bbox_inches="tight"); plt.close(fig)

# Figure 4: qualitative themes.
themes = sorted(theme_counts, key=theme_counts.get)
fig, ax = plt.subplots(figsize=(9.2, 5.7)); ax.barh(themes, [theme_counts[t] for t in themes], color="#6f8e87")
ax.set(xlabel="納入的短文本數", title="少量質性文本分析的主題分布"); ax.grid(axis="x", alpha=.2); fig.tight_layout()
fig.savefig(FIGURES / "figure-4-qualitative-themes.png", bbox_inches="tight"); plt.close(fig)

# Figure 7: source-quality sensitivity.
scenarios = [main] + sensitivities
names = [s["label"] for s in scenarios]
diffs = [s["portfolio_score"]["mean_difference_integrated_minus_traditional"] for s in scenarios]
ns = [s["n"] for s in scenarios]
fig, ax = plt.subplots(figsize=(9.2, 4.8)); bars=ax.barh(names, diffs, color=["#2f6f62", "#829d96", "#aab9b4"])
for b, d, n in zip(bars, diffs, ns): ax.text(d+.08, b.get_y()+b.get_height()/2, f"{d:.2f}（n={n}）", va="center")
ax.set(xlabel="整合型－傳統型的平均分數差", title="來源品質與證據年份敏感度分析"); ax.grid(axis="x", alpha=.2); fig.tight_layout()
fig.savefig(FIGURES / "figure-7-sensitivity-analysis.png", bbox_inches="tight"); plt.close(fig)

print(json.dumps({"n": len(rows), "portfolio": results["portfolio_score"], "grades": results["source_quality"]["grades"]}, ensure_ascii=False, indent=2))