Merge pull request #369 from VincentRouvreau/collapse_edges_stree_python

collapse edges for python simplex tree

6 files changed, 173 insertions, 38 deletions

diff --git a/src/python/doc/examples.rst b/src/python/doc/examples.rst
index a42227e3..76e5d4c7 100644
--- a/src/python/doc/examples.rst
+++ b/src/python/doc/examples.rst
@@ -7,27 +7,30 @@ Examples
 
 .. only:: builder_html
 
-    * :download:`rips_complex_from_points_example.py <../example/rips_complex_from_points_example.py>`
+    * :download:`alpha_complex_diagram_persistence_from_off_file_example.py <../example/alpha_complex_diagram_persistence_from_off_file_example.py>`
     * :download:`alpha_complex_from_points_example.py <../example/alpha_complex_from_points_example.py>`
-    * :download:`simplex_tree_example.py <../example/simplex_tree_example.py>`
     * :download:`alpha_rips_persistence_bottleneck_distance.py <../example/alpha_rips_persistence_bottleneck_distance.py>`
-    * :download:`tangential_complex_plain_homology_from_off_file_example.py <../example/tangential_complex_plain_homology_from_off_file_example.py>`
-    * :download:`alpha_complex_diagram_persistence_from_off_file_example.py <../example/alpha_complex_diagram_persistence_from_off_file_example.py>`
-    * :download:`periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py <../example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py>`
     * :download:`bottleneck_basic_example.py <../example/bottleneck_basic_example.py>`
-    * :download:`gudhi_graphical_tools_example.py <../example/gudhi_graphical_tools_example.py>`
-    * :download:`plot_simplex_tree_dim012.py <../example/plot_simplex_tree_dim012.py>`
-    * :download:`plot_rips_complex.py <../example/plot_rips_complex.py>`
-    * :download:`plot_alpha_complex.py <../example/plot_alpha_complex.py>`
-    * :download:`witness_complex_from_nearest_landmark_table.py <../example/witness_complex_from_nearest_landmark_table.py>`
+    * :download:`coordinate_graph_induced_complex.py <../example/coordinate_graph_induced_complex.py>`
+    * :download:`diagram_vectorizations_distances_kernels.py <../example/diagram_vectorizations_distances_kernels.py>`
     * :download:`euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py <../example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py>`
     * :download:`euclidean_witness_complex_diagram_persistence_from_off_file_example.py <../example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py>`
-    * :download:`rips_complex_diagram_persistence_from_off_file_example.py <../example/rips_complex_diagram_persistence_from_off_file_example.py>`
+    * :download:`functional_graph_induced_complex.py <../example/functional_graph_induced_complex.py>`
+    * :download:`gudhi_graphical_tools_example.py <../example/gudhi_graphical_tools_example.py>`
+    * :download:`nerve_of_a_covering.py <../example/nerve_of_a_covering.py>`
+    * :download:`periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py <../example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py>`
+    * :download:`plot_alpha_complex.py <../example/plot_alpha_complex.py>`
+    * :download:`plot_rips_complex.py <../example/plot_rips_complex.py>`
+    * :download:`plot_simplex_tree_dim012.py <../example/plot_simplex_tree_dim012.py>`
+    * :download:`random_cubical_complex_persistence_example.py <../example/random_cubical_complex_persistence_example.py>`
+    * :download:`rips_complex_diagram_persistence_from_correlation_matrix_file_example.py <../example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py>`
     * :download:`rips_complex_diagram_persistence_from_distance_matrix_file_example.py <../example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py>`
+    * :download:`rips_complex_diagram_persistence_from_off_file_example.py <../example/rips_complex_diagram_persistence_from_off_file_example.py>`
+    * :download:`rips_complex_edge_collapse_example.py <../example/rips_complex_edge_collapse_example.py>`
+    * :download:`rips_complex_from_points_example.py <../example/rips_complex_from_points_example.py>`
     * :download:`rips_persistence_diagram.py <../example/rips_persistence_diagram.py>`
+    * :download:`simplex_tree_example.py <../example/simplex_tree_example.py>`
     * :download:`sparse_rips_persistence_diagram.py <../example/sparse_rips_persistence_diagram.py>`
-    * :download:`random_cubical_complex_persistence_example.py <../example/random_cubical_complex_persistence_example.py>`
-    * :download:`coordinate_graph_induced_complex.py <../example/coordinate_graph_induced_complex.py>`
-    * :download:`functional_graph_induced_complex.py <../example/functional_graph_induced_complex.py>`
+    * :download:`tangential_complex_plain_homology_from_off_file_example.py <../example/tangential_complex_plain_homology_from_off_file_example.py>`
     * :download:`voronoi_graph_induced_complex.py <../example/voronoi_graph_induced_complex.py>`
