#!/usr/bin/env python3
"""
ANSUR II 一手人体测量数据分析。

数据源：2012 US Army Anthropometric Survey (ANSUR II)，NSRDEC 发布。
  - 官方声明 cleared for UNLIMITED PUBLIC RELEASE（美国政府 public domain）
  - 个体级原始数据：4,082 男 + 1,986 女 = 6,068 人，93 项直接测量
  - 技术报告：NATICK/TR-15/007
  - 镜像：https://github.com/pinarismiguzel/ANSUR-II

⚠️ 单位口径（本项目已在此栽过 5 次，改脚本前先读）：
  ANSUR II 里凡带 "sitting" 的字段（elbowrestheight / eyeheightsitting）
  一律是【座面以上】的高度，不是从地面。
  从地面 = poplitealheight + <sitting 字段>。

用法：
  python scripts/analyze-ansur.py            # 下载(缓存) + 分析 + 写 data/ansur-stats.json
  python scripts/analyze-ansur.py --keep     # 同时保留原始 CSV 到 .cache/
"""

import csv
import datetime
import json
import os
import shutil
import statistics as st
import sys
import urllib.request

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
CACHE = os.path.join(ROOT, ".cache", "ansur")
OUT = os.path.join(ROOT, "data", "ansur-stats.json")

SOURCES = {
    "M": "https://raw.githubusercontent.com/pinarismiguzel/ANSUR-II/main/ANSUR%20II%20MALE%20Public.csv",
    "F": "https://raw.githubusercontent.com/pinarismiguzel/ANSUR-II/main/ANSUR%20II%20FEMALE%20Public.csv",
}


def fetch(sex):
    os.makedirs(CACHE, exist_ok=True)
    path = os.path.join(CACHE, f"{sex}.csv")
    if os.path.exists(path) and os.path.getsize(path) > 100_000:
        return path
    print(f"  下载 {sex} ...", file=sys.stderr)
    req = urllib.request.Request(SOURCES[sex], headers={"User-Agent": "yangmk-research/1.0"})
    with urllib.request.urlopen(req, timeout=180) as r, open(path, "wb") as f:
        f.write(r.read())
    return path


def load():
    rows = []
    for sex in ("M", "F"):
        # ANSUR II 导出为 latin-1，用 utf-8 会 UnicodeDecodeError
        with open(fetch(sex), encoding="latin-1") as fh:
            for r in csv.DictReader(fh):
                try:
                    rows.append({
                        "sex": sex,
                        "stature": int(r["stature"]),                  # 身高（赤脚）
                        "elbow_rest": int(r["elbowrestheight"]),      # 座面以上肘高
                        "popliteal": int(r["poplitealheight"]),       # 腘窝高 ≈ 座高
                        "eye_sit": int(r["eyeheightsitting"]),        # 座面以上眼高
                    })
                except (ValueError, KeyError, TypeError):
                    continue
    return rows


def pct(sorted_vals, p):
    return sorted_vals[min(int(p * len(sorted_vals)), len(sorted_vals) - 1)]


