"""独立可复现的中文科研演示：小样本预测模型不能只看 AUC。
仅使用固定种子、CPU 轻量计算；所有预处理/筛选/调参均在训练折内。
"""
from __future__ import annotations
import base64, hashlib, json, platform, sys, time, warnings
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import font_manager
from sklearn.compose import ColumnTransformer
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (average_precision_score, brier_score_loss, confusion_matrix,
                             precision_recall_curve, roc_auc_score, roc_curve)
from sklearn.model_selection import GridSearchCV, StratifiedKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.base import clone
from sklearn.calibration import calibration_curve
import sklearn

warnings.filterwarnings("ignore")
SEED = 20260911
ROOT = Path(__file__).resolve().parent
OUT = ROOT
FONT_URL = "https://raw.githubusercontent.com/notofonts/noto-cjk/main/Sans/OTF/SimplifiedChinese/NotoSansCJKsc-Regular.otf"
FONT_PATH = ROOT / "NotoSansCJKsc-Regular.otf"

def choose_font():
    # 可复现获取说明：仅从 Google Noto 官方维护的 GitHub 仓库下载此文件，
    # 然后将其放在本脚本同目录：
    # python3 -c "import urllib.request; urllib.request.urlretrieve(\
    # 'https://raw.githubusercontent.com/notofonts/noto-cjk/main/Sans/OTF/SimplifiedChinese/NotoSansCJKsc-Regular.otf','NotoSansCJKsc-Regular.otf')"
    # 本脚本不静默回退到可能缺少中文 glyph 的字体；缺文件即清晰失败。
    if not FONT_PATH.is_file() or FONT_PATH.stat().st_size == 0:
        raise FileNotFoundError(
            f"缺少官方中文字体：{FONT_PATH}。请从 {FONT_URL} 下载为该文件后重试；"
            "禁止静默回退字体，以免生成方框乱码。")
    font_manager.fontManager.addfont(str(FONT_PATH))
    prop = font_manager.FontProperties(fname=str(FONT_PATH))
    if "Noto Sans CJK SC" not in prop.get_name():
        raise RuntimeError(f"字体文件名称校验失败：{prop.get_name()}，期望 Noto Sans CJK SC")
    plt.rcParams["axes.unicode_minus"] = False
    return prop

def simulate():
    rng = np.random.default_rng(SEED); n = 240
    size = rng.lognormal(3.0, .42, n).round(1)
    team = rng.integers(2, 13, n)
    prep = np.clip(rng.normal(67, 16, n), 15, 100).round(1)
    resource = np.clip(rng.normal(60, 18, n), 10, 100).round(1)
    field = rng.choice(["基础研究", "工程技术", "社会科学", "计算方法"], n, p=[.27,.28,.25,.20])
    field_effect = pd.Series(field).map({"基础研究":0.15,"工程技术":0.25,"社会科学":-0.20,"计算方法":0.10}).to_numpy()
    # 明确的虚构机制：logit 由准备、资源、团队边际作用及领域偏移组成。
    logit = -1.05 + .035*(prep-60) + .025*(resource-60) + .10*(team-6) + .16*np.log(size) + field_effect + rng.normal(0,.75,n)
    p = 1/(1+np.exp(-logit)); y = rng.binomial(1,p)
    df = pd.DataFrame({"project_id":[f"P{i:03d}" for i in range(1,n+1)], "project_size":size,
      "team_size":team, "preparation_score":prep, "research_field":field, "resource_input":resource,
      "redundant_size_copy":size+rng.normal(0,.12,n), "redundant_resource_copy":resource+rng.normal(0,.8,n),
      "completion":y})
    # 教学专用：直接由结果派生，法律模型绝不使用；保留在原始 CSV 以便审计。
    df["outcome_derived_leakage"] = y
    return df

def build_models(n_features):
    num = ["project_size","team_size","preparation_score","resource_input","redundant_size_copy","redundant_resource_copy"]
    cat = ["research_field"]
    prep = ColumnTransformer([("num", Pipeline([("imp",SimpleImputer(strategy="median")),("scale",StandardScaler())]), num),
                              ("cat", Pipeline([("imp",SimpleImputer(strategy="most_frequent")),("oh",OneHotEncoder(handle_unknown="ignore"))]), cat)])
    common = [("pre", prep), ("select", SelectKBest(f_classif, k="all"))]
    models = {
      "简单基线": (Pipeline([("model", DummyClassifier(strategy="prior"))]), {}),
      "正则化Logistic回归": (Pipeline(common + [("model",LogisticRegression(max_iter=3000, solver="liblinear", random_state=SEED))]),
        {"select__k":[5,8,"all"],"model__C":[.1,1,10]}),
      "随机森林": (Pipeline(common + [("model",RandomForestClassifier(n_estimators=160, random_state=SEED, n_jobs=1, class_weight="balanced"))]),
        {"select__k":[5,8,"all"],"model__max_depth":[2,4,None],"model__min_samples_leaf":[3,7]})}
    return models, num, cat

