示例11 代表上循环与系数域

"""ex11_cocycles_coeff.py — 第 8 章 代表上循环与系数域

对应教程：tutorials/08-cocycles-coefficients.md。

演示：

• do_cocycles=True 提取代表上循环
• cocycle 数据格式 (m, d+1)：前 d 列顶点索引，最后列系数
• 从 H1 cocycle 提取支撑边，叠加到点云上可视化「洞」绕哪一圈
• coeff=2 vs coeff=17 对圆环持久图的影响（通常一致）

运行：
python ex11_cocycles_coeff.py
"""

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

from common import sample_circle

def cocycle_support_edges(cocycle_1d, prime):
"""提取 1-cocycle 非零边 (i, j)。"""
edges = []
for row in cocycle_1d:
i, j, val = int(row[0]), int(row[1]), int(row[2])
if val % prime != 0:
edges.append((min(i, j), max(i, j)))
return edges

def main():
data = sample_circle(n=100, noise=0.05, seed=8)
prime = 17

# ── 1. 开启 cocycles ──
result = ripser(data, maxdim=1, coeff=prime, do_cocycles=True)
dgms = result['dgms']
cocycles = result['cocycles']

print("=== dgms 与 cocycles 一一对应 ===")
for d in range(len(dgms)):
    print(f"  H{d}: {len(dgms[d])} 个特征 ↔ "
          f"{len(cocycles[d])} 个 cocycle")

# ── 2. 取最显著的 H1 特征的 cocycle ──
pers = dgms[1][:, 1] - dgms[1][:, 0]
pers_finite = pers.copy()
pers_finite[~np.isfinite(pers_finite)] = -np.inf
k_top = int(np.argmax(pers_finite))
top_dgm = dgms[1][k_top]
top_cyc = cocycles[1][k_top]
print(f"\n最显著 H1: birth={top_dgm[0]:.3f}  death={top_dgm[1]:.3f}")
print(f"  对应 cocycle shape={top_cyc.shape}")
print(f"  前 5 行 [i, j, value mod {prime}]:")
for row in top_cyc[:5]:
    print(f"    {int(row[0]):3d} {int(row[1]):3d} {int(row[2])}")

# ── 3. 提取支撑边并可视化 ──
edges = cocycle_support_edges(top_cyc, prime)
print(f"  非零边数: {len(edges)}")

fig, ax = plt.subplots(figsize=(7, 7))
ax.scatter(data[:, 0], data[:, 1], s=15, c='steelblue', zorder=2)
for i, j in edges:
    ax.plot([data[i, 0], data[j, 0]],
            [data[i, 1], data[j, 1]],
            'r-', lw=0.5, alpha=0.5, zorder=1)
ax.set_aspect('equal')
ax.set_title(f'H1 代表 cocycle 的支撑边（绕洞一圈）\n'
             f'prime={prime}, 非零边 {len(edges)} 条')
ax.grid(alpha=0.3)

# ── 4. coeff=2 vs coeff=17 对比 ──
print("\n=== coeff 对比（圆环通常一致）===")
for coeff in [2, 17]:
    r = ripser(data, maxdim=1, coeff=coeff)
    p = r['dgms'][1][:, 1] - r['dgms'][1][:, 0]
    p = p[np.isfinite(p)]
    print(f"  coeff={coeff:3d}: H1 点数={len(r['dgms'][1])}, "
          f"max persistence={p.max() if len(p) else 0:.4f}")

plt.show()
​

if name == "main":
main()