-    * :download:`nerve_of_a_covering.py <../example/nerve_of_a_covering.py>`
+    * :download:`witness_complex_from_nearest_landmark_table.py <../example/witness_complex_from_nearest_landmark_table.py>`
diff --git a/src/python/example/rips_complex_edge_collapse_example.py b/src/python/example/rips_complex_edge_collapse_example.py
new file mode 100755
index 00000000..b26eb9fc
--- /dev/null
+++ b/src/python/example/rips_complex_edge_collapse_example.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+import gudhi
+import matplotlib.pyplot as plt
+import time
+
+""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+    Author(s):       Vincent Rouvreau
+
+    Copyright (C) 2016 Inria
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+__author__ = "Vincent Rouvreau"
+__copyright__ = "Copyright (C) 2020 Inria"
+__license__ = "MIT"
+
+
+print("#####################################################################")
+print("RipsComplex (only the one-skeleton) creation from tore3D_300.off file")
+
+off_file = gudhi.__root_source_dir__ + '/data/points/tore3D_300.off'
+point_cloud = gudhi.read_points_from_off_file(off_file = off_file)
+rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=12.0)
+simplex_tree = rips_complex.create_simplex_tree(max_dimension=1)
+print('1. Rips complex is of dimension ', simplex_tree.dimension(), ' - ',
+      simplex_tree.num_simplices(), ' simplices - ',
+      simplex_tree.num_vertices(), ' vertices.')
+
+# Expansion of this one-skeleton would require a lot of memory. Let's collapse it
+start = time.process_time()
+simplex_tree.collapse_edges()
+print('2. Rips complex is of dimension ', simplex_tree.dimension(), ' - ',
+      simplex_tree.num_simplices(), ' simplices - ',
+      simplex_tree.num_vertices(), ' vertices.')
+simplex_tree.expansion(3)
+diag = simplex_tree.persistence()
+print("Collapse, expansion and persistence computation took ", time.process_time() - start, " sec.")
+
+# Use subplots to display diagram and density side by side
+fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
+gudhi.plot_persistence_diagram(diag, axes=axes[0])
+axes[0].set_title("Persistence after 1 collapse")
+
+# Collapse can be performed several times. Let's collapse it 3 times
+start = time.process_time()
+simplex_tree.collapse_edges(nb_iterations = 3)
+print('3. Rips complex is of dimension ', simplex_tree.dimension(), ' - ',
+      simplex_tree.num_simplices(), ' simplices - ',
+      simplex_tree.num_vertices(), ' vertices.')
+simplex_tree.expansion(3)
+diag = simplex_tree.persistence()
+print("Collapse, expansion and persistence computation took ", time.process_time() - start, " sec.")
+
+gudhi.plot_persistence_diagram(diag, axes=axes[1])
+axes[1].set_title("Persistence after 3 more collapses")
+
+# Plot the 2 persistence diagrams side to side to check the persistence is the same
+plt.show()
\ No newline at end of file
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index e748ac40..75e94e0b 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -57,6 +57,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         bool make_filtration_non_decreasing() nogil
         void compute_extended_filtration() nogil
         vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) nogil
+        Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil
         # Iterators over Simplex tree
         pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) nogil
         Simplex_tree_simplices_iterator get_simplices_iterator_begin() nogil
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 20e66d9f..dfb1d985 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -69,7 +69,7 @@ cdef class SimplexTree:
         this simplicial complex, or +infinity if it is not in the complex.
 
         :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
+        :type simplex: list of int
         :returns:  The simplicial complex filtration value.
         :rtype:  float
         """
@@ -80,7 +80,7 @@ cdef class SimplexTree:
         given N-simplex.
 
         :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
+        :type simplex: list of int
         :param filtration:  The new filtration value.
         :type filtration:  float
 
@@ -153,7 +153,7 @@ cdef class SimplexTree:
         """This function sets the dimension of the simplicial complex.
 
         :param dimension: The new dimension value.
-        :type dimension: int.
+        :type dimension: int
 
         .. note::
 
@@ -172,7 +172,7 @@ cdef class SimplexTree:
         complex or not.
 
         :param simplex: The N-simplex to find, represented by a list of vertex.
