示例07 预计算与稀疏距离矩阵

"""ex07_distance_matrix.py — 第 6 章 预计算与稀疏距离矩阵

对应教程：tutorials/06-sparse-greedy.md。

演示：

• pairwise_distances + distance_matrix=True：支持任意度量
• 稀疏 CSR 矩阵：只保留 ≤ thresh 的边，走 C++ 稀疏路径
• num_edges 随稀疏度变化

运行：
python ex07_distance_matrix.py
"""

import numpy as np
from scipy import sparse
from sklearn.metrics import pairwise_distances
from ripser import ripser

from common import sample_circle

def main():
data = sample_circle(n=200, noise=0.05, seed=2)

# ── 1. 预计算稠密距离矩阵：曼哈顿距离 ──
print("=== 预计算稠密距离矩阵（曼哈顿度量）===")
D_manhattan = pairwise_distances(data, metric='manhattan')
print(f"  D shape={D_manhattan.shape}, 对称: "
      f"{np.allclose(D_manhattan, D_manhattan.T)}")
result = ripser(D_manhattan, distance_matrix=True, maxdim=1, thresh=1.5)
pers1 = result['dgms'][1][:, 1] - result['dgms'][1][:, 0]
pers1 = pers1[np.isfinite(pers1)]
print(f"  H1: {len(result['dgms'][1])} 点, "
      f"max persistence={pers1.max() if len(pers1) else 0:.4f}, "
      f"edges={result['num_edges']}")

# ── 2. 欧氏距离对照 ──
D_euclid = pairwise_distances(data, metric='euclidean')
result2 = ripser(D_euclid, distance_matrix=True, maxdim=1, thresh=1.5)
pers2 = result2['dgms'][1][:, 1] - result2['dgms'][1][:, 0]
pers2 = pers2[np.isfinite(pers2)]
print(f"  欧氏 H1: {len(result2['dgms'][1])} 点, "
      f"max persistence={pers2.max() if len(pers2) else 0:.4f}")

# ── 3. 稀疏距离矩阵：只保留 ≤ thresh 的边 ──
print("\n=== 稀疏距离矩阵（只保留短边）===")
thresh = 1.0
D_sparse = sparse.csr_matrix(D_euclid * (D_euclid <= thresh))
nnz = D_sparse.nnz
total = D_euclid.size
print(f"  非零元素: {nnz}/{total} = {nnz/total:.1%}")

result3 = ripser(D_sparse, distance_matrix=True, maxdim=1, thresh=thresh)
print(f"  稀疏路径 edges={result3['num_edges']}")
pers3 = result3['dgms'][1][:, 1] - result3['dgms'][1][:, 0]
pers3 = pers3[np.isfinite(pers3)]
print(f"  稀疏 H1 max persistence="
      f"{pers3.max() if len(pers3) else 0:.4f}")

# ── 4. thresh 对 num_edges 的影响 ──
print("\n=== num_edges 随 thresh 变化 ===")
print(f"{'thresh':>8s} {'nnz':>8s} {'edges':>8s}")
for t in [0.3, 0.5, 0.8, 1.0, 1.5]:
    Ds = sparse.csr_matrix(D_euclid * (D_euclid <= t))
    r = ripser(Ds, distance_matrix=True, maxdim=1, thresh=t)
    print(f"{t:8.2f} {Ds.nnz:8d} {r['num_edges']:8d}")

# ── 5. 验证 n_perm + sparse 会报错 ──
print("\n=== 验证 n_perm + sparse 触发 ValueError ===")
try:
    ripser(D_sparse, distance_matrix=True, maxdim=1, n_perm=50)
    print("  未报错（可能与版本行为不同）")
except (ValueError, Exception) as e:
    print(f"  预期报错: {type(e).__name__}: {e}")
​

if name == "main":
main()