Possible bug in computing diagrams for rips filtration

[cdeepakroy]

I am just getting started with using the dionysus on some toy data.

Below i tried to sample n points from a circle perturbed by gaussian noise of a specified scale

import numpy as np

num_pts = 500
r = 10
noise_sigma = 0.05 * r

theta = np.random.rand(num_pts, 1) * 2 * np.pi

pts = r * np.hstack((np.cos(theta), np.sin(theta))) 
pts = pts + np.random.randn(num_pts, 2) * noise_sigma

plt.figure(figsize=(8, 8))
_ = plt.plot(pts[:, 0], pts[:, 1], 'k+')

Then i tried to compute persistence diagrams upto dimension 1 using rips filtration using the code below:

import dionysus as d

f = d.fill_rips(pts, 1, 3*r)
m = d.homology_persistence(f)
dgms = d.init_diagrams(m, f)
print dgms

Output:

Took 1.673817 seconds for 500 points
[Diagram with 500 points, Diagram with 124251 points]

Then i tried to plot the diagrams of dimension 0 and 1 using the following code:

import matplotlib.pyplot as plt

plt.figure(figsize=(8,8))
d.plot.plot_diagram(dgms[0], show=False)

plt.figure(figsize=(8,8))
d.plot.plot_diagram(dgms[1], show=False)

The diagram of dimension 0 looks as expected:

However the plot of diagram of dimension 1 produced the following error

ValueError                                Traceback (most recent call last)
<ipython-input-88-65ea1b07beec> in <module>()
      1 plt.figure(figsize=(8,8))
----> 2 d.plot.plot_diagram(dgms[1], show=False)

/home/cdeepakroy/work/Libs/dionysus/build/bindings/python/dionysus/plot.pyc in plot_diagram(dgm, show)
      7     min_birth = min(p.birth for p in dgm if p.birth != inf)
      8     max_birth = max(p.birth for p in dgm if p.birth != inf)
----> 9     min_death = min(p.death for p in dgm if p.death != inf)
     10     max_death = max(p.death for p in dgm if p.death != inf)
     11 

ValueError: min() arg is an empty sequence

<matplotlib.figure.Figure at 0x7f78f7e35c10>

The reason is the death value of all birth-death pairs in diagram of dimension 1 is infinity.

This seems to be a bug or did i misinterpret any of the parameters of the rips filtration.

[mrzv]

You need at least 2-skeleton to compute 1-dimensional homology. So the skeleton parameter to fill_rips() should be 2.

[mrzv]

By the way, when working with Rips complexes, you are much better off using cohomology. So in your snippet above, replace m = d.homology_persistence(f) with m = d.cohomology_persistence(f). The computation will be significantly faster.

[cdeepakroy]

Thanks! It works!