-        :type simplex: list of int.
+        :type simplex: list of int
         :returns:  true if the simplex was found, false otherwise.
         :rtype:  bool
         """
@@ -186,9 +186,9 @@ cdef class SimplexTree:
 
         :param simplex: The N-simplex to insert, represented by a list of
             vertex.
-        :type simplex: list of int.
+        :type simplex: list of int
         :param filtration: The filtration value of the simplex.
-        :type filtration: float.
+        :type filtration: float
         :returns:  true if the simplex was not yet in the complex, false
             otherwise (whatever its original filtration value).
         :rtype:  bool
@@ -228,7 +228,7 @@ cdef class SimplexTree:
         """This function returns a generator with the (simplices of the) skeleton of a maximum given dimension.
 
         :param dimension: The skeleton dimension value.
-        :type dimension: int.
+        :type dimension: int
         :returns:  The (simplices of the) skeleton of a maximum dimension.
         :rtype:  generator with tuples(simplex, filtration)
         """
@@ -243,7 +243,7 @@ cdef class SimplexTree:
         """This function returns the star of a given N-simplex.
 
         :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
+        :type simplex: list of int
         :returns:  The (simplices of the) star of a simplex.
         :rtype:  list of tuples(simplex, filtration)
         """
@@ -265,10 +265,10 @@ cdef class SimplexTree:
         given codimension.
 
         :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
+        :type simplex: list of int
         :param codimension: The codimension. If codimension = 0, all cofaces
             are returned (equivalent of get_star function)
-        :type codimension: int.
+        :type codimension: int
         :returns:  The (simplices of the) cofaces of a simplex
         :rtype:  list of tuples(simplex, filtration)
         """
@@ -290,7 +290,7 @@ cdef class SimplexTree:
         complex.
 
         :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
+        :type simplex: list of int
 
         .. note::
 
@@ -308,7 +308,7 @@ cdef class SimplexTree:
         """Prune above filtration value given as parameter.
 
         :param filtration: Maximum threshold value.
-        :type filtration: float.
+        :type filtration: float
         :returns: The filtration modification information.
         :rtype: bool
 
@@ -342,7 +342,7 @@ cdef class SimplexTree:
         1 when calling the method.
 
         :param max_dim: The maximal dimension.
-        :type max_dim: int.
+        :type max_dim: int
         """
         cdef int maxdim = max_dim
         with nogil:
@@ -383,12 +383,12 @@ cdef class SimplexTree:
 
         :param homology_coeff_field: The homology coefficient field. Must be a
             prime number. Default value is 11.
-        :type homology_coeff_field: int.
+        :type homology_coeff_field: int
         :param min_persistence: The minimum persistence value (i.e., the absolute value of the difference between the persistence diagram point coordinates) to take into
             account (strictly greater than min_persistence). Default value is
             0.0.
             Sets min_persistence to -1.0 to see all values.
-        :type min_persistence: float.
+        :type min_persistence: float
         :returns: A list of four persistence diagrams in the format described in :func:`persistence`. The first one is Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. See https://link.springer.com/article/10.1007/s10208-008-9027-z and/or section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes.
 
         .. note::
@@ -415,12 +415,12 @@ cdef class SimplexTree:
 
         :param homology_coeff_field: The homology coefficient field. Must be a
             prime number. Default value is 11.
-        :type homology_coeff_field: int.
+        :type homology_coeff_field: int
         :param min_persistence: The minimum persistence value to take into
             account (strictly greater than min_persistence). Default value is
             0.0.
             Set min_persistence to -1.0 to see all values.
-        :type min_persistence: float.
+        :type min_persistence: float
         :param persistence_dim_max: If true, the persistent homology for the
             maximal dimension in the complex is computed. If false, it is
             ignored. Default is false.
@@ -438,12 +438,12 @@ cdef class SimplexTree:
 
         :param homology_coeff_field: The homology coefficient field. Must be a
             prime number. Default value is 11.
-        :type homology_coeff_field: int.
+        :type homology_coeff_field: int
         :param min_persistence: The minimum persistence value to take into
             account (strictly greater than min_persistence). Default value is
             0.0.
             Sets min_persistence to -1.0 to see all values.
-        :type min_persistence: float.
+        :type min_persistence: float
         :param persistence_dim_max: If true, the persistent homology for the
             maximal dimension in the complex is computed. If false, it is
             ignored. Default is false.
@@ -478,10 +478,10 @@ cdef class SimplexTree:
 
         :param from_value: The persistence birth limit to be added in the
             numbers (persistent birth <= from_value).
-        :type from_value: float.
+        :type from_value: float
         :param to_value: The persistence death limit to be added in the
             numbers (persistent death > to_value).
-        :type to_value: float.
+        :type to_value: float
 
         :returns: The persistent Betti numbers ([B0, B1, ..., Bn]).
         :rtype:  list of int
