Remove comments and identify both contributions. Include what is used and add some std::. Use std::sqrt instead of std::pow(x,0.5)

4 files changed, 381 insertions, 303 deletions

diff --git a/src/Persistence_representations/example/persistence_heat_maps.cpp b/src/Persistence_representations/example/persistence_heat_maps.cpp
index 0beb1052..45208b68 100644
--- a/src/Persistence_representations/example/persistence_heat_maps.cpp
+++ b/src/Persistence_representations/example/persistence_heat_maps.cpp
@@ -2,9 +2,12 @@
  *    (Geometric Understanding in Higher Dimensions) is a generic C++
  *    library for computational topology.
  *
- *    Author(s):       Pawel Dlotko
+ *    Author(s):       Pawel Dlotko and Mathieu Carriere
  *
- *    Copyright (C) 2016 Inria
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modifications:
+ *      - 2018/04 MC: Add persistence heat maps computation
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -26,10 +29,14 @@
 #include <iostream>
 #include <vector>
 #include <utility>
-
-
-std::function<double(std::pair<double, double>, std::pair<double, double>)> Gaussian_function(double sigma){
-  return [=](std::pair<double, double> p, std::pair<double, double> q){return std::exp(  -( (p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second) )/(sigma)  );};
+#include <functional>
+#include <cmath>
+
+std::function<double(std::pair<double, double>, std::pair<double, double>)> Gaussian_function(double sigma) {
+  return [=](std::pair<double, double> p, std::pair<double, double> q) {
+    return std::exp(-((p.first - q.first) * (p.first - q.first) + (p.second - q.second) * (p.second - q.second)) /
+                    (sigma));
+  };
 }
 
 using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function;
@@ -82,14 +89,13 @@ int main(int argc, char** argv) {
   // to compute scalar product of hm1 and hm2:
   std::cout << "Scalar product is : " << hm1.compute_scalar_product(hm2) << std::endl;
 
-  // Mathieu's code ************************************************************************************************************
   Persistence_heat_maps hm1k(persistence1, Gaussian_function(1.0));
   Persistence_heat_maps hm2k(persistence2, Gaussian_function(1.0));
   Persistence_heat_maps hm1i(persistence1, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11);
   Persistence_heat_maps hm2i(persistence2, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11);
-  std::cout << "Scalar product computed with exact 2D kernel on grid is : " << hm1i.compute_scalar_product(hm2i) << std::endl;
-  std::cout << "Scalar product computed with exact 2D kernel is : "         << hm1k.compute_scalar_product(hm2k) << std::endl;
-  // ***************************************************************************************************************************
+  std::cout << "Scalar product computed with exact 2D kernel on grid is : " << hm1i.compute_scalar_product(hm2i)
+            << std::endl;
+  std::cout << "Scalar product computed with exact 2D kernel is : " << hm1k.compute_scalar_product(hm2k) << std::endl;
 
   return 0;
 }
diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
index db51cc14..a8458bda 100644
--- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
+++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
@@ -2,9 +2,12 @@
  *    (Geometric Understanding in Higher Dimensions) is a generic C++
  *    library for computational topology.
  *
- *    Author(s):       Pawel Dlotko
+ *    Author(s):       Pawel Dlotko and Mathieu Carriere
  *
- *    Copyright (C) 2016 Inria
+ *    Modifications:
+ *      - 2018/04 MC: Add discrete/non-discrete mechanism and non-discrete version
+ *
+ *    Copyright (C) 2019 Inria
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -44,7 +47,7 @@ namespace Persistence_representations {
 /**
  * This is a simple procedure to create n by n (or 2*pixel_radius times 2*pixel_radius cubical approximation of a
  *Gaussian kernel.
-**/
+ **/
 std::vector<std::vector<double> > create_Gaussian_filter(size_t pixel_radius, double sigma) {
   bool dbg = false;
   // we are computing the kernel mask to 2 standard deviations away from the center. We discretize it in a grid of a
@@ -74,7 +77,7 @@ std::vector<std::vector<double> > create_Gaussian_filter(size_t pixel_radius, do
     for (int y = -pixel_radius; y <= static_cast<int>(pixel_radius); y++) {
       double real_x = 2 * sigma * x / pixel_radius;
       double real_y = 2 * sigma * y / pixel_radius;
-      r = sqrt(real_x * real_x + real_y * real_y);
+      r = std::sqrt(real_x * real_x + real_y * real_y);
       kernel[x + pixel_radius][y + pixel_radius] = (exp(-(r * r) / sigma_sqr)) / (3.141592 * sigma_sqr);
       sum += kernel[x + pixel_radius][y + pixel_radius];
     }
@@ -100,18 +103,18 @@ std::vector<std::vector<double> > create_Gaussian_filter(size_t pixel_radius, do
 }
 
 /*
-* There are various options to scale the points depending on their location. One can for instance:
-* (1) do nothing (scale all of them with the weight 1), as in the function constant_function
-* (2) Scale them by the distance to the diagonal. This is implemented in function
-* (3) Scale them with the square of their distance to diagonal. This is implemented in function
-* (4) Scale them with
-*/
+ * There are various options to scale the points depending on their location. One can for instance:
+ * (1) do nothing (scale all of them with the weight 1), as in the function constant_function
+ * (2) Scale them by the distance to the diagonal. This is implemented in function
+ * (3) Scale them with the square of their distance to diagonal. This is implemented in function
+ * (4) Scale them with
+ */
 
 /**
  * This is one of a scaling functions used to weight points depending on their persistence and/or location in the
  *diagram.
  * This particular functionality is a function which always assign value 1 to a point in the diagram.
-**/
+ **/
 class constant_scaling_function {
  public:
   double operator()(const std::pair<double, double>& point_in_diagram) { return 1; }
@@ -121,13 +124,13 @@ class constant_scaling_function {
  * This is one of a scaling functions used to weight points depending on their persistence and/or location in the
  *diagram.
  * The scaling given by this function to a point (b,d) is Euclidean distance of (b,d) from diagonal.
-**/
+ **/
 class distance_from_diagonal_scaling {
  public:
   double operator()(const std::pair<double, double>& point_in_diagram) {
     // (point_in_diagram.first+point_in_diagram.second)/2.0
-    return sqrt(pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) +
-                pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2));
+    return std::sqrt(std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) +
+                     std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2));
   }
 };
 
@@ -135,12 +138,12 @@ class distance_from_diagonal_scaling {
  * This is one of a scaling functions used to weight points depending on their persistence and/or location in the
  *diagram.
  * The scaling given by this function to a point (b,d) is a square of Euclidean distance of (b,d) from diagonal.
-**/
+ **/
 class squared_distance_from_diagonal_scaling {
  public:
   double operator()(const std::pair<double, double>& point_in_diagram) {
-    return pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) +
-           pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2);
+    return std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) +
+           std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2);
   }
 };
 