def evaluate(df):
    legal_cols = ["project_size","team_size","preparation_score","research_field","resource_input","redundant_size_copy","redundant_resource_copy"]
    X, y = df[legal_cols], df.completion.to_numpy()
    outer = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)
    inner = StratifiedKFold(n_splits=3, shuffle=True, random_state=SEED+1)
    models, _, _ = build_models(len(legal_cols)); preds=[]; rows=[]; importances=[]
    for name,(est,grid) in models.items():
        oof=np.zeros(len(df)); fold_best=[]
        for fold,(tr,te) in enumerate(outer.split(X,y),1):
            if grid:
                search=GridSearchCV(est,grid,cv=inner,scoring="roc_auc",n_jobs=1,refit=True)
                search.fit(X.iloc[tr],y[tr]); fitted=search.best_estimator_; best=str(search.best_params_)
            else:
                fitted=clone(est).fit(X.iloc[tr],y[tr]); best="固定先验概率"
            oof[te]=fitted.predict_proba(X.iloc[te])[:,1]
            fold_best.append(best)
            if name != "简单基线":
                names=fitted.named_steps["pre"].get_feature_names_out()
                mask=fitted.named_steps["select"].get_support()
                m=fitted.named_steps["model"]
                imp=np.abs(m.coef_[0]) if hasattr(m,"coef_") else m.feature_importances_
                for fn,v in zip(np.asarray(names)[mask], imp): importances.append({"model":name,"feature":fn,"importance":float(v),"fold":fold})
        ap=average_precision_score(y,oof)
        rows.append({"model":name,"roc_auc":roc_auc_score(y,oof),"pr_auc":ap,"brier":brier_score_loss(y,oof),
                     "positive_rate_baseline":y.mean(),"threshold":.5,"accuracy":((oof>=.5)==y).mean(),
                     "sensitivity":(((oof>=.5)&(y==1)).sum()/max(1,(y==1).sum())),"specificity":(((oof<.5)&(y==0)).sum()/max(1,(y==0).sum())),"best_params_by_fold":" | ".join(fold_best)})
        preds.extend([{"project_id":df.project_id.iloc[i],"model":name,"fold":int(next(f for f,(_,te) in enumerate(outer.split(X,y),1) if i in te)),"observed":int(y[i]),"predicted_probability":float(oof[i]),"predicted_class_0_5":int(oof[i]>=.5)} for i in range(len(df))])
    # 独立错误示范：把结果派生列输入，不能进入合法模型选择。
    leakX=df[["project_size","team_size","preparation_score","research_field","resource_input","outcome_derived_leakage"]]
    leak_pre=ColumnTransformer([("num",SimpleImputer(strategy="median"),["project_size","team_size","preparation_score","resource_input","outcome_derived_leakage"]),("cat",OneHotEncoder(handle_unknown="ignore"),["research_field"])])
    leak=Pipeline([("pre",leak_pre),("model",LogisticRegression(max_iter=1000,random_state=SEED))]).fit(leakX,y); lp=leak.predict_proba(leakX)[:,1]
    leakrow={"model":"错误示范（结果派生变量泄漏）","roc_auc":roc_auc_score(y,lp),"pr_auc":average_precision_score(y,lp),"brier":brier_score_loss(y,lp),"positive_rate_baseline":y.mean(),"threshold":.5,"accuracy":((lp>=.5)==y).mean(),"sensitivity":1.0,"specificity":1.0,"best_params_by_fold":"不参与合法选择；仅教学"}
    rows.append(leakrow)
    return pd.DataFrame(rows),pd.DataFrame(preds),pd.DataFrame(importances),leakX,lp

