Unexpected homology of boundary of a cube

[mbr085]

When a regular grid of the boundary of the 3-cube is fed into ripser, the expected 2-dimensional homology class does not appear.

Sample code: (Python interface)

from ripser import ripser, plot_dgms
n=3
segment = [np.linspace(0,1,20)]
endpoints = [np.linspace(0,1,2)]
face = segment * (n - 1) + endpoints
corners = []
for k in range(n):
corners.extend(itertools.product(*(face[k:] + face[:k])))
coords = np.array(corners)

rips = ripser(coords, maxdim=2, thresh=.4)
plot_dgms(rips['dgms'])

[ubauer]

Hi Morten, this is a duplicate of #16. In fact, surprisingly, in this example there is no nontrivial homology in degree 2, so the computation is correct.

[ubauer]

Hi Morten,

thanks for pointing this out. Can you attach a point cloud to run it in C++ Ripser?

Best, Uli

[mbr085]

Hi Uli,

yes I can. First, I checked that GUDHI behaves differently, and indeed it does. Below I paste in the python code I used in a jupyter notebook. In the next comment I will paste in the point cloud.

import gudhi
import numpy as np
import itertools
from ripser import ripser, plot_dgms

def cube_surface(dimension=2, num_segment_points=3):
segment = [np.linspace(0, 1, num_segment_points)]
endpoints = [np.linspace(0,1, dimension)]
face = segment * dimension + endpoints
corners = []
for k in range(dimension + 1):
corners.extend(itertools.product(*(face[k:] + face[:k])))
return np.array(corners)

rips_complex = gudhi.RipsComplex(points=cube_surface(dimension=2, num_segment_points=3),
max_edge_length=12.0)

simplex_tree = rips_complex.create_simplex_tree(max_dimension=2)
result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' +
repr(simplex_tree.num_simplices()) + ' simplices - ' +
repr(simplex_tree.num_vertices()) + ' vertices.'
print(result_str)

diag = simplex_tree.persistence(homology_coeff_field=2, min_persistence=0)
pplot = gudhi.plot_persistence_diagram(diag)

coords = cube_surface(dimension=2, num_segment_points=3)
rips = ripser(coords, maxdim=2)

plot_dgms(rips['dgms'])

[mbr085]

[[0.0, 0.0, 0.0],
[0.0, 0.0, 1.0],
[0.0, 0.5, 0.0],
[0.0, 0.5, 1.0],
[0.0, 1.0, 0.0],
[0.0, 1.0, 1.0],
[0.5, 0.0, 0.0],
[0.5, 0.0, 1.0],
[0.5, 0.5, 0.0],
[0.5, 0.5, 1.0],
[0.5, 1.0, 0.0],
[0.5, 1.0, 1.0],
[1.0, 0.0, 0.0],
[1.0, 0.0, 1.0],
[1.0, 0.5, 0.0],
[1.0, 0.5, 1.0],
[1.0, 1.0, 0.0],
[1.0, 1.0, 1.0],
[0.0, 0.0, 0.0],
[0.0, 0.0, 0.5],
[0.0, 0.0, 1.0],
[0.0, 1.0, 0.0],
[0.0, 1.0, 0.5],
[0.0, 1.0, 1.0],
[0.5, 0.0, 0.0],
[0.5, 0.0, 0.5],
[0.5, 0.0, 1.0],
[0.5, 1.0, 0.0],
[0.5, 1.0, 0.5],
[0.5, 1.0, 1.0],
[1.0, 0.0, 0.0],
[1.0, 0.0, 0.5],
[1.0, 0.0, 1.0],
[1.0, 1.0, 0.0],
[1.0, 1.0, 0.5],
[1.0, 1.0, 1.0],
[0.0, 0.0, 0.0],
[0.0, 0.0, 0.5],
[0.0, 0.0, 1.0],
[0.0, 0.5, 0.0],
[0.0, 0.5, 0.5],
[0.0, 0.5, 1.0],
[0.0, 1.0, 0.0],
[0.0, 1.0, 0.5],
[0.0, 1.0, 1.0],
[1.0, 0.0, 0.0],
[1.0, 0.0, 0.5],
[1.0, 0.0, 1.0],
[1.0, 0.5, 0.0],
[1.0, 0.5, 0.5],
[1.0, 0.5, 1.0],
[1.0, 1.0, 0.0],
[1.0, 1.0, 0.5],
[1.0, 1.0, 1.0]]

[mbr085]

This time I used a 2 by 2 grid on each face of the cube. It suffices to obtain the unexpected behaviour. I hope this will be easy to fix.
Best regards,
Morten

[ubauer]

Hi Morten,

I ran Ripser on this data set. The class in degree 2 shows up, but its interval is [0.5,1), so your threshold 0.4 is too small to see it. Guhdi might use a different scaling (radius instead of diameter as the parameter).

[mbr085]

I failed to reproduce when I had no threshold parameter. The threshold parameter is the culpit.

This succeeds:
coords = cube_surface(dimension=2, num_segment_points=5)
rips = ripser(coords, maxdim=2)
rips['dgms'][2]

array([[0.35355338, 1. ]])

This fails:
coords = cube_surface(dimension=2, num_segment_points=5)
rips = ripser(coords, maxdim=2, thresh=1.2)
rips['dgms'][2]

array([[0.35355338, 0.5 ],
[0.5 , 0.559017 ],
[0.559017 , 0.61237246],
[0.61237246, 0.70710677],
[0.70710677, 0.75 ],
[0.75 , 0.79056942],
[0.79056942, 0.82915622],
[0.82915622, 0.86602539],
[0.86602539, 0.90138781],
[0.90138781, 0.93541437],
[0.93541437, 1. ]])

[mbr085]

Here I have added more grid points. I will see if I can reproduce the behaviour with 3 times 3 grids.

[mbr085]

Indeed, on the data I sent you, setting thresh=1.2 yields the following classes:
array([[0.70710677, 0.86602539],
[0.86602539, 1. ]])

With no threshold it yields the class
array([[0.70710677, 1. ]])

[mbr085]

Now I tested this data in live ripser, and there I got the expected result. With and without threshold.

[mbr085]

Ripser from command line also gives the expected output.
So - it seems to be a problem with the python bindings. I guess this is good news.

[ubauer]

Could you please file a report at https://github.com/scikit-tda/ripser.py/issues?