@@ -498,7 +498,7 @@ cdef class SimplexTree:
         complex in a specific dimension.
 
         :param dimension: The specific dimension.
-        :type dimension: int.
+        :type dimension: int
         :returns: The persistence intervals.
         :rtype:  numpy array of dimension 2
 
@@ -527,7 +527,7 @@ cdef class SimplexTree:
         complex in a user given file name.
 
         :param persistence_file: Name of the file.
-        :type persistence_file: string.
+        :type persistence_file: string
 
         :note: intervals_in_dim function requires
             :func:`compute_persistence`
@@ -581,3 +581,23 @@ cdef class SimplexTree:
             infinite0 = np_array(next(l))
             infinites = [np_array(d).reshape(-1,2) for d in l]
         return (normal0, normals, infinite0, infinites)
+
+    def collapse_edges(self, nb_iterations = 1):
+        """Assuming the simplex tree is a 1-skeleton graph, this method collapse edges (simplices of higher dimension
+        are ignored) and resets the simplex tree from the remaining edges.
+        A good candidate is to build a simplex tree on top of a :class:`~gudhi.RipsComplex` of dimension 1 before
+        collapsing edges
+        (cf. :download:`rips_complex_edge_collapse_example.py <../example/rips_complex_edge_collapse_example.py>`).
+        For implementation details, please refer to :cite:`edgecollapsesocg2020`.
+
+        :param nb_iterations: The number of edge collapse iterations to perform. Default is 1.
+        :type nb_iterations: int
+        """
+        # Backup old pointer
+        cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr()
+        cdef int nb_iter = nb_iterations
+        with nogil:
+            # New pointer is a new collapsed simplex tree
+            self.thisptr = <intptr_t>(ptr.collapse_edges(nb_iter))
+            # Delete old pointer
+            del ptr
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 56d7c41d..e288a8cf 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -15,10 +15,13 @@
 #include <gudhi/distance_functions.h>
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Points_off_io.h>
+#include <gudhi/Flag_complex_edge_collapser.h>
 
 #include <iostream>
 #include <vector>
 #include <utility>  // std::pair
+#include <tuple>
+#include <iterator>  // for std::distance
 
 namespace Gudhi {
 
@@ -157,6 +160,34 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return new_dgm;
   }
 
+  Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) {
+    using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>;
+    std::vector<Filtered_edge> edges;
+    for (Simplex_handle sh : Base::skeleton_simplex_range(1)) {
+      if (Base::dimension(sh) == 1) {
+        typename Base::Simplex_vertex_range rg = Base::simplex_vertex_range(sh);
+        auto vit = rg.begin();
+        Vertex_handle v = *vit;
+        Vertex_handle w = *++vit;
+        edges.emplace_back(v, w, Base::filtration(sh));
+      }
+    }
+
+    for (int iteration = 0; iteration < nb_collapse_iteration; iteration++) {
+      edges = Gudhi::collapse::flag_complex_collapse_edges(edges);
+    }
+    Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface();
+    // Copy the original 0-skeleton
+    for (Simplex_handle sh : Base::skeleton_simplex_range(0)) {
+      collapsed_stree_ptr->insert({*(Base::simplex_vertex_range(sh).begin())}, Base::filtration(sh));
+    }
+    // Insert remaining edges
+    for (auto remaining_edge : edges) {
+      collapsed_stree_ptr->insert({std::get<0>(remaining_edge), std::get<1>(remaining_edge)}, std::get<2>(remaining_edge));
+    }
+    return collapsed_stree_ptr;
+  }
+
   // Iterator over the simplex tree
   Complex_simplex_iterator get_simplices_iterator_begin() {
     // this specific case works because the range is just a pair of iterators - won't work if range was a vector
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 2137d822..83be0602 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -340,3 +340,21 @@ def test_simplices_iterator():
         assert st.find(simplex[0]) == True
         print("filtration is: ", simplex[1])
         assert st.filtration(simplex[0]) == simplex[1]
+
+def test_collapse_edges():
+    st = SimplexTree()
+    
+    assert st.insert([0, 1], filtration=1.0) == True
+    assert st.insert([1, 2], filtration=1.0) == True
+    assert st.insert([2, 3], filtration=1.0) == True
+    assert st.insert([0, 3], filtration=1.0) == True
+    assert st.insert([0, 2], filtration=2.0) == True
+    assert st.insert([1, 3], filtration=2.0) == True
+
+    assert st.num_simplices() == 10
+
+    st.collapse_edges()
+    assert st.num_simplices() == 9
+    assert st.find([1, 3]) == False
+    for simplex in st.get_skeleton(0): 
+        assert simplex[1] == 1.