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

Former-commit-id: 12f32a1c8ca31e7e0a40e1c3502e2a3d810d5bc5

9 files changed, 158 insertions, 634 deletions

diff --git a/src/Persistence_representations/include/gudhi/Betti_sequence.h b/src/Persistence_representations/include/gudhi/Betti_sequence.h
deleted file mode 100644
index 57c52ad2..00000000
--- a/src/Persistence_representations/include/gudhi/Betti_sequence.h
+++ /dev/null
@@ -1,95 +0,0 @@
-/*    This file is part of the Gudhi Library. The Gudhi library
- *    (Geometric Understanding in Higher Dimensions) is a generic C++
- *    library for computational topology.
- *
- *    Author(s):       Mathieu Carriere
- *
- *    Copyright (C) 2018  INRIA (France)
- *
- *    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
- *    the Free Software Foundation, either version 3 of the License, or
- *    (at your option) any later version.
- *
- *    This program is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef BETTI_SEQUENCE_H_
-#define BETTI_SEQUENCE_H_
-
-// gudhi include
-#include <gudhi/read_persistence_from_file.h>
-#include <gudhi/common_persistence_representations.h>
-#include <gudhi/Debug_utils.h>
-
-// standard include
-#include <cmath>
-#include <iostream>
-#include <vector>
-#include <limits>
-#include <fstream>
-#include <sstream>
-#include <algorithm>
-#include <string>
-#include <utility>
-#include <functional>
-
-namespace Gudhi {
-namespace Persistence_representations {
-
-/**
- * \class Betti_sequence gudhi/Betti_sequence.h
- * \brief A class implementing Betti sequences
- *
- * \ingroup Persistence_representations
- *
- * \details
-**/
-
-class Betti_sequence {
-
-  protected:
-   Persistence_diagram diagram;
-   int res_x, nb_cv;
-   double min_x, max_x;
-
-  public:
-
-   /** \brief Betti_sequence constructor.
-    * \ingroup Betti_sequence
-    *
-    * @param[in] _diagram    persistence diagram.
-    * @param[in] _min_x      minimum value of samples.
-    * @param[in] _max_x      maximum value of samples.
-    * @param[in] _res_x      number of samples.
-    *
-    */
-   Betti_sequence(const Persistence_diagram & _diagram, double _min_x = 0.0, double _max_x = 1.0, int _res_x = 10){diagram = _diagram; min_x = _min_x; max_x = _max_x; res_x = _res_x;}
-
-   /** \brief Computes the Betti sequences of a diagram.
-    * \ingroup Betti_sequence
-    *
-    */
-   std::vector<int> vectorize() const {
-     int num_pts = diagram.size(); double step = (max_x - min_x)/(res_x - 1);
-     std::vector<int> bs(res_x); for(int i = 0; i < res_x; i++)  bs[i] = 0;
-     for(int j = 0; j < num_pts; j++){
-       double px = diagram[j].first; double py = diagram[j].second;
-       int first = std::ceil((px-min_x)/step); int last = std::ceil((py-min_x)/step);
-       for(int i = first; i < last; i++)  bs[i] += 1;
-     }
-
-     return bs;
-   }
-
-}; // class Betti_sequence
-}  // namespace Persistence_representations
-}  // namespace Gudhi
-
-#endif  // BETTI_SEQUENCE_H_
diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
index 35e51e63..63c6e239 100644
--- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
+++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
@@ -245,6 +245,20 @@ class Persistence_heat_maps {
                         unsigned dimension = std::numeric_limits<unsigned>::max());
 
   /**
+   * Construction that takes as inputs (1) the diagram, (2) grid parameters (min, max and number of samples for x and y axes), and (3) a universal kernel on the plane used
+   * to turn the diagram into a function.
+  **/
+  Persistence_heat_maps(const Persistence_diagram & interval, 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, const Kernel & kernel = Gaussian_kernel(1.0));
+
+  /**
+   * Construction that takes as inputs (1) the diagram and (2) a universal kernel on the plane used
+   * to turn the diagram into a function. Note that this construction is infinite dimensional so
+   * only compute_scalar_product() method is valid after calling this constructor.
+  **/
+  Persistence_heat_maps(const Persistence_diagram & interval, const Kernel & kernel = Gaussian_kernel(1.0));
+
+  /**
    * 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.
@@ -512,15 +526,27 @@ class Persistence_heat_maps {
   size_t number_of_functions_for_projections_to_reals;
   void construct(const std::vector<std::pair<double, double> >& intervals_,
                  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());
 
+  void construct_image_from_exact_universal_kernel(const Persistence_diagram & interval,
+                                                   size_t number_of_x_pixels = 10, size_t number_of_y_pixels = 10,
+                                                   double min_x = 0, double max_x = 1, double min_y = 0, double max_y = 1, const Kernel & kernel = Gaussian_kernel(1.0));
+  void construct_kernel_from_exact_universal_kernel(const Persistence_diagram & interval, const Kernel & kernel = Gaussian_kernel(1.0));
+
   void set_up_parameters_for_basic_classes() {
     this->number_of_functions_for_vectorization = 1;
     this->number_of_functions_for_projections_to_reals = 1;
   }
 
+  // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false)
+  bool discrete = true;
+
+  // PWGK
+  Kernel k;
+  Persistence_diagram d;
+  std::vector<double> weights;
+
   // data
   Scalling_of_kernels f;
   bool erase_below_diagonal;
@@ -529,6 +555,59 @@ class Persistence_heat_maps {
   std::vector<std::vector<double> > heat_map;
 };
 
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::construct_image_from_exact_universal_kernel(const Persistence_diagram & diagram,
+                                                                                             size_t number_of_x_pixels, size_t number_of_y_pixels,
+                                                                                             double min_x, double max_x,
+                                                                                             double min_y, double max_y, const Kernel & kernel) {
+
+  this->discrete = true; Scalling_of_kernels f; this->f = f; this->min_ = min_x; this->max_ = max_x;
+  for(size_t i = 0; i < number_of_y_pixels; i++)  this->heat_map.emplace_back();
+  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 = diagram.size();
+
+  for(size_t i = 0; i < number_of_y_pixels; i++){
+    double y = min_y + i*step_y;
+    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++){
+        double px = diagram[k].first; double py = diagram[k].second; std::pair<double, double> diagram_point(px,py);
+        pixel_value += this->f(diagram_point) * kernel(diagram_point, grid_point);
+      }
+      this->heat_map[i].push_back(pixel_value);
+
+    }
+  }
+
+}
+
+
+template <typename Scalling_of_kernels>
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(const Persistence_diagram & diagram,
+                                                                  size_t number_of_x_pixels, size_t number_of_y_pixels,
+                                                                  double min_x, double max_x,
+                                                                  double min_y, double max_y, const Kernel & kernel) {
+  this->construct_image_from_exact_universal_kernel(diagram, number_of_x_pixels, number_of_y_pixels, min_x, max_x, min_y, max_y, kernel);
+  this->set_up_parameters_for_basic_classes();
+}
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::construct_kernel_from_exact_universal_kernel(const Persistence_diagram & diagram, const Kernel & kernel){
+  this->discrete = false; Scalling_of_kernels f; this->f = f; this->k = kernel; this->d = diagram;
+  int num_pts = this->d.size();
+  for (int i = 0; i < num_pts; i++)  this->weights.push_back(this->f(this->d[i]));
+}
+
+
+template <typename Scalling_of_kernels>
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(const Persistence_diagram& diagram, const Kernel & kernel) {
+  this->construct_kernel_from_exact_universal_kernel(diagram, kernel);
+  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>
 void Persistence_heat_maps<Scalling_of_kernels>::construct(const std::vector<std::pair<double, double> >& intervals_,
@@ -826,13 +905,16 @@ 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;}
+
   // convert this->heat_map into one large vector:
   size_t size_of_result = 0;
   for (size_t i = 0; i != this->heat_map.size(); ++i) {
     size_of_result += this->heat_map[i].size();
   }
 
-  std::vector<double> result;
   result.reserve(size_of_result);
 
   for (size_t i = 0; i != this->heat_map.size(); ++i) {
@@ -846,34 +928,39 @@ 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 {
-  // 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 "
-                 "them. The program will now terminate";
-    throw "The persistence images are of non compatible sizes. The program will now terminate";
-  }
+  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 "
+                   "them. The program will now terminate";
+      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);
+    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);
+        }
       }
-    }
-  } 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]);
+    } 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]);
+          }
         }
       }
     }
+    return distance;
+  } else {
+
+    return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) -2 * this->compute_scalar_product(second));
+
   }
-  return distance;
 }
 
 template <typename Scalling_of_kernels>
@@ -895,22 +982,37 @@ 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 {
-  // 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 "
-                 "them. The program will now terminate";
-    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 scalar product:
-  double scalar_prod = 0;
-  for (size_t i = 0; i != this->heat_map.size(); ++i) {
-    for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
-      scalar_prod += this->heat_map[i][j] * second.heat_map[i][j];
+  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 "
+                   "them. The program will now terminate";
+      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 scalar product:
+    double scalar_prod = 0;
+    for (size_t i = 0; i != this->heat_map.size(); ++i) {
+      for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
+        scalar_prod += this->heat_map[i][j] * second.heat_map[i][j];
+      }
+    }
+    return scalar_prod;
   }
-  return scalar_prod;
+
+  else{
+    GUDHI_CHECK(this->approx != second.approx || this->f != second.f, std::invalid_argument("Error: different values for representations"));
+
+    int num_pts1 = this->d.size(); int num_pts2 = second.d.size(); double kernel_val = 0;
+    for(int i = 0; i < num_pts1; i++)
+      for(int j = 0; j < num_pts2; j++)
+        kernel_val += this->weights[i] * second.weights[j] * this->k(this->d[i], second.d[j]);
+    return kernel_val;
+  }
+
+
 }
 
 }  // namespace Persistence_representations
diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps_exact.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps_exact.h
deleted file mode 100644
index 7c5b2fdc..00000000
--- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps_exact.h
+++ /dev/null
@@ -1,125 +0,0 @@
-/*    This file is part of the Gudhi Library. The Gudhi library
- *    (Geometric Understanding in Higher Dimensions) is a generic C++
- *    library for computational topology.
- *
- *    Author(s):       Mathieu Carriere
- *
- *    Copyright (C) 2018  INRIA (France)
- *
- *    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
- *    the Free Software Foundation, either version 3 of the License, or
- *    (at your option) any later version.
- *
- *    This program is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef PERSISTENCE_HEAT_MAPS_EXACT_H_
-#define PERSISTENCE_HEAT_MAPS_EXACT_H_
-
-// 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>
-
-// standard include
-#include <cmath>
-#include <iostream>
-#include <vector>
-#include <limits>
-#include <fstream>
-#include <sstream>
-#include <algorithm>
-#include <string>
-#include <utility>
-#include <functional>
-
-namespace Gudhi {
-namespace Persistence_representations {
-
-/**
- * \class Persistence_heat_maps_exact gudhi/Persistence_heat_maps_exact.h
- * \brief A class implementing exact persistence heat maps.
- *
- * \ingroup Persistence_representations
- *
- * \details
- *
- * In this class, we propose a way to approximate persistence heat maps, or persistence surfaces, by centering weighted Gaussians on each point of the persistence diagram, and evaluating these (exact) weighted Gaussian functions
- * on the pixels of a 2D grid. Note that this scheme is different from the one proposed in Persistence_heat_maps, which first maps the points of the diagram to a 2D grid, and then evaluates the (approximate) weighted Gaussian functions.
- * Hence, the difference is that we do not modify the diagram in this implementation, but the code can be slower to run.
-**/
-
-class Persistence_heat_maps_exact {
-
- protected:
-  Persistence_diagram diagram;
-  int res_x, res_y;
-  double min_x, max_x, min_y, max_y;
-  Weight weight;
-  double sigma;
-
- public:
-
-  /** \brief Persistence_heat_maps_exact constructor.
-   * \ingroup Persistence_heat_maps_exact
-   *
-   * @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_heat_maps_exact(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;
-  }
-
-  /** \brief Computes the persistence image of a diagram.
-   * \ingroup Persistence_heat_maps_exact
-   *
-   */
-  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 - 1); double step_y = (max_y - min_y)/(res_y - 1);
-
-    int num_pts = diagram.size();
-
-    for(int i = 0; i < res_y; i++){
-      double y = min_y + i*step_y;
-      for(int j = 0; j < res_x; j++){
-        double x = min_x + j*step_x;
-
-        double pixel_value = 0;
-        for(int k = 0; k < num_pts; k++){
-          double px = diagram[k].first; double py = diagram[k].second;
-          pixel_value += weight(std::pair<double,double>(px,py)) * std::exp(   -((x-px)*(x-px) + (y-(py-px))*(y-(py-px))) / (2*sigma*sigma)   ) / (sigma*std::sqrt(2*pi));
-        }
-        im[i].push_back(pixel_value);
-
-      }
-    }
-
-    return im;
-
-  }
-
-
-
-
-}; // class Persistence_heat_maps_exact
-}  // namespace Persistence_representations
-}  // namespace Gudhi
-
-#endif  // PERSISTENCE_HEAT_MAPS_EXACT_H_
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 fd8a181c..db0e362a 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 * (grid_interval_begin + grid_interval_end));
+    size_t grid_interval_midpoint = (size_t)(0.5 * (p[int_no].first + p[int_no].second) - grid_min + 1);
 
     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 = dx;
