Add cocycles feature

gph/bindings/ripser_bindings.cpp

@@ -39,15 +39,13 @@ to_numpy_barcodes(std::vector<barcodes_t> barcodes)
 }
 
 #if defined USE_COEFFICIENTS
-PYBIND11_MODULE(gph_ripser_coeff, m)
-{
+PYBIND11_MODULE(gph_ripser_coeff, m) {
 #else
-PYBIND11_MODULE(gph_ripser, m)
-{
+PYBIND11_MODULE(gph_ripser, m) {
 #endif
 
-    using namespace pybind11::literals;
-    m.doc() = "Ripser python interface";
+  using namespace pybind11::literals;
+  m.doc() = "Ripser python interface";
 
     py::class_<flagPersGen>(m, "flagPersGen", py::module_local())
         .def_readwrite("finite_0", &flagPersGen::finite_0)
@@ -63,14 +61,15 @@ PYBIND11_MODULE(gph_ripser, m)
              [&](ripserResults& res) {
                  return to_numpy_barcodes(res.births_and_deaths_by_dim);
              })
+        .def_readwrite("cocycles_by_dim", &ripserResults::cocycles_by_dim)
         .def_readwrite("flag_persistence_generators_by_dim",
                        &ripserResults::flag_persistence_generators);
 
     m.def(
         "rips_dm",
         [](py::array_t<value_t>& D, py::array_t<value_t>& diag, int modulus,
            int dim_max, float threshold, int num_threads,
-           bool return_generators) {
+           bool return_generators, const bool do_cocycles) {
             // Setup distance matrix and figure out threshold
             auto D_ = static_cast<value_t*>(D.request().ptr);
             std::vector<value_t> distances(D_, D_ + D.size());
@@ -84,36 +83,38 @@ PYBIND11_MODULE(gph_ripser, m)
 
             ripserResults res;
             ripser<compressed_lower_distance_matrix> r(
-                std::move(dist), dim_max, threshold, modulus, num_threads,
-                return_generators);
+                std::move(dist), dim_max, threshold, modulus,
+                num_threads, return_generators, do_cocycles);
             r.compute_barcodes();
             r.copy_results(res);
             return res;
         },
-        "D"_a, "diag"_a, "modulus"_a, "dim_max"_a, "threshold"_a,
-        "num_threads"_a, "return_generators"_a, "ripser distance matrix");
+        "D"_a, "N"_a, "modulus"_a, "dim_max"_a, "threshold"_a, "num_threads"_a,
+        "return_generators"_a, "do_cocycles"_a, "ripser distance matrix");
 
     m.def(
         "rips_dm_sparse",
         [](py::array_t<index_t>& I, py::array_t<index_t>& J,
            py::array_t<value_t>& V, int NEdges, int N, int modulus, int dim_max,
-           float threshold, int num_threads, bool return_generators) {
+           float threshold, int num_threads, bool return_generators, 
+           const bool do_cocycles) {
             auto I_ = static_cast<index_t*>(I.request().ptr);
             auto J_ = static_cast<index_t*>(J.request().ptr);
             auto V_ = static_cast<value_t*>(V.request().ptr);
 
             // Setup distance matrix and figure out threshold
             ripser<sparse_distance_matrix> r(
                 sparse_distance_matrix(I_, J_, V_, NEdges, N, threshold),
-                dim_max, threshold, modulus, num_threads, return_generators);
+                dim_max, threshold, modulus, num_threads,
+                return_generators, do_cocycles);
             r.compute_barcodes();
 
             ripserResults res;
             r.copy_results(res);
             return res;
         },
         "I"_a, "J"_a, "V"_a, "NEdges"_a, "N"_a, "modulus"_a, "dim_max"_a,
-        "threshold"_a, "num_threads"_a, "return_generators"_a,
+        "threshold"_a, "num_threads"_a, "return_generators"_a, "do_cocycles"_a,
         "ripser sparse distance matrix");
 
     m.def("get_max_coefficient_field_supported",

gph/python/ripser_interface.py

@@ -35,26 +35,31 @@
 
 
 def _compute_ph_vr_dense(DParam, diagonal, maxHomDim, thresh=-1, coeff=2,
-                         n_threads=1, return_generators=False):
+                         n_threads=1, return_generators=False,
+                         return_cocycles=False):
     if coeff == 2:
-        ret = gph_ripser.rips_dm(DParam, diagonal, coeff, maxHomDim, thresh,
-                                 n_threads, return_generators)
+        ret = gph_ripser.rips_dm(DParam, diagonal, coeff,
+                                 maxHomDim, thresh, n_threads,
+                                 return_generators, return_cocycles)
     else:
-        ret = gph_ripser_coeff.rips_dm(DParam, diagonal, coeff, maxHomDim,
-                                       thresh, n_threads, return_generators)
+        ret = gph_ripser_coeff.rips_dm(DParam, diagonal, coeff,
+                                       maxHomDim, thresh, n_threads,
+                                       return_generators, return_cocycles)
     return ret
 
 
 def _compute_ph_vr_sparse(I, J, V, N, maxHomDim, thresh=-1, coeff=2,
-                          n_threads=1, return_generators=False):
+                          n_threads=1, return_generators=False,
+                          return_cocycles=False):
     if coeff == 2:
         ret = gph_ripser.rips_dm_sparse(I, J, V, I.size, N, coeff,
                                         maxHomDim, thresh, n_threads,
-                                        return_generators)
+                                        return_generators, return_cocycles)
     else:
         ret = gph_ripser_coeff.rips_dm_sparse(I, J, V, I.size, N, coeff,
                                               maxHomDim, thresh, n_threads,
-                                              return_generators)
+                                              return_generators,
+                                              return_cocycles)
     return ret
 
 