def histogram(values, bin_mm, start, end):
    """把一组 mm 值切成等宽 bin，返回 (counts, start, bin_mm)。
    页面用 SVG 画，所以 bin 数不能太多 —— 30 个左右刚好能看出形状。"""
    nbins = (end - start) // bin_mm
    counts = [0] * nbins
    for v in values:
        i = int((v - start) // bin_mm)
        if 0 <= i < nbins:
            counts[i] += 1
    return counts, start, bin_mm


def coverage_for(values, fixed_list, tol=20):
    """给定一组桌面高，算每个固定桌高下「太高 / 合适 / 太矮」的占比。"""
    n = len(values)
    out = {}
    for fixed in fixed_list:
        out[str(fixed)] = {
            "too_high_pct": round(100 * sum(1 for v in values if v < fixed - tol) / n),
            "fits_pm20mm_pct": round(100 * sum(1 for v in values if abs(v - fixed) <= tol) / n),
            "too_low_pct": round(100 * sum(1 for v in values if v > fixed + tol) / n),
        }
    return out


def write_derived_csv(rows, path):
    """导出派生数据集供他人复算。
    只保留 4 个测量值 + 性别，不含任何身份标识 —— 原始个体数据本身
    已是匿名测量，这里再收一层，避免再分发时的任何隐私顾虑。"""
    os.makedirs(os.path.dirname(path), exist_ok=True)
    with open(path, "w", newline="") as f:
        w = csv.writer(f)
        w.writerow([
            "# Derived from ANSUR II (2012), NATICK/TR-15/007. Public release.",
            "Units: mm. All heights measured from the floor, barefoot.",
        ])
        w.writerow(["sex", "stature_mm", "desk_target_mm", "chair_target_mm", "eye_height_mm"])
        for r in rows:
            w.writerow([r["sex"], r["stature"], r["desk"], r["chair"], r["eye"]])


def main():
    rows = load()
    n = len(rows)
    print(f"N = {n}")

    # 全部换算成「从地面」
    for r in rows:
        r["desk"] = r["popliteal"] + r["elbow_rest"]     # 桌面高目标
        r["chair"] = r["popliteal"]                      # 座高目标
        r["eye"] = r["popliteal"] + r["eye_sit"]         # 坐眼高

    out = {
        "n": n,
        "source": "ANSUR II (2012), NATICK/TR-15/007",
        "unit_mm": True,
        "generated": datetime.date.today().isoformat(),
    }

    # 1) 三系数验证：真实中位比值 vs 站点当前使用值
    coeff = {}
    for key, ours, label in [
        ("chair", 0.25, "chair height (popliteal)"),
        ("desk", 0.40, "desk height (popliteal + elbow rest)"),
        ("eye", 0.70, "sitting eye height"),
    ]:
        v = sorted(r[key] / r["stature"] for r in rows)
        coeff[key] = {
            "label": label,
            "site_uses": ours,
            "true_median": round(st.median(v), 4),
            "p5": round(pct(v, 0.05), 4),
            "p95": round(pct(v, 0.95), 4),
            "delta_mm_at_175cm": round((st.median(v) - ours) * 1750),
        }
    out["coefficients"] = coeff

    # 2) 0.40 x 身高 的预测误差分布
    err = sorted(0.40 * r["stature"] - r["desk"] for r in rows)
    out["formula_040_error"] = {
        "median_mm": round(st.median(err)),
        "mean_mm": round(st.mean(err)),
        "sd_mm": round(st.pstdev(err)),
        "p5_mm": round(pct(err, 0.05)),
        "p95_mm": round(pct(err, 0.95)),
        "within_10mm_pct": round(100 * sum(1 for e in err if abs(e) <= 10) / n, 1),
        "within_20mm_pct": round(100 * sum(1 for e in err if abs(e) <= 20) / n, 1),
        "within_30mm_pct": round(100 * sum(1 for e in err if abs(e) <= 30) / n, 1),
        "within_50mm_pct": round(100 * sum(1 for e in err if abs(e) <= 50) / n, 1),
    }

    # 2b) 误差直方图（副图用：0.40 规则到底偏多少）
    err_hist, err_start, err_bin = histogram(err, 10, -60, 130)
    out["error_histogram"] = {
        "bin_mm": err_bin, "start_mm": err_start,
        "counts": err_hist,
        "median_mm": round(st.median(err)),
    }

    # 3) 身高能解释多少桌高差异
    xs = [r["stature"] for r in rows]
    ys = [r["desk"] for r in rows]
    mx, my = st.mean(xs), st.mean(ys)
    corr = (sum((a - mx) * (b - my) for a, b in zip(xs, ys)) / n) / (st.pstdev(xs) * st.pstdev(ys))
    out["stature_vs_desk_r2"] = round(corr ** 2, 3)

    # 4) 真实桌高分布
    d = sorted(r["desk"] for r in rows)
    out["desk_distribution"] = {
        "p5": pct(d, 0.05), "p25": pct(d, 0.25),
        "median": round(st.median(d)),
        "p75": pct(d, 0.75), "p95": pct(d, 0.95),
    }

    # 5) 主直方图（SVG 分布图的数据源）
    hist, hist_start, hist_bin = histogram(d, 10, 500, 820)
    out["desk_histogram"] = {"bin_mm": hist_bin, "start_mm": hist_start, "counts": hist}

    # 6) 市面固定桌高覆盖率
    #    ⚠️ 两个版本都要：赤脚值会被挑刺（"你没算鞋"），
    #    穿鞋 +25mm 才是能对外站得住的版本。页面以穿鞋版为主。
    SHOD_MM = 25
    shod = sorted(v + SHOD_MM for v in d)
    out["shoe_correction_mm"] = SHOD_MM
    out["desk_distribution_shod"] = {
        "p5": pct(shod, 0.05), "p25": pct(shod, 0.25),
        "median": round(st.median(shod)),
        "p75": pct(shod, 0.75), "p95": pct(shod, 0.95),
    }
    out["fixed_desk_coverage_barefoot"] = coverage_for(d, (710, 730, 740, 750))
    out["fixed_desk_coverage_shod"] = coverage_for(shod, (710, 730, 740, 750))
    coverage = out["fixed_desk_coverage_barefoot"]

    # 7) 分身高段表（页面图表用）
    bands = []
    for lo in range(1500, 1951, 50):
        sub = [r for r in rows if lo <= r["stature"] < lo + 50]
        if len(sub) < 30:
            continue
        t = sorted(r["desk"] for r in sub)
        c = sorted(r["chair"] for r in sub)
        bands.append({
            "band_cm": f"{lo // 10}-{lo // 10 + 5}",
            "n": len(sub),
            "desk_p5": pct(t, 0.05), "desk_median": round(st.median(t)), "desk_p95": pct(t, 0.95),
            "chair_median": round(st.median(c)),
            "formula_040": round(0.40 * st.median([r["stature"] for r in sub])),
            "formula_error_mm": round(0.40 * st.median([r["stature"] for r in sub]) - st.median(t)),
        })
    out["height_bands"] = bands

    # 8) 性别差异
    by_sex = {}
    for sex in ("M", "F"):
        s = [r for r in rows if r["sex"] == sex]
        by_sex[sex] = {
            "n": len(s),
            "stature_median": round(st.median([r["stature"] for r in s])),
            "chair_coef": round(st.median([r["chair"] / r["stature"] for r in s]), 4),
            "desk_coef": round(st.median([r["desk"] / r["stature"] for r in s]), 4),
            "eye_coef": round(st.median([r["eye"] / r["stature"] for r in s]), 4),
        }
    out["by_sex"] = by_sex

    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    with open(OUT, "w") as f:
        json.dump(out, f, indent=2)
    print(f"\n写入 {OUT}\n")

    # 打印摘要
    print("=== 系数真实值 vs 站点用值 ===")
    for k, v in coeff.items():
        print(f"  {k:6s} 站点={v['site_uses']:.3f}  真实中位={v['true_median']:.4f}  "
              f"175cm处偏差={v['delta_mm_at_175cm']:+d}mm")
    e = out["formula_040_error"]
    print(f"\n=== 0.40x身高 误差：中位 {e['median_mm']:+d}mm，"
          f"仅 {e['within_30mm_pct']}% 的人在 ±30mm 内 ===")
    print(f"=== 身高只能解释 {out['stature_vs_desk_r2']*100:.0f}% 的桌高差异 ===")
    dd = out["desk_distribution"]
    ds = out["desk_distribution_shod"]
    print(f"=== 真实桌高(赤脚)：P5={dd['p5']} 中位={dd['median']} P95={dd['p95']} ===")
    print(f"=== 穿鞋 +{out['shoe_correction_mm']}mm：P5={ds['p5']} 中位={ds['median']} P95={ds['p95']} ===")
    print("\n=== 固定桌覆盖率（赤脚 / 穿鞋）===")
    for fixed in ("710", "730", "740", "750"):
        b = coverage[fixed]
        s = out["fixed_desk_coverage_shod"][fixed]
        print(f"    {fixed}mm：太高 {b['too_high_pct']}% / {s['too_high_pct']}%　"
              f"合适 {b['fits_pm20mm_pct']}% / {s['fits_pm20mm_pct']}%")

    # 派生 CSV：一份留在 data/（页面构建时读），一份放进 public/ 供下载。
    # 页面说「数据可复算」就得真的给得出数据，否则又是空口。
    derived = os.path.join(ROOT, "data", "ansur-derived.csv")
    write_derived_csv(rows, derived)
    pub = os.path.join(ROOT, "public", "data", "ansur-derived.csv")
    os.makedirs(os.path.dirname(pub), exist_ok=True)
    shutil.copyfile(derived, pub)
    print(f"\n派生数据集：{derived} → {pub} ({os.path.getsize(pub) // 1024} KB)")

    # 脚本自己也放进 public/，页面上的「下载分析脚本」链接才不会和代码脱节。
    # 每次跑分析自动同步，避免手改脚本后忘了复制。
    script_pub = os.path.join(ROOT, "public", "data", "analyze-ansur.py")
    shutil.copyfile(os.path.abspath(__file__), script_pub)
    print(f"分析脚本：{script_pub} ({os.path.getsize(script_pub) // 1024} KB)")


if __name__ == "__main__":
    main()