+    double landscape_value = grid_min + dx * (grid_interval_begin + 1) - p[int_no].first;
     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,6 +1030,8 @@ 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/Persistence_landscape_on_grid_exact.h b/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid_exact.h
deleted file mode 100644
index 52f24195..00000000
--- a/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid_exact.h
+++ /dev/null
@@ -1,108 +0,0 @@
-/*    This file is part of the Gudhi Library. The Gudhi library
- *    (Geometric Understanding in Higher Dimensions) is a generic C++
- *    library for computational topology.
- *
- *    Author(s):       Mathieu Carriere
- *
- *    Copyright (C) 2018  INRIA (France)
- *
- *    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
- *    the Free Software Foundation, either version 3 of the License, or
- *    (at your option) any later version.
- *
- *    This program is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef LANDSCAPE_H_
-#define LANDSCAPE_H_
-
-// gudhi include
-#include <gudhi/read_persistence_from_file.h>
-#include <gudhi/common_persistence_representations.h>
-#include <gudhi/Debug_utils.h>
-
-// standard include
-#include <cmath>
-#include <iostream>
-#include <vector>
-#include <limits>
-#include <fstream>
-#include <sstream>
-#include <algorithm>
-#include <string>
-#include <utility>
-#include <functional>
-
-namespace Gudhi {
-namespace Persistence_representations {
-
-/**
- * \class Persistence_landscape_on_grid_exact gudhi/Persistence_landscape_on_grid_exact.h
- * \brief A class implementing exact persistence landscapes by approximating them on a collection of grid points
- *
- * \ingroup Persistence_representations
- *
- * \details
- * In this class, we propose a way to approximate landscapes by sampling the x-axis of the persistence diagram and evaluating the (exact) landscape functions on the sample projections onto the diagonal. Note that this is a different approximation scheme 
- * from the one proposed in Persistence_landscape_on_grid, which puts a grid on the diagonal, maps the persistence intervals on this grid and computes the (approximate) landscape functions on the samples. 
- * Hence, the difference is that we do not modify the diagram in this implementation, but the code can be slower to run.
-**/
-
-class Persistence_landscape_on_grid_exact {
-
-  protected:
-   Persistence_diagram diagram;
-   int res_x, nb_ls;
-   double min_x, max_x;
-
-  public:
-
-   /** \brief Persistence_landscape_on_grid_exact constructor.
-    * \ingroup Persistence_landscape_on_grid_exact
-    *
-    * @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.
-    *
-    */
-   Persistence_landscape_on_grid_exact(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 approximation of a diagram.
-    * \ingroup Persistence_landscape_on_grid_exact
-    *
-    */
-   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 - 1);
-
-     std::vector<std::vector<double> > ls_t; for(int i = 0; i < res_x; i++)  ls_t.emplace_back();
-     for(int j = 0; j < num_pts; j++){
-       double px = diagram[j].first; double py = diagram[j].second; double mid = (px+py)/2;
-       int first = std::ceil((px-min_x)/step); int middle = std::ceil((mid-min_x)/step); int last = std::ceil((py-min_x)/step); double x = min_x + first*step;
-       for(int i = first; i < middle; i++){  double value = std::sqrt(2)*(x-px); ls_t[i].push_back(value); x += step;  }
-       for(int i = middle; i < last; i++){   double value = std::sqrt(2)*(py-x); ls_t[i].push_back(value); x += step;  }
-     }
-
-     for(int i = 0; i < res_x; i++){
-       std::sort(ls_t[i].begin(), ls_t[i].end(), [](const double & a, const double & b){return a > b;});
-       int nb_events_i = ls_t[i].size();
-       for (int j = 0; j < nb_ls; j++){  if(j < nb_events_i)  ls[j].push_back(ls_t[i][j]);  else  ls[j].push_back(0);  }
-     }
-
-     return ls;
-   }
-
-}; // class Persistence_landscape_on_grid_exact
-}  // namespace Persistence_representations
-}  // namespace Gudhi
-
-#endif  // LANDSCAPE_H_
diff --git a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h
deleted file mode 100644
index 9ef47bf1..00000000
--- a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h
+++ /dev/null
@@ -1,182 +0,0 @@
-/*    This file is part of the Gudhi Library. The Gudhi library
- *    (Geometric Understanding in Higher Dimensions) is a generic C++
- *    library for computational topology.
- *
- *    Author(s):       Mathieu Carriere
- *
- *    Copyright (C) 2018  INRIA (France)
- *
- *    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
- *    the Free Software Foundation, either version 3 of the License, or
- *    (at your option) any later version.
- *
- *    This program is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef PERSISTENCE_WEIGHTED_GAUSSIAN_H_
-#define PERSISTENCE_WEIGHTED_GAUSSIAN_H_
-
-// gudhi include
-#include <gudhi/read_persistence_from_file.h>
-#include <gudhi/common_persistence_representations.h>
-#include <gudhi/Weight_functions.h>
-
-// standard include
-#include <cmath>
-#include <iostream>
-#include <vector>
-#include <limits>
-#include <fstream>
-#include <sstream>
-#include <algorithm>
-#include <string>
-#include <utility>
-#include <functional>
-#include <random>
-
-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 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
- * 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$,
- *
- * Usually, the weight function is chosen to be an arctan function of the distance of the point to the diagonal:
- * \f$w(p) = {\rm arctan}(C\,|y-x|^\alpha)\f$, for some parameters \f$C,\alpha >0\f$.
- * Then, their scalar product in this space is computed:
- *
- * \f$ k(D_1,D_2)=\langle\Phi(D_1),\Phi(D_2)\rangle
- * \,=\,\sum_{p\in D_1}\,\sum_{q\in D_2}\,w(p)\,w(q)\,{\rm exp}\left(-\frac{\|p-q\|_2^2}{2\sigma^2}\right).\f$
- *
- * Note that one may apply a second Gaussian kernel to their distance in this space and still get a kernel.
- *
- * It follows that the computation time is \f$O(n^2)\f$ where \f$n\f$ is the number of points
- * 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 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 see \cite Kusano_Fukumizu_Hiraoka_PWGK
- * and \cite Reininghaus_Huber_ALL_PSSK .
- *
-**/
-class Persistence_weighted_gaussian{
-
- protected:
-    Persistence_diagram diagram;
-    Weight weight;
-    double sigma;
-    int approx;
-
- public:
-
-  /** \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] _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(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.
-  // **********************************
-
-  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;
-      for(int j = 0; j < md; j++){
-        double x = diag[j].first; double y = diag[j].second;
-        d1 += weight(diag[j])*cos(x*zx + y*zy);
-        d2 += weight(diag[j])*sin(x*zx + y*zy);
-      }
-      b.emplace_back(d1,d2);
-    }
-    return b;
-  }
-
-  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());
-      double zx = distrib(e1); double zy = distrib(e2);
-      z.emplace_back(zx/sigma,zy/sigma);
-    }
-    return z;
-  }
-
-
-
-  // **********************************
-  // Scalar product + distance.
-  // **********************************
-
-  /** \brief Evaluation of the kernel on a pair of diagrams.
-   * \ingroup Persistence_weighted_gaussian
-   *
-   * @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(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"));
-    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;
-      for(int i = 0; i < num_pts1; i++)
-        for(int j = 0; j < num_pts2; j++)
-          k += this->weight(diagram1[i])*this->weight(diagram2[j])*exp(-((diagram1[i].first  - diagram2[j].first)  *  (diagram1[i].first  - diagram2[j].first) +
-                                                                         (diagram1[i].second - diagram2[j].second) *  (diagram1[i].second - diagram2[j].second))
-                                                                        /(2*this->sigma*this->sigma));
-      return k;
-    }
-    else{
-      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;
-      return d/this->approx;
-    }
-  }
-
-  /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel.
-   * \ingroup Persistence_weighted_gaussian
-   *
-   * @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(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);
-  }
-
-
-}; // 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 d8ed0d98..8c92ab54 100644
--- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
+++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
@@ -62,7 +62,7 @@ namespace Persistence_representations {
  * in the diagrams, or approximated 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$
- * 
+ *
  * For more details, please see \cite pmlr-v70-carriere17a .
  *
 **/
@@ -80,7 +80,7 @@ class Sliced_Wasserstein {
   void build_rep(){
 
     if(approx > 0){
-      
+
       double step = pi/this->approx;
       int n = diagram.size();
 
@@ -188,7 +188,7 @@ class Sliced_Wasserstein {
    * \ingroup Sliced_Wasserstein
    *
    * @pre       approx attribute needs to be the same for both instances.
-   * @param[in] second other instance of class Sliced_Wasserstein. 
+   * @param[in] second other instance of class Sliced_Wasserstein.
    *
    *
    */
diff --git a/src/Persistence_representations/include/gudhi/Weight_functions.h b/src/Persistence_representations/include/gudhi/Weight_functions.h
deleted file mode 100644
index 78de406d..00000000
--- a/src/Persistence_representations/include/gudhi/Weight_functions.h
+++ /dev/null
@@ -1,81 +0,0 @@
-/*    This file is part of the Gudhi Library. The Gudhi library
- *    (Geometric Understanding in Higher Dimensions) is a generic C++
- *    library for computational topology.
- *
- *    Author(s):       Mathieu Carriere
- *
- *    Copyright (C) 2018  INRIA (France)
- *
- *    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
- *    the Free Software Foundation, either version 3 of the License, or
- *    (at your option) any later version.
- *
- *    This program is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef WEIGHT_FUNCTIONS_H_
-#define WEIGHT_FUNCTIONS_H_
-
-// gudhi include
-#include <gudhi/read_persistence_from_file.h>
-#include <gudhi/common_persistence_representations.h>
-
-// standard include
-#include <cmath>
-#include <iostream>
-#include <vector>
-#include <limits>
-#include <fstream>
-#include <sstream>
-#include <algorithm>
-#include <string>
-#include <utility>
-#include <functional>
-
-namespace Gudhi {
-namespace Persistence_representations {
-
-/** \fn  static double pss_weight(std::pair<double,double> p)
- * \brief Persistence Scale Space kernel weight function.
- * \ingroup Persistence_representations
- *
- * @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;}
-
-/** \fn static double linear_weight(std::pair<double,double> p)
- * \brief Linear weight function.
- * \ingroup Persistence_representations
- *
- * @param[in] p point in 2D.
- */
-static double linear_weight(std::pair<double,double> p)  {return std::abs(p.second - p.first);}
-
-/** \fn static double const_weight(std::pair<double,double> p)
- * \brief Constant weight function.
- * \ingroup Persistence_representations
- *
- * @param[in] p point in 2D.
- */
-static double const_weight(std::pair<double,double> p)   {return 1;}
-
-/** \fn static std::function<double (std::pair<double,double>) > arctan_weight(double C, double alpha)
- * \brief Returns the arctan weight function with parameters C and alpha.
- * \ingroup Persistence_representations
- *
- * @param[in] C     positive constant.
- * @param[in] alpha positive power.
- */
-static std::function<double (std::pair<double,double>) > arctan_weight(double C, double alpha)  {return [=](std::pair<double,double> p){return C * atan(std::pow(std::abs(p.second - p.first), alpha));};}
-
-}  // namespace Persistence_representations
-}  // namespace Gudhi
-
-#endif  // WEIGHT_FUNCTIONS_H_
diff --git a/src/Persistence_representations/include/gudhi/common_persistence_representations.h b/src/Persistence_representations/include/gudhi/common_persistence_representations.h
index 539eee60..024c99ec 100644
--- a/src/Persistence_representations/include/gudhi/common_persistence_representations.h
+++ b/src/Persistence_representations/include/gudhi/common_persistence_representations.h
@@ -40,12 +40,23 @@ static constexpr double pi = boost::math::constants::pi<double>();
 /**
  * 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> >;
+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>) >;
+using Weight = std::function<double (std::pair<double, double>) >;
+using Kernel = std::function<double (std::pair<double, double>, std::pair<double, double> )>;
+
+Kernel Gaussian_kernel(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*sigma)   );};
+}
+
+Kernel polynomial_kernel(double c, double d){
+  return [=](std::pair<double, double> p, std::pair<double, double> q){return std::pow( p.first*q.first + p.second*q.second + c, d);};
+}
+
+
 // double epsi = std::numeric_limits<double>::epsilon();
 double epsi = 0.000005;