kanshan
b0f7833560
- Extract shared fingerprint/GP/EI into platform/decision/bo_core.py - Replace BO stub with propose_next_screen_batch(method=gp_ei) - Add scripts/41 multi-round screen→gate→repro→scaleup→BO loop - Keep scripts/40 as rounds=1 backward-compat wrapper - Add scripts/42_ledger_query.py for hits/scaleup/fails/rounds - Unit tests for EI ranking; docs + README updated
186 lines
9.1 KiB
Python
186 lines
9.1 KiB
Python
#!/usr/bin/env python3
|
||
"""04_bayesian_opt.py — 高斯过程 + 期望改进 (EI) 模拟钙钛矿主动学习闭环
|
||
场景: 在 GGA 钙钛矿全池 (卤化物+氧化物) 中寻找带隙 = 1.4 eV 的候选 (光伏最优带隙附近)
|
||
目标: reward = -|band_gap - 1.4 eV| (越大越好, 0 表示精确命中)
|
||
实现: GP (Matern ν=2.5 + WhiteKernel) 拟合带隙回归 → EI 采集函数批量推荐
|
||
对比: Random Search 基线 (同等实验预算)
|
||
闭环类比: 每轮 12 个推荐 = 机器人单轮合成的 96 孔板中的一板; "实验" = DFT 计算结果
|
||
|
||
指纹 / GP+EI 核心已抽至 platform/decision/bo_core.py，供 Nature 多轮闭环复用。
|
||
"""
|
||
import json
|
||
import sys
|
||
import warnings
|
||
from pathlib import Path
|
||
|
||
import numpy as np
|
||
import pandas as pd
|
||
from sklearn.model_selection import GroupShuffleSplit
|
||
|
||
warnings.filterwarnings("ignore")
|
||
|
||
ROOT = Path(__file__).resolve().parents[1]
|
||
# Shared fingerprint / GP+EI — path-inject decision/ (avoid stdlib `platform` clash)
|
||
_DEC = ROOT / "platform" / "decision"
|
||
if str(_DEC) not in sys.path:
|
||
sys.path.insert(0, str(_DEC))
|
||
from bo_core import ( # noqa: E402
|
||
DEFAULT_TARGET_GAP,
|
||
expected_improvement,
|
||
fingerprint,
|
||
fit_gp_bandgap,
|
||
reward_from_gap,
|
||
)
|
||
|
||
AGG = ROOT / "data" / "nomad_perovskite_aggregated.csv"
|
||
TARGET_GAP = DEFAULT_TARGET_GAP # 目标带隙 (eV), 接近 S-Q 最优
|
||
INIT, BATCH, ROUNDS = 15, 12, 22
|
||
|
||
|
||
def main():
|
||
agg = pd.read_csv(AGG)
|
||
pool_df = agg[(agg.functional_class == "gga")
|
||
& (agg.band_gap_eV > 0.05) & (agg.band_gap_eV < 5.0)].copy()
|
||
print(f"GGA 钙钛矿全池 (0.05 < Eg < 5 eV): {len(pool_df)}")
|
||
|
||
X_all = np.stack([fingerprint(f) for f in pool_df.formula])
|
||
y_all = pool_df.band_gap_eV.values # GP 拟合目标: 带隙回归
|
||
r_all = reward_from_gap(y_all, TARGET_GAP) # reward (越大越好)
|
||
formulas = pool_df.formula.values
|
||
|
||
# 划分: 训练池 (可"合成/计算"的候选) vs 未来池 (从未观测, 用于理论上界)
|
||
gss = GroupShuffleSplit(n_splits=1, test_size=0.3, random_state=42)
|
||
pool_idx, future_idx = next(gss.split(X_all, y_all, groups=formulas))
|
||
X_pool, y_pool, r_pool = X_all[pool_idx], y_all[pool_idx], r_all[pool_idx]
|
||
formulas_pool = formulas[pool_idx]
|
||
print(f" 训练池 {len(X_pool)} | 未来未探索 {len(future_idx)}")
|
||
|
||
rng = np.random.RandomState(42)
|
||
init_idx = rng.choice(len(X_pool), INIT, replace=False)
|
||
known_mask = np.zeros(len(X_pool), dtype=bool); known_mask[init_idx] = True
|
||
known_x, known_y = X_pool[init_idx].copy(), y_pool[init_idx].copy()
|
||
|
||
bo_best_r = [float(np.max(reward_from_gap(known_y, TARGET_GAP)))]
|
||
rand_best_r = [bo_best_r[0]]
|
||
bo_best_bg = [float(known_y[np.argmax(reward_from_gap(known_y, TARGET_GAP))])]
|
||
rand_best_bg = [bo_best_bg[0]]
|
||
bo_pick_gap_traj = [] # 每轮 BO 推荐样本的带隙 (用于可视化探索→收敛)
|
||
|
||
for r in range(ROUNDS):
|
||
# ---------- BO: GP 带隙回归 + EI (shared bo_core) ----------
|
||
gp = fit_gp_bandgap(known_x, known_y, random_state=42)
|
||
cand_x = X_pool[~known_mask]
|
||
mu, sigma = gp.predict(cand_x, return_std=True)
|
||
best_r = float(np.max(reward_from_gap(known_y, TARGET_GAP)))
|
||
ei = expected_improvement(mu, sigma, best_r, target_gap=TARGET_GAP)
|
||
pick = np.argsort(ei)[-BATCH:]
|
||
new_y = y_pool[~known_mask][pick]
|
||
new_idx = np.where(~known_mask)[0][pick]
|
||
known_mask[new_idx] = True
|
||
known_x = np.vstack([known_x, X_pool[new_idx]])
|
||
known_y = np.concatenate([known_y, new_y])
|
||
bo_pick_gap_traj.append(new_y.tolist())
|
||
cur = reward_from_gap(known_y, TARGET_GAP)
|
||
bo_best_r.append(float(np.max(cur)))
|
||
bo_best_bg.append(float(known_y[np.argmax(cur)]))
|
||
|
||
# ---------- Random 基线: 等预算随机采样 ----------
|
||
avail = np.where(~known_mask)[0]
|
||
rand_pick = rng.choice(avail, min(BATCH, len(avail)), replace=False) if len(avail) else np.array([], int)
|
||
if len(rand_pick):
|
||
ry = r_pool[rand_pick]
|
||
if float(np.max(ry)) > rand_best_r[-1]:
|
||
rand_best_r.append(float(np.max(ry)))
|
||
rand_best_bg.append(float(y_pool[rand_pick[np.argmax(ry)]]))
|
||
else:
|
||
rand_best_r.append(rand_best_r[-1]); rand_best_bg.append(rand_best_bg[-1])
|
||
else:
|
||
rand_best_r.append(rand_best_r[-1]); rand_best_bg.append(rand_best_bg[-1])
|
||
|
||
bo_final_bg, rand_final_bg = bo_best_bg[-1], rand_best_bg[-1]
|
||
n_bo_hit = next((i for i, v in enumerate(bo_best_r) if v >= -0.01), None)
|
||
n_rand_hit = next((i for i, v in enumerate(rand_best_r) if v >= -0.01), None)
|
||
print(f"\nBO 最终带隙 {bo_final_bg:.3f} eV (|ΔEg| = {abs(bo_final_bg-TARGET_GAP):.4f})")
|
||
print(f"Rand 最终带隙 {rand_final_bg:.3f} eV (|ΔEg| = {abs(rand_final_bg-TARGET_GAP):.4f})")
|
||
if n_bo_hit: print(f"BO 第 {n_bo_hit} 轮即达 |ΔEg| ≤ 0.01 eV (共 {INIT + n_bo_hit*BATCH} 次实验)")
|
||
if n_rand_hit: print(f"Rand 第 {n_rand_hit} 轮才达 |ΔEg| ≤ 0.01 eV")
|
||
else: print("Random 在全部预算内未达 |ΔEg| ≤ 0.01 eV")
|
||
|
||
# ============ 图 7: BO 收敛曲线 ============
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
plt.rcParams["font.sans-serif"] = ["Noto Sans CJK JP", "DejaVu Sans"]
|
||
plt.rcParams["axes.unicode_minus"] = False
|
||
|
||
rounds = np.arange(len(bo_best_r))
|
||
fig, axes = plt.subplots(1, 2, figsize=(12.2, 4.6))
|
||
|
||
# 左图: 距目标带隙的偏差收敛 (log 轴)
|
||
ax = axes[0]
|
||
ax.semilogy(rounds, np.abs(np.array(bo_best_bg) - TARGET_GAP), "o-",
|
||
color="#7b4fd0", lw=2, ms=5, label="贝叶斯优化 (GP + EI)")
|
||
ax.semilogy(rounds, np.abs(np.array(rand_best_bg) - TARGET_GAP), "s--",
|
||
color="#888", lw=1.5, ms=4, label="随机搜索基线")
|
||
if n_bo_hit:
|
||
ax.axvline(n_bo_hit, color="#7b4fd0", ls=":", lw=1.2, alpha=0.7)
|
||
ax.annotate(f"第 {n_bo_hit} 轮命中\n(|ΔEg|≤0.01 eV)", xy=(n_bo_hit, abs(bo_best_bg[n_bo_hit]-TARGET_GAP)),
|
||
xytext=(n_bo_hit + 1.5, abs(bo_best_bg[n_bo_hit]-TARGET_GAP) * 8),
|
||
fontsize=9, color="#5a35a0",
|
||
arrowprops=dict(arrowstyle="->", color="#5a35a0", lw=1))
|
||
if n_rand_hit:
|
||
ax.annotate(f"随机: 第 {n_rand_hit} 轮\n(|ΔEg|={abs(rand_best_bg[n_rand_hit]-TARGET_GAP):.3f})",
|
||
xy=(n_rand_hit, abs(rand_best_bg[n_rand_hit]-TARGET_GAP)),
|
||
xytext=(n_rand_hit * 0.55, 0.002), fontsize=9, color="#555",
|
||
arrowprops=dict(arrowstyle="->", color="#555", lw=1))
|
||
ax.set_xlabel("闭环轮次 (每轮 12 个推荐 = 1 块 96 孔板)")
|
||
ax.set_ylabel("当前最优候选距目标带隙 |ΔEg| (eV)")
|
||
ax.set_title(f"主动学习收敛速度: 寻找带隙 = {TARGET_GAP} eV 候选\n(GGA 全池 {len(X_pool)} 个钙钛矿, GP 回归 + 期望改进)")
|
||
ax.legend(loc="upper right", fontsize=9)
|
||
ax.grid(alpha=0.3, which="both")
|
||
ax.text(0.02, 0.05, f"初始 {INIT} 个随机样本 → 共 {INIT + ROUNDS*BATCH} 次『实验』",
|
||
transform=ax.transAxes, fontsize=8.5,
|
||
bbox=dict(facecolor="white", alpha=0.7, edgecolor="#999"))
|
||
|
||
# 右图: BO 每轮推荐样本的带隙轨迹 (探索→收敛)
|
||
ax = axes[1]
|
||
for i, gaps in enumerate(bo_pick_gap_traj, start=1):
|
||
ax.scatter(np.full(len(gaps), i), gaps, s=18, color="#7b4fd0", alpha=0.45, zorder=3)
|
||
ax.axhline(TARGET_GAP, color="#d62728", ls="--", lw=1.4, label=f"目标带隙 {TARGET_GAP} eV")
|
||
ax.plot(rounds[1:], bo_best_bg[1:], color="#2ca02c", lw=1.8, marker=".", ms=6,
|
||
label="BO 当前最优候选", zorder=4)
|
||
ax.set_xlabel("闭环轮次")
|
||
ax.set_ylabel("推荐候选的 DFT 带隙 (eV)")
|
||
ax.set_title("BO 推荐轨迹: 前期大范围探索 → 后期在目标附近精修\n(紫色点 = 每轮 GP+EI 推荐的 12 个组分)")
|
||
ax.legend(loc="upper right", fontsize=9)
|
||
ax.grid(alpha=0.3)
|
||
fig.tight_layout()
|
||
fig.savefig(ROOT / "figures/fig7_bayesian_opt.png", bbox_inches="tight", dpi=150)
|
||
plt.close(fig)
|
||
|
||
# ============ 输出 JSON ============
|
||
out = {
|
||
"scenario": f"GGA 全池寻找带隙 = {TARGET_GAP} eV 的最佳材料",
|
||
"description": "BO 目标: 最小化 |带隙 - 1.4 eV|; 真实场景: 寻找『光伏最优带隙』候选",
|
||
"pool_size": int(len(X_pool)),
|
||
"future_size": int(len(future_idx)),
|
||
"init": INIT, "batch": BATCH, "rounds": ROUNDS,
|
||
"total_experiments": INIT + ROUNDS * BATCH,
|
||
"bo_final_reward": round(bo_best_r[-1], 4),
|
||
"rand_final_reward": round(rand_best_r[-1], 4),
|
||
"bo_final_bg": round(bo_final_bg, 3),
|
||
"rand_final_bg": round(rand_final_bg, 3),
|
||
"bo_round_to_001": n_bo_hit,
|
||
"rand_round_to_001": n_rand_hit,
|
||
"history_bo_reward": [round(v, 4) for v in bo_best_r],
|
||
"history_rand_reward": [round(v, 4) for v in rand_best_r],
|
||
"bo_core": "platform/decision/bo_core.py",
|
||
}
|
||
with open(ROOT / "output/bayesian_opt.json", "w") as f:
|
||
json.dump(out, f, ensure_ascii=False, indent=2)
|
||
print("\n已写入 output/bayesian_opt.json, 已重绘 figures/fig7_bayesian_opt.png")
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|