kanshan
63b72b1880
- Point all scripts at repo data/figures/output/magpie via Path(__file__) - Skip 01_process when raw JSONL absent but cleaned CSVs exist - Vendor Magpie elemental .table files; harden 03 loader (54 feats) - Tighten oxide: sole anion must be O (exclude mixed chalcogenide/pnictide; Si treated as B-site / oxide-perovskite-like, not anion) - Regenerate figures, ML/BO metrics, report.html end-to-end
366 lines
15 KiB
Python
366 lines
15 KiB
Python
#!/usr/bin/env python3
|
|
"""35_embed_visualize.py — 把钙钛矿化学空间投影到 2D 验证 denovo 覆盖度
|
|
思路: 历史池 23k 化学式作为背景, 模拟 denovo 候选作为前景, 一起投到 UMAP/PCA
|
|
看 denovo 是否落到历史池未覆盖的"空白区"
|
|
输入: data/nomad_perovskite_aggregated.csv (历史池)
|
|
data/denovo_candidates.csv (可选, denovo 输出格式: formula, predicted_eg, ...)
|
|
输出: figures/fig_d1~d4_*.png + output/embed_coords.npz
|
|
"""
|
|
import warnings
|
|
from pathlib import Path
|
|
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
import pandas as pd
|
|
import re
|
|
from sklearn.decomposition import PCA
|
|
from sklearn.preprocessing import StandardScaler
|
|
|
|
warnings.filterwarnings("ignore")
|
|
|
|
try:
|
|
import umap
|
|
HAS_UMAP = True
|
|
except ImportError:
|
|
HAS_UMAP = False
|
|
print("⚠️ UMAP 未装, 仅用 PCA")
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
AGG = ROOT / "data" / "nomad_perovskite_aggregated.csv"
|
|
FIG = ROOT / "figures"
|
|
OUT = ROOT / "output"
|
|
FIG.mkdir(exist_ok=True); OUT.mkdir(exist_ok=True)
|
|
MAGPIE_DIR = ROOT / "magpie"
|
|
|
|
# ============ Magpie 元素特征 (取上次 03_ml_model 已缓存的 65 维) ============
|
|
# 注意: Magpie 表文件是「按原子序数顺序」排列的纯数值 (H=1, He=2, Li=3, ...),
|
|
# 不是「元素名=值」格式。
|
|
magpie_tables = {}
|
|
for p in MAGPIE_DIR.glob("*.table"):
|
|
name = p.stem
|
|
arr = []
|
|
with open(p) as f:
|
|
for line in f:
|
|
line = line.strip()
|
|
if not line: continue
|
|
try:
|
|
arr.append(float(line.split()[0]))
|
|
except (ValueError, IndexError):
|
|
continue
|
|
if arr:
|
|
magpie_tables[name] = np.array(arr, dtype=float)
|
|
|
|
print(f"Magpie 元素表: {len(magpie_tables)} 个, 每表 {len(next(iter(magpie_tables.values())))} 个元素")
|
|
|
|
# ============ 元素词典 ============
|
|
Z2SYM = ["H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Ar",
|
|
"K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
|
|
"Rb","Sr","Y","Zr","Nb","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I","Xe",
|
|
"Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu",
|
|
"Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn"]
|
|
SYM2IDX = {s: i for i, s in enumerate(Z2SYM)}
|
|
|
|
FORMULA_RE = re.compile(r"([A-Z][a-z]?)(\d*\.?\d*)")
|
|
|
|
|
|
def parse_formula(s):
|
|
out = {}
|
|
for el, num in FORMULA_RE.findall(str(s)):
|
|
n = float(num) if num else 1.0
|
|
out[el] = out.get(el, 0.0) + n
|
|
return out
|
|
|
|
|
|
def magpie_features(formula):
|
|
comp = parse_formula(formula)
|
|
feats = []
|
|
total = sum(comp.values())
|
|
if total == 0: total = 1.0
|
|
for tbl_name, tbl in magpie_tables.items():
|
|
vals = []
|
|
for el, n in comp.items():
|
|
if el in SYM2IDX and SYM2IDX[el] < len(tbl):
|
|
v = tbl[SYM2IDX[el]]
|
|
vals.append((v, n / total))
|
|
if vals:
|
|
weighted = [v * w for v, w in vals]
|
|
feats.extend([np.nansum(weighted), np.nanmax(weighted) - np.nanmin(weighted),
|
|
np.nanstd(weighted)])
|
|
else:
|
|
feats.extend([np.nan, np.nan, np.nan])
|
|
v = np.array(feats, dtype=float)
|
|
return np.nan_to_num(v, nan=0.0)
|
|
|
|
|
|
# Oxide = O is sole anion (exclude silicide-like / mixed-anion without O-only anion).
|
|
# Si treated as B-site (oxide perovskite-like), not anion.
|
|
ANION_ELEMS = {"O", "F", "Cl", "Br", "I", "S", "Se", "Te", "N", "P", "As", "Sb", "H", "C"}
|
|
def is_halide(f):
|
|
return bool({"I","Br","Cl","F"} & set(parse_formula(f)))
|
|
def is_oxide(f):
|
|
els = set(parse_formula(f))
|
|
if "O" not in els or is_halide(f):
|
|
return False
|
|
return (els & ANION_ELEMS) == {"O"}
|
|
|
|
|
|
# ============ 加载历史池 ============
|
|
print("\n加载历史池...")
|
|
agg = pd.read_csv(AGG)
|
|
hist = agg[agg.functional_class == "gga"].drop_duplicates("formula").reset_index(drop=True)
|
|
print(f" 历史池 GGA 唯一化学式: {len(hist)}")
|
|
|
|
# 分类标签
|
|
hist["category"] = hist.formula.map(lambda f: "halide" if is_halide(f) else ("oxide" if is_oxide(f) else "other"))
|
|
|
|
X_hist = np.stack([magpie_features(f) for f in hist.formula])
|
|
print(f" 特征矩阵: {X_hist.shape}")
|
|
|
|
# ============ 模拟 denovo 候选 (切 10% 历史作"代用", 真候选上来直接替换) ============
|
|
np.random.seed(42)
|
|
denovo_idx = np.random.choice(len(hist), size=len(hist)//10, replace=False)
|
|
denovo_idx_set = set(denovo_idx)
|
|
denovo = hist.iloc[denovo_idx].copy().reset_index(drop=True)
|
|
denovo["source"] = "simulated" # 上线真 denovo 时改 "denovo"
|
|
hist_labeled = hist.copy()
|
|
hist_labeled["source"] = "historical"
|
|
|
|
print(f" 模拟 denovo 候选: {len(denovo)}")
|
|
|
|
# ============ 合并 + PCA / UMAP ============
|
|
X_all = np.vstack([X_hist, magpie_features_batch := np.stack([magpie_features(f) for f in denovo.formula])])
|
|
y_source = np.array(["hist"]*len(hist) + ["denovo"]*len(denovo))
|
|
y_category_hist = hist.category.values
|
|
y_category_all = np.concatenate([y_category_hist, ["denovo"]*len(denovo)])
|
|
|
|
scaler = StandardScaler().fit(X_all)
|
|
Xs = scaler.transform(X_all)
|
|
|
|
# PCA
|
|
print("\nPCA 降维...")
|
|
pca = PCA(n_components=2, random_state=42).fit(Xs[:len(hist)]) # 只在历史池上训PCA
|
|
X_pca = pca.transform(Xs)
|
|
explained = pca.explained_variance_ratio_
|
|
print(f" PCA 前2主成分解释方差: {explained[0]:.3f}, {explained[1]:.3f} (合计 {sum(explained):.3f})")
|
|
|
|
# UMAP
|
|
if HAS_UMAP:
|
|
print("UMAP 降维 (历史池 fit, denovo transform)...")
|
|
um = umap.UMAP(n_neighbors=15, min_dist=0.1, n_components=2,
|
|
metric="euclidean", random_state=42)
|
|
um.fit(Xs[:len(hist)])
|
|
X_umap = um.transform(Xs)
|
|
else:
|
|
X_umap = None
|
|
|
|
# 存坐标
|
|
if X_umap is not None:
|
|
np.savez(OUT / "embed_coords.npz",
|
|
hist_pca=X_pca[:len(hist)], hist_umap=X_umap[:len(hist)],
|
|
denovo_pca=X_pca[len(hist):], denovo_umap=X_umap[len(hist):],
|
|
explained_var=explained,
|
|
hist_formula=hist.formula.values, denovo_formula=denovo.formula.values,
|
|
hist_category=hist.category.values)
|
|
|
|
# ============ 图 D1: UMAP 化学空间 (主图) ============
|
|
plt.rcParams["font.sans-serif"] = ["Noto Sans CJK JP", "DejaVu Sans"]
|
|
plt.rcParams["axes.unicode_minus"] = False
|
|
|
|
fig, ax = plt.subplots(figsize=(9, 7))
|
|
if X_umap is not None:
|
|
# 历史池: 浅灰背景
|
|
ax.scatter(X_umap[:len(hist), 0], X_umap[:len(hist), 1],
|
|
s=3, c="#cccccc", alpha=0.4, label="历史池 (n=23k)", zorder=1)
|
|
# 按类别着色 (小样本, 只画halide+oxide)
|
|
for cat, color, lbl in [("halide","#e74c3c","卤化物"),
|
|
("oxide","#3498db","氧化物")]:
|
|
m = (y_category_hist == cat)
|
|
ax.scatter(X_umap[:len(hist)][m, 0], X_umap[:len(hist)][m, 1],
|
|
s=6, c=color, alpha=0.6, label=f"{lbl} (历史池内)", zorder=2)
|
|
# denovo 候选: 紫色星号
|
|
ax.scatter(X_umap[len(hist):, 0], X_umap[len(hist):, 1],
|
|
s=70, c="#7b4fd0", marker="*", edgecolors="white", linewidth=1.2,
|
|
label=f"模拟 denovo 候选 (n={len(denovo)})", zorder=5)
|
|
ax.set_xlabel("UMAP-1")
|
|
ax.set_ylabel("UMAP-2")
|
|
else:
|
|
ax.text(0.5, 0.5, "UMAP 未安装", ha="center")
|
|
ax.set_title("钙钛矿化学空间 UMAP 投影\n(Magpie 65 维特征, denovo 候选 vs 历史池覆盖度)")
|
|
ax.legend(loc="best", fontsize=10, framealpha=0.9)
|
|
ax.grid(alpha=0.2)
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d1_chem_space_umap.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d1_chem_space_umap.png")
|
|
|
|
# ============ 图 D2: PCA baseline 对照 ============
|
|
fig, ax = plt.subplots(figsize=(9, 7))
|
|
ax.scatter(X_pca[:len(hist), 0], X_pca[:len(hist), 1],
|
|
s=3, c="#cccccc", alpha=0.4, label="历史池", zorder=1)
|
|
m_h = (y_category_hist == "halide")
|
|
m_o = (y_category_hist == "oxide")
|
|
ax.scatter(X_pca[:len(hist)][m_h, 0], X_pca[:len(hist)][m_h, 1],
|
|
s=6, c="#e74c3c", alpha=0.6, label="卤化物 (历史)", zorder=2)
|
|
ax.scatter(X_pca[:len(hist)][m_o, 0], X_pca[:len(hist)][m_o, 1],
|
|
s=6, c="#3498db", alpha=0.6, label="氧化物 (历史)", zorder=2)
|
|
ax.scatter(X_pca[len(hist):, 0], X_pca[len(hist):, 1],
|
|
s=70, c="#7b4fd0", marker="*", edgecolors="white", linewidth=1.2,
|
|
label=f"模拟 denovo (n={len(denovo)})", zorder=5)
|
|
ax.set_xlabel(f"PC1 ({explained[0]*100:.1f}%)")
|
|
ax.set_ylabel(f"PC2 ({explained[1]*100:.1f}%)")
|
|
ax.set_title(f"PCA 对照 (前2主成分累计解释 {sum(explained)*100:.1f}%)\n注: 解释方差偏低, 印证 UMAP 非线性投影的必要性")
|
|
ax.legend(loc="best", fontsize=10, framealpha=0.9)
|
|
ax.grid(alpha=0.2)
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d2_pca_baseline.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d2_pca_baseline.png")
|
|
|
|
# ============ 图 D3: 性质空间 (带隙色) ============
|
|
fig, axes = plt.subplots(1, 2, figsize=(15, 6.5))
|
|
for ax_i, (X_proj, title) in enumerate([(X_pca, "PCA"), (X_umap, "UMAP")]):
|
|
ax = axes[ax_i]
|
|
ax.scatter(X_proj[:len(hist), 0], X_proj[:len(hist), 1],
|
|
c=hist.band_gap_eV.values[:len(hist)], cmap="viridis",
|
|
s=4, alpha=0.35, vmin=0, vmax=5, zorder=1)
|
|
sc = ax.scatter(X_proj[len(hist):, 0], X_proj[len(hist):, 1],
|
|
c=denovo.band_gap_eV.values, cmap="viridis",
|
|
s=120, marker="*", edgecolors="white", linewidth=1.4,
|
|
vmin=0, vmax=5, zorder=5)
|
|
# 光伏窗口圈
|
|
m_pv = (denovo.band_gap_eV.values >= 1.1) & (denovo.band_gap_eV.values <= 1.7)
|
|
if m_pv.any():
|
|
ax.scatter(X_proj[len(hist):][m_pv, 0], X_proj[len(hist):][m_pv, 1],
|
|
s=400, facecolors="none", edgecolors="#d62728", linewidth=2.2,
|
|
label=f"denovo 光伏窗口候选 ({m_pv.sum()})", zorder=4)
|
|
ax.set_xlabel(f"{title}-1")
|
|
ax.set_ylabel(f"{title}-2")
|
|
ax.set_title(f"带隙着色 (历史池 vs denovo)\n红圈 = 光伏窗口 [1.1, 1.7] eV 候选")
|
|
ax.grid(alpha=0.2)
|
|
if ax_i == 1:
|
|
plt.colorbar(sc, ax=ax, label="DFT 带隙 (eV)")
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d3_property_space.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d3_property_space.png")
|
|
|
|
# ============ 图 D4: 覆盖度指标 (4 子图) ============
|
|
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
|
|
|
|
# D4a: B 位元素分布对比
|
|
def b_site(formula):
|
|
comp = parse_formula(formula)
|
|
# 启发式: 中等原子序数金属 (除 X 位卤素/O)
|
|
b_candidates = {k: v for k, v in comp.items()
|
|
if k not in ["I","Br","Cl","F","O","H","N","S"]
|
|
and 13 <= SYM2IDX.get(k, -1) + 1 <= 83}
|
|
return sorted(b_candidates.keys())
|
|
|
|
b_hist = [b_site(f) for f in hist.formula]
|
|
b_denovo = [b_site(f) for f in denovo.formula]
|
|
from collections import Counter
|
|
top_hist = Counter([b for bs in b_hist for b in bs]).most_common(10)
|
|
top_denovo = Counter([b for bs in b_denovo for b in bs]).most_common(10)
|
|
ax = axes[0, 0]
|
|
els = [t[0] for t in top_hist]
|
|
y_h = [t[1] for t in top_hist]
|
|
y_d = [Counter([b for bs in b_denovo for b in bs]).get(e, 0) for e in els]
|
|
x = np.arange(len(els))
|
|
w = 0.4
|
|
ax.bar(x - w/2, y_h, w, color="#888", label="历史池")
|
|
ax.bar(x + w/2, y_d, w, color="#7b4fd0", label="denovo")
|
|
ax.set_xticks(x); ax.set_xticklabels(els)
|
|
ax.set_ylabel("出现次数"); ax.set_title("(a) B 位元素分布对比")
|
|
ax.set_yscale("log"); ax.legend()
|
|
|
|
# D4b: Voronoi 体积（历史池外覆盖比例）
|
|
from scipy.spatial import Voronoi, voronoi_plot_2d
|
|
ax = axes[0, 1]
|
|
if X_umap is not None and len(hist) > 4:
|
|
# 在历史池UMAP上画Voronoi，看denovo落点
|
|
sample_size = min(2000, len(hist))
|
|
rng = np.random.RandomState(42)
|
|
s_idx = rng.choice(len(hist), sample_size, replace=False)
|
|
pts_hist = X_umap[s_idx]
|
|
pts_denovo = X_umap[len(hist):]
|
|
try:
|
|
vor = Voronoi(pts_hist)
|
|
# 简化: 只画denovo点在历史池Voronoi单元内的占比
|
|
from scipy.spatial import cKDTree
|
|
tree = cKDTree(pts_hist)
|
|
dist, _ = tree.query(pts_denovo, k=1)
|
|
median_dist = float(np.median(dist))
|
|
max_dist = float(np.percentile(dist, 90))
|
|
# "远离历史池" 阈值 = 第90百分位
|
|
outlier = (dist > max_dist).sum()
|
|
ax.scatter(pts_hist[:, 0], pts_hist[:, 1], s=1, c="#888", alpha=0.3, label="历史池采样")
|
|
ax.scatter(pts_denovo[:, 0], pts_denovo[:, 1], s=80, c="#7b4fd0", marker="*",
|
|
edgecolors="white", linewidth=1, label="denovo")
|
|
ax.set_title(f"(b) 化学空间覆盖度\ndenovo 到历史池最近邻: 中位 {median_dist:.2f}, P90 {max_dist:.2f}")
|
|
ax.set_xlabel("UMAP-1"); ax.set_ylabel("UMAP-2")
|
|
ax.grid(alpha=0.2); ax.legend()
|
|
except Exception as e:
|
|
ax.text(0.5, 0.5, f"Voronoi err: {e}", ha="center")
|
|
else:
|
|
ax.text(0.5, 0.5, "无 UMAP", ha="center")
|
|
|
|
# D4c: Pairwise 距离分布（多样性）
|
|
ax = axes[1, 0]
|
|
def pairwise_dist(X, n=200):
|
|
rng = np.random.RandomState(42)
|
|
if len(X) > n: X = X[rng.choice(len(X), n, replace=False)]
|
|
from scipy.spatial.distance import pdist
|
|
return pdist(X, metric="euclidean")
|
|
d_hist_full = pairwise_dist(Xs[:len(hist)], n=500)
|
|
d_denovo_full = pairwise_dist(Xs[len(hist):], n=min(500, len(denovo)))
|
|
ax.hist(d_hist_full, bins=40, alpha=0.5, color="#888", label=f"历史池内部 (μ={d_hist_full.mean():.2f})")
|
|
ax.hist(d_denovo_full, bins=40, alpha=0.6, color="#7b4fd0", label=f"denovo 内部 (μ={d_denovo_full.mean():.2f})")
|
|
ax.set_xlabel("Magpie 空间 pairwise 欧氏距离")
|
|
ax.set_ylabel("频次")
|
|
ax.set_title("(c) 多样性: 池内/候选内化学距离分布")
|
|
ax.legend()
|
|
|
|
# D4d: 候选与历史池的距离分布
|
|
ax = axes[1, 1]
|
|
from scipy.spatial.distance import cdist
|
|
if len(denovo) <= 500:
|
|
sample = denovo
|
|
else:
|
|
sample = denovo.sample(500, random_state=42)
|
|
X_d = np.stack([magpie_features(f) for f in sample.formula])
|
|
Xd_s = scaler.transform(X_d)
|
|
dist_min = cdist(Xd_s, Xs[:len(hist)]).min(axis=1)
|
|
ax.hist(dist_min, bins=30, color="#7b4fd0", edgecolor="white")
|
|
ax.axvline(np.percentile(dist_min, 50), color="green", ls="--",
|
|
label=f"P50 = {np.percentile(dist_min, 50):.2f}")
|
|
ax.axvline(np.percentile(dist_min, 90), color="orange", ls="--",
|
|
label=f"P90 = {np.percentile(dist_min, 90):.2f}")
|
|
ax.set_xlabel("denovo 候选到历史池的最小距离 (Magpie 标准化空间)")
|
|
ax.set_ylabel("候选数")
|
|
ax.set_title("(d) denovo 到历史池的覆盖距离\n候选越靠右 → 越可能是『历史外』新发现")
|
|
ax.legend()
|
|
|
|
fig.suptitle("Denovo 覆盖度评估四联图 (UMAP 投影 + B 位分布 + 多样性 + 覆盖距离)",
|
|
fontsize=13, y=1.005)
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d4_coverage_metrics.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d4_coverage_metrics.png")
|
|
|
|
# ============ 摘要输出 ============
|
|
print("\n" + "="*60)
|
|
print("化学空间可视化摘要")
|
|
print("="*60)
|
|
print(f"历史池 GGA 唯一化学式: {len(hist)}")
|
|
print(f"模拟 denovo 候选: {len(denovo)} (切10%历史池作代用)")
|
|
print(f"PCA 前2主成分累计解释: {sum(explained)*100:.1f}%")
|
|
print(f"denovo→历史池最近邻距离 P50/P90: "
|
|
f"{np.percentile(dist_min, 50):.3f} / {np.percentile(dist_min, 90):.3f}")
|
|
print(f"\n生成图:")
|
|
print(f" fig_d1_chem_space_umap.png 主图: UMAP 化学空间")
|
|
print(f" fig_d2_pca_baseline.png 对照: 线性 PCA 解释方差不足, 印证 UMAP 必要性")
|
|
print(f" fig_d3_property_space.png 带隙着色 + 光伏窗口候选圈")
|
|
print(f" fig_d4_coverage_metrics.png 覆盖度四联图")
|
|
print(f"\n坐标已存: output/embed_coords.npz (后续接真 denovo 直接复用)") |