@@ -310,7 +315,8 @@ def _is_prime_and_larger_than_2(x, N):
 def ripser_parallel(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
                     metric_params={}, nearest_neighbors_params={},
                     weights=None, weight_params=None, collapse_edges=False,
-                    n_threads=1, return_generators=False):
+                    n_threads=1, return_generators=False,
+                    return_cocycles=False):
     """Compute persistence diagrams from an input dense array or sparse matrix.
 
     If `X` represents a point cloud, a distance matrix will be internally
@@ -441,6 +447,10 @@ def ripser_parallel(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
         the return dictionary under the key `gens`. Cannot be ``True`` if
         `collapse_edges` is also ``True``.
 
+    return_cocycles: bool, optional, default: ``False``
+        Computed cocycles will be available in the `cocycles` value
+        of the return dictionary.
+
     Returns
     -------
     A dictionary holding the results of the computation. Keys and values are as
@@ -472,6 +482,15 @@ def ripser_parallel(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
                 Essential simplices corresponding to infinite bars in
                 dimensions 1 to maxdim, with 2 vertices (edge) for each birth.
 
+        'cocycles': list (size maxdim) of list of ndarray, optional
+            For each dimension less than `maxdim`, a list of representative
+            cocycles. Each representative cocycle in dimension ``d`` is
+            represented as a ndarray of ``(k, d + 1)`` elements. Each non zero value
+            of the cocycle is laid out in a row, first the 'd' indices of the
+            simplex and then the value of the cocycle on the simplex. The
+            indices of the simplex reference the original point cloud.
+            'cocycles' are only available if 'return_cocycles' is True.
+
     Notes
     -----
     The C++ backend and Python API for the computation of Vietoris–Rips
@@ -498,7 +517,7 @@ def ripser_parallel(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
            Python", *Journal of Open Source Software*, **3**(29), 2021;
            `DOI: 10.21105/joss.00925
            <https://doi.org/10.21105/joss.00925>`_.
-   
+
     .. [2] U. Bauer, "Ripser: efficient computation of Vietoris–Rips
            persistence barcodes", *J Appl. and Comput. Topology*, **5**, pp.
            391–423, 2021; `DOI: 10.1007/s41468-021-00071-5
@@ -609,14 +628,16 @@ def ripser_parallel(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
             thresh,
             coeff,
             n_threads,
-            return_generators
+            return_generators,
+            return_cocycles
             )
     else:
         # Only consider upper diagonal
         diagonal = np.diagonal(dm).astype(np.float32)
         DParam = squareform(dm, checks=False).astype(np.float32)
-        res = _compute_ph_vr_dense(DParam, diagonal, maxdim, thresh,
-                                   coeff, n_threads, return_generators)
+        res = _compute_ph_vr_dense(DParam, diagonal, maxdim, thresh, coeff,
+                                   n_threads, return_generators,
+                                   return_cocycles)
 
     # Unwrap persistence diagrams
     # Barcodes must match the inner type of C++ core filtration value.
@@ -645,4 +666,20 @@ def ripser_parallel(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
             ]
         ret['gens'] = (finite_0, finite_higher, essential_0, essential_higher)
 
+    # Unwrap cocycles
+    if return_cocycles:
+        nb_dim = len(res.cocycles_by_dim)
+        cocycles = []
+        for dim in range(nb_dim):
+            cocycles.append([])
+            for j in range(len(res.cocycles_by_dim[dim])):
+                ccl = res.cocycles_by_dim[dim][j]
+                n = int(len(ccl) / (dim + 2))
+                ccl = np.reshape(np.array(ccl, dtype=np.int64), [n, dim + 2])
+                ccl[:, -1] = np.mod(ccl[:, -1], coeff)
+                cocycles[dim].append(ccl)
+        cocycles = np.asarray(cocycles, dtype=object)
+
+        ret["cocycles"] = cocycles
+
     return ret

gph/python/test/test_ripser.py

@@ -436,7 +436,6 @@ def test_unsupported_coefficient():
                coeff=gph_ripser.get_max_coefficient_field_supported()+1)
 
 
-@settings(deadline=500)
 @given(dm_dense=get_dense_distance_matrices())
 def test_non_0_diagonal_internal_representation(dm_dense):
     """Checks that, when passing a full distance matrix with non-zero values in
@@ -476,3 +475,35 @@ def test_infinite_deaths_always_essential():
 
     # With this example no finite generators in dimension 1 shall be found
     assert len(gens_fin_dim1) == 0
+
+
+def test_cocycles_in_persistence_pair():
+    """Test for a simple cocycle in a persistence pair
+       The data and the results are taken from cocycles notebook of
+       ripser.py package"""
+    X = np.array([[0.05546534, -1.03551662],
+                  [-0.97905569,  0.02129887],
+                  [-0.66797382, -0.99523678],
+                  [-0.95970867, -0.0622854],
+                  [-0.61541581, -1.02874641],
+                  [0.32049913,  0.96618272],
+                  [0.49975854,  1.02644011],
+                  [1.10357502,  0.36892592],
+                  [0.79793975, -0.50982899],
+                  [-0.15508265,  1.0381044],
+                  [0.96130068, -0.04503272],
+                  [1.11553304,  0.36782575]])
+    expected = np.array([[], [np.array([[9, 1, 1],
+                                        [9, 3, 1],
+                                        [5, 1, 1],
+                                        [5, 3, 1],
+                                        [6, 1, 1],
+                                        [8, 3, 1]])
+                              ]], dtype=object
+                        )
+    cocycles = ripser(X, coeff=17, return_cocycles=True)['cocycles']
+
+    assert np.array_equal(expected.shape, cocycles.shape)
+
+    for dim in range(2):
+        assert np.array_equal(expected[dim], cocycles[dim])