Improve ripser bindings performance (#501)

* Update ripser & collapser bindings to allow pass by reference with numpy

* Update making the distance matrix triangular more memory friendly

* Remove unnecessary dict inherited from ripser.py cython

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>

gtda/externals/bindings/collapser_bindings.cpp

@@ -6,14 +6,15 @@
 #include <gudhi/Flag_complex_edge_collapser.h>
 
 #include <pybind11/eigen.h>
+#include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
 namespace py = pybind11;
 
 /* GUDHI Collapser required types */
 using Filtration_value = float;
-using Vertex_handle = uint32_t;
+using Vertex_handle = int32_t;
 using Filtered_edge =
     std::tuple<Vertex_handle, Vertex_handle, Filtration_value>;
 using Filtered_edge_list = std::vector<Filtered_edge>;
@@ -29,7 +30,7 @@ using Filtration_values = std::vector<Filtration_value>;
 using COO_data = std::tuple<Row_idx, Col_idx, Filtration_values>;
 
 /* Dense matrix input types */
-using Distance_matrix = std::vector<std::vector<Filtration_value>>;
+using Distance_matrix_np = py::array_t<Filtration_value>;
 
 /* constants */
 const Filtration_value filtration_max =
@@ -86,13 +87,18 @@ PYBIND11_MODULE(gtda_collapser, m) {
         "collapses edges while preserving the persistent homology");
 
   m.def("flag_complex_collapse_edges_coo",
-        [](Row_idx& row, Col_idx& col, Filtration_values& data,
+        [](py::array_t<Vertex_handle>& row_, py::array_t<Vertex_handle>& col_,
+           py::array_t<Filtration_value>& data_,
            Filtration_value thresh = filtration_max) {
           Filtered_edge_list graph;
 
+          Vertex_handle* row = (Vertex_handle*)row_.request().ptr;
+          Vertex_handle* col = (Vertex_handle*)col_.request().ptr;
+          Filtration_value* data = (Filtration_value*)data_.request().ptr;
+
           /* Convert from COO input format to Filtered_edge_list */
           /* Applying threshold to the input data */
-          int size = data.size();
+          int size = data_.size();
           for (size_t k = 0; k < size; ++k)
             if (data[k] <= thresh)
               graph.push_back(Filtered_edge(row[k], col[k], data[k]));
@@ -108,15 +114,15 @@ PYBIND11_MODULE(gtda_collapser, m) {
         "collapses edges while preserving the persistent homology");
 
   m.def("flag_complex_collapse_edges_dense",
-        [](Distance_matrix& dm, Filtration_value thresh = filtration_max) {
+        [](Distance_matrix_np& dm, Filtration_value thresh = filtration_max) {
           Filtered_edge_list graph;
 
           /* Convert from dense format to Filtered edge list */
           /* Applying threshold to the input data */
-          for (size_t i = 0; i < dm.size(); i++)
-            for (size_t j = 0; j < dm[i].size(); j++)
-              if (j > i && (dm[i][j] <= thresh))
-                graph.push_back(Filtered_edge(i, j, dm[i][j]));
+          for (size_t i = 0; i < dm.shape(0); i++)
+            for (size_t j = 0; j < dm.shape(1); j++)
+              if (j > i && (*(dm.data(i, j)) <= thresh))
+                graph.push_back(Filtered_edge(i, j, *(dm.data(i, j))));
 
           /* Start collapser */
           auto vec_triples =

gtda/externals/bindings/ripser_bindings.cpp

@@ -7,8 +7,10 @@
 #include <ripser.cpp>
 
 // PYBIND11
+#include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
 
@@ -32,21 +34,22 @@ PYBIND11_MODULE(gtda_ripser, m) {
       .def_readwrite("num_edges", &ripserResults::num_edges);
 
   m.def("rips_dm",
-        [](std::vector<float> D, int N, int modulus, int dim_max,
+        [](py::array_t<float>& D, int N, int modulus, int dim_max,
            float threshold, int do_cocycles) {
-          ripserResults ret =
-              rips_dm(&D[0], N, modulus, dim_max, threshold, do_cocycles);
+          ripserResults ret = rips_dm((float*)D.request().ptr, N, modulus,
+                                      dim_max, threshold, do_cocycles);
           return ret;
         },
         "D"_a, "N"_a, "modulus"_a, "dim_max"_a, "threshold"_a, "do_cocycles"_a,
         "ripser distance matrix");
   m.def("rips_dm_sparse",
-        [](std::vector<int> I, std::vector<int> J, std::vector<float> V,
+        [](py::array_t<int>& I, py::array_t<int>& J, py::array_t<float>& V,
            int NEdges, int N, int modulus, int dim_max, float threshold,
            int do_cocycles) {
           ripserResults ret =
-              rips_dm_sparse(&I[0], &J[0], &V[0], NEdges, N, modulus, dim_max,
-                             threshold, do_cocycles);
+              rips_dm_sparse((int*)I.request().ptr, (int*)J.request().ptr,
+                             (float*)V.request().ptr, NEdges, N, modulus,
+                             dim_max, threshold, do_cocycles);
           return ret;
         },
         "I"_a, "J"_a, "V"_a, "NEdges"_a, "N"_a, "modulus"_a, "dim_max"_a,

gtda/externals/python/ripser_interface.py

@@ -21,10 +21,7 @@ def DRFDM(DParam, maxHomDim, thresh=-1, coeff=2, do_cocycles=0):
     else:
         ret = gtda_ripser_coeff.rips_dm(DParam, DParam.shape[0], coeff,
                                         maxHomDim, thresh, do_cocycles)
-    ret_rips = {}
-    ret_rips.update({"births_and_deaths_by_dim": ret.births_and_deaths_by_dim})
-    ret_rips.update({"num_edges": ret.num_edges})
-    return ret_rips
+    return ret
 
 
 def DRFDMSparse(I, J, V, N, maxHomDim, thresh=-1, coeff=2, do_cocycles=0):
@@ -34,10 +31,7 @@ def DRFDMSparse(I, J, V, N, maxHomDim, thresh=-1, coeff=2, do_cocycles=0):
     else:
         ret = gtda_ripser_coeff.rips_dm_sparse(I, J, V, I.size, N, coeff,
                                                maxHomDim, thresh, do_cocycles)
-    ret_rips = {}
-    ret_rips.update({"births_and_deaths_by_dim": ret.births_and_deaths_by_dim})
-    ret_rips.update({"num_edges": ret.num_edges})
-    return ret_rips
+    return ret
 
 
 def dpoint2pointcloud(X, i, metric):
@@ -295,19 +289,19 @@ def ripser(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
             )
     else:
         # Only consider strict upper diagonal
-        idx = np.triu_indices(n_points, 1)
-        DParam = dm[idx].astype(np.float32)
+        DParam = dm[np.invert(np.tri(n_points, k=0, dtype=np.bool))].astype(
+            np.float32).flatten()
         res = DRFDM(DParam, maxdim, thresh, coeff)
 
     # Unwrap persistence diagrams
-    dgms = res["births_and_deaths_by_dim"]
+    dgms = res.births_and_deaths_by_dim
     for dim in range(len(dgms)):
         N = int(len(dgms[dim]) / 2)
         dgms[dim] = np.reshape(np.array(dgms[dim]), [N, 2])
 
     ret = {
         "dgms": dgms,
-        "num_edges": res["num_edges"],
+        "num_edges": res.num_edges,
         "dperm2all": dperm2all,
         "idx_perm": idx_perm,
         "r_cover": r_cover,