13 files changed, 564 insertions, 110 deletions

diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h
new file mode 100644
index 00000000..41eb72c1
--- /dev/null
+++ b/src/python/include/Alpha_complex_factory.h
@@ -0,0 +1,156 @@
+/*    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) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef INCLUDE_ALPHA_COMPLEX_FACTORY_H_
+#define INCLUDE_ALPHA_COMPLEX_FACTORY_H_
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Alpha_complex.h>
+#include <gudhi/Alpha_complex_3d.h>
+#include <gudhi/Alpha_complex_options.h>
+#include <CGAL/Epeck_d.h>
+#include <CGAL/Epick_d.h>
+
+#include <boost/range/adaptor/transformed.hpp>
+
+#include "Simplex_tree_interface.h"
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <memory>  // for std::unique_ptr
+
+namespace Gudhi {
+
+namespace alpha_complex {
+
+// template Functor that transforms a CGAL point to a vector of double as expected by cython
+template<typename CgalPointType, bool Weighted>
+struct Point_cgal_to_cython;
+
+// Specialized Unweighted Functor
+template<typename CgalPointType>
+struct Point_cgal_to_cython<CgalPointType, false> {
+  std::vector<double> operator()(CgalPointType const& point) const 
+  {
+    std::vector<double> vd;
+    vd.reserve(point.dimension());
+    for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++)
+      vd.push_back(CGAL::to_double(*coord));
+    return vd;
+  }
+};
+
+// Specialized Weighted Functor
+template<typename CgalPointType>
+struct Point_cgal_to_cython<CgalPointType, true> {
+  std::vector<double> operator()(CgalPointType const& weighted_point) const 
+  {
+    const auto& point = weighted_point.point();
+    return Point_cgal_to_cython<decltype(point), false>()(point);
+  }
+};
+
+// Function that transforms a cython point (aka. a vector of double) to a CGAL point
+template <typename CgalPointType>
+static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) {
+  return CgalPointType(vec.size(), vec.begin(), vec.end());
+}
+
+class Abstract_alpha_complex {
+ public:
+  virtual std::vector<double> get_point(int vh) = 0;
+
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) = 0;
+  
+  virtual std::size_t num_vertices() const = 0;
+
+  virtual ~Abstract_alpha_complex() = default;
+};
+
+template <bool Weighted = false>
+class Exact_alpha_complex_dD final : public Abstract_alpha_complex {
+ private:
+  using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+  using Bare_point = typename Kernel::Point_d;
+  using Point = std::conditional_t<Weighted, typename Kernel::Weighted_point_d,
+                                             typename Kernel::Point_d>;
+
+ public:
+  Exact_alpha_complex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
+    : exact_version_(exact_version),
+      alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>)) {
+  }
+
+  Exact_alpha_complex_dD(const std::vector<std::vector<double>>& points,
+                           const std::vector<double>& weights, bool exact_version)
+    : exact_version_(exact_version),
+      alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>), weights) {
+  }
+
+  virtual std::vector<double> get_point(int vh) override {
+    // Can be a Weighted or a Bare point in function of Weighted
+    return Point_cgal_to_cython<Point, Weighted>()(alpha_complex_.get_point(vh));
+  }
+
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) override {
+    return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value);
+  }
+
+  virtual std::size_t num_vertices() const override {
+    return alpha_complex_.num_vertices();
+  }
+
+ private:
+  bool exact_version_;
+  Alpha_complex<Kernel, Weighted> alpha_complex_;
+};
+
+template <bool Weighted = false>
+class Inexact_alpha_complex_dD final : public Abstract_alpha_complex {
+ private:
+  using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+  using Bare_point = typename Kernel::Point_d;
+  using Point = std::conditional_t<Weighted, typename Kernel::Weighted_point_d,
+                                             typename Kernel::Point_d>;
+
+ public:
+  Inexact_alpha_complex_dD(const std::vector<std::vector<double>>& points)
+    : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>)) {
+  }
+
+  Inexact_alpha_complex_dD(const std::vector<std::vector<double>>& points, const std::vector<double>& weights)
+    : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>), weights) {
+  }
+
+  virtual std::vector<double> get_point(int vh) override {
+    // Can be a Weighted or a Bare point in function of Weighted
+    return Point_cgal_to_cython<Point, Weighted>()(alpha_complex_.get_point(vh));
+  }
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) override {
+    return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, false, default_filtration_value);
+  }
+
+  virtual std::size_t num_vertices() const override {
+    return alpha_complex_.num_vertices();
+  }
+
+ private:
+  Alpha_complex<Kernel, Weighted> alpha_complex_;
+};
+
+}  // namespace alpha_complex
+
+}  // namespace Gudhi
+
+#endif  // INCLUDE_ALPHA_COMPLEX_FACTORY_H_
diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h
index b3553d32..469b91ce 100644
--- a/src/python/include/Alpha_complex_interface.h
+++ b/src/python/include/Alpha_complex_interface.h
@@ -11,56 +11,64 @@
 #ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_H_
 #define INCLUDE_ALPHA_COMPLEX_INTERFACE_H_
 
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Alpha_complex.h>
-#include <CGAL/Epick_d.h>
+#include "Alpha_complex_factory.h"
+#include <gudhi/Alpha_complex_options.h>
 
 #include "Simplex_tree_interface.h"
 
 #include <iostream>
 #include <vector>
 #include <string>
+#include <memory>  // for std::unique_ptr
 
 namespace Gudhi {
 
 namespace alpha_complex {
 
 class Alpha_complex_interface {
-  using Dynamic_kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
-  using Point_d = Dynamic_kernel::Point_d;
-
  public:
-  Alpha_complex_interface(const std::vector<std::vector<double>>& points) {
-    alpha_complex_ = new Alpha_complex<Dynamic_kernel>(points);
+  Alpha_complex_interface(const std::vector<std::vector<double>>& points,
+                          const std::vector<double>& weights,
+                          bool fast_version, bool exact_version) {
+    const bool weighted = (weights.size() > 0);
+    if (fast_version) {
+      if (weighted) {
+        alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD<true>>(points, weights);
+      } else {
+        alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD<false>>(points);
+      }
+    } else {
+      if (weighted) {
+        alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD<true>>(points, weights, exact_version);
+      } else {
+        alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD<false>>(points, exact_version);
+      }
+    }
   }
 
-  Alpha_complex_interface(const std::string& off_file_name, bool from_file = true) {
-    alpha_complex_ = new Alpha_complex<Dynamic_kernel>(off_file_name);
+  std::vector<double> get_point(int vh) {
+    return alpha_ptr_->get_point(vh);
   }
 
-  ~Alpha_complex_interface() {
-    delete alpha_complex_;
+  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                           bool default_filtration_value) {
+    // Nothing to be done in case of an empty point set
+    if (alpha_ptr_->num_vertices() > 0)
+      alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value);
   }
 
-  std::vector<double> get_point(int vh) {
-    std::vector<double> vd;
-    try {
-      Point_d ph = alpha_complex_->get_point(vh);
-      for (auto coord = ph.cartesian_begin(); coord < ph.cartesian_end(); coord++)
-        vd.push_back(*coord);
-    } catch (std::out_of_range const&) {
-      // std::out_of_range is thrown in case not found. Other exceptions must be re-thrown
-    }
-    return vd;
+  static void set_float_relative_precision(double precision) {
+    // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h
+    CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::set_relative_precision_of_to_double(precision);
   }
 
-  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) {
-    alpha_complex_->create_complex(*simplex_tree, max_alpha_square);
-    simplex_tree->initialize_filtration();
+  static double get_float_relative_precision() {
+    // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h
+    return CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::get_relative_precision_of_to_double();
   }
 
  private:
-  Alpha_complex<Dynamic_kernel>* alpha_complex_;
+  std::unique_ptr<Abstract_alpha_complex> alpha_ptr_;
 };
 
 }  // namespace alpha_complex
diff --git a/src/python/include/Euclidean_strong_witness_complex_interface.h b/src/python/include/Euclidean_strong_witness_complex_interface.h
index c1c72737..f94c51ef 100644
--- a/src/python/include/Euclidean_strong_witness_complex_interface.h
+++ b/src/python/include/Euclidean_strong_witness_complex_interface.h
@@ -50,12 +50,10 @@ class Euclidean_strong_witness_complex_interface {
   void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square,
                            std::size_t limit_dimension) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension);
-    simplex_tree->initialize_filtration();
   }
 
   void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square);
-    simplex_tree->initialize_filtration();
   }
 
   std::vector<double> get_point(unsigned vh) {
diff --git a/src/python/include/Euclidean_witness_complex_interface.h b/src/python/include/Euclidean_witness_complex_interface.h
index 5d7dbdc2..4411ae79 100644
--- a/src/python/include/Euclidean_witness_complex_interface.h
+++ b/src/python/include/Euclidean_witness_complex_interface.h
@@ -49,12 +49,10 @@ class Euclidean_witness_complex_interface {
 
   void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square, std::size_t limit_dimension) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension);
-    simplex_tree->initialize_filtration();
   }
 
   void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square);
-    simplex_tree->initialize_filtration();
   }
 
   std::vector<double> get_point(unsigned vh) {
diff --git a/src/python/include/Nerve_gic_interface.h b/src/python/include/Nerve_gic_interface.h
index 5e7f8ae6..ab14c318 100644
--- a/src/python/include/Nerve_gic_interface.h
+++ b/src/python/include/Nerve_gic_interface.h
@@ -29,7 +29,6 @@ class Nerve_gic_interface : public Cover_complex<std::vector<double>> {
  public:
   void create_simplex_tree(Simplex_tree_interface<>* simplex_tree) {
     create_complex(*simplex_tree);
-    simplex_tree->initialize_filtration();
   }
   void set_cover_from_Euclidean_Voronoi(int m) {
     set_cover_from_Voronoi(Gudhi::Euclidean_distance(), m);
diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h
index 8c79e6f3..945378a0 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -12,10 +12,14 @@
 #define INCLUDE_PERSISTENT_COHOMOLOGY_INTERFACE_H_
 
 #include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Simplex_tree.h>  // for Extended_simplex_type
 
+
+#include <cstdlib>
 #include <vector>
 #include <utility>  // for std::pair
 #include <algorithm>  // for sort
+#include <unordered_map>
 
 namespace Gudhi {
 
@@ -23,82 +27,242 @@ template<class FilteredComplex>
 class Persistent_cohomology_interface : public
 persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp> {
  private:
+   typedef persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp> Base;
   /*
    * 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:
-  Persistent_cohomology_interface(FilteredComplex* stptr)
-      : persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp>(*stptr),
-      stptr_(stptr) { }
-
-  Persistent_cohomology_interface(FilteredComplex* stptr, bool persistence_dim_max)
-      : persistent_cohomology::Persistent_cohomology<FilteredComplex,
-          persistent_cohomology::Field_Zp>(*stptr, persistence_dim_max),
+  Persistent_cohomology_interface(FilteredComplex* stptr, bool persistence_dim_max=false)
+      : Base(*stptr, persistence_dim_max),
         stptr_(stptr) { }
 
-  std::vector<std::pair<int, std::pair<double, double>>> get_persistence(int homology_coeff_field,
-                                                                         double min_persistence) {
-    persistent_cohomology::Persistent_cohomology<FilteredComplex,
-      persistent_cohomology::Field_Zp>::init_coefficients(homology_coeff_field);
-    persistent_cohomology::Persistent_cohomology<FilteredComplex,
-      persistent_cohomology::Field_Zp>::compute_persistent_cohomology(min_persistence);
-
-    // Custom sort and output persistence
-    cmp_intervals_by_dim_then_length cmp(stptr_);
-    auto persistent_pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex,
-      persistent_cohomology::Field_Zp>::get_persistent_pairs();
-    std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
+  // TODO: move to the constructors?
+  void compute_persistence(int homology_coeff_field, double min_persistence) {
+    Base::init_coefficients(homology_coeff_field);
+    Base::compute_persistent_cohomology(min_persistence);
+  }
 
+  std::vector<std::pair<int, std::pair<double, double>>> get_persistence() {
     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 = persistent_cohomology::Persistent_cohomology<FilteredComplex,
+  // This function computes the top-dimensional cofaces associated to the positive and negative
+  // simplices of a cubical complex. The output format is a vector of vectors of three integers,
+  // which are [homological dimension, index of top-dimensional coface of positive simplex,
+  // index of top-dimensional coface of negative simplex]. If the topological feature is essential,
+  // then the index of top-dimensional coface of negative simplex is arbitrarily set to -1.
+  std::vector<std::vector<int>> cofaces_of_cubical_persistence_pairs() {
+
+    // Warning: this function is meant to be used with CubicalComplex only!!
+
+    auto&& pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex,
       persistent_cohomology::Field_Zp>::get_persistent_pairs();
 
+    // Gather all top-dimensional cells and store their simplex handles
+    std::vector<std::size_t> max_splx;
+    for (auto splx : stptr_->top_dimensional_cells_range())
+      max_splx.push_back(splx);
+    // Sort these simplex handles and compute the ordering function
+    // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data
+    std::unordered_map<std::size_t, int> order;
+    //std::sort(max_splx.begin(), max_splx.end());
+    for (unsigned int i = 0; i < max_splx.size(); i++)  order.emplace(max_splx[i], i);
+
+    std::vector<std::vector<int>> persistence_pairs;
+    for (auto pair : pairs) {
+      int h = stptr_->dimension(get<0>(pair));
+      // Recursively get the top-dimensional cell / coface associated to the persistence generator
+      std::size_t face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair));
+      // Retrieve the index of the corresponding top-dimensional cell in the input data
+      int splx0 = order[face0];
+
+      int splx1 = -1;
+      if (get<1>(pair) != stptr_->null_simplex()){
+        // Recursively get the top-dimensional cell / coface associated to the persistence generator
+        std::size_t face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair));
+        // Retrieve the index of the corresponding top-dimensional cell in the input data
+        splx1 = order[face1];
+      }
+      persistence_pairs.push_back({ h, splx0, splx1 });
+    }
+    return persistence_pairs;
+  }
+
+  std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_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;
+  }
+
+  using Filtration_value = typename FilteredComplex::Filtration_value;
+  using Birth_death = std::pair<Filtration_value, Filtration_value>;
+  using Persistence_subdiagrams = std::vector<std::vector<std::pair<int, Birth_death>>>;
+
+  Persistence_subdiagrams compute_extended_persistence_subdiagrams(Filtration_value min_persistence){
+    Persistence_subdiagrams pers_subs(4);
+    auto const& persistent_pairs = Base::get_persistent_pairs();
+    for (auto pair : persistent_pairs) {
+      std::pair<Filtration_value, Extended_simplex_type> px = stptr_->decode_extended_filtration(stptr_->filtration(get<0>(pair)),
+                                                                                                        stptr_->efd);
+      std::pair<Filtration_value, Extended_simplex_type> py = stptr_->decode_extended_filtration(stptr_->filtration(get<1>(pair)),
+                                                                                                        stptr_->efd);
+      std::pair<int, Birth_death> pd_point = std::make_pair(stptr_->dimension(get<0>(pair)),
+                                                            std::make_pair(px.first, py.first));
+      if(std::abs(px.first - py.first) > min_persistence){
+        //Ordinary
+        if (px.second == Extended_simplex_type::UP && py.second == Extended_simplex_type::UP){
+          pers_subs[0].push_back(pd_point);
+        }
+        // Relative
+        else if (px.second == Extended_simplex_type::DOWN && py.second == Extended_simplex_type::DOWN){
+          pers_subs[1].push_back(pd_point);
+        }
+        else{
+          // Extended+
+          if (px.first < py.first){
+            pers_subs[2].push_back(pd_point);
+          }
+          //Extended-
+          else{
+            pers_subs[3].push_back(pd_point);
+          }
+        }
+      }
+    }
+    return pers_subs;
+  }
+
  private:
   // A copy
   FilteredComplex* stptr_;
diff --git a/src/python/include/Rips_complex_interface.h b/src/python/include/Rips_complex_interface.h
index a66b0e5b..d98b0226 100644
--- a/src/python/include/Rips_complex_interface.h
+++ b/src/python/include/Rips_complex_interface.h
@@ -53,7 +53,6 @@ class Rips_complex_interface {
       rips_complex_->create_complex(*simplex_tree, dim_max);
     else
       sparse_rips_complex_->create_complex(*simplex_tree, dim_max);
-    simplex_tree->initialize_filtration();
   }
 
  private:
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 06f31341..0317ea39 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -15,12 +15,13 @@
 #include <gudhi/distance_functions.h>
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Points_off_io.h>
-
-#include "Persistent_cohomology_interface.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 {
 
@@ -33,22 +34,60 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   using Simplex_handle = typename Base::Simplex_handle;
   using Insertion_result = typename std::pair<Simplex_handle, bool>;
   using Simplex = std::vector<Vertex_handle>;
-  using Filtered_simplices = std::vector<std::pair<Simplex, Filtration_value>>;
+  using Simplex_and_filtration = std::pair<Simplex, Filtration_value>;
+  using Filtered_simplices = std::vector<Simplex_and_filtration>;
+  using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator;
+  using Complex_simplex_iterator = typename Base::Complex_simplex_iterator;
+  using Extended_filtration_data = typename Base::Extended_filtration_data;
+  using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator;
+  using Siblings = typename Base::Siblings;
+  using Node = typename Base::Node;
+  typedef bool (*blocker_func_t)(Simplex simplex, void *user_data);
 
  public:
-  bool find_simplex(const Simplex& vh) {
-    return (Base::find(vh) != Base::null_simplex());
+
+  Extended_filtration_data efd;
+  
+  bool find_simplex(const Simplex& simplex) {
+    return (Base::find(simplex) != Base::null_simplex());
   }
 
-  void assign_simplex_filtration(const Simplex& vh, Filtration_value filtration) {
-    Base::assign_filtration(Base::find(vh), filtration);
+  void assign_simplex_filtration(const Simplex& simplex, Filtration_value filtration) {
+    Base::assign_filtration(Base::find(simplex), filtration);
+    Base::clear_filtration();
   }
 
   bool insert(const Simplex& simplex, Filtration_value filtration = 0) {
     Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration);
+    if (result.first != Base::null_simplex())
+      Base::clear_filtration();
     return (result.second);
   }
 
+  void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) {
+    // We could delegate to insert_graph, but wrapping the matrix in a graph interface is too much work,
+    // and this is a bit more efficient.
+    auto& rm = this->root()->members_;
+    for(int i=0; i<n; ++i) {
+      char* p = reinterpret_cast<char*>(filtrations) + i * stride0;
+      double fv = *reinterpret_cast<double*>(p + i * stride1);
+      if(fv > max_filtration) continue;
+      auto sh = rm.emplace_hint(rm.end(), i, Node(this->root(), fv));
+      Siblings* children = nullptr;
+      // Should we make a first pass to count the number of edges so we can reserve the right space?
+      for(int j=i+1; j<n; ++j) {
+        double fe = *reinterpret_cast<double*>(p + j * stride1);
+        if(fe > max_filtration) continue;
+        if(!children) {
+          children = new Siblings(this->root(), i);
+          sh->second.assign_children(children);
+        }
+        children->members().emplace_hint(children->members().end(), j, Node(children, fe));
+      }
+    }
+
+  }
+
   // Do not interface this function, only used in alpha complex interface for complex creation
   bool insert_simplex(const Simplex& simplex, Filtration_value filtration = 0) {
     Insertion_result result = Base::insert_simplex(simplex, filtration);
@@ -79,32 +118,15 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
 
   void remove_maximal_simplex(const Simplex& simplex) {
     Base::remove_maximal_simplex(Base::find(simplex));
-    Base::initialize_filtration();
+    Base::clear_filtration();
   }
 
-  Filtered_simplices get_filtration() {
-    Base::initialize_filtration();
-    Filtered_simplices filtrations;
-    for (auto f_simplex : Base::filtration_simplex_range()) {
-      Simplex simplex;
-      for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
-        simplex.insert(simplex.begin(), vertex);
-      }
-      filtrations.push_back(std::make_pair(simplex, Base::filtration(f_simplex)));
+  Simplex_and_filtration get_simplex_and_filtration(Simplex_handle f_simplex) {
+    Simplex simplex;
+    for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
+      simplex.insert(simplex.begin(), vertex);
     }
-    return filtrations;
-  }
-
-  Filtered_simplices get_skeleton(int dimension) {
-    Filtered_simplices skeletons;
-    for (auto f_simplex : Base::skeleton_simplex_range(dimension)) {
-      Simplex simplex;
-      for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
-        simplex.insert(simplex.begin(), vertex);
-      }
-      skeletons.push_back(std::make_pair(simplex, Base::filtration(f_simplex)));
-    }
-    return skeletons;
+    return std::make_pair(std::move(simplex), Base::filtration(f_simplex));
   }
 
   Filtered_simplices get_star(const Simplex& simplex) {
@@ -131,9 +153,85 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return cofaces;
   }
 
-  void create_persistence(Gudhi::Persistent_cohomology_interface<Base>* pcoh) {
-    Base::initialize_filtration();
-    pcoh = new Gudhi::Persistent_cohomology_interface<Base>(*this);
+  void compute_extended_filtration() {
+    this->efd = this->extend_filtration();
+    return;
+  }
+
+  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(std::move(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;
+  }
+
+  void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) {
+    Base::expansion_with_blockers(dimension, [&](Simplex_handle sh){
+      Simplex simplex(Base::simplex_vertex_range(sh).begin(), Base::simplex_vertex_range(sh).end());
+      return user_func(simplex, user_data);
+    });
+  }
+
+  // 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
+    return Base::complex_simplex_range().begin();
+  }
+
+  Complex_simplex_iterator get_simplices_iterator_end() {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::complex_simplex_range().end();
+  }
+
+  typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_begin() {
+    // Base::initialize_filtration(); already performed in filtration_simplex_range
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::filtration_simplex_range().begin();
+  }
+
+  typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_end() {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::filtration_simplex_range().end();
+  }
+
+  Skeleton_simplex_iterator get_skeleton_iterator_begin(int dimension) {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::skeleton_simplex_range(dimension).begin();
+  }
+
+  Skeleton_simplex_iterator get_skeleton_iterator_end(int dimension) {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::skeleton_simplex_range(dimension).end();
+  }
+
+  std::pair<Boundary_simplex_iterator, Boundary_simplex_iterator> get_boundary_iterators(const Simplex& simplex) {
+    auto bd_sh = Base::find(simplex);
+    if (bd_sh == Base::null_simplex())
+      throw std::runtime_error("simplex not found - cannot find boundaries");
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    auto boundary_srange = Base::boundary_simplex_range(bd_sh);
+    return std::make_pair(boundary_srange.begin(), boundary_srange.end());
   }
 };
 
diff --git a/src/python/include/Strong_witness_complex_interface.h b/src/python/include/Strong_witness_complex_interface.h
index cda5b514..e9ab0c7b 100644
--- a/src/python/include/Strong_witness_complex_interface.h
+++ b/src/python/include/Strong_witness_complex_interface.h
@@ -41,13 +41,11 @@ class Strong_witness_complex_interface {
   void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double  max_alpha_square,
                            std::size_t limit_dimension) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension);
-    simplex_tree->initialize_filtration();
   }
 
   void create_simplex_tree(Simplex_tree_interface<>* simplex_tree,
                            double  max_alpha_square) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square);
-    simplex_tree->initialize_filtration();
   }
 
  private:
diff --git a/src/python/include/Subsampling_interface.h b/src/python/include/Subsampling_interface.h
index cdda851f..6aee7231 100644
--- a/src/python/include/Subsampling_interface.h
+++ b/src/python/include/Subsampling_interface.h
@@ -11,6 +11,7 @@
 #ifndef INCLUDE_SUBSAMPLING_INTERFACE_H_
 #define INCLUDE_SUBSAMPLING_INTERFACE_H_
 
+#include <gudhi/distance_functions.h>
 #include <gudhi/choose_n_farthest_points.h>
 #include <gudhi/pick_n_random_points.h>
 #include <gudhi/sparsify_point_set.h>
@@ -27,14 +28,13 @@ namespace subsampling {
 
 using Subsampling_dynamic_kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
 using Subsampling_point_d = Subsampling_dynamic_kernel::Point_d;
-using Subsampling_ft = Subsampling_dynamic_kernel::FT;
 
 // ------ choose_n_farthest_points ------
 std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector<std::vector<double>>& points,
                                                                unsigned nb_points) {
   std::vector<std::vector<double>> landmarks;
-  Subsampling_dynamic_kernel k;
-  choose_n_farthest_points(k, points, nb_points, random_starting_point, std::back_inserter(landmarks));
+  choose_n_farthest_points(Euclidean_distance(), points, nb_points,
+                           random_starting_point, std::back_inserter(landmarks));
 
   return landmarks;
 }
@@ -42,8 +42,8 @@ std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector
 std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector<std::vector<double>>& points,
                                                                unsigned nb_points, unsigned starting_point) {
   std::vector<std::vector<double>> landmarks;
-  Subsampling_dynamic_kernel k;
-  choose_n_farthest_points(k, points, nb_points, starting_point, std::back_inserter(landmarks));
+  choose_n_farthest_points(Euclidean_distance(), points, nb_points,
+                           starting_point, std::back_inserter(landmarks));
 
   return landmarks;
 }
diff --git a/src/python/include/Tangential_complex_interface.h b/src/python/include/Tangential_complex_interface.h
index 698226cc..b1afce94 100644
--- a/src/python/include/Tangential_complex_interface.h
+++ b/src/python/include/Tangential_complex_interface.h
@@ -90,7 +90,6 @@ class Tangential_complex_interface {
 
   void create_simplex_tree(Simplex_tree<>* simplex_tree) {
     tangential_complex_->create_complex<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>(*simplex_tree);
-    simplex_tree->initialize_filtration();
   }
 
   void set_max_squared_edge_length(double max_squared_edge_length) {
diff --git a/src/python/include/Witness_complex_interface.h b/src/python/include/Witness_complex_interface.h
index 45e14253..76947e53 100644
--- a/src/python/include/Witness_complex_interface.h
+++ b/src/python/include/Witness_complex_interface.h
@@ -41,13 +41,11 @@ class Witness_complex_interface {
   void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double  max_alpha_square,
                            std::size_t limit_dimension) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension);
-    simplex_tree->initialize_filtration();
   }
 
   void create_simplex_tree(Simplex_tree_interface<>* simplex_tree,
                            double  max_alpha_square) {
     witness_complex_->create_complex(*simplex_tree, max_alpha_square);
-    simplex_tree->initialize_filtration();
   }
 
  private:
diff --git a/src/python/include/pybind11_diagram_utils.h b/src/python/include/pybind11_diagram_utils.h
new file mode 100644
index 00000000..2d5194f4
--- /dev/null
+++ b/src/python/include/pybind11_diagram_utils.h
@@ -0,0 +1,39 @@
+/*    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
+ */
+
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+#include <boost/range/counting_range.hpp>
+#include <boost/range/adaptor/transformed.hpp>
+
+namespace py = pybind11;
+typedef py::array_t<double> Dgm;
+
+// Get m[i,0] and m[i,1] as a pair
+static auto pairify(void* p, py::ssize_t h, py::ssize_t w) {
+  return [=](py::ssize_t i){
+    char* birth = (char*)p + i * h;
+    char* death = birth + w;
+    return std::make_pair(*(double*)birth, *(double*)death);
+  };
+}
+
+inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
+  py::buffer_info buf = dgm.request();
+  // shape (n,2) or (0) for empty
+  if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
+    throw std::runtime_error("Diagram must be an array of size n x 2");
+  // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
+  py::ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
+  auto cnt = boost::counting_range<py::ssize_t>(0, buf.shape[0]);
+  return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
+  // Be careful that the returned range cannot contain references to dead temporaries.
+}