def plot_all(df, metrics, preds, lp, font):
    colors={"简单基线":"#9aa6b2","正则化Logistic回归":"#176b87","随机森林":"#e07a5f"}
    legal=preds.model.unique(); y=df.completion.to_numpy()
    fig,ax=plt.subplots(1,2,figsize=(13,5));
    for name in legal:
        p=preds[preds.model==name].sort_values("project_id").predicted_probability
        fpr,tpr,_=roc_curve(y,p); prec,rec,_=precision_recall_curve(y,p)
        ax[0].plot(fpr,tpr,label=f"{name} (AUC={roc_auc_score(y,p):.3f})",color=colors[name])
        ax[1].plot(rec,prec,label=f"{name} (AP={average_precision_score(y,p):.3f})",color=colors[name])
    ax[0].plot([0,1],[0,1],"--",color="#888",label="随机基线"); ax[1].axhline(y.mean(),ls="--",color="#888",label=f"阳性率基线={y.mean():.3f}")
    ax[0].set_xlabel("假阳性率",fontproperties=font); ax[0].set_ylabel("真阳性率",fontproperties=font); ax[0].set_title("ROC 曲线",fontproperties=font)
    ax[1].set_xlabel("召回率",fontproperties=font); ax[1].set_ylabel("精确率",fontproperties=font); ax[1].set_title("PR 曲线",fontproperties=font)
    for a in ax:a.legend(prop=font,fontsize=8);a.grid(alpha=.2)
    fig.tight_layout(); fig.savefig(OUT/"figure-01-roc-pr.png",dpi=180);plt.close(fig)
    fig,ax=plt.subplots(figsize=(7,5)); ax.plot([0,1],[0,1],"--",color="#888",label="完美校准参考")
    for name in legal:
        p=preds[preds.model==name].sort_values("project_id").predicted_probability; frac,mean=calibration_curve(y,p,n_bins=8,strategy="quantile")
        ax.plot(mean,frac,"o-",label=name,color=colors[name])
    ax.set_xlabel("平均预测概率",fontproperties=font); ax.set_ylabel("实际阳性比例",fontproperties=font); ax.set_title("折外概率校准图",fontproperties=font);ax.legend(prop=font);ax.grid(alpha=.2);fig.tight_layout();fig.savefig(OUT/"figure-02-calibration.png",dpi=180);plt.close(fig)
    fig,axs=plt.subplots(2,2,figsize=(11,8)); legalrows=metrics[metrics.model.isin(legal)]
    for j,(col,title) in enumerate([("roc_auc","ROC-AUC"),("pr_auc","PR-AUC"),("brier","Brier（越低越好）"),("accuracy","准确率（阈值0.5）")]):
        a=axs.flat[j]; a.bar(legalrows.model,legalrows[col],color=[colors[x] for x in legalrows.model]);a.set_title(title,fontproperties=font);a.tick_params(axis="x",rotation=18);[t.set_fontproperties(font) for t in a.get_xticklabels()];a.grid(axis="y",alpha=.2)
    fig.suptitle("合法模型指标比较（折外预测）",fontweight="bold",fontproperties=font);fig.tight_layout();fig.savefig(OUT/"figure-03-metrics-confusion.png",dpi=180);plt.close(fig)
    fig,ax=plt.subplots(1,2,figsize=(12,5)); normal=metrics[metrics.model.isin(legal)];
    for col,title in [("roc_auc","ROC-AUC"),("brier","Brier（越低越好）")]:
        vals=list(normal[col])+[metrics.loc[metrics.model.str.startswith("错误示范"),col].iloc[0]]; labels=list(normal.model)+["错误示范（泄漏）"]; a=ax[0 if col=="roc_auc" else 1];a.bar(labels,vals,color=[colors[x] for x in normal.model]+["#b2182b"]);a.set_title(title,fontproperties=font);a.tick_params(axis="x",rotation=20);[t.set_fontproperties(font) for t in a.get_xticklabels()];a.grid(axis="y",alpha=.2)
    fig.suptitle("合法方案与泄漏方案对照\n不正确做法，仅用于教学",color="#b2182b",fontweight="bold",fontproperties=font);fig.tight_layout();fig.savefig(OUT/"figure-04-leakage.png",dpi=180);plt.close(fig)

