示例08 贪心置换子采样

"""ex08_greedy_perm.py — 第 6 章 贪心置换子采样（最远点采样）

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

演示：

• n_perm 子采样的最远点采样
• idx_perm / r_cover / dperm2all 字段
• 不同 n_perm 下的速度与精度（与全量的瓶颈距离）

运行：
python ex08_greedy_perm.py
"""

import time
import numpy as np
from ripser import ripser
from persim import bottleneck

from common import sample_sphere

def main():
# 用一个稍大的点云放大子采样的收益
data = sample_sphere(n=1500, noise=0.02, seed=4)
print(f"原始点云 n={len(data)}")

# ── 1. 全量计算（基准）──
t0 = time.perf_counter()
r_full = ripser(data, maxdim=1, thresh=1.5)
t_full = time.perf_counter() - t0
print(f"\n全量: time={t_full:.2f}s, edges={r_full['num_edges']}, "
      f"H1 点={len(r_full['dgms'][1])}")

# ── 2. 不同 n_perm 子采样 ──
print(f"\n{'n_perm':>8s} {'time':>8s} {'r_cover':>10s} "
      f"{'edges':>8s} {'bottleneck_H1':>14s}")
for n_perm in [200, 400, 800]:
    t0 = time.perf_counter()
    r = ripser(data, maxdim=1, n_perm=n_perm)
    dt = time.perf_counter() - t0

    # 子采样索引与覆盖半径
    idx = r['idx_perm']
    r_cover = r['r_cover']
    subsample = data[idx]

    # 与全量 H1 的瓶颈距离（衡量近似误差）
    try:
        d, _ = bottleneck(r['dgms'][1], r_full['dgms'][1])
    except Exception:
        d = float('nan')
    print(f"{n_perm:8d} {dt:8.2f}s {r_cover:10.4f} "
          f"{r['num_edges']:8d} {d:14.4f}")

# ── 3. dperm2all：子集到全集的距离 ──
r = ripser(data, maxdim=1, n_perm=400)
D_sub = r['dperm2all']
print(f"\ndperm2all shape={D_sub.shape}  (n_perm × n_total)")
print("  含义：子采样点到全部点的距离，可用于回射或加权。")

# ── 4. 子采样点可视化 ──
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(data[:, 0], data[:, 1], data[:, 2], s=2, alpha=0.2,
           c='lightgray', label='全集')
sub = data[r['idx_perm']]
ax.scatter(sub[:, 0], sub[:, 1], sub[:, 2], s=20, c='red',
           label=f'贪心子采样 (n={len(sub)})')
ax.set_title(f'最远点采样  r_cover={r["r_cover"]:.4f}')
ax.legend()
plt.show()
​

if name == "main":
main()