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

15 files changed, 303 insertions, 127 deletions

diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index aec8c1b1..a7372f19 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -17,8 +17,7 @@
 // to construct Alpha_complex from a OFF file of points
 #include <gudhi/Points_off_io.h>
 
-#include <stdlib.h>
-#include <math.h>  // isnan, fmax
+#include <cmath>  // isnan, fmax
 #include <memory>  // for std::unique_ptr
 #include <cstddef>  // for std::size_t
 
@@ -45,6 +44,7 @@
 #include <utility>  // std::pair
 #include <stdexcept>
 #include <numeric>  // for std::iota
+#include <algorithm>  // for std::sort
 
 // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
 #if CGAL_VERSION_NR < 1041101000
@@ -101,13 +101,17 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val
  */
 template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>, bool Weighted = false>
 class Alpha_complex {
+ private:
+  // Vertex_handle internal type (required by triangulation_ and vertices_).
+  using Internal_vertex_handle = std::ptrdiff_t;
+
  public:
   /** \brief Geometric traits class that provides the geometric types and predicates needed by the triangulations.*/
   using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>;
 
   // Add an int in TDS to save point index in the structure
   using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension,
-                                                 CGAL::Triangulation_vertex<Geom_traits, std::ptrdiff_t>,
+                                                 CGAL::Triangulation_vertex<Geom_traits, Internal_vertex_handle>,
                                                  CGAL::Triangulation_full_cell<Geom_traits> >;
 
   /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/
@@ -132,9 +136,6 @@ class Alpha_complex {
   // Vertex_iterator type from CGAL.
   using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator;
 
-  // size_type type from CGAL.
-  using size_type = typename Triangulation::size_type;
-
   // Structure to switch from simplex tree vertex handle to CGAL vertex iterator.
   using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >;
 
@@ -146,6 +147,10 @@ class Alpha_complex {
   std::unique_ptr<Triangulation> triangulation_;
   /** \brief Kernel for triangulation_ functions access.*/
   A_kernel_d kernel_;
+  /** \brief Vertices to be inserted first by the create_complex method to avoid quadratic complexity.
+   * It isn't just [0, n) if some points have multiplicity (only one copy appears in the complex).
+  */
+  std::vector<Internal_vertex_handle> vertices_;
 
   /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/
   std::vector<Sphere> cache_, old_cache_;
@@ -257,11 +262,11 @@ class Alpha_complex {
       std::vector<Point_d> point_cloud(first, last);
 
       // Creates a vector {0, 1, ..., N-1}
-      std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
-                                          boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
+      std::vector<Internal_vertex_handle> indices(boost::counting_iterator<Internal_vertex_handle>(0),
+                                                  boost::counting_iterator<Internal_vertex_handle>(point_cloud.size()));
 
       using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator,
-                                                              CGAL::Identity_property_map<std::ptrdiff_t>>;
+                                                              CGAL::Identity_property_map<Internal_vertex_handle>>;
       using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d<Geom_traits, Point_property_map>;
 
       CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
@@ -279,6 +284,9 @@ class Alpha_complex {
       // structure to retrieve CGAL points from vertex handle - one vertex handle per point.
       // Needs to be constructed before as vertex handles arrives in no particular order.
       vertex_handle_to_iterator_.resize(point_cloud.size());
+      // List of sorted unique vertices in the triangulation. We take advantage of the existing loop to construct it
+      // Vertices list avoids quadratic complexity with the Simplex_tree. We should not fill it up with Toplex_map e.g.
+      vertices_.reserve(triangulation_->number_of_vertices());
       // Loop on triangulation vertices list
       for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
         if (!triangulation_->is_infinite(*vit)) {
@@ -286,8 +294,10 @@ class Alpha_complex {
           std::clog << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl;
 #endif  // DEBUG_TRACES
           vertex_handle_to_iterator_[vit->data()] = vit;
+          vertices_.push_back(vit->data());
         }
       }
+      std::sort(vertices_.begin(), vertices_.end());
       // --------------------------------------------------------------------------------------------
     }
   }
@@ -384,12 +394,21 @@ class Alpha_complex {
     // --------------------------------------------------------------------------------------------
     // Simplex_tree construction from loop on triangulation finite full cells list
     if (num_vertices() > 0) {
+      std::vector<Vertex_handle> one_vertex(1);
+      for (auto vertex : vertices_) {
+#ifdef DEBUG_TRACES
+        std::clog << "SimplicialComplex insertion " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        one_vertex[0] = vertex;
+        complex.insert_simplex_and_subfaces(one_vertex, std::numeric_limits<double>::quiet_NaN());
+      }
+
       for (auto cit = triangulation_->finite_full_cells_begin();
            cit != triangulation_->finite_full_cells_end();
            ++cit) {
         Vector_vertex vertexVector;
 #ifdef DEBUG_TRACES
-        std::clog << "Simplex_tree insertion ";
+        std::clog << "SimplicialComplex insertion ";
 #endif  // DEBUG_TRACES
         for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
           if (*vit != nullptr) {
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/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/Tangential_complex/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h
index cc424810..56a24af0 100644
--- a/src/Tangential_complex/include/gudhi/Tangential_complex.h
+++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h
@@ -36,7 +36,6 @@
 #include <Eigen/Eigen>
 #include <Eigen/src/Core/util/Macros.h>  // for EIGEN_VERSION_AT_LEAST
 
-#include <boost/optional.hpp>
 #include <boost/iterator/transform_iterator.hpp>
 #include <boost/range/adaptor/transformed.hpp>
 #include <boost/range/counting_range.hpp>
@@ -56,6 +55,7 @@
 #include <cmath>  // for std::sqrt
 #include <string>
 #include <cstddef>  // for std::size_t
+#include <optional>
 
 #ifdef GUDHI_USE_TBB
 #include <tbb/parallel_for.h>
@@ -994,7 +994,7 @@ class Tangential_complex {
     // circumspheres of the star of "center_vertex"
     // If th the m_max_squared_edge_length is set the maximal radius of the "star sphere"
     // is at most square root of m_max_squared_edge_length
-    boost::optional<FT> squared_star_sphere_radius_plus_margin = m_max_squared_edge_length;
+    std::optional<FT> squared_star_sphere_radius_plus_margin = m_max_squared_edge_length;
 
     // Insert points until we find a point which is outside "star sphere"
     for (auto nn_it = ins_range.begin(); nn_it != ins_range.end(); ++nn_it) {
@@ -1036,7 +1036,7 @@ class Tangential_complex {
 
           // Let's recompute squared_star_sphere_radius_plus_margin
           if (triangulation.current_dimension() >= tangent_space_dim) {
-            squared_star_sphere_radius_plus_margin = boost::none;
+            squared_star_sphere_radius_plus_margin = std::nullopt;
             // Get the incident cells and look for the biggest circumsphere
             std::vector<Tr_full_cell_handle> incident_cells;
             triangulation.incident_full_cells(center_vertex, std::back_inserter(incident_cells));
@@ -1044,7 +1044,7 @@ class Tangential_complex {
                  cit != incident_cells.end(); ++cit) {
               Tr_full_cell_handle cell = *cit;
               if (triangulation.is_infinite(cell)) {
-                squared_star_sphere_radius_plus_margin = boost::none;
+                squared_star_sphere_radius_plus_margin = std::nullopt;
                 break;
               } else {
                 // Note that this uses the perturbed point since it uses
@@ -2030,7 +2030,7 @@ class Tangential_complex {
   // and their center vertex
   Stars_container m_stars;
   std::vector<FT> m_squared_star_spheres_radii_incl_margin;
-  boost::optional<FT> m_max_squared_edge_length;
+  std::optional<FT> m_max_squared_edge_length;
 
 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
   Points m_points_for_tse;
diff --git a/src/cmake/modules/GUDHI_compilation_flags.cmake b/src/cmake/modules/GUDHI_compilation_flags.cmake
index e2d3d872..b43ccf73 100644
--- a/src/cmake/modules/GUDHI_compilation_flags.cmake
+++ b/src/cmake/modules/GUDHI_compilation_flags.cmake
@@ -12,6 +12,7 @@ macro(add_cxx_compiler_flag _flag)
 endmacro()
 
 set (CMAKE_CXX_STANDARD 17)
+# This number needs to be changed in python/CMakeLists.txt at the same time
 
 enable_testing()
 
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 643c3c91..d5584475 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -129,9 +129,10 @@ if(PYTHONINTERP_FOUND)
 
     # Gudhi and CGAL compilation option
     if(MSVC)
+      set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'/std:c++17', ")
       set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'/fp:strict', ")
     else(MSVC)
-      set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-std=c++14', ")
+      set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-std=c++17', ")
     endif(MSVC)
     if(CMAKE_COMPILER_IS_GNUCXX)
       set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-frounding-math', ")
@@ -246,8 +247,8 @@ if(PYTHONINTERP_FOUND)
 
     # Specific for Mac
     if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
-        set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-mmacosx-version-min=10.12', ")
-        set(GUDHI_PYTHON_EXTRA_LINK_ARGS "${GUDHI_PYTHON_EXTRA_LINK_ARGS}'-mmacosx-version-min=10.12', ")
+        set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-mmacosx-version-min=10.14', ")
+        set(GUDHI_PYTHON_EXTRA_LINK_ARGS "${GUDHI_PYTHON_EXTRA_LINK_ARGS}'-mmacosx-version-min=10.14', ")
     endif(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
 
     # Loop on INCLUDE_DIRECTORIES PROPERTY
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 50ddabfe..d5dfefd7 100644
--- a/src/python/doc/installation.rst
+++ b/src/python/doc/installation.rst
@@ -399,7 +399,7 @@ The :doc:`persistence graphical tools </persistence_graphical_tools_user>` and
 mathematics, science, and engineering.
 
 :class:`~gudhi.point_cloud.knn.KNearestNeighbors` can use the Python package
-`SciPy <http://scipy.org>`_ as a backend if explicitly requested.
+`SciPy <http://scipy.org>`_ :math:`\geq` 1.6.0 as a backend if explicitly requested.
 
 TensorFlow
 ----------
diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py
index de5844f9..7dc83817 100644
--- a/src/python/gudhi/point_cloud/knn.py
+++ b/src/python/gudhi/point_cloud/knn.py
@@ -314,7 +314,9 @@ class KNearestNeighbors:
             return None
 
         if self.params["implementation"] == "ckdtree":
-            qargs = {key: val for key, val in self.params.items() if key in {"p", "eps", "n_jobs"}}
+            qargs = {key: val for key, val in self.params.items() if key in {"p", "eps"}}
+            # SciPy renamed n_jobs to workers
+            qargs["workers"] = self.params.get("workers") or self.params.get("n_jobs") or 1
             distances, neighbors = self.kdtree.query(X, k=self.k, **qargs)
             if k == 1:
                 # SciPy decided to squeeze the last dimension for k=1
diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py
index 69ff5e1e..a169aee8 100644
--- a/src/python/gudhi/representations/vector_methods.py
+++ b/src/python/gudhi/representations/vector_methods.py
@@ -85,7 +85,7 @@ class PersistenceImage(BaseEstimator, TransformerMixin):
 
             Xfit.append(image.flatten()[np.newaxis,:])
 
-        Xfit = np.concatenate(Xfit,0)
+        Xfit = np.concatenate(Xfit, 0)
 
         return Xfit
 
@@ -123,6 +123,15 @@ def _automatic_sample_range(sample_range, X, y):
                 pass
         return sample_range
 
+
+def _trim_on_edges(x, are_endpoints_nan):
+    if are_endpoints_nan[0]:
+        x = x[1:]
+    if are_endpoints_nan[1]:
+        x = x[:-1]
+    return x
+
+
 class Landscape(BaseEstimator, TransformerMixin):
     """
     This is a class for computing persistence landscapes from a list of persistence diagrams. A persistence landscape is a collection of 1D piecewise-linear functions computed from the rank function associated to the persistence diagram. These piecewise-linear functions are then sampled evenly on a given range and the corresponding vectors of samples are concatenated and returned. See http://jmlr.org/papers/v16/bubenik15a.html for more details.
@@ -149,6 +158,8 @@ class Landscape(BaseEstimator, TransformerMixin):
             y (n x 1 array): persistence diagram labels (unused).
         """
         self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
+        self.im_range = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution)
+        self.im_range = _trim_on_edges(self.im_range, self.nan_in_range)
         return self
 
     def transform(self, X):
@@ -161,53 +172,26 @@ class Landscape(BaseEstimator, TransformerMixin):
         Returns:
             numpy array with shape (number of diagrams) x (number of samples = **num_landscapes** x **resolution**): output persistence landscapes.
         """
-        num_diag, Xfit = len(X), []
-        x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution)
-        step_x = x_values[1] - x_values[0]
-
-        for i in range(num_diag):
-
-            diagram, num_pts_in_diag = X[i], X[i].shape[0]
 
-            ls = np.zeros([self.num_landscapes, self.new_resolution])
+        Xfit = []
+        x_values = self.im_range
+        for diag in X:
+            midpoints, heights = (diag[:, 0] + diag[:, 1]) / 2., (diag[:, 1] - diag[:, 0]) / 2.
+            tent_functions = np.maximum(heights[None, :] - np.abs(x_values[:, None] - midpoints[None, :]), 0)
+            n_points = diag.shape[0]
+            # Complete the array with zeros to get the right number of landscapes
+            if self.num_landscapes > n_points:
+                tent_functions = np.concatenate(
+                    [tent_functions, np.zeros((tent_functions.shape[0], self.num_landscapes-n_points))],
+                    axis=1
+                )
+            tent_functions.partition(tent_functions.shape[1]-self.num_landscapes, axis=1)
+            landscapes = np.sort(tent_functions[:, -self.num_landscapes:], axis=1)[:, ::-1].T
 
-            events = []
-            for j in range(self.new_resolution):
-                events.append([])
+            landscapes = np.sqrt(2) * np.ravel(landscapes)
+            Xfit.append(landscapes)
 
-            for j in range(num_pts_in_diag):
-                [px,py] = diagram[j,:2]
-                min_idx = np.clip(np.ceil((px          - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution)
-                mid_idx = np.clip(np.ceil((0.5*(py+px) - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution)
-                max_idx = np.clip(np.ceil((py          - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution)
-
-                if min_idx < self.new_resolution and max_idx > 0:
-
-                    landscape_value = self.sample_range[0] + min_idx * step_x - px
-                    for k in range(min_idx, mid_idx):
-                        events[k].append(landscape_value)
-                        landscape_value += step_x
-
-                    landscape_value = py - self.sample_range[0] - mid_idx * step_x
-                    for k in range(mid_idx, max_idx):
-                        events[k].append(landscape_value)
-                        landscape_value -= step_x
-
-            for j in range(self.new_resolution):
-                events[j].sort(reverse=True)
-                for k in range( min(self.num_landscapes, len(events[j])) ):
-                    ls[k,j] = events[j][k]
-
-            if self.nan_in_range[0]:
-                ls = ls[:,1:]
-            if self.nan_in_range[1]:
-                ls = ls[:,:-1]
-            ls = np.sqrt(2)*np.reshape(ls,[1,-1])
-            Xfit.append(ls)
-
-        Xfit = np.concatenate(Xfit,0)
-
-        return Xfit
+        return np.stack(Xfit, axis=0)
 
     def __call__(self, diag):
         """
@@ -219,7 +203,7 @@ class Landscape(BaseEstimator, TransformerMixin):
         Returns:
             numpy array with shape (number of samples = **num_landscapes** x **resolution**): output persistence landscape.
         """
-        return self.fit_transform([diag])[0,:]
+        return self.fit_transform([diag])[0, :]
 
 class Silhouette(BaseEstimator, TransformerMixin):
     """
@@ -235,6 +219,8 @@ class Silhouette(BaseEstimator, TransformerMixin):
             sample_range ([double, double]): minimum and maximum for the weighted average domain, of the form [x_min, x_max] (default [numpy.nan, numpy.nan]). It is the interval on which samples will be drawn evenly. If one of the values is numpy.nan, it can be computed from the persistence diagrams with the fit() method.
         """
         self.weight, self.resolution, self.sample_range = weight, resolution, sample_range
+        self.nan_in_range = np.isnan(np.array(self.sample_range))
+        self.new_resolution = self.resolution + self.nan_in_range.sum()
 
     def fit(self, X, y=None):
         """
@@ -245,6 +231,8 @@ class Silhouette(BaseEstimator, TransformerMixin):
             y (n x 1 array): persistence diagram labels (unused).
         """
         self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
+        self.im_range = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution)
+        self.im_range = _trim_on_edges(self.im_range, self.nan_in_range)
         return self
 
     def transform(self, X):
@@ -257,44 +245,19 @@ class Silhouette(BaseEstimator, TransformerMixin):
         Returns:
             numpy array with shape (number of diagrams) x (**resolution**): output persistence silhouettes.
         """
-        num_diag, Xfit = len(X), []
-        x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.resolution)
-        step_x = x_values[1] - x_values[0]
-
-        for i in range(num_diag):
-
-            diagram, num_pts_in_diag = X[i], X[i].shape[0]
+        Xfit = []
+        x_values = self.im_range
 
-            sh, weights = np.zeros(self.resolution), np.zeros(num_pts_in_diag)
-            for j in range(num_pts_in_diag):
-                weights[j] = self.weight(diagram[j,:])
+        for diag in X:
+            midpoints, heights = (diag[:, 0] + diag[:, 1]) / 2., (diag[:, 1] - diag[:, 0]) / 2.
+            weights = np.array([self.weight(pt) for pt in diag])
             total_weight = np.sum(weights)
 
-            for j in range(num_pts_in_diag):
-
-                [px,py] = diagram[j,:2]
-                weight  = weights[j] / total_weight
-                min_idx = np.clip(np.ceil((px          - self.sample_range[0]) / step_x).astype(int), 0, self.resolution)
-                mid_idx = np.clip(np.ceil((0.5*(py+px) - self.sample_range[0]) / step_x).astype(int), 0, self.resolution)
-                max_idx = np.clip(np.ceil((py          - self.sample_range[0]) / step_x).astype(int), 0, self.resolution)
-
-                if min_idx < self.resolution and max_idx > 0:
-
-                    silhouette_value = self.sample_range[0] + min_idx * step_x - px
-                    for k in range(min_idx, mid_idx):
-                        sh[k] += weight * silhouette_value
-                        silhouette_value += step_x
-
-                    silhouette_value = py - self.sample_range[0] - mid_idx * step_x
-                    for k in range(mid_idx, max_idx):
-                        sh[k] += weight * silhouette_value
-                        silhouette_value -= step_x
-
-            Xfit.append(np.reshape(np.sqrt(2) * sh, [1,-1]))
-
-        Xfit = np.concatenate(Xfit, 0)
+            tent_functions = np.maximum(heights[None, :] - np.abs(x_values[:, None] - midpoints[None, :]), 0)
+            silhouette = np.sum(weights[None, :] / total_weight * tent_functions, axis=1)
+            Xfit.append(silhouette * np.sqrt(2))
 
-        return Xfit 
+        return np.stack(Xfit, axis=0)
 
     def __call__(self, diag):
         """
diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py
index 4a455bb6..58caab21 100755
--- a/src/python/test/test_representations.py
+++ b/src/python/test/test_representations.py
@@ -187,3 +187,67 @@ def test_kernel_empty_diagrams():
 #    PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1.)(empty_diag, empty_diag)
 #    PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1., kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))(empty_diag, empty_diag)
 
+
+def test_silhouette_permutation_invariance():
+    dgm = _n_diags(1)[0]
+    dgm_permuted = dgm[np.random.permutation(dgm.shape[0]).astype(int)]
+    random_resolution = random.randint(50, 100) * 10
+    slt = Silhouette(resolution=random_resolution, weight=pow(2))
+
+    assert np.all(np.isclose(slt(dgm), slt(dgm_permuted)))
+
+
+def test_silhouette_multiplication_invariance():
+    dgm = _n_diags(1)[0]
+    n_repetitions = np.random.randint(2, high=10)
+    dgm_augmented = np.repeat(dgm, repeats=n_repetitions, axis=0)
+
+    random_resolution = random.randint(50, 100) * 10
+    slt = Silhouette(resolution=random_resolution, weight=pow(2))
+    assert np.all(np.isclose(slt(dgm), slt(dgm_augmented)))
+
+
+def test_silhouette_numeric():
+    dgm = np.array([[2., 3.], [5., 6.]])
+    slt = Silhouette(resolution=9, weight=pow(1), sample_range=[2., 6.])
+    #slt.fit([dgm])
+    # x_values = array([2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6.])
+
+    expected_silhouette = np.array([0., 0.5, 0., 0., 0., 0., 0., 0.5, 0.])/np.sqrt(2)
+    output_silhouette = slt(dgm)
+    assert np.all(np.isclose(output_silhouette, expected_silhouette))
+
+
+def test_landscape_small_persistence_invariance():
+    dgm = np.array([[2., 6.], [2., 5.], [3., 7.]])
+    small_persistence_pts = np.random.rand(10, 2)
+    small_persistence_pts[:, 1] += small_persistence_pts[:, 0]
+    small_persistence_pts += np.min(dgm)
+    dgm_augmented = np.concatenate([dgm, small_persistence_pts], axis=0)
+
+    lds = Landscape(num_landscapes=2, resolution=5)
+    lds_dgm, lds_dgm_augmented = lds(dgm), lds(dgm_augmented)
+
+    assert np.all(np.isclose(lds_dgm, lds_dgm_augmented))
+
+
+def test_landscape_numeric():
+    dgm = np.array([[2., 6.], [3., 5.]])
+    lds_ref = np.array([
+        0., 0.5, 1., 1.5, 2., 1.5, 1., 0.5, 0.,  # tent of [2, 6]
+        0., 0., 0., 0.5, 1., 0.5, 0., 0., 0.,
+        0., 0., 0., 0., 0., 0., 0., 0., 0.,
+        0., 0., 0., 0., 0., 0., 0., 0., 0.,
+    ])
+    lds_ref *= np.sqrt(2)
+    lds = Landscape(num_landscapes=4, resolution=9, sample_range=[2., 6.])
+    lds_dgm = lds(dgm)
+    assert np.all(np.isclose(lds_dgm, lds_ref))
+
+
+def test_landscape_nan_range():
+    dgm = np.array([[2., 6.], [3., 5.]])
+    lds = Landscape(num_landscapes=2, resolution=9, sample_range=[np.nan, 6.])
+    lds_dgm = lds(dgm)
+    assert (lds.sample_range[0] == 2) & (lds.sample_range[1] == 6)
+    assert lds.new_resolution == 10