git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3387 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 3fe2ae4af0c7cadf507fc5148c05dcf664c5e151

9 files changed, 242 insertions, 255 deletions

diff --git a/src/Persistence_representations/doc/Persistence_representations_doc.h b/src/Persistence_representations/doc/Persistence_representations_doc.h
index 6d4cc96c..ca283017 100644
--- a/src/Persistence_representations/doc/Persistence_representations_doc.h
+++ b/src/Persistence_representations/doc/Persistence_representations_doc.h
@@ -24,7 +24,6 @@
 #define DOC_GUDHI_STAT_H_
 
 namespace Gudhi {
-
 namespace Persistence_representations {
 
 /**  \defgroup Persistence_representations Persistence representations
@@ -254,11 +253,11 @@ namespace Persistence_representations {
 
 
 
-\section sec_persistence_kernels Kernels on Persistence Diagrams
+\section sec_persistence_kernels Kernels on persistence diagrams
  <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Sliced_Wasserstein <br>
  <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_weighted_gaussian <br>
 
- Kernels for Persistence Diagrams can be regarded as infinite-dimensional vectorizations. More specifically,
+ Kernels for persistence diagrams can be regarded as infinite-dimensional vectorizations. More specifically,
  they are similarity functions whose evaluations on pairs of persistence diagrams equals the scalar products
  between images of these pairs under a map \f$\Phi\f$ taking values in a specific (possibly non Euclidean) Hilbert space \f$k(D_i, D_j) = \langle \Phi(D_i),\Phi(D_j)\rangle\f$.
  Reciprocally, classical results of learning theory ensure that such a \f$\Phi\f$ exists for a given similarity function \f$k\f$ if and only if \f$k\f$ is <i>positive semi-definite</i>.
diff --git a/src/Persistence_representations/example/persistence_image.cpp b/src/Persistence_representations/example/persistence_image.cpp
index dfa469d4..cdce3bbf 100644
--- a/src/Persistence_representations/example/persistence_image.cpp
+++ b/src/Persistence_representations/example/persistence_image.cpp
@@ -40,7 +40,7 @@ int main(int argc, char** argv) {
   persistence.push_back(std::make_pair(0, 4));
   persistence.push_back(std::make_pair(3, 8));
 
-  double min_x = 0.0; double max_x = 10.0; int res_x = 100; double min_y = 0.0; double max_y = 10.0; int res_y = 100; double sigma = 1.0; Weight weight = Gudhi::Persistence_representations::Persistence_weighted_gaussian::linear_weight;
+  double min_x = 0.0; double max_x = 10.0; int res_x = 100; double min_y = 0.0; double max_y = 10.0; int res_y = 100; double sigma = 1.0; Weight weight = Gudhi::Persistence_representations::linear_weight;
 
   PI pim(persistence, min_x, max_x, res_x, min_y, max_y, res_y, weight, sigma);
   std::vector<std::vector<double> > P = pim.vectorize();
diff --git a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp
index 234f6323..db60755f 100644
--- a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp
+++ b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp
@@ -26,13 +26,11 @@
 #include <vector>
 #include <utility>
 
-using PD = std::vector<std::pair<double,double> >;
 using PWG = Gudhi::Persistence_representations::Persistence_weighted_gaussian;
 
 int main(int argc, char** argv) {
 
-  std::vector<std::pair<double, double> > persistence1;
-  std::vector<std::pair<double, double> > persistence2;
+  Persistence_diagram persistence1, persistence2;
 
   persistence1.push_back(std::make_pair(1, 2));
   persistence1.push_back(std::make_pair(6, 8));
@@ -48,11 +46,11 @@ int main(int argc, char** argv) {
   double tau = 1;
   int m = 10000;
 
-  PWG PWG1(persistence1, sigma, m, PWG::arctan_weight(1,1));
-  PWG PWG2(persistence2, sigma, m, PWG::arctan_weight(1,1));
+  PWG PWG1(persistence1, sigma, m, Gudhi::Persistence_representations::arctan_weight(1,1));
+  PWG PWG2(persistence2, sigma, m, Gudhi::Persistence_representations::arctan_weight(1,1));
 
-  PWG PWGex1(persistence1, sigma, -1, PWG::arctan_weight(1,1));
-  PWG PWGex2(persistence2, sigma, -1, PWG::arctan_weight(1,1));
+  PWG PWGex1(persistence1, sigma, -1, Gudhi::Persistence_representations::arctan_weight(1,1));
+  PWG PWGex2(persistence2, sigma, -1, Gudhi::Persistence_representations::arctan_weight(1,1));
 
 
   // Linear PWG
@@ -82,14 +80,14 @@ int main(int argc, char** argv) {
 
   // PSS
 
-  PD pd1 = persistence1; int numpts = persistence1.size();    for(int i = 0; i < numpts; i++)  pd1.emplace_back(persistence1[i].second,persistence1[i].first);
-  PD pd2 = persistence2;     numpts = persistence2.size();    for(int i = 0; i < numpts; i++)  pd2.emplace_back(persistence2[i].second,persistence2[i].first);
+  Persistence_diagram pd1 = persistence1; int numpts = persistence1.size();    for(int i = 0; i < numpts; i++)  pd1.emplace_back(persistence1[i].second,persistence1[i].first);
+  Persistence_diagram pd2 = persistence2;     numpts = persistence2.size();    for(int i = 0; i < numpts; i++)  pd2.emplace_back(persistence2[i].second,persistence2[i].first);
 
-  PWG pwg1(pd1, 2*std::sqrt(sigma), m, PWG::pss_weight);
-  PWG pwg2(pd2, 2*std::sqrt(sigma), m, PWG::pss_weight);
+  PWG pwg1(pd1, 2*std::sqrt(sigma), m, Gudhi::Persistence_representations::pss_weight);
+  PWG pwg2(pd2, 2*std::sqrt(sigma), m, Gudhi::Persistence_representations::pss_weight);
 
-  PWG pwgex1(pd1, 2*std::sqrt(sigma), -1, PWG::pss_weight);
-  PWG pwgex2(pd2, 2*std::sqrt(sigma), -1, PWG::pss_weight);
+  PWG pwgex1(pd1, 2*std::sqrt(sigma), -1, Gudhi::Persistence_representations::pss_weight);
+  PWG pwgex2(pd2, 2*std::sqrt(sigma), -1, Gudhi::Persistence_representations::pss_weight);
 
   std::cout << "Approx PSS kernel: " << pwg1.compute_scalar_product  (pwg2) / (16*Gudhi::Persistence_representations::pi*sigma) << std::endl;
   std::cout << "Exact  PSS kernel: " << pwgex1.compute_scalar_product  (pwgex2) / (16*Gudhi::Persistence_representations::pi*sigma) << std::endl;
diff --git a/src/Persistence_representations/example/sliced_wasserstein.cpp b/src/Persistence_representations/example/sliced_wasserstein.cpp
index f1aeea5c..d37cb23c 100644
--- a/src/Persistence_representations/example/sliced_wasserstein.cpp
+++ b/src/Persistence_representations/example/sliced_wasserstein.cpp
@@ -30,8 +30,7 @@ using SW = Gudhi::Persistence_representations::Sliced_Wasserstein;
 
 int main(int argc, char** argv) {
 
-  std::vector<std::pair<double, double> > persistence1;
-  std::vector<std::pair<double, double> > persistence2;
+  Persistence_diagram persistence1, persistence2;
 
   persistence1.push_back(std::make_pair(1, 2));
   persistence1.push_back(std::make_pair(6, 8));
diff --git a/src/Persistence_representations/include/gudhi/Landscape.h b/src/Persistence_representations/include/gudhi/Landscape.h
index d6608a57..bbbca36b 100644
--- a/src/Persistence_representations/include/gudhi/Landscape.h
+++ b/src/Persistence_representations/include/gudhi/Landscape.h
@@ -40,64 +40,69 @@
 #include <utility>
 #include <functional>
 
-using PD = std::vector<std::pair<double,double> >;
-
 namespace Gudhi {
 namespace Persistence_representations {
 
 /**
  * \class Landscape gudhi/Landscape.h
- * \brief A class implementing the Landscapes.
+ * \brief A class implementing landscapes.
  *
  * \ingroup Persistence_representations
  *
  * \details
  *
+ * The landscape is a way to turn a persistence diagram into \f$L^2\f$ functions. Roughly, the idea is to see the boundaries of the rank functions as scalar functions taking values on the diagonal.
+ * See \cite bubenik_landscapes_2015 for more details. Here we provide a way to approximate such functions by computing their values on a set of samples.
+ *
 **/
 
 class Landscape {
 
- protected:
-    PD diagram;
-    int res_x, nb_ls;
-    double min_x, max_x;
+  protected:
+   Persistence_diagram diagram;
+   int res_x, nb_ls;
+   double min_x, max_x;
 
   public:
 
-  /** \brief Landscape constructor.
-   * \ingroup Landscape
-   *
-   */
-  Landscape(PD _diagram, int _nb_ls = 5, double _min_x = 0.0, double _max_x = 1.0, int _res_x = 10){diagram = _diagram; nb_ls = _nb_ls; min_x = _min_x; max_x = _max_x; res_x = _res_x;}
-
-  /** \brief Computes the landscape of a diagram.
-   * \ingroup Landscape
-   *
-   */
-  std::vector<std::vector<double> > vectorize() {
-    std::vector<std::vector<double> >  ls; for(int i = 0; i < nb_ls; i++)  ls.emplace_back();
-    int num_pts = diagram.size(); double step = (max_x - min_x)/res_x;
-
-    for(int i = 0; i < res_x; i++){
-      double x = min_x + i*step; double t = x / std::sqrt(2); std::vector<double> events;
-      for(int j = 0; j < num_pts; j++){
-        double px = diagram[j].first; double py = diagram[j].second;
-        if(t >= px && t <= py){  if(t >= (px+py)/2)  events.push_back(std::sqrt(2)*(py-t)); else  events.push_back(std::sqrt(2)*(t-px));  }
-      }
-
-      std::sort(events.begin(), events.end(), [](const double & a, const double & b){return a > b;}); int nb_events = events.size();
-      for (int j = 0; j < nb_ls; j++){  if(j < nb_events)  ls[j].push_back(events[j]);  else  ls[j].push_back(0);  }
-    }
-
-    return ls;
-  }
-
-
-
-
-};
-
-}  // namespace Landscape
+   /** \brief Landscape constructor.
+    * \ingroup Landscape
+    *
+    * @param[in] _diagram    persistence diagram.
+    * @param[in] _nb_ls      number of landscape functions.
+    * @param[in] _min_x      minimum value of samples.
+    * @param[in] _max_x      maximum value of samples.
+    * @param[in] _res_x      number of samples.
+    *
+    */
+   Landscape(const Persistence_diagram & _diagram, int _nb_ls = 5, double _min_x = 0.0, double _max_x = 1.0, int _res_x = 10){diagram = _diagram; nb_ls = _nb_ls; min_x = _min_x; max_x = _max_x; res_x = _res_x;}
+
+   /** \brief Computes the landscape of a diagram.
+    * \ingroup Landscape
+    *
+    */
+   std::vector<std::vector<double> > vectorize() const {
+     std::vector<std::vector<double> >  ls; for(int i = 0; i < nb_ls; i++)  ls.emplace_back();
+     int num_pts = diagram.size(); double step = (max_x - min_x)/res_x;
+
+     for(int i = 0; i < res_x; i++){
+       double x = min_x + i*step; double t = x / std::sqrt(2); std::vector<double> events;
+       for(int j = 0; j < num_pts; j++){
+         double px = diagram[j].first; double py = diagram[j].second;
+         if(t >= px && t <= py){  if(t >= (px+py)/2)  events.push_back(std::sqrt(2)*(py-t)); else  events.push_back(std::sqrt(2)*(t-px));  }
+       }
+
+       std::sort(events.begin(), events.end(), [](const double & a, const double & b){return a > b;}); int nb_events = events.size();
+       for (int j = 0; j < nb_ls; j++){  if(j < nb_events)  ls[j].push_back(events[j]);  else  ls[j].push_back(0);  }
+     }
+     return ls;
+   }
+
+
+
+
+}; // class Landscape
+}  // namespace Persistence_representations
 }  // namespace Gudhi
 
 #endif  // LANDSCAPE_H_