@@ -148,7 +151,7 @@ class squared_distance_from_diagonal_scaling {
  * This is one of a scaling functions used to weight points depending on their persistence and/or location in the
  *diagram.
  * The scaling given by this function to a point (b,d) is an arctan of a persistence of a point (i.e. arctan( b-d ).
-**/
+ **/
 class arc_tan_of_persistence_of_point {
  public:
   double operator()(const std::pair<double, double>& point_in_diagram) {
@@ -162,32 +165,24 @@ class arc_tan_of_persistence_of_point {
  * This scaling function do not only depend on a point (p,d) in the diagram, but it depends on the whole diagram.
  * The longest persistence pair get a scaling 1. Any other pair get a scaling belong to [0,1], which is proportional
  * to the persistence of that pair.
-**/
+ **/
 class weight_by_setting_maximal_interval_to_have_length_one {
  public:
-  weight_by_setting_maximal_interval_to_have_length_one(double len) : letngth_of_maximal_interval(len) {}
+  weight_by_setting_maximal_interval_to_have_length_one(double len) : length_of_maximal_interval(len) {}
   double operator()(const std::pair<double, double>& point_in_diagram) {
-    return (point_in_diagram.second - point_in_diagram.first) / this->letngth_of_maximal_interval;
+    return (point_in_diagram.second - point_in_diagram.first) / this->length_of_maximal_interval;
   }
 
  private:
-  double letngth_of_maximal_interval;
+  double length_of_maximal_interval;
 };
 
-// Mathieu's code ***********************************************************************************************************************
-//class Gaussian_kernel_2D {
-// public:
-//  double operator()(const std::pair<double, double>& p, const std::pair<double, double>& q) const {
-//    return (1.0/(std::sqrt(2*pi))) * std::exp(-((p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second)) / 2); }
-//};
-// **************************************************************************************************************************************
-
 /**
  * \class Persistence_heat_maps Persistence_heat_maps.h gudhi/Persistence_heat_maps.h
  * \brief A class implementing persistence heat maps.
  *
  * \ingroup Persistence_representations
-**/
+ **/
 
 // This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances,
 // Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product
@@ -197,7 +192,7 @@ class Persistence_heat_maps {
   /**
    * The default constructor. A scaling function from the diagonal is set up to a constant function. The image is not
    *erased below the diagonal. The Gaussian have diameter 5.
-  **/
+   **/
   Persistence_heat_maps() {
     Scalling_of_kernels f;
     this->f = f;
@@ -218,7 +213,7 @@ class Persistence_heat_maps {
    *std::numeric_limits<double>::max(), in which case the program compute the values based on the data,
    * (6) a max x and y value of points that are to be taken into account. By default it is set to
    *std::numeric_limits<double>::max(), in which case the program compute the values based on the data.
-  **/
+   **/
   Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
                         std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1),
                         bool erase_below_diagonal = false, size_t number_of_pixels = 1000,
@@ -226,12 +221,12 @@ class Persistence_heat_maps {
                         double max_ = std::numeric_limits<double>::max());
 
   /**
-       * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by
-    *default), a scaling function (constant by default), a boolean value which determines if the area of image below
-    *diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each
-    *direction (set to 1000 by default) and min and max values of images (both set to std::numeric_limits<double>::max()
-    *by default. If this is the case, the program will pick the right values based on the data).
-      **/
+   * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by
+   *default), a scaling function (constant by default), a boolean value which determines if the area of image below
+   *diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each
+   *direction (set to 1000 by default) and min and max values of images (both set to std::numeric_limits<double>::max()
+   *by default. If this is the case, the program will pick the right values based on the data).
+   **/
   /**
    * Construction that takes at the input the following parameters:
    * (1) A name of a file with persistence intervals. The file should be readable by the function
@@ -245,48 +240,45 @@ class Persistence_heat_maps {
    *std::numeric_limits<double>::max(), in which case the program compute the values based on the data,
    * (6) a max x and y value of points that are to be taken into account. By default it is set to
    *std::numeric_limits<double>::max(), in which case the program compute the values based on the data.
-  **/
+   **/
   Persistence_heat_maps(const char* filename, std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1),
                         bool erase_below_diagonal = false, size_t number_of_pixels = 1000,
                         double min_ = std::numeric_limits<double>::max(),
                         double max_ = std::numeric_limits<double>::max(),
                         unsigned dimension = std::numeric_limits<unsigned>::max());
 
-// Mathieu's code ***********************************************************************************************************************
   /**
-   * Construction that takes as inputs (1) the diagram, and (2) the grid parameters (min, max and number of samples for x and y axes)
-  **/
+   * Construction that takes as inputs (1) the diagram, and (2) the grid parameters (min, max and number of samples for
+   * x and y axes)
+   **/
   Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
                         const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel,
-                        size_t number_of_x_pixels, size_t number_of_y_pixels,
-                        double min_x = 0, double max_x = 1, double min_y = 0, double max_y = 1);
+                        size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x = 0, double max_x = 1,
+                        double min_y = 0, double max_y = 1);
 
   /**
    * Construction that takes only the diagram as input (weight and 2D kernel are template parameters)
-  **/
+   **/
   Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
                         const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel);
-// **************************************************************************************************************************************
-
-
 
   /**
    * Compute a mean value of a collection of heat maps and store it in the current object. Note that all the persistence
    *maps send in a vector to this procedure need to have the same parameters.
    * If this is not the case, the program will throw an exception.
-  **/
+   **/
   void compute_mean(const std::vector<Persistence_heat_maps*>& maps);
 
   /**
    * Compute a median value of a collection of heat maps and store it in the current object. Note that all the
    *persistence maps send in a vector to this procedure need to have the same parameters.
    * If this is not the case, the program will throw an exception.
-  **/
+   **/
   void compute_median(const std::vector<Persistence_heat_maps*>& maps);
 
   /**
    * Compute a percentage of active (i.e) values above the cutoff of a collection of heat maps.
-  **/
+   **/
   void compute_percentage_of_active(const std::vector<Persistence_heat_maps*>& maps, size_t cutoff = 1);
 
   // put to file subroutine
@@ -294,18 +286,18 @@ class Persistence_heat_maps {
    * The function outputs the persistence image to a text file. The format as follow:
    * In the first line, the values min and max of the image are stored
    * In the next lines, we have the persistence images in a form of a bitmap image.
-  **/
+   **/
   void print_to_file(const char* filename) const;
 
   /**
    * A function that load a heat map from file to the current object (and erase whatever was stored in the current
    *object before).
-  **/
+   **/
   void load_from_file(const char* filename);
 
   /**
    * The procedure checks if min_, max_ and this->heat_maps sizes are the same.
-  **/
+   **/
   inline bool check_if_the_same(const Persistence_heat_maps& second) const {
     bool dbg = false;
     if (this->heat_map.size() != second.heat_map.size()) {
@@ -328,17 +320,17 @@ class Persistence_heat_maps {
 
   /**
    * Return minimal range value of persistent image.
-  **/
+   **/
   inline double get_min() const { return this->min_; }
 
   /**
    * Return maximal range value of persistent image.
-  **/
+   **/
   inline double get_max() const { return this->max_; }
 
   /**
    * Operator == to check if to persistence heat maps are the same.
-  **/
+   **/
   bool operator==(const Persistence_heat_maps& rhs) const {
     bool dbg = false;
     if (!this->check_if_the_same(rhs)) {
@@ -361,12 +353,12 @@ class Persistence_heat_maps {
 
   /**
    * Operator != to check if to persistence heat maps are different.
-  **/
+   **/
   bool operator!=(const Persistence_heat_maps& rhs) const { return !((*this) == rhs); }
 
   /**
    * A function to generate a gnuplot script to visualize the persistent image.
-  **/
+   **/
   void plot(const char* filename) const;
 
   template <typename Operation_type>
@@ -396,7 +388,7 @@ class Persistence_heat_maps {
   /**
    * Multiplication of Persistence_heat_maps by scalar (so that all values of the heat map gets multiplied by that
    *scalar).
-  **/
+   **/
   Persistence_heat_maps multiply_by_scalar(double scalar) const {
     Persistence_heat_maps result;
     result.min_ = this->min_;
@@ -415,56 +407,56 @@ class Persistence_heat_maps {
 
   /**
    * This function computes a sum of two objects of a type Persistence_heat_maps.
-  **/
+   **/
   friend Persistence_heat_maps operator+(const Persistence_heat_maps& first, const Persistence_heat_maps& second) {
     return operation_on_pair_of_heat_maps(first, second, std::plus<double>());
   }
   /**
-* This function computes a difference of two objects of a type Persistence_heat_maps.
-**/
+   * This function computes a difference of two objects of a type Persistence_heat_maps.
+   **/
   friend Persistence_heat_maps operator-(const Persistence_heat_maps& first, const Persistence_heat_maps& second) {
     return operation_on_pair_of_heat_maps(first, second, std::minus<double>());
   }
   /**
-* This function computes a product of an object of a type Persistence_heat_maps with real number.
-**/
+   * This function computes a product of an object of a type Persistence_heat_maps with real number.
+   **/
   friend Persistence_heat_maps operator*(double scalar, const Persistence_heat_maps& A) {
     return A.multiply_by_scalar(scalar);
   }
   /**
-* This function computes a product of an object of a type Persistence_heat_maps with real number.
-**/
+   * This function computes a product of an object of a type Persistence_heat_maps with real number.
+   **/
   friend Persistence_heat_maps operator*(const Persistence_heat_maps& A, double scalar) {
     return A.multiply_by_scalar(scalar);
   }
   /**
-* This function computes a product of an object of a type Persistence_heat_maps with real number.
-**/
+   * This function computes a product of an object of a type Persistence_heat_maps with real number.
+   **/
   Persistence_heat_maps operator*(double scalar) { return this->multiply_by_scalar(scalar); }
   /**
    * += operator for Persistence_heat_maps.
-  **/
+   **/
   Persistence_heat_maps operator+=(const Persistence_heat_maps& rhs) {
     *this = *this + rhs;
     return *this;
   }
   /**
-       * -= operator for Persistence_heat_maps.
-      **/
+   * -= operator for Persistence_heat_maps.
+   **/
   Persistence_heat_maps operator-=(const Persistence_heat_maps& rhs) {
     *this = *this - rhs;
     return *this;
   }
   /**
-       * *= operator for Persistence_heat_maps.
-      **/
+   * *= operator for Persistence_heat_maps.
+   **/
   Persistence_heat_maps operator*=(double x) {
     *this = *this * x;
     return *this;
   }
   /**
-       * /= operator for Persistence_heat_maps.
-      **/
+   * /= operator for Persistence_heat_maps.
+   **/
   Persistence_heat_maps operator/=(double x) {
     if (x == 0) throw("In operator /=, division by 0. Program terminated.");
     *this = *this * (1 / x);
@@ -474,14 +466,14 @@ class Persistence_heat_maps {
   // Implementations of functions for various concepts.
 
   /**
-       * This function produce a vector of doubles based on a persistence heat map. It is required in a concept
+   * This function produce a vector of doubles based on a persistence heat map. It is required in a concept
    * Vectorized_topological_data
-      */
+   */
   std::vector<double> vectorize(int number_of_function) const;
   /**
-       * This function return the number of functions that allows vectorization of persistence heat map. It is required
-    *in a concept Vectorized_topological_data.
-       **/
+   * This function return the number of functions that allows vectorization of persistence heat map. It is required
+   *in a concept Vectorized_topological_data.
+   **/
   size_t number_of_vectorize_functions() const { return this->number_of_functions_for_vectorization; }
 
   /**
@@ -490,45 +482,45 @@ class Persistence_heat_maps {
    * At the moment this function is not tested, since it is quite likely to be changed in the future. Given this, when
    *using it, keep in mind that it
    * will be most likely changed in the next versions.
-  **/
+   **/
   double project_to_R(int number_of_function) const;
   /**
    * The function gives the number of possible projections to R. This function is required by the
    *Real_valued_topological_data concept.
-  **/
+   **/
   size_t number_of_projections_to_R() const { return this->number_of_functions_for_projections_to_reals; }
 
   /**
-  * A function to compute distance between persistence heat maps.
-  * The parameter of this function is a const reference to an object of a class Persistence_heat_maps.
-  * This function is required in Topological_data_with_distances concept.
-* For max norm distance, set power to std::numeric_limits<double>::max()
-  **/
+   * A function to compute distance between persistence heat maps.
+   * The parameter of this function is a const reference to an object of a class Persistence_heat_maps.
+   * This function is required in Topological_data_with_distances concept.
+   * For max norm distance, set power to std::numeric_limits<double>::max()
+   **/
   double distance(const Persistence_heat_maps& second_, double power = 1) const;
 
   /**
    * A function to compute averaged persistence heat map, based on vector of persistence heat maps.
    * This function is required by Topological_data_with_averages concept.
-  **/
+   **/
   void compute_average(const std::vector<Persistence_heat_maps*>& to_average);
 
   /**
-  * A function to compute scalar product of persistence heat maps.
-  * The parameter of this function is a const reference to an object of a class Persistence_heat_maps.
-  * This function is required in Topological_data_with_scalar_product concept.
-  **/
+   * A function to compute scalar product of persistence heat maps.
+   * The parameter of this function is a const reference to an object of a class Persistence_heat_maps.
+   * This function is required in Topological_data_with_scalar_product concept.
+   **/
   double compute_scalar_product(const Persistence_heat_maps& second_) const;
 
   // end of implementation of functions needed for concepts.
 
   /**
    * The x-range of the persistence heat map.
-  **/
+   **/
   std::pair<double, double> get_x_range() const { return std::make_pair(this->min_, this->max_); }
 
   /**
    * The y-range of the persistence heat map.
-  **/
+   **/
   std::pair<double, double> get_y_range() const { return this->get_x_range(); }
 
  protected:
@@ -546,12 +538,11 @@ class Persistence_heat_maps {
     this->number_of_functions_for_projections_to_reals = 1;
   }
 
-// Mathieu's code ***********************************************************************************************************************
-  bool discrete = true; // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false)
-  std::function<double(std::pair<double, double>, std::pair<double,double>)> kernel;
+  // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false)
+  bool discrete = true;
+  std::function<double(std::pair<double, double>, std::pair<double, double>)> kernel;
   std::vector<std::pair<double, double> > interval;
   std::vector<double> weights;
-// **************************************************************************************************************************************
 
   // data
   Scalling_of_kernels f;
@@ -561,27 +552,27 @@ class Persistence_heat_maps {
   std::vector<std::vector<double> > heat_map;
 };
 
-
-
-// Mathieu's code ***********************************************************************************************************************
 template <typename Scalling_of_kernels>
-Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
-                                                                  const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel,
-                                                                  size_t number_of_x_pixels, size_t number_of_y_pixels,
-                                                                  double min_x, double max_x,
-                                                                  double min_y, double max_y) {
-
-  this->discrete = true; this->min_ = min_x; this->max_ = max_x;
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(
+    const std::vector<std::pair<double, double> >& interval,
+    const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel,
+    size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x, double max_x, double min_y, double max_y) {
+  this->discrete = true;
+  this->min_ = min_x;
+  this->max_ = max_x;
   this->heat_map.resize(number_of_y_pixels);
-  double step_x = (max_x - min_x)/(number_of_x_pixels - 1); double step_y = (max_y - min_y)/(number_of_y_pixels - 1);
+  double step_x = (max_x - min_x) / (number_of_x_pixels - 1);
+  double step_y = (max_y - min_y) / (number_of_y_pixels - 1);
 
   int num_pts = interval.size();
-  for(size_t i = 0; i < number_of_y_pixels; i++){
-    double y = min_y + i*step_y; this->heat_map[i].reserve(number_of_x_pixels);
-    for(size_t j = 0; j < number_of_x_pixels; j++){
-      double x = min_x + j*step_x;
-      std::pair<double, double> grid_point(x,y); double pixel_value = 0;
-      for(int k = 0; k < num_pts; k++)  pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point);
+  for (size_t i = 0; i < number_of_y_pixels; i++) {
+    double y = min_y + i * step_y;
+    this->heat_map[i].reserve(number_of_x_pixels);
+    for (size_t j = 0; j < number_of_x_pixels; j++) {
+      double x = min_x + j * step_x;
+      std::pair<double, double> grid_point(x, y);
+      double pixel_value = 0;
+      for (int k = 0; k < num_pts; k++) pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point);
       this->heat_map[i].push_back(pixel_value);
     }
   }
@@ -589,16 +580,16 @@ Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(const std::vec
 }
 
 template <typename Scalling_of_kernels>
-Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
-                                                                  const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel) {
-  this->discrete = false; this->interval = interval; this->kernel = kernel;
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(
+    const std::vector<std::pair<double, double> >& interval,
+    const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel) {
+  this->discrete = false;
+  this->interval = interval;
+  this->kernel = kernel;
   int num_pts = this->interval.size();
-  for (int i = 0; i < num_pts; i++)  this->weights.push_back(this->f(this->interval[i]));
+  for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->interval[i]));
   this->set_up_parameters_for_basic_classes();
 }
-// **************************************************************************************************************************************
-
-
 
 // if min_ == max_, then the program is requested to set up the values itself based on persistence intervals
 template <typename Scalling_of_kernels>
@@ -897,9 +888,11 @@ void Persistence_heat_maps<Scalling_of_kernels>::load_from_file(const char* file
 // Concretizations of virtual methods:
 template <typename Scalling_of_kernels>
 std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize(int number_of_function) const {
-
   std::vector<double> result;
-  if(!discrete){std::cout << "No vectorize method in case of infinite dimensional vectorization" << std::endl; return result;}
+  if (!discrete) {
+    std::cout << "No vectorize method in case of infinite dimensional vectorization" << std::endl;
+    return result;
+  }
 
   // convert this->heat_map into one large vector:
   size_t size_of_result = 0;
@@ -920,7 +913,7 @@ std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize(int nu
 
 template <typename Scalling_of_kernels>
 double Persistence_heat_maps<Scalling_of_kernels>::distance(const Persistence_heat_maps& second, double power) const {
-  if(this->discrete){
+  if (this->discrete) {
     // first we need to check if (*this) and second are defined on the same domain and have the same dimensions:
     if (!this->check_if_the_same(second)) {
       std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between "
@@ -928,30 +921,30 @@ double Persistence_heat_maps<Scalling_of_kernels>::distance(const Persistence_he
       throw "The persistence images are of non compatible sizes. The program will now terminate";
     }
 
-    // if we are here, we know that the two persistence images are defined on the same domain, so we can start computing their distances:
+    // if we are here, we know that the two persistence images are defined on the same domain, so we can start
+    // computing their distances:
 
     double distance = 0;
     if (power < std::numeric_limits<double>::max()) {
       for (size_t i = 0; i != this->heat_map.size(); ++i) {
         for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
-          distance += pow(fabs(this->heat_map[i][j] - second.heat_map[i][j]), power);
+          distance += std::pow(std::fabs(this->heat_map[i][j] - second.heat_map[i][j]), power);
         }
       }
     } else {
       // in this case, we compute max norm distance
       for (size_t i = 0; i != this->heat_map.size(); ++i) {
         for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
-          if (distance < fabs(this->heat_map[i][j] - second.heat_map[i][j])) {
-            distance = fabs(this->heat_map[i][j] - second.heat_map[i][j]);
+          if (distance < std::fabs(this->heat_map[i][j] - second.heat_map[i][j])) {
+            distance = std::fabs(this->heat_map[i][j] - second.heat_map[i][j]);
           }
         }
       }
     }
     return distance;
   } else {
-
-    return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) -2 * this->compute_scalar_product(second));
-
+    return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) -
+                     2 * this->compute_scalar_product(second));
   }
 }
 
@@ -974,8 +967,7 @@ void Persistence_heat_maps<Scalling_of_kernels>::compute_average(
 
 template <typename Scalling_of_kernels>
 double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product(const Persistence_heat_maps& second) const {
-
-  if(discrete){
+  if (discrete) {
     // first we need to check if (*this) and second are defined on the same domain and have the same dimensions:
     if (!this->check_if_the_same(second)) {
       std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between "
@@ -992,22 +984,19 @@ double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product(const
       }
     }
     return scalar_prod;
-  }
-
-// Mathieu's code ***********************************************************************************************************************
-  else{
-    int num_pts1 = this->interval.size(); int num_pts2 = second.interval.size(); double kernel_val = 0;
-    for(int i = 0; i < num_pts1; i++){
+  } else {
+    int num_pts1 = this->interval.size();
+    int num_pts2 = second.interval.size();
+    double kernel_val = 0;
+    for (int i = 0; i < num_pts1; i++) {
       std::pair<double, double> pi = this->interval[i];
-      for(int j = 0; j < num_pts2; j++){
+      for (int j = 0; j < num_pts2; j++) {
         std::pair<double, double> pj = second.interval[j];
         kernel_val += this->weights[i] * second.weights[j] * this->kernel(pi, pj);
       }
     }
     return kernel_val;
   }
-// **************************************************************************************************************************************
-
 }
 
 }  // namespace Persistence_representations
diff --git a/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid.h b/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid.h
index db0e362a..fd8a181c 100644
--- a/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid.h
+++ b/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid.h
@@ -986,7 +986,7 @@ void Persistence_landscape_on_grid::set_up_values_of_landscapes(const std::vecto
   for (size_t int_no = 0; int_no != p.size(); ++int_no) {
     size_t grid_interval_begin = (p[int_no].first - grid_min_) / dx;
     size_t grid_interval_end = (p[int_no].second - grid_min_) / dx;
-    size_t grid_interval_midpoint = (size_t)(0.5 * (p[int_no].first + p[int_no].second) - grid_min + 1);
+    size_t grid_interval_midpoint = (size_t)(0.5 * (grid_interval_begin + grid_interval_end));
 
     if (dbg) {
       std::cerr << "Considering an interval : " << p[int_no].first << "," << p[int_no].second << std::endl;
@@ -996,7 +996,7 @@ void Persistence_landscape_on_grid::set_up_values_of_landscapes(const std::vecto
       std::cerr << "grid_interval_midpoint : " << grid_interval_midpoint << std::endl;
     }
 
-    double landscape_value = grid_min + dx * (grid_interval_begin + 1) - p[int_no].first;
+    double landscape_value = dx;
     for (size_t i = grid_interval_begin + 1; i < grid_interval_midpoint; ++i) {
       if (dbg) {
         std::cerr << "Adding landscape value (going up) for a point : " << i << " equal : " << landscape_value
@@ -1030,8 +1030,6 @@ void Persistence_landscape_on_grid::set_up_values_of_landscapes(const std::vecto
       }
       landscape_value += dx;
     }
-
-    landscape_value = p[int_no].second - grid_min - dx * grid_interval_midpoint;
     for (size_t i = grid_interval_midpoint; i <= grid_interval_end; ++i) {
       if (landscape_value > 0) {
         if (number_of_levels != std::numeric_limits<unsigned>::max()) {
diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
index 5d0f4a5d..18165c5f 100644
--- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
+++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Mathieu Carriere
  *
- *    Copyright (C) 2018  INRIA (France)
+ *    Copyright (C) 2018  Inria
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -28,6 +28,13 @@
 #include <gudhi/common_persistence_representations.h>
 #include <gudhi/Debug_utils.h>
 
+#include <vector>     // for std::vector<>
+#include <utility>    // for std::pair<>, std::move
+#include <algorithm>  // for std::sort, std::max, std::merge
+#include <cmath>      // for std::abs, std::sqrt
+#include <stdexcept>  // for std::invalid_argument
+#include <random>     // for std::random_device
+
 namespace Gudhi {
 namespace Persistence_representations {
 
@@ -39,31 +46,33 @@ namespace Persistence_representations {
  *
  * \details
  * In this class, we compute infinite-dimensional representations of persistence diagrams by using the
- * Sliced Wasserstein kernel (see \ref sec_persistence_kernels for more details on kernels). We recall that infinite-dimensional
- * representations are defined implicitly, so only scalar products and distances are available for the representations defined in this class.
- * The Sliced Wasserstein kernel is defined as a Gaussian-like kernel between persistence diagrams, where the distance used for
- * comparison is the Sliced Wasserstein distance \f$SW\f$ between persistence diagrams, defined as the integral of the 1-norm
- * between the sorted projections of the diagrams onto all lines passing through the origin:
+ * Sliced Wasserstein kernel (see \ref sec_persistence_kernels for more details on kernels). We recall that
+ * infinite-dimensional representations are defined implicitly, so only scalar products and distances are available for
+ * the representations defined in this class.
+ * The Sliced Wasserstein kernel is defined as a Gaussian-like kernel between persistence diagrams, where the distance
+ * used for comparison is the Sliced Wasserstein distance \f$SW\f$ between persistence diagrams, defined as the
+ * integral of the 1-norm between the sorted projections of the diagrams onto all lines passing through the origin:
  *
- * \f$ SW(D_1,D_2)=\int_{\theta\in\mathbb{S}}\,\|\pi_\theta(D_1\cup\pi_\Delta(D_2))-\pi_\theta(D_2\cup\pi_\Delta(D_1))\|_1{\rm d}\theta\f$,
+ * \f$ SW(D_1,D_2)=\int_{\theta\in\mathbb{S}}\,\|\pi_\theta(D_1\cup\pi_\Delta(D_2))-\pi_\theta(D_2\cup\pi_\Delta(D_1))\
+ * |_1{\rm d}\theta\f$,
  *
- * where \f$\pi_\theta\f$ is the projection onto the line defined with angle \f$\theta\f$ in the unit circle \f$\mathbb{S}\f$,
- * and \f$\pi_\Delta\f$ is the projection onto the diagonal.
- * Assuming that the diagrams are in general position (i.e. there is no collinear triple), the integral can be computed exactly in \f$O(n^2{\rm log}(n))\f$ time, where \f$n\f$ is the number of points
- * in the diagrams. We provide two approximations of the integral: one in which we slightly perturb the diagram points so that they are in general position,
- * and another in which we approximate the integral by sampling \f$N\f$ lines in the circle in \f$O(Nn{\rm log}(n))\f$ time. The Sliced Wasserstein Kernel is then computed as:
+ * where \f$\pi_\theta\f$ is the projection onto the line defined with angle \f$\theta\f$ in the unit circle
+ * \f$\mathbb{S}\f$, and \f$\pi_\Delta\f$ is the projection onto the diagonal.
+ * Assuming that the diagrams are in general position (i.e. there is no collinear triple), the integral can be computed
+ * exactly in \f$O(n^2{\rm log}(n))\f$ time, where \f$n\f$ is the number of points in the diagrams. We provide two
+ * approximations of the integral: one in which we slightly perturb the diagram points so that they are in general
+ * position, and another in which we approximate the integral by sampling \f$N\f$ lines in the circle in
+ * \f$O(Nn{\rm log}(n))\f$ time. The Sliced Wasserstein Kernel is then computed as:
  *
  * \f$ k(D_1,D_2) = {\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right).\f$
  *
  * The first method is usually much more accurate but also
  * much slower. For more details, please see \cite pmlr-v70-carriere17a .
  *
-**/
+ **/
 
 class Sliced_Wasserstein {
-
  protected:
-
   Persistence_diagram diagram;
   int approx;
   double sigma;
@@ -73,27 +82,26 @@ class Sliced_Wasserstein {
   // Utils.
   // **********************************
 
-  void build_rep(){
-
-    if(approx > 0){
-
-      double step = pi/this->approx;
+  void build_rep() {
+    if (approx > 0) {
+      double step = pi / this->approx;
       int n = diagram.size();
 
-      for (int i = 0; i < this->approx; i++){
-        std::vector<double> l,l_diag;
-        for (int j = 0; j < n; j++){
-
-          double px = diagram[j].first; double py = diagram[j].second;
-          double proj_diag = (px+py)/2;
+      for (int i = 0; i < this->approx; i++) {
+        std::vector<double> l, l_diag;
+        for (int j = 0; j < n; j++) {
+          double px = diagram[j].first;
+          double py = diagram[j].second;
+          double proj_diag = (px + py) / 2;
 
-          l.push_back         (   px         * cos(-pi/2+i*step) + py         * sin(-pi/2+i*step)   );
-          l_diag.push_back    (   proj_diag  * cos(-pi/2+i*step) + proj_diag  * sin(-pi/2+i*step)   );
+          l.push_back(px * cos(-pi / 2 + i * step) + py * sin(-pi / 2 + i * step));
+          l_diag.push_back(proj_diag * cos(-pi / 2 + i * step) + proj_diag * sin(-pi / 2 + i * step));
         }
 
-        std::sort(l.begin(), l.end()); std::sort(l_diag.begin(), l_diag.end());
-        projections.push_back(std::move(l)); projections_diagonal.push_back(std::move(l_diag));
-
+        std::sort(l.begin(), l.end());
+        std::sort(l_diag.begin(), l_diag.end());
+        projections.push_back(std::move(l));
+        projections_diagonal.push_back(std::move(l_diag));
       }
 
       diagram.clear();
@@ -101,179 +109,254 @@ class Sliced_Wasserstein {
   }
 
   // Compute the angle formed by two points of a PD
-  double compute_angle(const Persistence_diagram & diag, int i, int j) const {
-    if(diag[i].second == diag[j].second) return pi/2; else return atan((diag[j].first-diag[i].first)/(diag[i].second-diag[j].second));
+  double compute_angle(const Persistence_diagram& diag, int i, int j) const {
+    if (diag[i].second == diag[j].second)
+      return pi / 2;
+    else
+      return atan((diag[j].first - diag[i].first) / (diag[i].second - diag[j].second));
   }
 
-  // Compute the integral of |cos()| between alpha and beta, valid only if alpha is in [-pi,pi] and beta-alpha is in [0,pi]
+  // Compute the integral of |cos()| between alpha and beta, valid only if alpha is in [-pi,pi] and beta-alpha is in
+  // [0,pi]
   double compute_int_cos(double alpha, double beta) const {
     double res = 0;
-    if (alpha >= 0 && alpha <= pi){
-      if (cos(alpha) >= 0){
-        if(pi/2 <= beta){res = 2-sin(alpha)-sin(beta);}
-        else{res = sin(beta)-sin(alpha);}
-      }
-      else{
-        if(1.5*pi <= beta){res = 2+sin(alpha)+sin(beta);}
-        else{res = sin(alpha)-sin(beta);}
+    if (alpha >= 0 && alpha <= pi) {
+      if (cos(alpha) >= 0) {
+        if (pi / 2 <= beta) {
+          res = 2 - sin(alpha) - sin(beta);
+        } else {
+          res = sin(beta) - sin(alpha);
+        }
+      } else {
+        if (1.5 * pi <= beta) {
+          res = 2 + sin(alpha) + sin(beta);
+        } else {
+          res = sin(alpha) - sin(beta);
+        }
       }
     }
-    if (alpha >= -pi && alpha <= 0){
-      if (cos(alpha) <= 0){
-        if(-pi/2 <= beta){res = 2+sin(alpha)+sin(beta);}
-        else{res = sin(alpha)-sin(beta);}
-      }
-      else{
-        if(pi/2 <= beta){res = 2-sin(alpha)-sin(beta);}
-        else{res = sin(beta)-sin(alpha);}
+    if (alpha >= -pi && alpha <= 0) {
+      if (cos(alpha) <= 0) {
+        if (-pi / 2 <= beta) {
+          res = 2 + sin(alpha) + sin(beta);
+        } else {
+          res = sin(alpha) - sin(beta);
+        }
+      } else {
+        if (pi / 2 <= beta) {
+          res = 2 - sin(alpha) - sin(beta);
+        } else {
+          res = sin(beta) - sin(alpha);
+        }
       }
     }
     return res;
   }
 
-  double compute_int(double theta1, double theta2, int p, int q, const Persistence_diagram & diag1, const Persistence_diagram & diag2) const {
-    double norm = std::sqrt(  (diag1[p].first-diag2[q].first)*(diag1[p].first-diag2[q].first) + (diag1[p].second-diag2[q].second)*(diag1[p].second-diag2[q].second)  );
+  double compute_int(double theta1, double theta2, int p, int q, const Persistence_diagram& diag1,
+                     const Persistence_diagram& diag2) const {
+    double norm = std::sqrt((diag1[p].first - diag2[q].first) * (diag1[p].first - diag2[q].first) +
+                            (diag1[p].second - diag2[q].second) * (diag1[p].second - diag2[q].second));
     double angle1;
-    if (diag1[p].first == diag2[q].first) angle1 = theta1 - pi/2; else angle1 = theta1 - atan((diag1[p].second-diag2[q].second)/(diag1[p].first-diag2[q].first));
+    if (diag1[p].first == diag2[q].first)
+      angle1 = theta1 - pi / 2;
+    else
+      angle1 = theta1 - atan((diag1[p].second - diag2[q].second) / (diag1[p].first - diag2[q].first));
     double angle2 = angle1 + theta2 - theta1;
-    double integral = compute_int_cos(angle1,angle2);
-    return norm*integral;
+    double integral = compute_int_cos(angle1, angle2);
+    return norm * integral;
   }
 
   // Evaluation of the Sliced Wasserstein Distance between a pair of diagrams.
-  double compute_sliced_wasserstein_distance(const Sliced_Wasserstein & second) const {
-
-    GUDHI_CHECK(this->approx == second.approx, std::invalid_argument("Error: different approx values for representations"));
+  double compute_sliced_wasserstein_distance(const Sliced_Wasserstein& second) const {
+    GUDHI_CHECK(this->approx == second.approx,
+                std::invalid_argument("Error: different approx values for representations"));
 
-    Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram; double sw = 0;
-
-    if(this->approx == -1){
+    Persistence_diagram diagram1 = this->diagram;
+    Persistence_diagram diagram2 = second.diagram;
+    double sw = 0;
 
+    if (this->approx == -1) {
       // Add projections onto diagonal.
-      int n1, n2; n1 = diagram1.size(); n2 = diagram2.size();
-      double min_ordinate = std::numeric_limits<double>::max();    double min_abscissa = std::numeric_limits<double>::max();
-      double max_ordinate = std::numeric_limits<double>::lowest(); double max_abscissa = std::numeric_limits<double>::lowest();
-      for (int i = 0; i < n2; i++){
-        min_ordinate = std::min(min_ordinate, diagram2[i].second); min_abscissa = std::min(min_abscissa, diagram2[i].first);
-        max_ordinate = std::max(max_ordinate, diagram2[i].second); max_abscissa = std::max(max_abscissa, diagram2[i].first);
-        diagram1.emplace_back(  (diagram2[i].first+diagram2[i].second)/2, (diagram2[i].first+diagram2[i].second)/2  );
+      int n1, n2;
+      n1 = diagram1.size();
+      n2 = diagram2.size();
+      double min_ordinate = std::numeric_limits<double>::max();
+      double min_abscissa = std::numeric_limits<double>::max();
+      double max_ordinate = std::numeric_limits<double>::lowest();
+      double max_abscissa = std::numeric_limits<double>::lowest();
+      for (int i = 0; i < n2; i++) {
+        min_ordinate = std::min(min_ordinate, diagram2[i].second);
+        min_abscissa = std::min(min_abscissa, diagram2[i].first);
+        max_ordinate = std::max(max_ordinate, diagram2[i].second);
+        max_abscissa = std::max(max_abscissa, diagram2[i].first);
+        diagram1.emplace_back((diagram2[i].first + diagram2[i].second) / 2,
+                              (diagram2[i].first + diagram2[i].second) / 2);
       }
-      for (int i = 0; i < n1; i++){
-        min_ordinate = std::min(min_ordinate, diagram1[i].second); min_abscissa = std::min(min_abscissa, diagram1[i].first);
-        max_ordinate = std::max(max_ordinate, diagram1[i].second); max_abscissa = std::max(max_abscissa, diagram1[i].first);
-        diagram2.emplace_back(  (diagram1[i].first+diagram1[i].second)/2, (diagram1[i].first+diagram1[i].second)/2  );
+      for (int i = 0; i < n1; i++) {
+        min_ordinate = std::min(min_ordinate, diagram1[i].second);
+        min_abscissa = std::min(min_abscissa, diagram1[i].first);
+        max_ordinate = std::max(max_ordinate, diagram1[i].second);
+        max_abscissa = std::max(max_abscissa, diagram1[i].first);
+        diagram2.emplace_back((diagram1[i].first + diagram1[i].second) / 2,
+                              (diagram1[i].first + diagram1[i].second) / 2);
       }
       int num_pts_dgm = diagram1.size();
 
       // Slightly perturb the points so that the PDs are in generic positions.
       double epsilon = 0.0001;
-      double thresh_y = (max_ordinate-min_ordinate) * epsilon; double thresh_x = (max_abscissa-min_abscissa) * epsilon;
-      std::random_device rd; std::default_random_engine re(rd()); std::uniform_real_distribution<double> uni(-1,1);
-      for (int i = 0; i < num_pts_dgm; i++){
+      double thresh_y = (max_ordinate - min_ordinate) * epsilon;
+      double thresh_x = (max_abscissa - min_abscissa) * epsilon;
+      std::random_device rd;
+      std::default_random_engine re(rd());
+      std::uniform_real_distribution<double> uni(-1, 1);
+      for (int i = 0; i < num_pts_dgm; i++) {
         double u = uni(re);
-        diagram1[i].first += u*thresh_x; diagram1[i].second += u*thresh_y;
-        diagram2[i].first += u*thresh_x; diagram2[i].second += u*thresh_y;
+        diagram1[i].first += u * thresh_x;
+        diagram1[i].second += u * thresh_y;
+        diagram2[i].first += u * thresh_x;
+        diagram2[i].second += u * thresh_y;
       }
 
       // Compute all angles in both PDs.
-      std::vector<std::pair<double, std::pair<int,int> > > angles1, angles2;
-      for (int i = 0; i < num_pts_dgm; i++){
-        for (int j = i+1; j < num_pts_dgm; j++){
-          double theta1 = compute_angle(diagram1,i,j); double theta2 = compute_angle(diagram2,i,j);
-          angles1.emplace_back(theta1, std::pair<int,int>(i,j));
-          angles2.emplace_back(theta2, std::pair<int,int>(i,j));
+      std::vector<std::pair<double, std::pair<int, int> > > angles1, angles2;
+      for (int i = 0; i < num_pts_dgm; i++) {
+        for (int j = i + 1; j < num_pts_dgm; j++) {
+          double theta1 = compute_angle(diagram1, i, j);
+          double theta2 = compute_angle(diagram2, i, j);
+          angles1.emplace_back(theta1, std::pair<int, int>(i, j));
+          angles2.emplace_back(theta2, std::pair<int, int>(i, j));
         }
       }
 
       // Sort angles.
-      std::sort(angles1.begin(), angles1.end(), [](const std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);});
-      std::sort(angles2.begin(), angles2.end(), [](const std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);});
+      std::sort(angles1.begin(), angles1.end(),
+                [](const std::pair<double, std::pair<int, int> >& p1,
+                   const std::pair<double, std::pair<int, int> >& p2) { return (p1.first < p2.first); });
+      std::sort(angles2.begin(), angles2.end(),
+                [](const std::pair<double, std::pair<int, int> >& p1,
+                   const std::pair<double, std::pair<int, int> >& p2) { return (p1.first < p2.first); });
 
       // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2).
       std::vector<int> orderp1, orderp2;
-      for (int i = 0; i < num_pts_dgm; i++){ orderp1.push_back(i); orderp2.push_back(i); }
-      std::sort( orderp1.begin(), orderp1.end(), [&](int i, int j){ if(diagram1[i].second != diagram1[j].second) return (diagram1[i].second < diagram1[j].second); else return (diagram1[i].first > diagram1[j].first); } );
-      std::sort( orderp2.begin(), orderp2.end(), [&](int i, int j){ if(diagram2[i].second != diagram2[j].second) return (diagram2[i].second < diagram2[j].second); else return (diagram2[i].first > diagram2[j].first); } );
+      for (int i = 0; i < num_pts_dgm; i++) {
+        orderp1.push_back(i);
+        orderp2.push_back(i);
+      }
+      std::sort(orderp1.begin(), orderp1.end(), [&](int i, int j) {
+        if (diagram1[i].second != diagram1[j].second)
+          return (diagram1[i].second < diagram1[j].second);
+        else
+          return (diagram1[i].first > diagram1[j].first);
+      });
+      std::sort(orderp2.begin(), orderp2.end(), [&](int i, int j) {
+        if (diagram2[i].second != diagram2[j].second)
+          return (diagram2[i].second < diagram2[j].second);
+        else
+          return (diagram2[i].first > diagram2[j].first);
+      });
 
       // Find the inverses of the orders.
-      std::vector<int> order1(num_pts_dgm); std::vector<int> order2(num_pts_dgm);
-      for(int i = 0; i < num_pts_dgm; i++){ order1[orderp1[i]] = i; order2[orderp2[i]] = i; }
+      std::vector<int> order1(num_pts_dgm);
+      std::vector<int> order2(num_pts_dgm);
+      for (int i = 0; i < num_pts_dgm; i++) {
+        order1[orderp1[i]] = i;
+        order2[orderp2[i]] = i;
+      }
 
       // Record all inversions of points in the orders as theta varies along the positive half-disk.
-      std::vector<std::vector<std::pair<int,double> > > anglePerm1(num_pts_dgm);
-      std::vector<std::vector<std::pair<int,double> > > anglePerm2(num_pts_dgm);
+      std::vector<std::vector<std::pair<int, double> > > anglePerm1(num_pts_dgm);
+      std::vector<std::vector<std::pair<int, double> > > anglePerm2(num_pts_dgm);
 
       int m1 = angles1.size();
-      for (int i = 0; i < m1; i++){
-        double theta = angles1[i].first; int p = angles1[i].second.first; int q = angles1[i].second.second;
-        anglePerm1[order1[p]].emplace_back(p,theta);
-        anglePerm1[order1[q]].emplace_back(q,theta);
-        int a = order1[p]; int b = order1[q]; order1[p] = b; order1[q] = a;
+      for (int i = 0; i < m1; i++) {
+        double theta = angles1[i].first;
+        int p = angles1[i].second.first;
+        int q = angles1[i].second.second;
+        anglePerm1[order1[p]].emplace_back(p, theta);
+        anglePerm1[order1[q]].emplace_back(q, theta);
+        int a = order1[p];
+        int b = order1[q];
+        order1[p] = b;
+        order1[q] = a;
       }
 
       int m2 = angles2.size();
-      for (int i = 0; i < m2; i++){
-        double theta = angles2[i].first; int p = angles2[i].second.first; int q = angles2[i].second.second;
-        anglePerm2[order2[p]].emplace_back(p,theta);
-        anglePerm2[order2[q]].emplace_back(q,theta);
-        int a = order2[p]; int b = order2[q]; order2[p] = b; order2[q] = a;
+      for (int i = 0; i < m2; i++) {
+        double theta = angles2[i].first;
+        int p = angles2[i].second.first;
+        int q = angles2[i].second.second;
+        anglePerm2[order2[p]].emplace_back(p, theta);
+        anglePerm2[order2[q]].emplace_back(q, theta);
+        int a = order2[p];
+        int b = order2[q];
+        order2[p] = b;
+        order2[q] = a;
       }
 
-      for (int i = 0; i < num_pts_dgm; i++){
-        anglePerm1[order1[i]].emplace_back(i,pi/2);
-        anglePerm2[order2[i]].emplace_back(i,pi/2);
+      for (int i = 0; i < num_pts_dgm; i++) {
+        anglePerm1[order1[i]].emplace_back(i, pi / 2);
+        anglePerm2[order2[i]].emplace_back(i, pi / 2);
       }
 
       // Compute the SW distance with the list of inversions.
-      for (int i = 0; i < num_pts_dgm; i++){
-        std::vector<std::pair<int,double> > u,v; u = anglePerm1[i]; v = anglePerm2[i];
-        double theta1, theta2; theta1 = -pi/2;
-        unsigned int ku, kv; ku = 0; kv = 0; theta2 = std::min(u[ku].second,v[kv].second);
-        while(theta1 != pi/2){
-          if(diagram1[u[ku].first].first != diagram2[v[kv].first].first || diagram1[u[ku].first].second != diagram2[v[kv].first].second)
-            if(theta1 != theta2)
-              sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2);
+      for (int i = 0; i < num_pts_dgm; i++) {
+        std::vector<std::pair<int, double> > u, v;
+        u = anglePerm1[i];
+        v = anglePerm2[i];
+        double theta1, theta2;
+        theta1 = -pi / 2;
+        unsigned int ku, kv;
+        ku = 0;
+        kv = 0;
+        theta2 = std::min(u[ku].second, v[kv].second);
+        while (theta1 != pi / 2) {
+          if (diagram1[u[ku].first].first != diagram2[v[kv].first].first ||
+              diagram1[u[ku].first].second != diagram2[v[kv].first].second)
+            if (theta1 != theta2) sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2);
           theta1 = theta2;
-          if (  (theta2 == u[ku].second)  &&  ku < u.size()-1  )  ku++;
-          if (  (theta2 == v[kv].second)  &&  kv < v.size()-1  )  kv++;
+          if ((theta2 == u[ku].second) && ku < u.size() - 1) ku++;
+          if ((theta2 == v[kv].second) && kv < v.size() - 1) kv++;
           theta2 = std::min(u[ku].second, v[kv].second);
         }
       }
-    }
-
-
-    else{
-
-      double step = pi/this->approx; std::vector<double> v1, v2;
-      for (int i = 0; i < this->approx; i++){
-
-        v1.clear(); v2.clear();
-        std::merge(this->projections[i].begin(), this->projections[i].end(), second.projections_diagonal[i].begin(), second.projections_diagonal[i].end(), std::back_inserter(v1));
-        std::merge(second.projections[i].begin(), second.projections[i].end(), this->projections_diagonal[i].begin(), this->projections_diagonal[i].end(), std::back_inserter(v2));
-
-        int n = v1.size(); double f = 0;
-        for (int j = 0; j < n; j++)  f += std::abs(v1[j] - v2[j]);
-        sw += f*step;
-
+    } else {
+      double step = pi / this->approx;
+      std::vector<double> v1, v2;
+      for (int i = 0; i < this->approx; i++) {
+        v1.clear();
+        v2.clear();
+        std::merge(this->projections[i].begin(), this->projections[i].end(), second.projections_diagonal[i].begin(),
+                   second.projections_diagonal[i].end(), std::back_inserter(v1));
+        std::merge(second.projections[i].begin(), second.projections[i].end(), this->projections_diagonal[i].begin(),
+                   this->projections_diagonal[i].end(), std::back_inserter(v2));
+
+        int n = v1.size();
+        double f = 0;
+        for (int j = 0; j < n; j++) f += std::abs(v1[j] - v2[j]);
+        sw += f * step;
       }
     }
 
-    return sw/pi;
+    return sw / pi;
   }
 
  public:
-
   /** \brief Sliced Wasserstein kernel constructor.
    * \implements Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_scalar_product
    * \ingroup Sliced_Wasserstein
    *
    * @param[in] _diagram  persistence diagram.
    * @param[in] _sigma    bandwidth parameter.
-   * @param[in] _approx   number of directions used to approximate the integral in the Sliced Wasserstein distance, set to -1 for random perturbation. If positive, then projections of the diagram
-   *                      points on all directions are stored in memory to reduce computation time.
+   * @param[in] _approx   number of directions used to approximate the integral in the Sliced Wasserstein distance, set
+   *                      to -1 for random perturbation. If positive, then projections of the diagram points on all
+   *                      directions are stored in memory to reduce computation time.
    *
    */
-  Sliced_Wasserstein(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 10):diagram(_diagram), approx(_approx), sigma(_sigma) {build_rep();}
+  Sliced_Wasserstein(const Persistence_diagram& _diagram, double _sigma = 1.0, int _approx = 10)
+      : diagram(_diagram), approx(_approx), sigma(_sigma) {
+    build_rep();
+  }
 
   /** \brief Evaluation of the kernel on a pair of diagrams.
    * \ingroup Sliced_Wasserstein
@@ -282,9 +365,10 @@ class Sliced_Wasserstein {
    * @param[in] second other instance of class Sliced_Wasserstein.
    *
    */
-  double compute_scalar_product(const Sliced_Wasserstein & second) const {
-    GUDHI_CHECK(this->sigma == second.sigma, std::invalid_argument("Error: different sigma values for representations"));
-    return std::exp(-compute_sliced_wasserstein_distance(second)/(2*this->sigma*this->sigma));
+  double compute_scalar_product(const Sliced_Wasserstein& second) const {
+    GUDHI_CHECK(this->sigma == second.sigma,
+                std::invalid_argument("Error: different sigma values for representations"));
+    return std::exp(-compute_sliced_wasserstein_distance(second) / (2 * this->sigma * this->sigma));
   }
 
   /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel.
@@ -294,13 +378,14 @@ class Sliced_Wasserstein {
    * @param[in] second  other instance of class Sliced_Wasserstein.
    *
    */
-  double distance(const Sliced_Wasserstein & second) const {
-    GUDHI_CHECK(this->sigma == second.sigma, std::invalid_argument("Error: different sigma values for representations"));
-    return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second),  0.5);
+  double distance(const Sliced_Wasserstein& second) const {
+    GUDHI_CHECK(this->sigma == second.sigma,
+                std::invalid_argument("Error: different sigma values for representations"));
+    return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) -
+                     2 * this->compute_scalar_product(second));
   }
 
-
-}; // class Sliced_Wasserstein
+};  // class Sliced_Wasserstein
 }  // namespace Persistence_representations
 }  // namespace Gudhi