def write_reports(df,metrics,preds,imp,font,start):
    elapsed=time.perf_counter()-start; pos=int(df.completion.sum())
    metric_show=metrics[["model","roc_auc","pr_auc","brier","accuracy","sensitivity","specificity"]].round(4).to_string(index=False)
    cm=[]
    for _,r in metrics[metrics.model.isin(["简单基线","正则化Logistic回归","随机森林"])].iterrows():
        p=preds[preds.model==r.model].sort_values("project_id").predicted_probability; cm.append(f"- **{r.model}**：\n\n  ```text\n  {confusion_matrix(df.completion,(p>=.5).astype(int)).tolist()}\n  ```（行=真实类别0/1，列=预测类别0/1）")
    top=imp.groupby(["model","feature"],as_index=False).importance.mean().sort_values(["model","importance"],ascending=[True,False]).groupby("model").head(8).round(4).to_string(index=False)
    font_name = font.get_name()
    body=f'''# 小样本预测模型为什么不能只看 AUC

> **独立演示项目；数据完全虚构，不代表真实人群证据，也不用于临床或其他现实决策。**

## 1. 研究问题与模拟机制

面向科研研究生，模拟“研究项目是否按期完成”的二分类预测任务。固定种子 `{SEED}` 生成 240 个项目；变量包括项目规模、团队规模、前期准备、研究领域、资源投入和两项冗余特征。先生成特征，再按如下虚构机制生成概率：`logit(p) = -1.05 + 0.035*(准备-60) + 0.025*(资源-60) + 0.10*(团队-6) + 0.16*log(规模) + 领域偏移 + N(0,0.75)`，最后按 Bernoulli(p) 生成结果。该机制只是软件演示设定，不能解释为变量的因果作用。

原始数据还含有 `outcome_derived_leakage`，它直接等于结果，是**结果派生变量**，仅用于下文隔离的错误示范，未进入任何合法模型的选择或训练。

## 2. 验证设计与模型

- 外层 5 折、内层 3 折，均为分层、打乱、固定种子；保留每个样本的外层折外概率。
- 合法模型：简单先验基线、带标准化/独热编码/单变量筛选与 C 调参的正则化 Logistic 回归、带同样折内预处理/筛选与深度和叶节点调参的随机森林。
- `GridSearchCV` 只在外层训练折中进行；预处理、特征筛选和调参全部在训练折内。阈值预先固定为 0.5，未用折外结果选阈值。
- 报告 ROC-AUC、PR-AUC、Brier 分数；PR-AUC 的阳性率基线为 `{df.completion.mean():.4f}`（阳性数 {pos}/240）。折间 SD 没有被冒充为置信区间；这里报告的是一次固定模拟的折外汇总，未给出置信区间。

## 3. 结果

{metric_show}

图 1 展示 ROC/PR，图 2 展示校准，图 3 展示合法模型指标。AUC 衡量排序区分度，不衡量概率是否可信；因此必须同时看 PR、Brier、校准和阈值下混淆矩阵。预测概率只表示本模拟数据生成条件下的模型输出，不是个体真实概率。

### 阈值 0.5 混淆矩阵

{chr(10).join(cm)}

### 特征重要性限制

下表是外层折内拟合后系数绝对值/树特征重要性的平均摘要，仅用于描述模型内部利用的关联信号；冗余特征可能分摊重要性，重要性不等于稳定性，更不等于因果作用。

{top}

## 4. 明确隔离的错误示范

图 4 和指标表中的“错误示范（结果派生变量泄漏）”将 `outcome_derived_leakage` 作为输入，且在同一数据上展示，属于 **不正确做法，仅用于教学**；它不得与合法方案比较后用于模型选择。由于输入直接携带结果，表现会虚高，正是数据泄漏审计需要拦截的情形。

## 5. 复现、耗时与审计

- 实际验证命令：`python3 analysis.py`
- 实际执行耗时：`{elapsed:.2f}` 秒（CPU 轻量分析；模型线程数设为 1）。
- 首轮视觉验收失败：初版 PNG 的中文标题、坐标轴和图例显示为方框，不能发布。本轮已改用 Google Noto 官方字体 `{font_name}`，来源：{FONT_URL}；通过 `FontProperties(fname=...)` 显式指定所有标题、坐标轴、图例和注释，并设置负号正常显示。
- 字体验证：字体文件存在且非空；字体名称解析为 `{font_name}`；4 张 PNG 均重新渲染、尺寸和非空检查通过。该记录说明代码执行成功不等同于视觉无乱码，仍需人工查看图像。
- 依赖版本见 `dependencies.txt`；机器可读摘要见 `result-summary.json`。
- 输出：`simulated-project-data.csv`、`simulated-project-data.xlsx`、`oof-predictions.csv`、`metrics.csv`、`feature-importance.csv`、`data-dictionary.csv`、四张 PNG、`report.md`、`report.html`、`analysis.py`。
- 审计结论：结果派生列被明确排除在合法特征列表外；无外部数据、无密钥/环境变量读取、无大型模型、无系统服务安装；所有随机过程固定种子。

## 6. 解读边界

这是一次小样本、完全虚构的可复现实验。性能受生成机制、样本量和随机种子影响；不能外推到真实科研团队或人群。若用于真实研究，应预先定义目标、时间切点、缺失处理、外部验证和不确定性分析，并由领域专家审查。
'''
    (OUT/"report.md").write_text(body,encoding="utf-8")
    # 自包含 HTML，便于离线独立阅读
    html=body.replace("\n","<br>\n")
    for img in ["figure-01-roc-pr.png","figure-02-calibration.png","figure-03-metrics-confusion.png","figure-04-leakage.png"]:
        b=base64.b64encode((OUT/img).read_bytes()).decode(); html += f'<br><h3>{img}</h3><img style="max-width:100%" src="data:image/png;base64,{b}">'
    (OUT/"report.html").write_text('<!doctype html><meta charset="utf-8"><title>小样本预测演示</title><style>body{font-family:Arial,sans-serif;max-width:1100px;margin:2em auto;line-height:1.6}table{border-collapse:collapse}pre{background:#f5f5f5;padding:1em}</style><article>'+html+'</article>',encoding="utf-8")

