示例04 VR 复形随半径演化

"""ex04_vr_progression.py — 第 3 章 Vietoris-Rips 复形随半径增长

对应教程：tutorials/03-vietoris-rips.md。

直觉化展示 VR 过滤：随半径 ε 增大，单纯形（边、三角形）依次加入，
连通分量合并、洞出现又消失。这里通过在不同 thresh 下计算，
观察 H0/H1 点数随 ε 的变化（Betti 曲线雏形）。

运行：
python ex04_vr_progression.py
"""

import numpy as np
from ripser import ripser
import matplotlib.pyplot as plt

from common import sample_circle

def main():
data = sample_circle(n=100, noise=0.05, seed=3)

# 一组递增的 thresh，对应过滤的不同「快照」
thresholds = np.linspace(0.2, 2.0, 15)

betti0 = []   # 每个 thresh 下「存活」的 H0 数
betti1 = []
for eps in thresholds:
    # thresh=eps 等价于过滤到 ε=eps 时停
    dgm = ripser(data, maxdim=1, thresh=float(eps))['dgms']
    # Betti_k = birth ≤ eps < death 的特征数
    b0 = np.sum((dgm[0][:, 0] <= eps) &
                (dgm[0][:, 1] > eps))
    b1 = np.sum((dgm[1][:, 0] <= eps) &
                (dgm[1][:, 1] > eps))
    betti0.append(b0)
    betti1.append(b1)

betti0 = np.array(betti0)
betti1 = np.array(betti1)

print("=== Betti 数随 thresh 变化 ===")
print(f"{'ε':>6s} {'β0':>4s} {'β1':>4s}")
for eps, b0, b1 in zip(thresholds, betti0, betti1):
    print(f"{eps:6.2f} {b0:4d} {b1:4d}")

# ── 可视化 ──
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

axes[0].scatter(data[:, 0], data[:, 1], s=12, alpha=0.7)
axes[0].set_aspect('equal')
axes[0].set_title('圆环点云')
axes[0].grid(alpha=0.3)

axes[1].plot(thresholds, betti0, 'b-o', ms=4, label='β0 (连通分量)')
axes[1].plot(thresholds, betti1, 'r-s', ms=4, label='β1 (洞)')
axes[1].set_xlabel('过滤半径 ε')
axes[1].set_ylabel('Betti 数')
axes[1].set_title('Betti 曲线（VR 过滤随 ε 演化）')
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

print("\n解读：")
print("  - β0 从 100（每点独立）随 ε 增大降到 1（全连通）")
print("  - β1 在某个 ε 区间为 1（圆环的洞存活），之后被填平归 0")
​

if name == "main":
main()