diff --git a/src/Persistence_representations/include/gudhi/Persistence_image.h b/src/Persistence_representations/include/gudhi/Persistence_image.h
index 6c9f75b7..76b34d8d 100644
--- a/src/Persistence_representations/include/gudhi/Persistence_image.h
+++ b/src/Persistence_representations/include/gudhi/Persistence_image.h
@@ -26,8 +26,8 @@
 // gudhi include
 #include <gudhi/read_persistence_from_file.h>
 #include <gudhi/common_persistence_representations.h>
+#include <gudhi/Weight_functions.h>
 #include <gudhi/Debug_utils.h>
-#include <gudhi/Persistence_weighted_gaussian.h>
 
 // standard include
 #include <cmath>
@@ -41,39 +41,49 @@
 #include <utility>
 #include <functional>
 
-using PD = std::vector<std::pair<double,double> >;
-using Weight = std::function<double (std::pair<double,double>) >;
-
 namespace Gudhi {
 namespace Persistence_representations {
 
 /**
  * \class Persistence_image gudhi/Persistence_image.h
- * \brief A class implementing the Persistence Images.
+ * \brief A class implementing the persistence images.
  *
  * \ingroup Persistence_representations
  *
  * \details
  *
+ * Persistence images are a way to build images from persistence diagrams. Roughly, the idea is to center Gaussians on each diagram point, with a weight that usually depends on
+ * the distance to the diagonal, so that the diagram is turned into a function, and then to discretize the plane into pixels, and integrate this function on each pixel.
+ * See \cite Persistence_Images_2017 for more details.
+ *
 **/
 
 class Persistence_image {
 
  protected:
-    PD diagram;
-    int res_x, res_y;
-    double min_x, max_x, min_y, max_y;
-    Weight weight;
-    double sigma;
+  Persistence_diagram diagram;
+  int res_x, res_y;
+  double min_x, max_x, min_y, max_y;
+  Weight weight;
+  double sigma;
 
  public:
 
   /** \brief Persistence Image constructor.
    * \ingroup Persistence_image
    *
+   * @param[in] _diagram   persistence diagram.
+   * @param[in] _min_x     minimum value of pixel abscissa.
+   * @param[in] _max_x     maximum value of pixel abscissa.
+   * @param[in] _res_x     number of pixels for the x-direction.
+   * @param[in] _min_y     minimum value of pixel ordinate.
+   * @param[in] _max_y     maximum value of pixel ordinate.
+   * @param[in] _res_y     number of pixels for the y-direction.
+   * @param[in] _weight    weight function for the Gaussians.
+   * @param[in] _sigma     bandwidth parameter for the Gaussians.
+   *
    */
-  Persistence_image(PD _diagram, double _min_x = 0.0, double _max_x = 1.0, int _res_x = 10, double _min_y = 0.0, double _max_y = 1.0, int _res_y = 10,
-                    Weight _weight = Gudhi::Persistence_representations::Persistence_weighted_gaussian::arctan_weight(1,1), double _sigma = 1.0){
+  Persistence_image(const Persistence_diagram & _diagram, double _min_x = 0.0, double _max_x = 1.0, int _res_x = 10, double _min_y = 0.0, double _max_y = 1.0, int _res_y = 10, const Weight & _weight = arctan_weight(1,1), double _sigma = 1.0){
       diagram = _diagram; min_x = _min_x; max_x = _max_x; res_x = _res_x; min_y = _min_y; max_y = _max_y; res_y = _res_y, weight = _weight; sigma = _sigma;
   }
 
@@ -81,7 +91,7 @@ class Persistence_image {
    * \ingroup Persistence_image
    *
    */
-  std::vector<std::vector<double> > vectorize() {
+  std::vector<std::vector<double> > vectorize() const {
     std::vector<std::vector<double> > im; for(int i = 0; i < res_y; i++)  im.emplace_back();
     double step_x = (max_x - min_x)/res_x; double step_y = (max_y - min_y)/res_y;
 
@@ -109,9 +119,8 @@ class Persistence_image {
 
 
 
-};
-
-}  // namespace Persistence_image
+}; // class Persistence_image
+}  // namespace Persistence_representations
 }  // namespace Gudhi
 
 #endif  // PERSISTENCE_IMAGE_H_
diff --git a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h
index 9a63fccd..76c43e65 100644
--- a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h
+++ b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h
@@ -26,6 +26,7 @@
 // gudhi include
 #include <gudhi/read_persistence_from_file.h>
 #include <gudhi/common_persistence_representations.h>
+#include <gudhi/Weight_functions.h>
 
 // standard include
 #include <cmath>
@@ -39,19 +40,16 @@
 #include <utility>
 #include <functional>
 
-using PD = std::vector<std::pair<double,double> >;
-using Weight = std::function<double (std::pair<double,double>) >;
-
 namespace Gudhi {
 namespace Persistence_representations {
 /**
  * \class Persistence_weighted_gaussian gudhi/Persistence_weighted_gaussian.h
- * \brief A class implementing the Persistence Weighted Gaussian Kernel and a specific case of it called the Persistence Scale Space Kernel.
+ * \brief A class implementing the Persistence Weighted Gaussian kernel and a specific case thereof called the Persistence Scale Space kernel.
  *
  * \ingroup Persistence_representations
  *
  * \details
- * The Persistence Weighted Gaussian Kernel is built with Gaussian Kernel Mean Embedding, meaning that each persistence diagram is first
+ * The Persistence Weighted Gaussian kernel is built with Gaussian Kernel Mean Embedding, meaning that each persistence diagram is first
  * sent to the Hilbert space of a Gaussian kernel with bandwidth parameter \f$\sigma >0\f$ using a weighted mean embedding \f$\Phi\f$:
  *
  * \f$ \Phi\,:\,D\,\rightarrow\,\sum_{p\in D}\,w(p)\,{\rm exp}\left(-\frac{\|p-\cdot\|_2^2}{2\sigma^2}\right) \f$,
@@ -69,59 +67,41 @@ namespace Persistence_representations {
  * in the diagrams. This time can be improved by computing approximations of the kernel
  * with \f$m\f$ Fourier features \cite Rahimi07randomfeatures. In that case, the computation time becomes \f$O(mn)\f$.
  *
- * The Persistence Scale Space Kernel is a Persistence Weighted Gaussian Kernel between modified diagrams:
+ * The Persistence Scale Space kernel is a Persistence Weighted Gaussian kernel between modified diagrams:
  * the symmetric of each point with respect to the diagonal is first added in each diagram, and then the weight function
  * is set to be +1 if the point is above the diagonal and -1 otherwise.
  * 
- * For more details, please consult <i>Persistence Weighted Kernel for Topological Data Analysis</i>\cite Kusano_Fukumizu_Hiraoka_PWGK 
- * and <i>A Stable Multi-Scale Kernel for Topological Machine Learning</i>\cite Reininghaus_Huber_ALL_PSSK . 
- * It implements the following concepts: Topological_data_with_distances, Topological_data_with_scalar_product.
+ * For more details, please see \cite Kusano_Fukumizu_Hiraoka_PWGK
+ * and \cite Reininghaus_Huber_ALL_PSSK .
  *
 **/
 class Persistence_weighted_gaussian{
 
  protected:
-    PD diagram;
+    Persistence_diagram diagram;
     Weight weight;
     double sigma;
     int approx;
 
  public:
 
-  /** \brief Persistence Weighted Gaussian Kernel constructor.
+  /** \brief Persistence Weighted Gaussian kernel constructor.
    * \ingroup Persistence_weighted_gaussian
    *
    * @param[in] _diagram       persistence diagram.
-   * @param[in] _sigma         bandwidth parameter of the Gaussian Kernel used for the Kernel Mean Embedding of the diagrams.
+   * @param[in] _sigma         bandwidth parameter of the Gaussian kernel used for the Kernel Mean Embedding of the diagrams.
    * @param[in] _approx        number of random Fourier features in case of approximate computation, set to -1 for exact computation.
    * @param[in] _weight        weight function for the points in the diagrams.
    *
    */
-  Persistence_weighted_gaussian(PD _diagram, double _sigma = 1.0, int _approx = 1000, Weight _weight = arctan_weight(1,1)){diagram = _diagram; sigma = _sigma; approx = _approx; weight = _weight;}
-  
-  PD get_diagram() const {return this->diagram;}
-  double get_sigma() const {return this->sigma;}
-  int get_approx() const {return this->approx;}
-  Weight get_weight() const {return this->weight;}
+  Persistence_weighted_gaussian(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 1000, const Weight & _weight = arctan_weight(1,1)){diagram = _diagram; sigma = _sigma; approx = _approx; weight = _weight;}
 
 
   // **********************************
   // Utils.
   // **********************************
 
-  /** \brief Specific weight of Persistence Scale Space Kernel.
-   * \ingroup Persistence_weighted_gaussian
-   *
-   * @param[in] p point in 2D.
-   *
-   */
-  static double pss_weight(std::pair<double,double> p)     {if(p.second > p.first)  return 1; else return -1;}
-  static double linear_weight(std::pair<double,double> p)  {return std::abs(p.second - p.first);}
-  static double const_weight(std::pair<double,double> p)   {return 1;}
-  static std::function<double (std::pair<double,double>) > arctan_weight(double C, double power)  {return [=](std::pair<double,double> p){return C * atan(std::pow(std::abs(p.second - p.first), power));};}
-
-
-  std::vector<std::pair<double,double> > Fourier_feat(PD diag, std::vector<std::pair<double,double> > z, Weight weight = arctan_weight(1,1)){
+  std::vector<std::pair<double,double> > Fourier_feat(const Persistence_diagram & diag, const std::vector<std::pair<double,double> > & z, const Weight & weight = arctan_weight(1,1)) const {
     int md = diag.size(); std::vector<std::pair<double,double> > b; int mz = z.size();
     for(int i = 0; i < mz; i++){
       double d1 = 0; double d2 = 0; double zx = z[i].first; double zy = z[i].second;
@@ -135,7 +115,7 @@ class Persistence_weighted_gaussian{
     return b;
   }
 
-  std::vector<std::pair<double,double> > random_Fourier(double sigma, int m = 1000){
+  std::vector<std::pair<double,double> > random_Fourier(double sigma, int m = 1000) const {
     std::normal_distribution<double> distrib(0,1); std::vector<std::pair<double,double> > z; std::random_device rd;
     for(int i = 0; i < m; i++){
       std::mt19937 e1(rd()); std::mt19937 e2(rd());
@@ -154,12 +134,14 @@ class Persistence_weighted_gaussian{
   /** \brief Evaluation of the kernel on a pair of diagrams.
    * \ingroup Persistence_weighted_gaussian
    *
-   * @param[in] second other instance of class Persistence_weighted_gaussian. Warning: sigma, approx and weight parameters need to be the same for both instances!!! 
+   * @pre       sigma, approx and weight attributes need to be the same for both instances.
+   * @param[in] second other instance of class Persistence_weighted_gaussian.
    *
    */
-  double compute_scalar_product(Persistence_weighted_gaussian second){
+  double compute_scalar_product(const Persistence_weighted_gaussian & second) const {
 
-    PD diagram1 = this->diagram; PD diagram2 = second.diagram;
+    GUDHI_CHECK(this->sigma != second.sigma || this->approx != second.approx || this->weight != second.weight, std::invalid_argument("Error: different values for representations"));
+    Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram;
 
     if(this->approx == -1){
       int num_pts1 = diagram1.size(); int num_pts2 = diagram2.size(); double k = 0;
@@ -171,7 +153,7 @@ class Persistence_weighted_gaussian{
       return k;
     }
     else{
-      std::vector<std::pair<double,double> > z =  random_Fourier(this->sigma, this->approx);
+      std::vector<std::pair<double,double> > z  = random_Fourier(this->sigma, this->approx);
       std::vector<std::pair<double,double> > b1 = Fourier_feat(diagram1,z,this->weight);
       std::vector<std::pair<double,double> > b2 = Fourier_feat(diagram2,z,this->weight);
       double d = 0; for(int i = 0; i < this->approx; i++) d += b1[i].first*b2[i].first + b1[i].second*b2[i].second;
@@ -182,20 +164,18 @@ class Persistence_weighted_gaussian{
   /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel.
    * \ingroup Persistence_weighted_gaussian
    *
-   * @param[in] second other instance of class Persistence_weighted_gaussian. Warning: sigma, approx and weight parameters need to be the same for both instances!!! 
+   * @pre       sigma, approx and weight attributes need to be the same for both instances.
+   * @param[in] second other instance of class Persistence_weighted_gaussian.
    *
    */
-  double distance(Persistence_weighted_gaussian second) {
-    if(this->sigma != second.get_sigma() || this->approx != second.get_approx()){
-      std::cout << "Error: different representations!" << std::endl; return 0;
-    }
-    else return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), 0.5);
+  double distance(const Persistence_weighted_gaussian & second) const {
+    GUDHI_CHECK(this->sigma != second.sigma || this->approx != second.approx || this->weight != second.weight, std::invalid_argument("Error: different values for representations"));
+    return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), 0.5);
   }
 
 
-};
-
-}  // namespace Persistence_weighted_gaussian
+}; // class Persistence_weighted_gaussian
+}  // namespace Persistence_representations
 }  // namespace Gudhi
 
 #endif  // PERSISTENCE_WEIGHTED_GAUSSIAN_H_
diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
index 6a9a607e..235918fe 100644
--- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
+++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
@@ -40,19 +40,17 @@
 #include <utility>
 #include <functional>
 
-using PD = std::vector<std::pair<double,double> >;
-
 namespace Gudhi {
 namespace Persistence_representations {
 
 /**
  * \class Sliced_Wasserstein gudhi/Sliced_Wasserstein.h
- * \brief A class implementing the Sliced Wasserstein Kernel.
+ * \brief A class implementing the Sliced Wasserstein kernel.
  *
  * \ingroup Persistence_representations
  *
  * \details
- * The Sliced Wasserstein Kernel is defined as a Gaussian-like Kernel between persistence diagrams, where the distance used for
+ * 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:
  *
@@ -65,15 +63,14 @@ namespace Persistence_representations {
  *
  * \f$ k(D_1,D_2) = {\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right).\f$
  * 
- * For more details, please consult <i>Sliced Wasserstein Kernel for Persistence Diagrams</i>\cite pmlr-v70-carriere17a .
- * It implements the following concepts: Topological_data_with_distances, Topological_data_with_scalar_product.
+ * For more details, please see \cite pmlr-v70-carriere17a .
  *
 **/
 
 class Sliced_Wasserstein {
 
  protected:
-    PD diagram;
+    Persistence_diagram diagram;
     int approx;
     double sigma;
     std::vector<std::vector<double> > projections, projections_diagonal;
@@ -107,7 +104,7 @@ class Sliced_Wasserstein {
 
   }
 
-  /** \brief Sliced Wasserstein Kernel constructor.
+  /** \brief Sliced Wasserstein kernel constructor.
    * \ingroup Sliced_Wasserstein
    *
    * @param[in] _diagram  persistence diagram.
@@ -115,21 +112,14 @@ class Sliced_Wasserstein {
    * @param[in] _approx number of directions used to approximate the integral in the Sliced Wasserstein distance, set to -1 for exact computation.
    *
    */
-  Sliced_Wasserstein(PD _diagram, double _sigma = 1.0, int _approx = 100){diagram = _diagram; approx = _approx; sigma = _sigma; build_rep();}
-  
-  PD get_diagram() const {return this->diagram;}
-  int get_approx() const {return this->approx;}
-  double get_sigma() const {return this->sigma;}
-
-  
-
+  Sliced_Wasserstein(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 100){diagram = _diagram; approx = _approx; sigma = _sigma; build_rep();}
 
   // **********************************
   // Utils.
   // **********************************
 
   // Compute the angle formed by two points of a PD
-  double compute_angle(PD diag, int i, int j){
+  double compute_angle(const Persistence_diagram & diag, int i, int j) const {
     std::pair<double,double> vect; double x1,y1, x2,y2;
     x1 = diag[i].first; y1 = diag[i].second;
     x2 = diag[j].first; y2 = diag[j].second;
@@ -150,7 +140,7 @@ class Sliced_Wasserstein {
   }
 
   // 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(const double & alpha, const double & beta){
+  double compute_int_cos(double alpha, double beta) const {
     double res = 0;
     if (alpha >= 0 && alpha <= pi){
       if (cos(alpha) >= 0){
@@ -175,13 +165,13 @@ class Sliced_Wasserstein {
     return res;
   }
 
-  double compute_int(const double & theta1, const double & theta2, const int & p, const int & q, const PD & PD1, const PD & PD2){
-    double norm = std::sqrt(  (PD1[p].first-PD2[q].first)*(PD1[p].first-PD2[q].first) + (PD1[p].second-PD2[q].second)*(PD1[p].second-PD2[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 (PD1[p].first > PD2[q].first)
-      angle1 = theta1 - asin( (PD1[p].second-PD2[q].second)/norm  );
+    if (diag1[p].first > diag2[q].first)
+      angle1 = theta1 - asin( (diag1[p].second-diag2[q].second)/norm  );
     else
-      angle1 = theta1 - asin( (PD2[q].second-PD1[p].second)/norm  );
+      angle1 = theta1 - asin( (diag2[q].second-diag1[p].second)/norm  );
     double angle2 = angle1 + theta2 - theta1;
     double integral = compute_int_cos(angle1,angle2);
     return norm*integral;
@@ -197,134 +187,133 @@ class Sliced_Wasserstein {
   /** \brief Evaluation of the Sliced Wasserstein Distance between a pair of diagrams.
    * \ingroup Sliced_Wasserstein
    *
+   * @pre       approx attribute needs to be the same for both instances.
    * @param[in] second other instance of class Sliced_Wasserstein. 
-   * For warning in red:
-   * @warning approx parameter needs to be the same for both instances.
+   *
    *
    */
-  double compute_sliced_wasserstein_distance(Sliced_Wasserstein second) {
+  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"));
+    GUDHI_CHECK(this->approx != second.approx, std::invalid_argument("Error: different approx values for representations"));
 
-      PD diagram1 = this->diagram; PD diagram2 = second.diagram; double sw = 0;
+    Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram; double sw = 0;
 
-      if(this->approx == -1){
+    if(this->approx == -1){
 
-        // Add projections onto diagonal.
-        int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double max_ordinate = std::numeric_limits<double>::lowest();
-        for (int i = 0; i < n2; i++){
-          max_ordinate = std::max(max_ordinate, diagram2[i].second);
-          diagram1.emplace_back(  (diagram2[i].first+diagram2[i].second)/2, (diagram2[i].first+diagram2[i].second)/2  );
-        }
-        for (int i = 0; i < n1; i++){
-          max_ordinate = std::max(max_ordinate, diagram1[i].second);
-          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.
-        int mag = 0; while(max_ordinate > 10){mag++; max_ordinate/=10;}
-        double thresh = pow(10,-5+mag);
-        srand(time(NULL));
-        for (int i = 0; i < num_pts_dgm; i++){
-          diagram1[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); diagram1[i].second += thresh*(1.0-2.0*rand()/RAND_MAX);
-          diagram2[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); diagram2[i].second += thresh*(1.0-2.0*rand()/RAND_MAX);
-        }
+      // Add projections onto diagonal.
+      int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double max_ordinate = std::numeric_limits<double>::lowest();
+      for (int i = 0; i < n2; i++){
+        max_ordinate = std::max(max_ordinate, diagram2[i].second);
+        diagram1.emplace_back(  (diagram2[i].first+diagram2[i].second)/2, (diagram2[i].first+diagram2[i].second)/2  );
+      }
+      for (int i = 0; i < n1; i++){
+        max_ordinate = std::max(max_ordinate, diagram1[i].second);
+        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.
+      int mag = 0; while(max_ordinate > 10){mag++; max_ordinate/=10;}
+      double thresh = pow(10,-5+mag);
+      srand(time(NULL));
+      for (int i = 0; i < num_pts_dgm; i++){
+        diagram1[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); diagram1[i].second += thresh*(1.0-2.0*rand()/RAND_MAX);
+        diagram2[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); diagram2[i].second += thresh*(1.0-2.0*rand()/RAND_MAX);
+      }
 
-        // 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));
-          }
+      // 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));
         }
+      }
 
-        // Sort angles.
-        std::sort(angles1.begin(), angles1.end(), [=](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(), [=](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); } );
-
-        // 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++)  for (int j = 0; j < num_pts_dgm; j++)  if(orderp1[j] == i){  order1[i] = j; break;  }
-        for(int i = 0; i < num_pts_dgm; i++)  for (int j = 0; j < num_pts_dgm; j++)  if(orderp2[j] == i){  order2[i] = j; break;  }
-
-        // 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);
-
-        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;
-        }
+      // Sort angles.
+      std::sort(angles1.begin(), angles1.end(), [=](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(), [=](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); } );
+
+      // 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++)  for (int j = 0; j < num_pts_dgm; j++)  if(orderp1[j] == i){  order1[i] = j; break;  }
+      for(int i = 0; i < num_pts_dgm; i++)  for (int j = 0; j < num_pts_dgm; j++)  if(orderp2[j] == i){  order2[i] = j; break;  }
+
+      // 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);
+
+      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;
+      }
 
-        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;
-        }
+      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 < 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);
-            theta1 = theta2;
-            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);
-          }
+      // 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);
+          theta1 = theta2;
+          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;
-
-        for (int i = 0; i < this->approx; i++){
+    else{
 
-          std::vector<double> v1; std::vector<double> l1 = this->projections[i]; std::vector<double> l1bis = second.projections_diagonal[i]; std::merge(l1.begin(), l1.end(), l1bis.begin(), l1bis.end(), std::back_inserter(v1));
-          std::vector<double> v2; std::vector<double> l2 = second.projections[i]; std::vector<double> l2bis = this->projections_diagonal[i]; std::merge(l2.begin(), l2.end(), l2bis.begin(), l2bis.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;
+      double step = pi/this->approx;
+      for (int i = 0; i < this->approx; i++){
 
-        }
+        std::vector<double> v1; std::vector<double> l1 = this->projections[i]; std::vector<double> l1bis = second.projections_diagonal[i]; std::merge(l1.begin(), l1.end(), l1bis.begin(), l1bis.end(), std::back_inserter(v1));
+        std::vector<double> v2; std::vector<double> l2 = second.projections[i]; std::vector<double> l2bis = this->projections_diagonal[i]; std::merge(l2.begin(), l2.end(), l2bis.begin(), l2bis.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;
   }
 
   /** \brief Evaluation of the kernel on a pair of diagrams.
    * \ingroup Sliced_Wasserstein
    *
-   * @param[in] second other instance of class Sliced_Wasserstein. Warning: sigma and approx parameters need to be the same for both instances!!! 
+   * @pre       approx and sigma attributes need to be the same for both instances.
+   * @param[in] second other instance of class Sliced_Wasserstein.
    *
    */
-  double compute_scalar_product(Sliced_Wasserstein second){
+  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));
   }
@@ -332,10 +321,11 @@ class Sliced_Wasserstein {
   /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel.
    * \ingroup Sliced_Wasserstein
    *
-   * @param[in] second  other instance of class Sliced_Wasserstein. Warning: sigma and approx parameters need to be the same for both instances!!! 
+   * @pre       approx and sigma attributes need to be the same for both instances.
+   * @param[in] second  other instance of class Sliced_Wasserstein.
    *
    */
-  double distance(Sliced_Wasserstein second) {
+  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);
   }
@@ -343,9 +333,8 @@ class Sliced_Wasserstein {
 
 
 
-};
-
-}  // namespace Sliced_Wasserstein
+}; // class Sliced_Wasserstein
+}  // namespace Persistence_representations
 }  // namespace Gudhi
 
 #endif  // SLICED_WASSERSTEIN_H_
diff --git a/src/Persistence_representations/include/gudhi/common_persistence_representations.h b/src/Persistence_representations/include/gudhi/common_persistence_representations.h
index 90f2626d..884fce58 100644
--- a/src/Persistence_representations/include/gudhi/common_persistence_representations.h
+++ b/src/Persistence_representations/include/gudhi/common_persistence_representations.h
@@ -28,17 +28,28 @@
 #include <cmath>
 #include <boost/math/constants/constants.hpp>
 
+/**
+ * In this module, we use the name Persistence_diagram for the representation of a diagram in a vector of pairs of two double.
+ */
+using Persistence_diagram = std::vector<std::pair<double,double> >;
+
+/**
+ * In this module, we use the name Weight for the representation of a function taking a pair of two double and returning a double.
+ */
+using Weight = std::function<double (std::pair<double,double>) >;
+
 namespace Gudhi {
 namespace Persistence_representations {
 // this file contain an implementation of some common procedures used in Persistence_representations.
 
 static constexpr double pi = boost::math::constants::pi<double>();
+
 // double epsi = std::numeric_limits<double>::epsilon();
 double epsi = 0.000005;
 
 /**
  *  A procedure used to compare doubles. Typically given two doubles A and B, comparing A == B is not good idea. In this
- *case, we use the procedure almostEqual with the epsi defined at
+ *  case, we use the procedure almostEqual with the epsi defined at
  *  the top of the file. Setting up the epsi gives the user a tolerance on what should be consider equal.
 **/
 inline bool almost_equal(double a, double b) {
@@ -55,8 +66,7 @@ double birth_plus_deaths(std::pair<double, double> a) { return a.first + a.secon
 
 // landscapes
 /**
- * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two
- *points.
+ * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two points.
 **/
 std::pair<double, double> compute_parameters_of_a_line(std::pair<double, double> p1, std::pair<double, double> p2) {
   double a = (p2.second - p1.second) / (p2.first - p1.first);
@@ -66,8 +76,7 @@ std::pair<double, double> compute_parameters_of_a_line(std::pair<double, double>
 
 // landscapes
 /**
- * This procedure given two points which lies on the opposite sides of x axis, compute x for which the line connecting
- *those two points crosses x axis.
+ * This procedure given two points which lies on the opposite sides of x axis, compute x for which the line connecting those two points crosses x axis.
 **/
 double find_zero_of_a_line_segment_between_those_two_points(std::pair<double, double> p1,
                                                             std::pair<double, double> p2) {
@@ -91,8 +100,7 @@ double find_zero_of_a_line_segment_between_those_two_points(std::pair<double, do
 // landscapes
 /**
  * This method provides a comparison of points that is used in construction of persistence landscapes. The ordering is
- *lexicographical for the first coordinate, and reverse-lexicographical for the
- * second coordinate.
+ * lexicographical for the first coordinate, and reverse-lexicographical for the second coordinate.
 **/
 bool compare_points_sorting(std::pair<double, double> f, std::pair<double, double> s) {
   if (f.first < s.first) {