def main():
    start=time.perf_counter(); font=choose_font(); df=simulate(); metrics,preds,imp,_,lp=evaluate(df)
    df.to_csv(OUT/"simulated-project-data.csv",index=False,encoding="utf-8-sig"); df.to_excel(OUT/"simulated-project-data.xlsx",index=False)
    preds.to_csv(OUT/"oof-predictions.csv",index=False,encoding="utf-8-sig");metrics.to_csv(OUT/"metrics.csv",index=False,encoding="utf-8-sig");imp.to_csv(OUT/"feature-importance.csv",index=False,encoding="utf-8-sig")
    pd.DataFrame([
      ["project_id","项目编号","标识符，不作为特征","字符串"],["project_size","项目规模","虚构连续数值","数值"],["team_size","团队规模","团队人数","整数"],["preparation_score","前期准备","0-100 虚构评分","数值"],["research_field","研究领域","四类虚构类别","分类"],["resource_input","资源投入","0-100 虚构评分","数值"],["redundant_size_copy","规模冗余特征","规模加小噪声","数值"],["redundant_resource_copy","资源冗余特征","资源加小噪声","数值"],["completion","按期完成","模拟二分类结果","0/1"],["outcome_derived_leakage","结果派生泄漏列","直接复制结果；仅错误示范","0/1"]],columns=["variable","中文名称","说明","类型"]).to_csv(OUT/"data-dictionary.csv",index=False,encoding="utf-8-sig")
    (OUT/"dependencies.txt").write_text(f"Python {sys.version.split()[0]}\nnumpy {np.__version__}\npandas {pd.__version__}\nscikit-learn {sklearn.__version__}\nmatplotlib {matplotlib.__version__}\nopenpyxl {__import__('openpyxl').__version__}\nfont {font}\n",encoding="utf-8")
    plot_all(df,metrics,preds,lp,font); write_reports(df,metrics,preds,imp,font,start)
    files=sorted([p for p in OUT.iterdir() if p.is_file() and p.name not in ["result-summary.json"]])
    image_checks={p.name:{"bytes":p.stat().st_size,"size_px":list(plt.imread(p).shape[:2][::-1]),"nonempty":p.stat().st_size>0} for p in sorted(OUT.glob("figure-*.png"))}
    summary={"seed":SEED,"n":240,"positive_rate":float(df.completion.mean()),"font":font.get_name(),"font_path":str(FONT_PATH),"font_source_url":FONT_URL,"font_fix_reason":"首轮视觉验收发现中文标题、坐标轴和图例为方框；本轮显式使用官方 Noto Sans CJK SC。","font_validation":{"exists":FONT_PATH.is_file(),"nonempty":FONT_PATH.stat().st_size>0,"resolved_name":font.get_name(),"no_silent_fallback":True,"manual_visual_review_required":True},"image_validation":image_checks,"elapsed_seconds":round(time.perf_counter()-start,3),"metrics":metrics.to_dict(orient="records"),"files":{p.name:hashlib.sha256(p.read_bytes()).hexdigest() for p in files}}
    (OUT/"result-summary.json").write_text(json.dumps(summary,ensure_ascii=False,indent=2),encoding="utf-8")
    print(json.dumps({"elapsed_seconds":summary["elapsed_seconds"],"positive_rate":summary["positive_rate"],"metrics":summary["metrics"]},ensure_ascii=False,indent=2))
if __name__=="__main__": main()