Python Ripser 获取图表点的顶点

Question

我正在使用 Python 包 ripser 进行持久性同源。我想利用它来帮助分割 2D 点云。

例如，我正在关注 Elizabeth Munch: Python Tutorial on Topological Data Analysis。在这里，我取 DoubleAnnulus 并增加两者之间的间隔：

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import persim
import ripser
import teaspoon.MakeData.DynSysLib.DynSysLib as DSL
import teaspoon.MakeData.PointCloud as makePtCloud
import teaspoon.TDA.Draw as Draw
from teaspoon.parameter_selection.MsPE import MsPE_tau
from teaspoon.SP.network import ordinal_partition_graph
from teaspoon.SP.network_tools import make_network
from teaspoon.TDA.PHN import PH_network


def DoubleAnnulus(r1=1, R1=2, r2=0.8, R2=1.3, xshift=3):
    P = makePtCloud.Annulus(r=r1, R=R1)
    Q = makePtCloud.Annulus(r=r2, R=R2)

    Q[:, 0] = Q[:, 0] + xshift

    return np.concatenate((P, Q))


def drawTDAtutorial(P, diagrams, R=2):
    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))

    plt.sca(axes[0])
    plt.title("Point Cloud")
    plt.scatter(P[:, 0], P[:, 1])

    plt.sca(axes[1])
    plt.title("0-dim Diagram")
    Draw.drawDgm(diagrams[0])

    plt.sca(axes[2])
    plt.title("1-dim Diagram")
    Draw.drawDgm(diagrams[1])
    plt.axis([0, R, 0, R])

    plt.show()


if __name__ == "__main__":

    P = DoubleAnnulus(r1=1, R1=2, r2=0.5, R2=1.3, xshift=6)
    plt.scatter(*zip(*P))
    plt.show()

    diagrams = ripser.ripser(P)["dgms"]

    drawTDAtutorial(P, diagrams)

如何抓取1-dim Diagram上最大的2个图表点对应的点云顶点？

Answer 1

K 均值聚类在这里会很好地工作（我们的数据相当凸，算法不需要标签）。由于持久性图上的关键特征与线 y = x 的 y 距离很大，我们可以生成与线 y = x 的 1-dim Diagram y 距离的直方图，并使用点数超过平均值 3 个标准偏差作为我们的 k 值。

代码

import pprint

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import persim
import ripser
import teaspoon.MakeData.DynSysLib.DynSysLib as DSL
import teaspoon.MakeData.PointCloud as makePtCloud
import teaspoon.TDA.Draw as Draw
from sklearn.cluster import KMeans
from teaspoon.parameter_selection.MsPE import MsPE_tau
from teaspoon.SP.network import ordinal_partition_graph
from teaspoon.SP.network_tools import make_network
from teaspoon.TDA.PHN import PH_network


def DoubleAnnulus(r1=1, R1=2, r2=0.8, R2=1.3, xshift=3):
    P = makePtCloud.Annulus(r=r1, R=R1)
    Q = makePtCloud.Annulus(r=r2, R=R2)

    Q[:, 0] = Q[:, 0] + xshift

    return np.concatenate((P, Q))


def drawTDAtutorial(P, diagrams, R=2):
    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))

    plt.sca(axes[0])
    plt.title("Point Cloud")
    plt.scatter(P[:, 0], P[:, 1])

    plt.sca(axes[1])
    plt.title("0-dim Diagram")
    Draw.drawDgm(diagrams[0])

    plt.sca(axes[2])
    plt.title("1-dim Diagram")
    Draw.drawDgm(diagrams[1])
    plt.axis([0, R, 0, R])

    plt.show()


if __name__ == "__main__":

    P = DoubleAnnulus(r1=1, R1=2, r2=0.5, R2=1.3, xshift=6)
    # plt.scatter(*zip(*P))
    # plt.show()

    result = ripser.ripser(P)

    diagrams = result["dgms"]
    # drawTDAtutorial(P, diagrams, R=3)

    loop_points = diagrams[1]  # H1 class
    y_heights = loop_points[:, 1] - loop_points[:, 0]
    features = y_heights[y_heights - np.mean(y_heights) >= 3 * np.std(y_heights)]
    k = len(features)

    print(f"Number of Features: {k}")

    kmeans = KMeans(n_clusters=k).fit(P)
    labels = kmeans.labels_

    points1 = P[labels == 0]
    points2 = P[labels == 1]

    fig, ax = plt.subplots()

    ax.scatter(*zip(*points1), color="red")
    ax.scatter(*zip(*points2), color="blue")

    plt.show()

输出

Number of Features: 2

Python Ripser 获取图表点的顶点

Python Ripser Get Vertices for Diagram Point

python

persistence

point-clouds

代码

输出