Merge branch 'master' into ST-batch-vertices

25 files changed, 614 insertions, 290 deletions

diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
index 4a6af3a4..29fabc6c 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
@@ -241,10 +241,16 @@ class Bitmap_cubical_complex : public T {
    **/
   class Filtration_simplex_range;
 
-  class Filtration_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> {
+  class Filtration_simplex_iterator {
     // Iterator over all simplices of the complex in the order of the indexing scheme.
     // 'value_type' must be 'Simplex_handle'.
    public:
+    typedef std::input_iterator_tag iterator_category;
+    typedef Simplex_handle value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef value_type* pointer;
+    typedef value_type reference;
+
     Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) {}
 
     Filtration_simplex_iterator() : b(NULL), position(0) {}
@@ -386,10 +392,16 @@ class Bitmap_cubical_complex : public T {
    **/
   class Skeleton_simplex_range;
 
-  class Skeleton_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> {
+  class Skeleton_simplex_iterator {
     // Iterator over all simplices of the complex in the order of the indexing scheme.
     // 'value_type' must be 'Simplex_handle'.
    public:
+    typedef std::input_iterator_tag iterator_category;
+    typedef Simplex_handle value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef value_type* pointer;
+    typedef value_type reference;
+
     Skeleton_simplex_iterator(Bitmap_cubical_complex* b, std::size_t d) : b(b), dimension(d) {
       if (globalDbg) {
         std::clog << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , std::size_t d )\n";
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
index bafe7981..2bf62f9b 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -251,8 +251,14 @@ class Bitmap_cubical_complex_base {
    * @brief Iterator through all cells in the complex (in order they appear in the structure -- i.e.
    * in lexicographical order).
    **/
-  class All_cells_iterator : std::iterator<std::input_iterator_tag, T> {
+  class All_cells_iterator {
    public:
+    typedef std::input_iterator_tag iterator_category;
+    typedef std::size_t value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef value_type* pointer;
+    typedef value_type reference;
+
     All_cells_iterator() { this->counter = 0; }
 
     All_cells_iterator operator++() {
@@ -355,8 +361,14 @@ class Bitmap_cubical_complex_base {
    * @brief Iterator through top dimensional cells of the complex. The cells appear in order they are stored
    * in the structure (i.e. in lexicographical order)
    **/
-  class Top_dimensional_cells_iterator : std::iterator<std::input_iterator_tag, T> {
+  class Top_dimensional_cells_iterator {
    public:
+    typedef std::input_iterator_tag iterator_category;
+    typedef std::size_t value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef value_type* pointer;
+    typedef value_type reference;
+
     Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) {
       this->counter = std::vector<std::size_t>(b.dimension());
       // std::fill( this->counter.begin() , this->counter.end() , 0 );
diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
index 33f03b9c..c1e10f8e 100644
--- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
@@ -20,6 +20,7 @@
 #include <vector>
 #include <algorithm>
 #include <limits>  // for numeric_limits
+#include <cmath>
 
 namespace Gudhi {
 
@@ -31,7 +32,7 @@ namespace persistence_diagram {
  * \ingroup bottleneck_distance
  */
 class Persistence_graph {
- public:
+public:
   /** \internal \brief Constructor taking 2 PersistenceDiagrams (concept) as parameters. */
   template<typename Persistence_diagram1, typename Persistence_diagram2>
   Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e);
@@ -58,7 +59,7 @@ class Persistence_graph {
   /** \internal \brief Returns the corresponding internal point */
   Internal_point get_v_point(int v_point_index) const;
 
- private:
+private:
   std::vector<Internal_point> u;
   std::vector<Internal_point> v;
   double b_alive;
@@ -67,30 +68,54 @@ class Persistence_graph {
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1,
                                      const Persistence_diagram2 &diag2, double e)
-    : u(), v(), b_alive(0.) {
+  : u(), v(), b_alive(0.) {
   std::vector<double> u_alive;
   std::vector<double> v_alive;
+  std::vector<double> u_nalive;
+  std::vector<double> v_nalive;
+  int u_inf = 0;
+  int v_inf = 0;
+  double inf = std::numeric_limits<double>::infinity();
+  double neginf = -inf;
+
   for (auto it = std::begin(diag1); it != std::end(diag1); ++it) {
-    if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
-      u_alive.push_back(std::get<0>(*it));
-    else if (std::get<1>(*it) - std::get<0>(*it) > e)
-      u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
+    if (std::get<0>(*it) != inf && std::get<1>(*it) != neginf){
+      if (std::get<0>(*it) == neginf && std::get<1>(*it) == inf)
+        u_inf++;
+      else if (std::get<0>(*it) == neginf)
+        u_nalive.push_back(std::get<1>(*it));
+      else if (std::get<1>(*it) == inf)
+        u_alive.push_back(std::get<0>(*it));
+      else if (std::get<1>(*it) - std::get<0>(*it) > e)
+        u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
+    }
   }
   for (auto it = std::begin(diag2); it != std::end(diag2); ++it) {
-    if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
-      v_alive.push_back(std::get<0>(*it));
-    else if (std::get<1>(*it) - std::get<0>(*it) > e)
-      v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
+    if (std::get<0>(*it) != inf && std::get<1>(*it) != neginf){
+      if (std::get<0>(*it) == neginf && std::get<1>(*it) == inf)
+        v_inf++;
+      else if (std::get<0>(*it) == neginf)
+        v_nalive.push_back(std::get<1>(*it));
+      else if (std::get<1>(*it) == inf)
+        v_alive.push_back(std::get<0>(*it));
+      else if (std::get<1>(*it) - std::get<0>(*it) > e)
+        v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
+    }
   }
   if (u.size() < v.size())
     swap(u, v);
-  std::sort(u_alive.begin(), u_alive.end());
-  std::sort(v_alive.begin(), v_alive.end());
-  if (u_alive.size() != v_alive.size()) {
+
+  if (u_alive.size() != v_alive.size() || u_nalive.size() != v_nalive.size() || u_inf != v_inf) {
     b_alive = std::numeric_limits<double>::infinity();
   } else {
+    std::sort(u_alive.begin(), u_alive.end());
+    std::sort(v_alive.begin(), v_alive.end());
+    std::sort(u_nalive.begin(), u_nalive.end());
+    std::sort(v_nalive.begin(), v_nalive.end());
     for (auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
       b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
+    for (auto it_u = u_nalive.cbegin(), it_v = v_nalive.cbegin(); it_u != u_nalive.cend(); ++it_u, ++it_v)
+      b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
   }
 }
 
@@ -104,12 +129,12 @@ inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const {
 
 inline int Persistence_graph::corresponding_point_in_u(int v_point_index) const {
   return on_the_v_diagonal(v_point_index) ?
-      v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
+        v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
 }
 
 inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const {
   return on_the_u_diagonal(u_point_index) ?
-      u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
+        u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
 }
 
 inline double Persistence_graph::distance(int u_point_index, int v_point_index) const {
diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
index 44141baa..9872f20c 100644
--- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
+++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
@@ -159,3 +159,81 @@ BOOST_AUTO_TEST_CASE(global) {
   BOOST_CHECK(bottleneck_distance(empty, empty) == 0);
   BOOST_CHECK(bottleneck_distance(empty, one) == 1);
 }
+
+BOOST_AUTO_TEST_CASE(neg_global) {
+  std::uniform_real_distribution<double> unif1(0., upper_bound);
+  std::uniform_real_distribution<double> unif2(upper_bound / 10000., upper_bound / 100.);
+  std::default_random_engine re;
+  std::vector< std::pair<double, double> > v1, v2;
+  for (int i = 0; i < n1; i++) {
+    double a = std::log(unif1(re));
+    double b = std::log(unif1(re));
+    double x = std::log(unif2(re));
+    double y = std::log(unif2(re));
+    v1.emplace_back(std::min(a, b), std::max(a, b));
+    v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y));
+    if (i % 5 == 0)
+      v1.emplace_back(std::min(a, b), std::min(a, b) + x);
+    if (i % 3 == 0)
+      v2.emplace_back(std::max(a, b), std::max(a, b) + y);
+  }
+  BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound / 100.);
+  BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound / 10000.) <= upper_bound / 100. + upper_bound / 10000.);
+  BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound / 10000.)) <= upper_bound / 10000.);
+
+  std::vector< std::pair<double, double> > empty;
+  std::vector< std::pair<double, double> > one = {{8, 10}};
+  BOOST_CHECK(bottleneck_distance(empty, empty) == 0);
+  BOOST_CHECK(bottleneck_distance(empty, one) == 1);
+}
+
+BOOST_AUTO_TEST_CASE(bottleneck_simple_test) {
+  std::vector< std::pair<double, double> > v1, v2;
+  double inf = std::numeric_limits<double>::infinity();
+  double neginf = -inf;
+  double b;
+
+  v1.emplace_back(9.6, 14.);
+  v2.emplace_back(9.5, 14.1);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK(b > 0.09 && b < 0.11);
+
+  v1.emplace_back(-34.974, -34.2);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK(b > 0.386 && b < 0.388);
+
+  v1.emplace_back(neginf, 3.7);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK_EQUAL(b, inf);
+
+  v2.emplace_back(neginf, 4.45);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK(b > 0.74 && b < 0.76);
+
+  v1.emplace_back(-60.6, 52.1);
+  v2.emplace_back(-61.5, 53.);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK(b > 0.89 && b < 0.91);
+
+  v1.emplace_back(3., inf);
+  v2.emplace_back(3.2, inf);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK(b > 0.89 && b < 0.91);
+
+  v1.emplace_back(neginf, inf);
+  v2.emplace_back(neginf, inf);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK(b > 0.89 && b < 0.91);
+
+  v2.emplace_back(6, inf);
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.);
+  BOOST_CHECK_EQUAL(b, inf);
+}
diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h
index 595fb64b..73093c07 100644
--- a/src/Cech_complex/doc/Intro_cech_complex.h
+++ b/src/Cech_complex/doc/Intro_cech_complex.h
@@ -17,7 +17,7 @@ namespace cech_complex {
 
 /**  \defgroup cech_complex Čech complex
  * 
- * \author    Vincent Rouvreau
+ * \author    Vincent Rouvreau, Hind montassif
  * 
  * @{
  * 
diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h
index 625f7c9c..dbdf5e93 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex.h
@@ -1,11 +1,12 @@
 /*    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
+ *    Author(s):       Vincent Rouvreau, Hind Montassif
  *
  *    Copyright (C) 2018 Inria
  *
  *    Modification(s):
  *      - YYYY/MM Author: Description of the modification
+ *      - 2022/02 Hind Montassif : Replace MiniBall with Sphere_circumradius
  */
 
 #ifndef CECH_COMPLEX_H_
diff --git a/src/GudhUI/view/Viewer.cpp b/src/GudhUI/view/Viewer.cpp
index 6b17c833..2c00f86f 100644
--- a/src/GudhUI/view/Viewer.cpp
+++ b/src/GudhUI/view/Viewer.cpp
@@ -31,7 +31,11 @@ void Viewer::set_bounding_box(const Point_3 & lower_left, const Point_3 & upper_
 }
 
 void Viewer::update_GL() {
+#if QGLVIEWER_VERSION >= 0x020700
+  this->update();
+#else
   this->updateGL();
+#endif
 }
 
 void Viewer::init_scene() {
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 32168d9e..4177a0b8 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -25,6 +25,7 @@
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/range/adaptor/reversed.hpp>
 #include <boost/range/adaptor/transformed.hpp>
+#include <boost/range/size.hpp>
 #include <boost/container/static_vector.hpp>
 
 #ifdef GUDHI_USE_TBB
@@ -703,10 +704,10 @@ class Simplex_tree {
     return true;
   }
 
- private:
-  /** \brief Inserts a simplex represented by a vector of vertex.
-   * @param[in]  simplex    vector of Vertex_handles, representing the vertices of the new simplex. The vector must be
-   * sorted by increasing vertex handle order.
+ protected:
+  /** \brief Inserts a simplex represented by a range of vertex.
+   * @param[in]  simplex    range of Vertex_handles, representing the vertices of the new simplex. The range must be
+   * sorted by increasing vertex handle order, and not empty.
    * @param[in]  filtration the filtration value assigned to the new simplex.
    * @return If the new simplex is inserted successfully (i.e. it was not in the
    * simplicial complex yet) the bool is set to true and the Simplex_handle is the handle assigned
@@ -718,12 +719,13 @@ class Simplex_tree {
    * null_simplex.
    * 
   */
-  std::pair<Simplex_handle, bool> insert_vertex_vector(const std::vector<Vertex_handle>& simplex,
+  template <class RandomVertexHandleRange = std::initializer_list<Vertex_handle>>
+  std::pair<Simplex_handle, bool> insert_simplex_raw(const RandomVertexHandleRange& simplex,
                                                      Filtration_value filtration) {
     Siblings * curr_sib = &root_;
     std::pair<Simplex_handle, bool> res_insert;
     auto vi = simplex.begin();
-    for (; vi != simplex.end() - 1; ++vi) {
+    for (; vi != std::prev(simplex.end()); ++vi) {
       GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
       res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
       if (!(has_children(res_insert.first))) {
@@ -744,9 +746,10 @@ class Simplex_tree {
       return std::pair<Simplex_handle, bool>(null_simplex(), false);
     }
     // otherwise the insertion has succeeded - size is a size_type
-    if (static_cast<int>(simplex.size()) - 1 > dimension_) {
+    int dim = static_cast<int>(boost::size(simplex)) - 1;
+    if (dim > dimension_) {
       // Update dimension if needed
-      dimension_ = static_cast<int>(simplex.size()) - 1;
+      dimension_ = dim;
     }
     return res_insert;
   }
@@ -787,7 +790,7 @@ class Simplex_tree {
     // Copy before sorting
     std::vector<Vertex_handle> copy(first, last);
     std::sort(std::begin(copy), std::end(copy));
-    return insert_vertex_vector(copy, filtration);
+    return insert_simplex_raw(copy, filtration);
   }
 
     /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
@@ -1610,7 +1613,7 @@ class Simplex_tree {
     Simplex_tree st_copy = *this;
 
     // Add point for coning the simplicial complex
-    this->insert_simplex({maxvert}, -3);
+    this->insert_simplex_raw({maxvert}, -3);
 
     // For each simplex
     std::vector<Vertex_handle> vr;
diff --git a/src/Spatial_searching/example/example_spatial_searching.cpp b/src/Spatial_searching/example/example_spatial_searching.cpp
index 8f9151fc..09c2dabf 100644
--- a/src/Spatial_searching/example/example_spatial_searching.cpp
+++ b/src/Spatial_searching/example/example_spatial_searching.cpp
@@ -25,7 +25,7 @@ int main(void) {
   // 10-nearest neighbor query
   std::clog << "10 nearest neighbors from points[20]:\n";
   auto knn_range = points_ds.k_nearest_neighbors(points[20], 10, true);
-  for (auto const& nghb : knn_range)
+  for (auto const nghb : knn_range)
     std::clog << nghb.first << " (sq. dist. = " << nghb.second << ")\n";
 
   // Incremental nearest neighbor query
@@ -38,7 +38,7 @@ int main(void) {
   // 10-furthest neighbor query
   std::clog << "10 furthest neighbors from points[20]:\n";
   auto kfn_range = points_ds.k_furthest_neighbors(points[20], 10, true);
-  for (auto const& nghb : kfn_range)
+  for (auto const nghb : kfn_range)
     std::clog << nghb.first << " (sq. dist. = " << nghb.second << ")\n";
 
   // Incremental furthest neighbor query
diff --git a/src/Spatial_searching/test/test_Kd_tree_search.cpp b/src/Spatial_searching/test/test_Kd_tree_search.cpp
index d6c6fba3..e9acfaa7 100644
--- a/src/Spatial_searching/test/test_Kd_tree_search.cpp
+++ b/src/Spatial_searching/test/test_Kd_tree_search.cpp
@@ -45,7 +45,7 @@ BOOST_AUTO_TEST_CASE(test_Kd_tree_search) {
 
   std::vector<std::size_t> knn_result;
   FT last_dist = -1.;
-  for (auto const& nghb : kns_range) {
+  for (auto const nghb : kns_range) {
     BOOST_CHECK(nghb.second > last_dist);
     knn_result.push_back(nghb.second);
     last_dist = nghb.second;
@@ -76,7 +76,7 @@ BOOST_AUTO_TEST_CASE(test_Kd_tree_search) {
 
   std::vector<std::size_t> kfn_result;
   last_dist = kfn_range.begin()->second;
-  for (auto const& nghb : kfn_range) {
+  for (auto const nghb : kfn_range) {
     BOOST_CHECK(nghb.second <= last_dist);
     kfn_result.push_back(nghb.second);
     last_dist = nghb.second;
diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md
index ce903405..9b7c2853 100644
--- a/src/common/doc/main_page.md
+++ b/src/common/doc/main_page.md
@@ -178,7 +178,7 @@
     The set of all simplices is filtered by the radius of their minimal enclosing ball.
     </td>
     <td width="15%">
-      <b>Author:</b> Vincent Rouvreau<br>
+      <b>Author:</b> Vincent Rouvreau, Hind Montassif<br>
       <b>Introduced in:</b> GUDHI 2.2.0<br>
       <b>Copyright:</b> MIT [(LGPL v3)](../../licensing/)<br>
       <b>Requires:</b> \ref cgal
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 8f8df138..32ec13bd 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -53,7 +53,7 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'datasets', ")
 
     # Cython modules
-    set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_reader', ")
+    set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_utils', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'cubical_complex', ")
@@ -152,7 +152,7 @@ if(PYTHONINTERP_FOUND)
       set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ")
     endif (EIGEN3_FOUND)
 
-    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_utils', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'rips_complex', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'cubical_complex', ")
@@ -546,6 +546,7 @@ if(PYTHONINTERP_FOUND)
 
     # Reader utils
     add_gudhi_py_test(test_reader_utils)
+    add_gudhi_py_test(test_off)
 
     # Wasserstein
     if(OT_FOUND)
diff --git a/src/python/doc/clustering.rst b/src/python/doc/clustering.rst
index c5a57d3c..62422682 100644
--- a/src/python/doc/clustering.rst
+++ b/src/python/doc/clustering.rst
@@ -17,9 +17,8 @@ As a by-product, we produce the persistence diagram of the merge tree of the ini
    :include-source:
 
    import gudhi
-   f = open(gudhi.__root_source_dir__ + '/data/points/spiral_2d.csv', 'r')
-   import numpy as np
-   data = np.loadtxt(f)
+   from gudhi.datasets.remote import fetch_spiral_2d
+   data = fetch_spiral_2d()
    import matplotlib.pyplot as plt
    plt.scatter(data[:,0],data[:,1],marker='.',s=1)
    plt.show()
diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst
index 276ac4e2..5491542f 100644
--- a/src/python/doc/installation.rst
+++ b/src/python/doc/installation.rst
@@ -150,7 +150,7 @@ You shall have something like:
     Cython version 0.29.25 
     Numpy version 1.21.4 
     Boost version 1.77.0
-    + Installed modules are: off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex;
+    + Installed modules are: off_utils;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex;
         persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex;
     + Missing modules are: bottleneck;nerve_gic;subsampling;tangential_complex;alpha_complex;euclidean_witness_complex;
         euclidean_strong_witness_complex;
@@ -188,7 +188,7 @@ A complete configuration would be :
     GMPXX_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmpxx.so
     MPFR_LIBRARIES = /usr/lib/x86_64-linux-gnu/libmpfr.so
     TBB version 9107 found and used
-    + Installed modules are: bottleneck;off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex;
+    + Installed modules are: bottleneck;off_utils;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex;
         persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex;nerve_gic;subsampling;
         tangential_complex;alpha_complex;euclidean_witness_complex;euclidean_strong_witness_complex;
     + Missing modules are: 
diff --git a/src/python/doc/point_cloud.rst b/src/python/doc/point_cloud.rst
index ffd8f85b..473b303f 100644
--- a/src/python/doc/point_cloud.rst
+++ b/src/python/doc/point_cloud.rst
@@ -13,6 +13,11 @@ File Readers
 
 .. autofunction:: gudhi.read_lower_triangular_matrix_from_csv_file
 
+File Writers
+------------
+
+.. autofunction:: gudhi.write_points_to_off_file
+
 Subsampling
 -----------
 
diff --git a/src/python/doc/representations_sum.inc b/src/python/doc/representations_sum.inc
index 4298aea9..9515f044 100644
--- a/src/python/doc/representations_sum.inc
+++ b/src/python/doc/representations_sum.inc
@@ -1,14 +1,14 @@
 .. table::
    :widths: 30 40 30
 
-   +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------+
-   | .. figure::                                                      | Vectorizations, distances and kernels that work on persistence | :Author: Mathieu Carrière, Martin Royer                     |
-   |      img/sklearn-tda.png                                         | diagrams, compatible with scikit-learn.                        |                                                             |
-   |                                                                  |                                                                | :Since: GUDHI 3.1.0                                         |
-   |                                                                  |                                                                |                                                             |
-   |                                                                  |                                                                | :License: MIT                                               |
-   |                                                                  |                                                                |                                                             |
-   |                                                                  |                                                                | :Requires: `Scikit-learn <installation.html#scikit-learn>`_ |
-   +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------+
-   | * :doc:`representations`                                                                                                                                                                        |
-   +------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------------+
+   +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------------------+
+   | .. figure::                                                      | Vectorizations, distances and kernels that work on persistence | :Author: Mathieu Carrière, Martin Royer, Gard Spreemann, Wojciech Reise |
+   |      img/sklearn-tda.png                                         | diagrams, compatible with scikit-learn.                        |                                                                         |
+   |                                                                  |                                                                | :Since: GUDHI 3.1.0                                                     |
+   |                                                                  |                                                                |                                                                         |
+   |                                                                  |                                                                | :License: MIT                                                           |
+   |                                                                  |                                                                |                                                                         |
+   |                                                                  |                                                                | :Requires: `Scikit-learn <installation.html#scikit-learn>`_             |
+   +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------------------+
+   | * :doc:`representations`                                                                                                                                                                                    |
+   +------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------------------------+
diff --git a/src/python/gudhi/off_reader.pyx b/src/python/gudhi/off_utils.pyx
index a3200704..9276c7b0 100644
--- a/src/python/gudhi/off_reader.pyx
+++ b/src/python/gudhi/off_utils.pyx
@@ -13,8 +13,10 @@ from __future__ import print_function
 from cython cimport numeric
 from libcpp.vector cimport vector
 from libcpp.string cimport string
+cimport cython
 import errno
 import os
+import numpy as np
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
@@ -24,7 +26,7 @@ cdef extern from "Off_reader_interface.h" namespace "Gudhi":
     vector[vector[double]] read_points_from_OFF_file(string off_file)
 
 def read_points_from_off_file(off_file=''):
-    """Read points from OFF file.
+    """Read points from an `OFF file <fileformats.html#off-file-format>`_.
 
     :param off_file: An OFF file style name.
     :type off_file: string
@@ -39,3 +41,22 @@ def read_points_from_off_file(off_file=''):
             raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT),
                                     off_file)
 
+@cython.embedsignature(True)
+def write_points_to_off_file(fname, points):
+    """Write points to an `OFF file <fileformats.html#off-file-format>`_.
+
+    A simple wrapper for `numpy.savetxt`.
+
+    :param fname: Name of the OFF file.
+    :type fname: str or file handle
+    :param points: Point coordinates.
+    :type points: numpy array of shape (n, dim)
+    """
+    points = np.array(points, copy=False)
+    assert len(points.shape) == 2
+    dim = points.shape[1]
+    if dim == 3:
+        head = 'OFF\n{} 0 0'.format(points.shape[0])
+    else:
+        head = 'nOFF\n{} {} 0 0'.format(dim, points.shape[0])
+    np.savetxt(fname, points, header=head, comments='')
diff --git a/src/python/gudhi/rips_complex.pyx b/src/python/gudhi/rips_complex.pyx
index c3470292..d748f91e 100644
--- a/src/python/gudhi/rips_complex.pyx
+++ b/src/python/gudhi/rips_complex.pyx
@@ -41,31 +41,30 @@ cdef class RipsComplex:
     cdef Rips_complex_interface thisref
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=None, distance_matrix=None,
+    def __init__(self, *, points=None, distance_matrix=None,
                  max_edge_length=float('inf'), sparse=None):
         """RipsComplex constructor.
 
-        :param max_edge_length: Rips value.
-        :type max_edge_length: float
-
         :param points: A list of points in d-Dimension.
-        :type points: list of list of float
+        :type points: List[List[float]]
 
         Or
 
         :param distance_matrix: A distance matrix (full square or lower
             triangular).
-        :type points: list of list of float
+        :type distance_matrix: List[List[float]]
 
         And in both cases
 
+        :param max_edge_length: Rips value.
+        :type max_edge_length: float
         :param sparse: If this is not None, it switches to building a sparse
             Rips and represents the approximation parameter epsilon.
         :type sparse: float
         """
 
     # The real cython constructor
-    def __cinit__(self, points=None, distance_matrix=None,
+    def __cinit__(self, *, points=None, distance_matrix=None,
                   max_edge_length=float('inf'), sparse=None):
         if sparse is not None:
           if distance_matrix is not None:
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 5642f82d..f86f1232 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -56,6 +56,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         int upper_bound_dimension() nogil
         bool find_simplex(vector[int] simplex) nogil
         bool insert(vector[int] simplex, double filtration) nogil
+        void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) nogil
         vector[pair[vector[int], double]] get_star(vector[int] simplex) nogil
         vector[pair[vector[int], double]] get_cofaces(vector[int] simplex, int dimension) nogil
         void expansion(int max_dim) nogil except +
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 05bfe22e..18215d2f 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -8,14 +8,23 @@
 #   - YYYY/MM Author: Description of the modification
 
 from cython.operator import dereference, preincrement
-from libc.stdint cimport intptr_t
+from libc.stdint cimport intptr_t, int32_t, int64_t
 import numpy as np
 cimport gudhi.simplex_tree
+cimport cython
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
 __license__ = "MIT"
 
+ctypedef fused some_int:
+    int32_t
+    int64_t
+
+ctypedef fused some_float:
+    float
+    double
+
 cdef bool callback(vector[int] simplex, void *blocker_func):
     return (<object>blocker_func)(simplex)
 
@@ -228,6 +237,87 @@ cdef class SimplexTree:
         """
         return self.get_ptr().insert(simplex, <double>filtration)
 
+    @staticmethod
+    @cython.boundscheck(False)
+    def create_from_array(filtrations, double max_filtration=np.inf):
+        """Creates a new, empty complex and inserts vertices and edges. The vertices are numbered from 0 to n-1, and
+        the filtration values are encoded in the array, with the diagonal representing the vertices. It is the
+        caller's responsibility to ensure that this defines a filtration, which can be achieved with either::
+
+            filtrations[np.diag_indices_from(filtrations)] = filtrations.min(axis=1)
+
+        or::
+
+            diag = filtrations.diagonal()
+            filtrations = np.fmax(np.fmax(filtrations, diag[:, None]), diag[None, :])
+
+        :param filtrations: the filtration values of the vertices and edges to insert. The matrix is assumed to be symmetric.
+        :type filtrations: numpy.ndarray of shape (n,n)
+        :param max_filtration: only insert vertices and edges with filtration values no larger than max_filtration
+        :type max_filtration: float
+        :returns: the new complex
+        :rtype: SimplexTree
+        """
+        # TODO: document which half of the matrix is actually read?
+        filtrations = np.asanyarray(filtrations, dtype=float)
+        cdef double[:,:] F = filtrations
+        ret = SimplexTree()
+        cdef int n = F.shape[0]
+        assert n == F.shape[1], 'create_from_array() expects a square array'
+        with nogil:
+            ret.get_ptr().insert_matrix(&F[0,0], n, F.strides[0], F.strides[1], max_filtration)
+        return ret
+
+    def insert_edges_from_coo_matrix(self, edges):
+        """Inserts edges given by a sparse matrix in `COOrdinate format
+        <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html>`_.
+        If an edge is repeated, the smallest filtration value is used. Missing entries are not inserted.
+        Diagonal entries are currently interpreted as vertices, although we do not guarantee this behavior
+        in the future, and this is only useful if you want to insert vertices with a smaller filtration value
+        than the smallest edge containing it, since vertices are implicitly inserted together with the edges.
+
+        :param edges: the edges to insert and their filtration values.
+        :type edges: scipy.sparse.coo_matrix of shape (n,n)
+
+        .. seealso:: :func:`insert_batch`
+        """
+        # TODO: optimize this?
+        for edge in zip(edges.row, edges.col, edges.data):
+            self.get_ptr().insert((edge[0], edge[1]), edge[2])
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def insert_batch(self, some_int[:,:] vertex_array, some_float[:] filtrations):
+        """Inserts k-simplices given by a sparse array in a format similar
+        to `torch.sparse <https://pytorch.org/docs/stable/sparse.html>`_.
+        The n-th simplex has vertices `vertex_array[0,n]`, ...,
+        `vertex_array[k,n]` and filtration value `filtrations[n]`.
+        If a simplex is repeated, the smallest filtration value is used.
+        Simplices with a repeated vertex are currently interpreted as lower
+        dimensional simplices, but we do not guarantee this behavior in the
+        future. Any time a simplex is inserted, its faces are inserted as well
+        if needed to preserve a simplicial complex.
+
+        :param vertex_array: the k-simplices to insert.
+        :type vertex_array: numpy.array of shape (k+1,n)
+        :param filtrations: the filtration values.
+        :type filtrations: numpy.array of shape (n,)
+        """
+        # This may be slow if we end up inserting vertices in a bad order (flat_map).
+        # We could first insert the vertices from np.unique(vertex_array), or leave it to the caller.
+        cdef Py_ssize_t k = vertex_array.shape[0]
+        cdef Py_ssize_t n = vertex_array.shape[1]
+        assert filtrations.shape[0] == n, 'inconsistent sizes for vertex_array and filtrations'
+        cdef Py_ssize_t i
+        cdef Py_ssize_t j
+        cdef vector[int] v
+        with nogil:
+            for i in range(n):
+                for j in range(k):
+                    v.push_back(vertex_array[j, i])
+                self.get_ptr().insert(v, filtrations[i])
+                v.clear()
+
     def get_simplices(self):
         """This function returns a generator with simplices and their given
         filtration values.
@@ -376,7 +466,7 @@ cdef class SimplexTree:
         """
         return self.get_ptr().prune_above_filtration(filtration)
 
-    def expansion(self, max_dim):
+    def expansion(self, max_dimension):
         """Expands the simplex tree containing only its one skeleton
         until dimension max_dim.
 
@@ -390,10 +480,10 @@ cdef class SimplexTree:
         The simplex tree must contain no simplex of dimension bigger than
         1 when calling the method.
 
-        :param max_dim: The maximal dimension.
-        :type max_dim: int
+        :param max_dimension: The maximal dimension.
+        :type max_dimension: int
         """
-        cdef int maxdim = max_dim
+        cdef int maxdim = max_dimension
         with nogil:
             self.get_ptr().expansion(maxdim)
 
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 3848c5ad..0317ea39 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -40,6 +40,8 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   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:
@@ -62,6 +64,30 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     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);
diff --git a/src/python/test/test_off.py b/src/python/test/test_off.py
new file mode 100644
index 00000000..aea1941b
--- /dev/null
+++ b/src/python/test/test_off.py
@@ -0,0 +1,21 @@
+""" 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) 2022 Inria
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+import gudhi as gd
+import numpy as np
+import pytest
+
+
+def test_off_rw():
+    for dim in range(2, 6):
+        X = np.random.rand(123, dim)
+        gd.write_points_to_off_file("rand.off", X)
+        Y = gd.read_points_from_off_file("rand.off")
+        assert Y == pytest.approx(X)
diff --git a/src/python/test/test_simplex_generators.py b/src/python/test/test_simplex_generators.py
index 8a9b4844..c567d4c1 100755
--- a/src/python/test/test_simplex_generators.py
+++ b/src/python/test/test_simplex_generators.py
@@ -14,7 +14,7 @@ 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)
+    r = gudhi.RipsComplex(points=pts, max_edge_length=4)
     st = r.create_simplex_tree(max_dimension=50)
     st.persistence()
     g = st.flag_persistence_generators()
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 54bafed5..2ccbfbf5 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -249,6 +249,7 @@ def test_make_filtration_non_decreasing():
     assert st.filtration([3, 4]) == 2.0
     assert st.filtration([4, 5]) == 2.0
 
+
 def test_extend_filtration():
 
     # Inserted simplex:
@@ -257,86 +258,87 @@ def test_extend_filtration():
     #     / \ /
     #    o   o
     #   /2\ /3
-    #  o   o        
-    #  1   0        
-
-    st = SimplexTree()                                                                                                                     
-    st.insert([0,2])
-    st.insert([1,2])
-    st.insert([0,3])
-    st.insert([2,5])
-    st.insert([3,4])
-    st.insert([3,5])                                                                                                         
-    st.assign_filtration([0], 1.)                                                                                                                
-    st.assign_filtration([1], 2.)                                                                                                                
-    st.assign_filtration([2], 3.)                                                                                                                
-    st.assign_filtration([3], 4.)                                                                                                                
-    st.assign_filtration([4], 5.)                                                                                                                
-    st.assign_filtration([5], 6.)                                                                                                                
-
-    assert list(st.get_filtration()) == [                                                                                                                         
-        ([0, 2], 0.0), 
-        ([1, 2], 0.0), 
-        ([0, 3], 0.0), 
-        ([3, 4], 0.0), 
-        ([2, 5], 0.0), 
-        ([3, 5], 0.0), 
-        ([0], 1.0), 
-        ([1], 2.0), 
-        ([2], 3.0), 
-        ([3], 4.0), 
-        ([4], 5.0), 
-        ([5], 6.0)
+    #  o   o
+    #  1   0
+
+    st = SimplexTree()
+    st.insert([0, 2])
+    st.insert([1, 2])
+    st.insert([0, 3])
+    st.insert([2, 5])
+    st.insert([3, 4])
+    st.insert([3, 5])
+    st.assign_filtration([0], 1.0)
+    st.assign_filtration([1], 2.0)
+    st.assign_filtration([2], 3.0)
+    st.assign_filtration([3], 4.0)
+    st.assign_filtration([4], 5.0)
+    st.assign_filtration([5], 6.0)
+
+    assert list(st.get_filtration()) == [
+        ([0, 2], 0.0),
+        ([1, 2], 0.0),
+        ([0, 3], 0.0),
+        ([3, 4], 0.0),
+        ([2, 5], 0.0),
+        ([3, 5], 0.0),
+        ([0], 1.0),
+        ([1], 2.0),
+        ([2], 3.0),
+        ([3], 4.0),
+        ([4], 5.0),
+        ([5], 6.0),
     ]
-        
+
     st.extend_filtration()
-    
-    assert list(st.get_filtration()) == [                                                                                                                         
-        ([6], -3.0), 
-        ([0], -2.0), 
-        ([1], -1.8), 
-        ([2], -1.6), 
-        ([0, 2], -1.6), 
-        ([1, 2], -1.6), 
-        ([3], -1.4), 
-        ([0, 3], -1.4), 
-        ([4], -1.2), 
-        ([3, 4], -1.2), 
-        ([5], -1.0), 
-        ([2, 5], -1.0), 
-        ([3, 5], -1.0), 
-        ([5, 6], 1.0), 
-        ([4, 6], 1.2), 
-        ([3, 6], 1.4), 
+
+    assert list(st.get_filtration()) == [
+        ([6], -3.0),
+        ([0], -2.0),
+        ([1], -1.8),
+        ([2], -1.6),
+        ([0, 2], -1.6),
+        ([1, 2], -1.6),
+        ([3], -1.4),
+        ([0, 3], -1.4),
+        ([4], -1.2),
+        ([3, 4], -1.2),
+        ([5], -1.0),
+        ([2, 5], -1.0),
+        ([3, 5], -1.0),
+        ([5, 6], 1.0),
+        ([4, 6], 1.2),
+        ([3, 6], 1.4),
         ([3, 4, 6], 1.4),
-        ([3, 5, 6], 1.4), 
-        ([2, 6], 1.6), 
-        ([2, 5, 6], 1.6), 
-        ([1, 6], 1.8), 
-        ([1, 2, 6], 1.8), 
-        ([0, 6], 2.0), 
-        ([0, 2, 6], 2.0), 
-        ([0, 3, 6], 2.0)
+        ([3, 5, 6], 1.4),
+        ([2, 6], 1.6),
+        ([2, 5, 6], 1.6),
+        ([1, 6], 1.8),
+        ([1, 2, 6], 1.8),
+        ([0, 6], 2.0),
+        ([0, 2, 6], 2.0),
+        ([0, 3, 6], 2.0),
     ]
 
-    dgms = st.extended_persistence(min_persistence=-1.)
+    dgms = st.extended_persistence(min_persistence=-1.0)
     assert len(dgms) == 4
     # Sort by (death-birth) descending - we are only interested in those with the longest life span
     for idx in range(4):
-        dgms[idx] = sorted(dgms[idx], key=lambda x:(-abs(x[1][0]-x[1][1])))
+        dgms[idx] = sorted(dgms[idx], key=lambda x: (-abs(x[1][0] - x[1][1])))
+
+    assert dgms[0][0][1][0] == pytest.approx(2.0)
+    assert dgms[0][0][1][1] == pytest.approx(3.0)
+    assert dgms[1][0][1][0] == pytest.approx(5.0)
+    assert dgms[1][0][1][1] == pytest.approx(4.0)
+    assert dgms[2][0][1][0] == pytest.approx(1.0)
+    assert dgms[2][0][1][1] == pytest.approx(6.0)
+    assert dgms[3][0][1][0] == pytest.approx(6.0)
+    assert dgms[3][0][1][1] == pytest.approx(1.0)
 
-    assert dgms[0][0][1][0] == pytest.approx(2.)
-    assert dgms[0][0][1][1] == pytest.approx(3.)
-    assert dgms[1][0][1][0] == pytest.approx(5.)
-    assert dgms[1][0][1][1] == pytest.approx(4.)
-    assert dgms[2][0][1][0] == pytest.approx(1.)
-    assert dgms[2][0][1][1] == pytest.approx(6.)
-    assert dgms[3][0][1][0] == pytest.approx(6.)
-    assert dgms[3][0][1][1] == pytest.approx(1.) 
 
 def test_simplices_iterator():
     st = SimplexTree()
-    
+
     assert st.insert([0, 1, 2], filtration=4.0) == True
     assert st.insert([2, 3, 4], filtration=2.0) == True
 
@@ -346,9 +348,10 @@ def test_simplices_iterator():
         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
@@ -360,31 +363,33 @@ def test_collapse_edges():
 
     st.collapse_edges()
     assert st.num_simplices() == 9
-    assert st.find([0, 2]) == False # [1, 3] would be fine as well
+    assert st.find([0, 2]) == False  # [1, 3] would be fine as well
     for simplex in st.get_skeleton(0):
-        assert simplex[1] == 1.
+        assert simplex[1] == 1.0
+
 
 def test_reset_filtration():
     st = SimplexTree()
-    
-    assert st.insert([0, 1, 2], 3.) == True
-    assert st.insert([0, 3], 2.) == True
-    assert st.insert([3, 4, 5], 3.) == True
-    assert st.insert([0, 1, 6, 7], 4.) == True
+
+    assert st.insert([0, 1, 2], 3.0) == True
+    assert st.insert([0, 3], 2.0) == True
+    assert st.insert([3, 4, 5], 3.0) == True
+    assert st.insert([0, 1, 6, 7], 4.0) == True
 
     # Guaranteed by construction
     for simplex in st.get_simplices():
-        assert st.filtration(simplex[0]) >= 2.
-    
+        assert st.filtration(simplex[0]) >= 2.0
+
     # dimension until 5 even if simplex tree is of dimension 3 to test the limits
     for dimension in range(5, -1, -1):
-        st.reset_filtration(0., dimension)
+        st.reset_filtration(0.0, dimension)
         for simplex in st.get_skeleton(3):
             print(simplex)
             if len(simplex[0]) < (dimension) + 1:
-                assert st.filtration(simplex[0]) >= 2.
+                assert st.filtration(simplex[0]) >= 2.0
             else:
-                assert st.filtration(simplex[0]) == 0.
+                assert st.filtration(simplex[0]) == 0.0
+
 
 def test_boundaries_iterator():
     st = SimplexTree()
@@ -400,16 +405,17 @@ def test_boundaries_iterator():
         list(st.get_boundaries([]))
 
     with pytest.raises(RuntimeError):
-        list(st.get_boundaries([0, 4])) # (0, 4) does not exist
+        list(st.get_boundaries([0, 4]))  # (0, 4) does not exist
 
     with pytest.raises(RuntimeError):
-        list(st.get_boundaries([6])) # (6) does not exist
+        list(st.get_boundaries([6]))  # (6) does not exist
+
 
 def test_persistence_intervals_in_dimension():
     # Here is our triangulation of a 2-torus - taken from https://dioscuri-tda.org/Paris_TDA_Tutorial_2021.html
     #   0-----3-----4-----0
     #   | \   | \   | \   | \   |
-    #   |   \ |   \ |    \|   \ | 
+    #   |   \ |   \ |    \|   \ |
     #   1-----8-----7-----1
     #   | \   | \   | \   | \   |
     #   |   \ |   \ |   \ |   \ |
@@ -418,50 +424,52 @@ def test_persistence_intervals_in_dimension():
     #   |   \ |   \ |   \ |   \ |
     #   0-----3-----4-----0
     st = SimplexTree()
-    st.insert([0,1,8])
-    st.insert([0,3,8])
-    st.insert([3,7,8])
-    st.insert([3,4,7])
-    st.insert([1,4,7])
-    st.insert([0,1,4])
-    st.insert([1,2,5])
-    st.insert([1,5,8])
-    st.insert([5,6,8])
-    st.insert([6,7,8])
-    st.insert([2,6,7])
-    st.insert([1,2,7])
-    st.insert([0,2,3])
-    st.insert([2,3,5])
-    st.insert([3,4,5])
-    st.insert([4,5,6])
-    st.insert([0,4,6])
-    st.insert([0,2,6])
+    st.insert([0, 1, 8])
+    st.insert([0, 3, 8])
+    st.insert([3, 7, 8])
+    st.insert([3, 4, 7])
+    st.insert([1, 4, 7])
+    st.insert([0, 1, 4])
+    st.insert([1, 2, 5])
+    st.insert([1, 5, 8])
+    st.insert([5, 6, 8])
+    st.insert([6, 7, 8])
+    st.insert([2, 6, 7])
+    st.insert([1, 2, 7])
+    st.insert([0, 2, 3])
+    st.insert([2, 3, 5])
+    st.insert([3, 4, 5])
+    st.insert([4, 5, 6])
+    st.insert([0, 4, 6])
+    st.insert([0, 2, 6])
     st.compute_persistence(persistence_dim_max=True)
-    
+
     H0 = st.persistence_intervals_in_dimension(0)
-    assert np.array_equal(H0, np.array([[ 0., float("inf")]]))
+    assert np.array_equal(H0, np.array([[0.0, float("inf")]]))
     H1 = st.persistence_intervals_in_dimension(1)
-    assert np.array_equal(H1, np.array([[ 0., float("inf")], [ 0., float("inf")]]))
+    assert np.array_equal(H1, np.array([[0.0, float("inf")], [0.0, float("inf")]]))
     H2 = st.persistence_intervals_in_dimension(2)
-    assert np.array_equal(H2, np.array([[ 0., float("inf")]]))
+    assert np.array_equal(H2, np.array([[0.0, float("inf")]]))
     # Test empty case
     assert st.persistence_intervals_in_dimension(3).shape == (0, 2)
 
+
 def test_equality_operator():
     st1 = SimplexTree()
     st2 = SimplexTree()
 
     assert st1 == st2
 
-    st1.insert([1,2,3], 4.)
+    st1.insert([1, 2, 3], 4.0)
     assert st1 != st2
 
-    st2.insert([1,2,3], 4.)
+    st2.insert([1, 2, 3], 4.0)
     assert st1 == st2
 
+
 def test_simplex_tree_deep_copy():
     st = SimplexTree()
-    st.insert([1, 2, 3], 0.)
+    st.insert([1, 2, 3], 0.0)
     # compute persistence only on the original
     st.compute_persistence()
 
@@ -480,14 +488,15 @@ def test_simplex_tree_deep_copy():
 
     for a_splx in a_filt_list:
         assert a_splx in st_filt_list
-    
+
     # test double free
     del st
     del st_copy
 
+
 def test_simplex_tree_deep_copy_constructor():
     st = SimplexTree()
-    st.insert([1, 2, 3], 0.)
+    st.insert([1, 2, 3], 0.0)
     # compute persistence only on the original
     st.compute_persistence()
 
@@ -506,56 +515,132 @@ def test_simplex_tree_deep_copy_constructor():
 
     for a_splx in a_filt_list:
         assert a_splx in st_filt_list
-    
+
     # test double free
     del st
     del st_copy
 
+
 def test_simplex_tree_constructor_exception():
     with pytest.raises(TypeError):
-        st = SimplexTree(other = "Construction from a string shall raise an exception")
+        st = SimplexTree(other="Construction from a string shall raise an exception")
+
+
+def test_create_from_array():
+    a = np.array([[1, 4, 13, 6], [4, 3, 11, 5], [13, 11, 10, 12], [6, 5, 12, 2]])
+    st = SimplexTree.create_from_array(a, max_filtration=5.0)
+    assert list(st.get_filtration()) == [([0], 1.0), ([3], 2.0), ([1], 3.0), ([0, 1], 4.0), ([1, 3], 5.0)]
+
+
+def test_insert_edges_from_coo_matrix():
+    try:
+        from scipy.sparse import coo_matrix
+        from scipy.spatial import cKDTree
+    except ImportError:
+        print("Skipping, no SciPy")
+        return
+
+    st = SimplexTree()
+    st.insert([1, 2, 7], 7)
+    row = np.array([2, 5, 3])
+    col = np.array([1, 4, 6])
+    dat = np.array([1, 2, 3])
+    edges = coo_matrix((dat, (row, col)))
+    st.insert_edges_from_coo_matrix(edges)
+    assert list(st.get_filtration()) == [
+        ([1], 1.0),
+        ([2], 1.0),
+        ([1, 2], 1.0),
+        ([4], 2.0),
+        ([5], 2.0),
+        ([4, 5], 2.0),
+        ([3], 3.0),
+        ([6], 3.0),
+        ([3, 6], 3.0),
+        ([7], 7.0),
+        ([1, 7], 7.0),
+        ([2, 7], 7.0),
+        ([1, 2, 7], 7.0),
+    ]
+
+    pts = np.random.rand(100, 2)
+    tree = cKDTree(pts)
+    edges = tree.sparse_distance_matrix(tree, max_distance=0.15, output_type="coo_matrix")
+    st = SimplexTree()
+    st.insert_edges_from_coo_matrix(edges)
+    assert 100 < st.num_simplices() < 1000
+
+
+def test_insert_batch():
+    st = SimplexTree()
+    # vertices
+    st.insert_batch(np.array([[6, 1, 5]]), np.array([-5.0, 2.0, -3.0]))
+    # triangles
+    st.insert_batch(np.array([[2, 10], [5, 0], [6, 11]]), np.array([4.0, 0.0]))
+    # edges
+    st.insert_batch(np.array([[1, 5], [2, 5]]), np.array([1.0, 3.0]))
+
+    assert list(st.get_filtration()) == [
+        ([6], -5.0),
+        ([5], -3.0),
+        ([0], 0.0),
+        ([10], 0.0),
+        ([0, 10], 0.0),
+        ([11], 0.0),
+        ([0, 11], 0.0),
+        ([10, 11], 0.0),
+        ([0, 10, 11], 0.0),
+        ([1], 1.0),
+        ([2], 1.0),
+        ([1, 2], 1.0),
+        ([2, 5], 4.0),
+        ([2, 6], 4.0),
+        ([5, 6], 4.0),
+        ([2, 5, 6], 4.0),
+    ]
+
 
 def test_expansion_with_blocker():
-    st=SimplexTree()
-    st.insert([0,1],0)
-    st.insert([0,2],1)
-    st.insert([0,3],2)
-    st.insert([1,2],3)
-    st.insert([1,3],4)
-    st.insert([2,3],5)
-    st.insert([2,4],6)
-    st.insert([3,6],7)
-    st.insert([4,5],8)
-    st.insert([4,6],9)
-    st.insert([5,6],10)
-    st.insert([6],10)
+    st = SimplexTree()
+    st.insert([0, 1], 0)
+    st.insert([0, 2], 1)
+    st.insert([0, 3], 2)
+    st.insert([1, 2], 3)
+    st.insert([1, 3], 4)
+    st.insert([2, 3], 5)
+    st.insert([2, 4], 6)
+    st.insert([3, 6], 7)
+    st.insert([4, 5], 8)
+    st.insert([4, 6], 9)
+    st.insert([5, 6], 10)
+    st.insert([6], 10)
 
     def blocker(simplex):
         try:
             # Block all simplices that contain vertex 6
             simplex.index(6)
-            print(simplex, ' is blocked')
+            print(simplex, " is blocked")
             return True
         except ValueError:
-            print(simplex, ' is accepted')
-            st.assign_filtration(simplex, st.filtration(simplex) + 1.)
+            print(simplex, " is accepted")
+            st.assign_filtration(simplex, st.filtration(simplex) + 1.0)
             return False
 
     st.expansion_with_blocker(2, blocker)
     assert st.num_simplices() == 22
     assert st.dimension() == 2
-    assert st.find([4,5,6]) == False
-    assert st.filtration([0,1,2]) == 4.
-    assert st.filtration([0,1,3]) == 5.
-    assert st.filtration([0,2,3]) == 6.
-    assert st.filtration([1,2,3]) == 6.
+    assert st.find([4, 5, 6]) == False
+    assert st.filtration([0, 1, 2]) == 4.0
+    assert st.filtration([0, 1, 3]) == 5.0
+    assert st.filtration([0, 2, 3]) == 6.0
+    assert st.filtration([1, 2, 3]) == 6.0
 
     st.expansion_with_blocker(3, blocker)
     assert st.num_simplices() == 23
     assert st.dimension() == 3
-    assert st.find([4,5,6]) == False
-    assert st.filtration([0,1,2]) == 4.
-    assert st.filtration([0,1,3]) == 5.
-    assert st.filtration([0,2,3]) == 6.
-    assert st.filtration([1,2,3]) == 6.
-    assert st.filtration([0,1,2,3]) == 7.
+    assert st.find([4, 5, 6]) == False
+    assert st.filtration([0, 1, 2]) == 4.0
+    assert st.filtration([0, 1, 3]) == 5.0
+    assert st.filtration([0, 2, 3]) == 6.0
+    assert st.filtration([1, 2, 3]) == 6.0
+    assert st.filtration([0, 1, 2, 3]) == 7.0
diff --git a/src/python/test/test_subsampling.py b/src/python/test/test_subsampling.py
index 3431f372..c1cb4e3f 100755
--- a/src/python/test/test_subsampling.py
+++ b/src/python/test/test_subsampling.py
@@ -16,17 +16,9 @@ __license__ = "MIT"
 
 
 def test_write_off_file_for_tests():
-    file = open("subsample.off", "w")
-    file.write("nOFF\n")
-    file.write("2 7 0 0\n")
-    file.write("1.0 1.0\n")
-    file.write("7.0 0.0\n")
-    file.write("4.0 6.0\n")
-    file.write("9.0 6.0\n")
-    file.write("0.0 14.0\n")
-    file.write("2.0 19.0\n")
-    file.write("9.0 17.0\n")
-    file.close()
+    gudhi.write_points_to_off_file(
+        "subsample.off", [[1.0, 1.0], [7.0, 0.0], [4.0, 6.0], [9.0, 6.0], [0.0, 14.0], [2.0, 19.0], [9.0, 17.0]]
+    )
 
 
 def test_simple_choose_n_farthest_points_with_a_starting_point():
@@ -34,54 +26,29 @@ def test_simple_choose_n_farthest_points_with_a_starting_point():
     i = 0
     for point in point_set:
         # The iteration starts with the given starting point
-        sub_set = gudhi.choose_n_farthest_points(
-            points=point_set, nb_points=1, starting_point=i
-        )
+        sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=1, starting_point=i)
         assert sub_set[0] == point_set[i]
         i = i + 1
 
     # The iteration finds then the farthest
-    sub_set = gudhi.choose_n_farthest_points(
-        points=point_set, nb_points=2, starting_point=1
-    )
+    sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=1)
     assert sub_set[1] == point_set[3]
-    sub_set = gudhi.choose_n_farthest_points(
-        points=point_set, nb_points=2, starting_point=3
-    )
+    sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=3)
     assert sub_set[1] == point_set[1]
-    sub_set = gudhi.choose_n_farthest_points(
-        points=point_set, nb_points=2, starting_point=0
-    )
+    sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=0)
     assert sub_set[1] == point_set[2]
-    sub_set = gudhi.choose_n_farthest_points(
-        points=point_set, nb_points=2, starting_point=2
-    )
+    sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=2)
     assert sub_set[1] == point_set[0]
 
     # Test the limits
-    assert (
-        gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=0) == []
-    )
-    assert (
-        gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=0) == []
-    )
-    assert (
-        gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=1) == []
-    )
-    assert (
-        gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=1) == []
-    )
+    assert gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=0) == []
+    assert gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=0) == []
+    assert gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=1) == []
+    assert gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=1) == []
 
     # From off file test
     for i in range(0, 7):
-        assert (
-            len(
-                gudhi.choose_n_farthest_points(
-                    off_file="subsample.off", nb_points=i, starting_point=i
-                )
-            )
-            == i
-        )
+        assert len(gudhi.choose_n_farthest_points(off_file="subsample.off", nb_points=i, starting_point=i)) == i
 
 
 def test_simple_choose_n_farthest_points_randomed():
@@ -104,10 +71,7 @@ def test_simple_choose_n_farthest_points_randomed():
 
     # From off file test
     for i in range(0, 7):
-        assert (
-            len(gudhi.choose_n_farthest_points(off_file="subsample.off", nb_points=i))
-            == i
-        )
+        assert len(gudhi.choose_n_farthest_points(off_file="subsample.off", nb_points=i)) == i
 
 
 def test_simple_pick_n_random_points():
@@ -130,9 +94,7 @@ def test_simple_pick_n_random_points():
 
     # From off file test
     for i in range(0, 7):
-        assert (
-            len(gudhi.pick_n_random_points(off_file="subsample.off", nb_points=i)) == i
-        )
+        assert len(gudhi.pick_n_random_points(off_file="subsample.off", nb_points=i)) == i
 
 
 def test_simple_sparsify_points():
@@ -152,31 +114,10 @@ def test_simple_sparsify_points():
     ]
     assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.001) == [[0, 1]]
 
-    assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=0.0))
-        == 7
-    )
-    assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=30.0))
-        == 5
-    )
-    assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.1))
-        == 4
-    )
-    assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=89.9))
-        == 3
-    )
-    assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=100.0))
-        == 2
-    )
-    assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=324.9))
-        == 2
-    )
-    assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=325.01))
-        == 1
-    )
+    assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=0.0)) == 7
+    assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=30.0)) == 5
+    assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.1)) == 4
+    assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=89.9)) == 3
+    assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=100.0)) == 2
+    assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=324.9)) == 2
+    assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=325.01)) == 1