Merge branch 'master' of https://github.com/GUDHI/gudhi-devel into wasserstein_representations

5 files changed, 222 insertions, 22 deletions

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 055d5b23..4771cef9 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -452,6 +452,7 @@ if(PYTHONINTERP_FOUND)
              ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/example/simplex_tree_example.py)
 
     add_gudhi_py_test(test_simplex_tree)
+    add_gudhi_py_test(test_simplex_generators)
 
     # Witness
     add_test(NAME witness_complex_from_nearest_landmark_table_py_test
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 5dea2449..1d4ed926 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -77,3 +77,5 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
         vector[pair[double,double]] intervals_in_dimension(int dimension)
         void write_output_diagram(string diagram_file_name) except +
         vector[pair[vector[int], vector[int]]] persistence_pairs()
+        pair[vector[vector[int]], vector[vector[int]]] lower_star_generators()
+        pair[vector[vector[int]], vector[vector[int]]] flag_generators()
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 9479118a..55115cca 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -9,6 +9,7 @@
 
 from cython.operator import dereference, preincrement
 from libc.stdint cimport intptr_t
+import numpy
 from numpy import array as np_array
 cimport simplex_tree
 
@@ -415,7 +416,7 @@ cdef class SimplexTree:
         :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.
+            Set min_persistence to -1.0 to see all values.
         :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
@@ -527,3 +528,49 @@ cdef class SimplexTree:
         """
         assert self.pcohptr != NULL, "compute_persistence() must be called before write_persistence_diagram()"
         self.pcohptr.write_output_diagram(persistence_file.encode('utf-8'))
+
+    def lower_star_persistence_generators(self):
+        """Assuming this is a lower-star filtration, this function returns the persistence pairs,
+        where each simplex is replaced with the vertex that gave it its filtration value.
+
+        :returns: First the regular persistence pairs, grouped by dimension, with one vertex per extremity,
+            and second the essential features, grouped by dimension, with one vertex each
+        :rtype: Tuple[List[numpy.array[int] of shape (n,2)], List[numpy.array[int] of shape (m,)]]
+
+        :note: lower_star_persistence_generators requires that `persistence()` be called first.
+        """
+        assert self.pcohptr != NULL, "lower_star_persistence_generators() requires that persistence() be called first."
+        gen = self.pcohptr.lower_star_generators()
+        normal = [np_array(d).reshape(-1,2) for d in gen.first]
+        infinite = [np_array(d) for d in gen.second]
+        return (normal, infinite)
+
+    def flag_persistence_generators(self):
+        """Assuming this is a flag complex, this function returns the persistence pairs,
+        where each simplex is replaced with the vertices of the edges that gave it its filtration value.
+
+        :returns: First the regular persistence pairs of dimension 0, with one vertex for birth and two for death;
+            then the other regular persistence pairs, grouped by dimension, with 2 vertices per extremity;
+            then the connected components, with one vertex each;
+            finally the other essential features, grouped by dimension, with 2 vertices for birth.
+        :rtype: Tuple[numpy.array[int] of shape (n,3), List[numpy.array[int] of shape (m,4)], numpy.array[int] of shape (l,), List[numpy.array[int] of shape (k,2)]]
+
+        :note: flag_persistence_generators requires that `persistence()` be called first.
+        """
+        assert self.pcohptr != NULL, "flag_persistence_generators() requires that persistence() be called first."
+        gen = self.pcohptr.flag_generators()
+        if len(gen.first) == 0:
+            normal0 = numpy.empty((0,3))
+            normals = []
+        else:
+            l = iter(gen.first)
+            normal0 = np_array(next(l)).reshape(-1,3)
+            normals = [np_array(d).reshape(-1,4) for d in l]
+        if len(gen.second) == 0:
+            infinite0 = numpy.empty(0)
+            infinites = []
+        else:
+            l = iter(gen.second)
+            infinite0 = np_array(next(l))
+            infinites = [np_array(d).reshape(-1,2) for d in l]
+        return (normal0, normals, infinite0, infinites)
diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h
index e2b69a52..0de9bd5c 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -28,18 +28,17 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
    * Compare two intervals by dimension, then by length.
    */
   struct cmp_intervals_by_dim_then_length {
-    explicit cmp_intervals_by_dim_then_length(FilteredComplex * sc)
-        : sc_(sc) { }
-
     template<typename Persistent_interval>
     bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
-      if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2)))
-        return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
-                > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
+      if (std::get<0>(p1) == std::get<0>(p2)) {
+        auto& i1 = std::get<1>(p1);
+        auto& i2 = std::get<1>(p2);
+        return std::get<1>(i1) - std::get<0>(i1) > std::get<1>(i2) - std::get<0>(i2);
+      }
       else
-        return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2)));
+        return (std::get<0>(p1) > std::get<0>(p2));
+        // Why does this sort by decreasing dimension?
     }
-    FilteredComplex* sc_;
   };
 
  public:
@@ -54,45 +53,132 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
   }
 
   std::vector<std::pair<int, std::pair<double, double>>> get_persistence() {
-    // Custom sort and output persistence
-    cmp_intervals_by_dim_then_length cmp(stptr_);
-    auto persistent_pairs = Base::get_persistent_pairs();
-    std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
-
     std::vector<std::pair<int, std::pair<double, double>>> persistence;
+    auto const& persistent_pairs = Base::get_persistent_pairs();
+    persistence.reserve(persistent_pairs.size());
     for (auto pair : persistent_pairs) {
-      persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)),
-                                           std::make_pair(stptr_->filtration(get<0>(pair)),
-                                                          stptr_->filtration(get<1>(pair)))));
+      persistence.emplace_back(stptr_->dimension(get<0>(pair)),
+                               std::make_pair(stptr_->filtration(get<0>(pair)),
+                                              stptr_->filtration(get<1>(pair))));
     }
+    // Custom sort and output persistence
+    cmp_intervals_by_dim_then_length cmp;
+    std::sort(std::begin(persistence), std::end(persistence), cmp);
     return persistence;
   }
 
   std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() {
-    auto pairs = Base::get_persistent_pairs();
-
     std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs;
+    auto const& pairs = Base::get_persistent_pairs();
     persistence_pairs.reserve(pairs.size());
+    std::vector<int> birth;
+    std::vector<int> death;
     for (auto pair : pairs) {
-      std::vector<int> birth;
+      birth.clear();
       if (get<0>(pair) != stptr_->null_simplex()) {
         for (auto vertex : stptr_->simplex_vertex_range(get<0>(pair))) {
           birth.push_back(vertex);
         }
       }
 
-      std::vector<int> death;
+      death.clear();
       if (get<1>(pair) != stptr_->null_simplex()) {
+        death.reserve(birth.size()+1);
         for (auto vertex : stptr_->simplex_vertex_range(get<1>(pair))) {
           death.push_back(vertex);
         }
       }
 
-      persistence_pairs.push_back(std::make_pair(birth, death));
+      persistence_pairs.emplace_back(birth, death);
     }
     return persistence_pairs;
   }
 
+  // TODO: (possibly at the python level)
+  // - an option to return only some of those vectors?
+  typedef std::pair<std::vector<std::vector<int>>, std::vector<std::vector<int>>> Generators;
+
+  Generators lower_star_generators() {
+    Generators out;
+    // diags[i] should be interpreted as vector<array<int,2>>
+    auto& diags = out.first;
+    // diagsinf[i] should be interpreted as vector<int>
+    auto& diagsinf = out.second;
+    for (auto pair : Base::get_persistent_pairs()) {
+      auto s = std::get<0>(pair);
+      auto t = std::get<1>(pair);
+      int dim = stptr_->dimension(s);
+      auto v = stptr_->vertex_with_same_filtration(s);
+      if(t == stptr_->null_simplex()) {
+        while(diagsinf.size() < dim+1) diagsinf.emplace_back();
+        diagsinf[dim].push_back(v);
+      } else {
+        while(diags.size() < dim+1) diags.emplace_back();
+        auto w = stptr_->vertex_with_same_filtration(t);
+        auto& d = diags[dim];
+        d.insert(d.end(), { v, w });
+      }
+    }
+    return out;
+  }
+
+  // An alternative, to avoid those different sizes, would be to "pad" vertex generator v as (v, v) or (v, -1). When using it as index, this corresponds to adding the vertex filtration values either on the diagonal of the distance matrix, or as an extra row or column.
+  // We could also merge the vectors for different dimensions into a single one, with an extra column for the dimension (converted to type double).
+  Generators flag_generators() {
+    Generators out;
+    // diags[0] should be interpreted as vector<array<int,3>> and other diags[i] as vector<array<int,4>>
+    auto& diags = out.first;
+    // diagsinf[0] should be interpreted as vector<int> and other diagsinf[i] as vector<array<int,2>>
+    auto& diagsinf = out.second;
+    for (auto pair : Base::get_persistent_pairs()) {
+      auto s = std::get<0>(pair);
+      auto t = std::get<1>(pair);
+      int dim = stptr_->dimension(s);
+      bool infinite = t == stptr_->null_simplex();
+      if(infinite) {
+        if(dim == 0) {
+          auto v = *std::begin(stptr_->simplex_vertex_range(s));
+          if(diagsinf.size()==0)diagsinf.emplace_back();
+          diagsinf[0].push_back(v);
+        } else {
+          auto e = stptr_->edge_with_same_filtration(s);
+          auto&& e_vertices = stptr_->simplex_vertex_range(e);
+          auto i = std::begin(e_vertices);
+          auto v1 = *i;
+          auto v2 = *++i;
+          GUDHI_CHECK(++i==std::end(e_vertices), "must be an edge");
+          while(diagsinf.size() < dim+1) diagsinf.emplace_back();
+          auto& d = diagsinf[dim];
+          d.insert(d.end(), { v1, v2 });
+        }
+      } else {
+        auto et = stptr_->edge_with_same_filtration(t);
+        auto&& et_vertices = stptr_->simplex_vertex_range(et);
+        auto it = std::begin(et_vertices);
+        auto w1 = *it;
+        auto w2 = *++it;
+        GUDHI_CHECK(++it==std::end(et_vertices), "must be an edge");
+        if(dim == 0) {
+          auto v = *std::begin(stptr_->simplex_vertex_range(s));
+          if(diags.size()==0)diags.emplace_back();
+          auto& d = diags[0];
+          d.insert(d.end(), { v, w1, w2 });
+        } else {
+          auto es = stptr_->edge_with_same_filtration(s);
+          auto&& es_vertices = stptr_->simplex_vertex_range(es);
+          auto is = std::begin(es_vertices);
+          auto v1 = *is;
+          auto v2 = *++is;
+          GUDHI_CHECK(++is==std::end(es_vertices), "must be an edge");
+          while(diags.size() < dim+1) diags.emplace_back();
+          auto& d = diags[dim];
+          d.insert(d.end(), { v1, v2, w1, w2 });
+        }
+      }
+    }
+    return out;
+  }
+
  private:
   // A copy
   FilteredComplex* stptr_;
diff --git a/src/python/test/test_simplex_generators.py b/src/python/test/test_simplex_generators.py
new file mode 100755
index 00000000..8a9b4844
--- /dev/null
+++ b/src/python/test/test_simplex_generators.py
@@ -0,0 +1,64 @@
+""" 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):       Marc Glisse
+
+    Copyright (C) 2020 Inria
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+import gudhi
+import numpy as np
+
+
+def test_flag_generators():
+    pts = np.array([[0, 0], [0, 1.01], [1, 0], [1.02, 1.03], [100, 0], [100, 3.01], [103, 0], [103.02, 3.03]])
+    r = gudhi.RipsComplex(pts, max_edge_length=4)
+    st = r.create_simplex_tree(max_dimension=50)
+    st.persistence()
+    g = st.flag_persistence_generators()
+    assert np.array_equal(g[0], [[2, 2, 0], [1, 1, 0], [3, 3, 1], [6, 6, 4], [5, 5, 4], [7, 7, 5]])
+    assert len(g[1]) == 1
+    assert np.array_equal(g[1][0], [[3, 2, 2, 1]])
+    assert np.array_equal(g[2], [0, 4])
+    assert len(g[3]) == 1
+    assert np.array_equal(g[3][0], [[7, 6]])
+    # Compare trivial cases (where the simplex is the generator) with persistence_pairs.
+    # This still makes assumptions on the order of vertices in a simplex and could be more robust.
+    pairs = st.persistence_pairs()
+    assert {tuple(i) for i in g[0]} == {(i[0][0],) + tuple(i[1]) for i in pairs if len(i[0]) == 1 and len(i[1]) != 0}
+    assert {(i[0], i[1]) for i in g[1][0]} == {tuple(i[0]) for i in pairs if len(i[0]) == 2 and len(i[1]) != 0}
+    assert set(g[2]) == {i[0][0] for i in pairs if len(i[0]) == 1 and len(i[1]) == 0}
+    assert {(i[0], i[1]) for i in g[3][0]} == {tuple(i[0]) for i in pairs if len(i[0]) == 2 and len(i[1]) == 0}
+
+
+def test_lower_star_generators():
+    st = gudhi.SimplexTree()
+    st.insert([0, 1, 2], -10)
+    st.insert([0, 3], -10)
+    st.insert([1, 3], -10)
+    st.assign_filtration([2], -1)
+    st.assign_filtration([3], 0)
+    st.assign_filtration([0], 1)
+    st.assign_filtration([1], 2)
+    st.make_filtration_non_decreasing()
+    st.persistence(min_persistence=-1)
+    g = st.lower_star_persistence_generators()
+    assert len(g[0]) == 2
+    assert np.array_equal(g[0][0], [[0, 0], [3, 0], [1, 1]])
+    assert np.array_equal(g[0][1], [[1, 1]])
+    assert len(g[1]) == 2
+    assert np.array_equal(g[1][0], [2])
+    assert np.array_equal(g[1][1], [1])
+
+
+def test_empty():
+    st = gudhi.SimplexTree()
+    st.persistence()
+    assert st.lower_star_persistence_generators() == ([], [])
+    g = st.flag_persistence_generators()
+    assert np.array_equal(g[0], np.empty((0, 3)))
+    assert g[1] == []
+    assert np.array_equal(g[2], [])
+    assert g[3] == []