示例10 时间序列 Lower-Star

"""ex10_lower_star_timeseries.py — 第 7 章 时间序列 Lower-Star 过滤

对应教程：tutorials/07-lower-star.md 7.4 节。

把一维时间序列视作路径图上的函数：
顶点 = 时刻，边 = 相邻时刻，边权 = max(f[i], f[j])
构造稀疏距离矩阵后用 ripser(distance_matrix=True) 计算。

观察：周期信号会产生规律的 H0 bar 模式。

运行：
python ex10_lower_star_timeseries.py
"""

import numpy as np
from scipy import sparse
from ripser import ripser
from persim import plot_diagrams, plot_bars
import matplotlib.pyplot as plt

def time_series_lower_star(signal):
"""构造时序 Lower-Star 稀疏距离矩阵。

顶点 i 的「距离」= signal[i]；边 (i, i+1) 的权 = max(signal[i], signal[i+1])。
ripser 把这视作距离矩阵，distance_matrix=True 触发下星语义。
"""
n = len(signal)
rows, cols, vals = [], [], []
# 顶点（自环权 = 自身函数值）
for i in range(n):
    rows.append(i); cols.append(i); vals.append(float(signal[i]))
# 相邻边
for i in range(n - 1):
    w = float(max(signal[i], signal[i + 1]))
    # 对称两条
    rows.extend([i, i + 1]); cols.extend([i + 1, i]); vals.extend([w, w])
return sparse.coo_matrix((vals, (rows, cols)), shape=(n, n))
​

def main():
n = 200
t = np.linspace(0, 4 * np.pi, n)
rng = np.random.default_rng(0)
signal = np.sin(t) + 0.1 * rng.standard_normal(n)

D = time_series_lower_star(signal)
result = ripser(D, distance_matrix=True, maxdim=0)
dgm = result['dgms'][0]

pers = dgm[:, 1] - dgm[:, 0]
pers = pers[np.isfinite(pers)]
print(f"时序 Lower-Star H0: {len(dgm)} 点, "
      f"max persistence={pers.max() if len(pers) else 0:.4f}")
print(f"  显著 H0 数（persistence > 1.0）: "
      f"{int(np.sum(pers > 1.0))}")
print("  解读：每个显著 H0 对应信号的一个「峰谷」结构。")

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

axes[0].plot(t, signal, lw=0.8)
axes[0].set_title('正弦时间序列（带噪声）')
axes[0].set_xlabel('t'); axes[0].set_ylabel('f(t)')
axes[0].grid(alpha=0.3)

plot_bars(dgm, ax=axes[1])
axes[1].set_title('H0 条形码（峰值结构）')

plt.tight_layout()
plt.show()
​

if name == "main":
main()