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

Former-commit-id: 12f32a1c8ca31e7e0a40e1c3502e2a3d810d5bc5

16 files changed, 217 insertions, 993 deletions

diff --git a/src/Persistence_representations/doc/Persistence_representations_doc.h b/src/Persistence_representations/doc/Persistence_representations_doc.h
index d0b02739..73800d0d 100644
--- a/src/Persistence_representations/doc/Persistence_representations_doc.h
+++ b/src/Persistence_representations/doc/Persistence_representations_doc.h
@@ -24,6 +24,7 @@
 #define DOC_GUDHI_STAT_H_
 
 namespace Gudhi {
+
 namespace Persistence_representations {
 
 /**  \defgroup Persistence_representations Persistence representations
@@ -128,35 +129,33 @@ namespace Persistence_representations {
  function \f$L : \mathbb{N} \times \mathbb{R} \to [0,\infty)\f$ of two
  variables, if we define \f$L(k,t) = \lambda_k(t)\f$.
 
- The detailed description of algorithms used to compute persistence landscapes can be found in \cite bubenik_dlotko_landscapes_2016.
- Note that this implementation provides exact representation of landscapes. That have many advantages, but also a few drawbacks. 
- For instance, as discussed in \cite bubenik_dlotko_landscapes_2016, the exact representation of landscape may be of quadratic size with respect
- to the input persistence diagram. It may therefore happen that, for very large diagrams, using this representation may be memory--prohibitive. 
- In such a case, there are two possible ways to proceed:
+ The detailed description of algorithms used to compute persistence landscapes can be found in
+ \cite bubenik_dlotko_landscapes_2016.
+ Note that this implementation provides exact representation of landscapes. That have many advantages, but also a few
+ drawbacks. For instance, as discussed
+ in \cite bubenik_dlotko_landscapes_2016, the exact representation of landscape may be of quadratic size with respect
+ to the input persistence diagram. It may therefore happen
+ that, for very large diagrams, using this representation may be memory--prohibitive. In such a case, there are two
+ possible ways to proceed:
 
- \li Use representation on a grid---see section \ref sec_landscapes_on_grid.
+ \li Use non exact representation on a grid described in the Section \ref sec_landscapes_on_grid.
  \li Compute just a number of initial nonzero landscapes. This option is available from C++ level as a last parameter of
  the constructor of persistence landscape (set by default to std::numeric_limits<size_t>::max()).
 
 
 
  \section sec_landscapes_on_grid Persistence Landscapes on a grid
-
  <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_landscape_on_grid <br>
- <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_landscape_on_grid_exact <br>
-
- Here, we provide alternative, not exact, representations of persistence landscapes defined in Section \ref sec_persistence_landscapes. 
- Unlike Section \ref sec_persistence_landscapes, we build representations of persistence landscapes by evaluating the landscape functions on a finite, equally distributed grid of points.
- We propose two different representations depending on whether the persistence intervals are also mapped on the grid (Persistence_landscape_on_grid) or not (Persistence_landscape_on_grid_exact).
- This makes a big difference since mapping the intervals on the grid makes the computation time smaller but only provides an approximation of the landscape values.
+ This is an alternative, not--exact, representation of persistence landscapes defined in the Section \ref
+ sec_persistence_landscapes. Unlike in the Section \ref sec_persistence_landscapes we build a
+ representation of persistence landscape by sampling its values on a finite, equally distributed grid of points.
+ Since, the persistence landscapes that originate from persistence diagrams have slope \f$1\f$ or \f$-1\f$, we have an
+ estimate of a region between the grid points where the landscape cab be located.
+ That allows to estimate an error make when performing various operations on landscape. Note that for average
+ landscapes the slope is in range \f$[-1,1]\f$ and similar estimate can be used.
 
- Since persistence landscapes originating from persistence diagrams have slope \f$1\f$ or \f$-1\f$, we have an
- estimate of a region between the grid points where the landscapes can be located.
- That allows to estimate an error made when performing various operations on landscapes. Note that for average
- landscapes the slope is in range \f$[-1,1]\f$ and similar estimates can be used.
-
- Due to the lack of rigorous description of the algorithms for these non rigorous representations of persistence
- landscapes in the literature, we provide a short discussion below.
+ Due to a lack of rigorous description of the algorithms to deal with this non--rigorous representation of persistence
+ landscapes in the literature, we are providing a short discussion of them in below.
 
  Let us assume that we want to compute persistence landscape on a interval \f$[x,y]\f$. Let us assume that we want to
  use \f$N\f$ grid points for that purpose.
@@ -168,11 +167,11 @@ namespace Persistence_representations {
  functions) on the i-th point of a grid, i.e. \f$x + i \frac{y-x}{N}\f$.
 
  When averaging two persistence landscapes represented by a grid we need to make sure that they are defined in a
- compatible grids, i.e. the intervals \f$[x,y]\f$ on which they are defined are
+ compatible grids. I.e. the intervals \f$[x,y]\f$ on which they are defined are
  the same, and the numbers of grid points \f$N\f$ are the same in both cases. If this is the case, we simply compute
- point-wise averages of the entries of the corresponding
- vectors (in this whole section we assume that if one vector of numbers is shorter than the other, we extend the shortest
- one with zeros so that they have the same length).
+ point-wise averages of the entries of corresponding
+ vectors (In this whole section we assume that if one vector of numbers is shorter than another, we extend the shorter
+ one with zeros so that they have the same length.)
 
  Computations of distances between two persistence landscapes on a grid is not much different than in the rigorous
  case. In this case, we sum up the distances between the same levels of
@@ -181,11 +180,11 @@ namespace Persistence_representations {
 
  Similarly as in case of distance, when computing the scalar product of two persistence landscapes on a grid, we sum up
  the scalar products of corresponding levels of landscapes. For each level,
- we assume that the persistence landscape on a grid between two grid points is approximated by a linear function.
- Therefore to compute the scalar product of two corresponding levels of landscapes,
+ we assume that the persistence landscape on a grid between two grid points is approximated by linear function.
+ Therefore to compute scalar product of two corresponding levels of landscapes,
  we sum up the integrals of products of line segments for every pair of constitutive grid points.
 
- Note that for these representations we need to specify a few parameters:
+ Note that for this representation we need to specify a few parameters:
 
  \li Begin and end point of a grid -- the interval \f$[x,y]\f$ (real numbers).
  \li Number of points in a grid (positive integer \f$N\f$).
@@ -194,33 +193,29 @@ namespace Persistence_representations {
  Note that the same representation is used in TDA R-package \cite Fasy_Kim_Lecci_Maria_tda.
 
  \section sec_persistence_heat_maps Persistence heat maps
-
  <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_heat_maps <br>
- <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_heat_maps_exact <br>
-
- This is a general class of discrete structures which are based on idea of placing a kernel in the points of persistence diagrams.
+ This is a general class of discrete structures which are based on idea of placing a kernel in the points of
+ persistence diagrams.
  This idea appeared in work by many authors over the last 15 years. As far as we know this idea was firstly described
  in the work of Bologna group in \cite Ferri_Frosini_comparision_sheme_1 and \cite Ferri_Frosini_comparision_sheme_2.
  Later it has been described by Colorado State University group in \cite Persistence_Images_2017. The presented paper
- in the first time provided a discussion of stability of this representation.
- Also, the same ideas are used in the construction of two recent kernels used for machine learning:
- \cite Kusano_Fukumizu_Hiraoka_PWGK and \cite Reininghaus_Huber_ALL_PSSK. Both the kernels use
- interesting ideas to ensure stability of the representations with respect to the 1-Wasserstein metric. In the kernel
+ in the first time provide a discussion of stability of the representation.
+ Also, the same ideas are used in construction of two recent kernels used for machine learning:
+ \cite Kusano_Fukumizu_Hiraoka_PWGK and \cite Reininghaus_Huber_ALL_PSSK. Both the kernel's construction uses
+ interesting ideas to ensure stability of the representation with respect to Wasserstein metric. In the kernel
  presented in \cite Kusano_Fukumizu_Hiraoka_PWGK, a scaling function is used to multiply the Gaussian kernel in the
- way that the points close to diagonal have low weights and consequently do not have a big influence on the resulting
+ way that the points close to diagonal got low weight and consequently do not have a big influence on the resulting
  distribution. In \cite Reininghaus_Huber_ALL_PSSK for every point \f$(b,d)\f$ two Gaussian kernels
  are added: first, with a weight 1 in a point \f$(b,d)\f$, and the second, with the weight -1 for a point \f$(b,d)\f$.
  In both cases, the representations are stable with respect to 1-Wasserstein distance.
 
- In Persistence_representations package, we currently implement a discretization of the distributions described above.
- The base of this implementation is a 2-dimensional array of pixels. To each pixel is assigned a real value which
- is the sum of the distribution values induced by each point of the persistence diagram. 
- As for Persistence_landscapes, we propose two different representations depending on whether the persistence intervals are also mapped on the pixels 
- (Persistence_heat_maps) or not (Persistence_heat_maps_exact).
- At the moment we compute the sum over the evaluations of the distributions on the pixel centers. It can be easily extended to any other function
- (like for instance the sum of the integrals of the distributions over the pixels).
+ In Persistence\_representations package we currently implement a discretization of the distributions described above.
+ The base of this implementation is 2-dimensional array of pixels. Each pixel have assigned a real value which
+ is a sum of values of distributions induced by each point of the persistence diagram. At the moment we compute the
+ sum of values on a center of a pixels. It can be easily extended to any other function
+ (like for instance sum of integrals of the intermediate distribution on a pixel).
 
- Concerning Persistence_heat_maps, the parameters that determine the structure are the following:
+ The parameters that determine the structure are the following:
 
  \li A positive integer k determining the size of the kernel we used (we always assume that the kernels are square).
  \li A filter: in practice a square matrix of a size \f$2k+1 \times 2k+1\f$. By default, this is a discretization of
@@ -232,7 +227,6 @@ namespace Persistence_representations {
  to diagonal are given then sometimes the kernel have support that reaches the region
  below the diagonal. If the value of this parameter is true, then the values below diagonal can be erased.
 
- Concerning Persistence_heat_maps_exact, only Gaussian kernels are implemented, so the parameters are the array of pixels, the weight functions for the Gaussians and the bandwidth of the Gaussians. 
 
  \section sec_persistence_vectors Persistence vectors
  <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Vector_distances_in_diagram <br>
@@ -256,11 +250,7 @@ namespace Persistence_representations {
  absolute value of differences between coordinates. A scalar product is a sum of products of
  values at the corresponding positions of two vectors.
 
-
-
-
-
-\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>
 
@@ -269,53 +259,26 @@ namespace Persistence_representations {
  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>.
  Kernels are designed for algorithms that can be <i>kernelized</i>, i.e., algorithms that only require to know scalar products between instances in order to run.
- Examples of such algorithms include Support Vector Machines, Principal Component Analysis and Ridge Regression. 
+ Examples of such algorithms include Support Vector Machines, Principal Component Analysis and Ridge Regression.
 
  There have been several attempts at defining kernels, i.e., positive semi-definite functions, between persistence diagrams within the last few years. We provide implementation
- for three of them:
-
- \li the <i>Persistence Scale Space Kernel</i>---see \cite Reininghaus_Huber_ALL_PSSK, which is the classical scalar product between \f$L^2\f$ functions, where persistence diagrams
- are turned into functions by centering and summing Gaussian functions over the diagram points and their symmetric counterparts w.r.t. the diagonal: \f$k(D_1,D_2)=\int \Phi(D_1)\Phi(D_2)\f$,
- where \f$\Phi(D)=\sum_{p\in D} {\rm exp}\left(-\frac{\|p-\cdot\|_2^2}{2\sigma^2}\right)\f$. 
- 
- \li the <i>Persistence Weighted Gaussian Kernel</i>---see \cite Kusano_Fukumizu_Hiraoka_PWGK, which is a slight generalization of the previous kernel, is the scalar product between
- weighted Kernel Mean Embeddings of persistence diagrams w.r.t. the Gaussian Kernel \f$k_G\f$ (with corresponding map \f$\Phi_G\f$) in \f$\mathbb{R}^2\f$: 
- \f$k(D_1,D_2)=\langle\sum_{p\in D_1} w(p)\Phi_G(p), \sum_{q\in D_2} w(q)\Phi_G(q)\rangle\f$
- 
- \li the  <i>Sliced Wasserstein Kernel</i>---see \cite pmlr-v70-carriere17a, which takes the form of a Gaussian kernel with a specific distance between persistence diagrams
-     called the <i>Sliced Wasserstein Distance</i>: \f$k(D_1,D_2)={\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right)\f$
+ for the <i>Sliced Wasserstein Kernel</i>---see \cite pmlr-v70-carriere17a, which takes the form of a Gaussian kernel with a specific distance between persistence diagrams
+ called the <i>Sliced Wasserstein Distance</i>: \f$k(D_1,D_2)={\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right)\f$. Other kernels such as the Persistence Weighted Gaussian Kernel or
+ the Persistence Scale Space Kernel are implemented in Persistence_heat_maps.
 
  When launching:
 
  \code $> ./Sliced_Wasserstein
  \endcode
- 
+
  the program output is:
- 
+
  \code $> Approx SW distance: 5.33648
  $> Exact SW distance: 5.33798
  $> Approx SW kernel: 0.0693743
  $> Exact  SW kernel: 0.0693224
  $> Distance induced by approx SW kernel: 1.36428
- $> Distance induced by exact  SW kernel: 1.3643 
- \endcode
-
-
- and when launching:
-
- \code $> ./Persistence_weighted_gaussian
- \endcode
- 
- the program output is:
- 
- \code $> Approx PWG kernel: 1.21509
- $> Exact  PWG kernel: 1.13628
- $> Distance induced by approx PWG kernel: 3.23354
- $> Distance induced by exact  PWG kernel: 3.25697
- $> Approx Gaussian PWG kernel: 0.0194222
- $> Exact  Gaussian PWG kernel: 0.0192524
- $> Approx PSS kernel: 0.134413
- $> Exact  PSS kernel: 0.133394
+ $> Distance induced by exact  SW kernel: 1.3643
  \endcode
 
  */
diff --git a/src/Persistence_representations/example/CMakeLists.txt b/src/Persistence_representations/example/CMakeLists.txt
index 9be22085..d236c3a6 100644
--- a/src/Persistence_representations/example/CMakeLists.txt
+++ b/src/Persistence_representations/example/CMakeLists.txt
@@ -31,23 +31,3 @@ add_executable ( Sliced_Wasserstein sliced_wasserstein.cpp )
 add_test(NAME Sliced_Wasserstein
     COMMAND $<TARGET_FILE:Sliced_Wasserstein>)
 install(TARGETS Sliced_Wasserstein DESTINATION bin)
-
-add_executable ( Persistence_weighted_gaussian persistence_weighted_gaussian.cpp )
-add_test(NAME Persistence_weighted_gaussian
-    COMMAND $<TARGET_FILE:Persistence_weighted_gaussian>)
-install(TARGETS Persistence_weighted_gaussian DESTINATION bin)
-
-add_executable ( Persistence_heat_maps_exact persistence_heat_maps_exact.cpp )
-add_test(NAME Persistence_heat_maps_exact
-    COMMAND $<TARGET_FILE:Persistence_heat_maps_exact>)
-install(TARGETS Persistence_heat_maps_exact DESTINATION bin)
-
-add_executable ( Persistence_landscape_on_grid_exact persistence_landscape_on_grid_exact.cpp )
-add_test(NAME Persistence_landscape_on_grid_exact
-    COMMAND $<TARGET_FILE:Persistence_landscape_on_grid_exact>)
-install(TARGETS Persistence_landscape_on_grid_exact DESTINATION bin)
-
-add_executable ( Betti_sequence betti_sequence.cpp )
-add_test(NAME Betti_sequence
-    COMMAND $<TARGET_FILE:Betti_sequence>)
-install(TARGETS Betti_sequence DESTINATION bin)
diff --git a/src/Persistence_representations/example/betti_sequence.cpp b/src/Persistence_representations/example/betti_sequence.cpp
deleted file mode 100644
index a422a822..00000000
--- a/src/Persistence_representations/example/betti_sequence.cpp
+++ /dev/null
@@ -1,49 +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/>.
- */
-
-#include <gudhi/Betti_sequence.h>
-
-#include <iostream>
-#include <vector>
-#include <utility>
-
-using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram;
-using BS = Gudhi::Persistence_representations::Betti_sequence;
-
-int main(int argc, char** argv) {
-
-  Persistence_diagram persistence;
-
-  persistence.push_back(std::make_pair(1, 2));
-  persistence.push_back(std::make_pair(6, 8));
-  persistence.push_back(std::make_pair(0, 4));
-  persistence.push_back(std::make_pair(3, 8));
-
-  double min_x = 0; double max_x = 8; int res_x = 1000;
-
-  BS bs(persistence, min_x, max_x, res_x);
-  std::vector<int> B = bs.vectorize();
-
-  for(int i = 0; i < res_x; i++)  std::cout << B[i] << ", ";
-
-  return 0;
-}
diff --git a/src/Persistence_representations/example/persistence_heat_maps.cpp b/src/Persistence_representations/example/persistence_heat_maps.cpp
index 323b57e9..f1791e97 100644
--- a/src/Persistence_representations/example/persistence_heat_maps.cpp
+++ b/src/Persistence_representations/example/persistence_heat_maps.cpp
@@ -21,6 +21,7 @@
  */
 
 #include <gudhi/Persistence_heat_maps.h>
+#include <gudhi/common_persistence_representations.h>
 
 #include <iostream>
 #include <vector>
@@ -76,5 +77,16 @@ 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;
 
+  Gudhi::Persistence_representations::Kernel k = Gudhi::Persistence_representations::Gaussian_kernel(1.0);
+
+  Persistence_heat_maps hm1k(persistence1, k);
+  Persistence_heat_maps hm2k(persistence2, k);
+
+  Persistence_heat_maps hm1i(persistence1, 20, 20, 0, 11, 0, 11, k);
+  Persistence_heat_maps hm2i(persistence2, 20, 20, 0, 11, 0, 11, k);
+
+  std::cout << "Scalar product computed with exact kernel is : " << hm1i.compute_scalar_product(hm2i) << std::endl;
+  std::cout << "Kernel value between PDs seen as functions is : " << hm1k.compute_scalar_product(hm2k) << std::endl;
+
   return 0;
 }
diff --git a/src/Persistence_representations/example/persistence_heat_maps_exact.cpp b/src/Persistence_representations/example/persistence_heat_maps_exact.cpp
deleted file mode 100644
index f15b710d..00000000
--- a/src/Persistence_representations/example/persistence_heat_maps_exact.cpp
+++ /dev/null
@@ -1,55 +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/>.
- */
-
-#include <gudhi/Persistence_heat_maps_exact.h>
-#include <gudhi/Weight_functions.h>
-
-#include <iostream>
-#include <vector>
-#include <utility>
-#include <string>
-
-using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram;
-using PI = Gudhi::Persistence_representations::Persistence_heat_maps_exact;
-using Weight = std::function<double (std::pair<double,double>) >;
-
-int main(int argc, char** argv) {
-
-  Persistence_diagram persistence;
-
-  persistence.push_back(std::make_pair(1, 2));
-  persistence.push_back(std::make_pair(6, 8));
-  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::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();
-
-  for(int i = 0; i < res_y; i++){
-    for(int j = 0; j < res_x; j++) std::cout << P[i][j] << " ";
-    std::cout << std::endl;
-  }
-
-  return 0;
-}
diff --git a/src/Persistence_representations/example/persistence_landscape_on_grid_exact.cpp b/src/Persistence_representations/example/persistence_landscape_on_grid_exact.cpp
deleted file mode 100644
index 9ce42649..00000000
--- a/src/Persistence_representations/example/persistence_landscape_on_grid_exact.cpp
+++ /dev/null
@@ -1,52 +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/>.
- */
-
-#include <gudhi/Persistence_landscape_on_grid_exact.h>
-
-#include <iostream>
-#include <vector>
-#include <utility>
-
-using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram;
-using LS = Gudhi::Persistence_representations::Persistence_landscape_on_grid_exact;
-
-int main(int argc, char** argv) {
-
-  Persistence_diagram persistence;
-
-  persistence.push_back(std::make_pair(1, 2));
-  persistence.push_back(std::make_pair(6, 8));
-  persistence.push_back(std::make_pair(0, 4));
-  persistence.push_back(std::make_pair(3, 8));
-
-  int nb_ls = 2; double min_x = 0; double max_x = 8; int res_x = 1000;
-
-  LS ls(persistence, nb_ls, min_x, max_x, res_x);
-  std::vector<std::vector<double> > L = ls.vectorize();
-
-  for(int i = 0; i < nb_ls; i++){
-    for(int j = 0; j < res_x; j++) std::cout << L[i][j] << ", ";
-    std::cout << std::endl;
-  }
-
-  return 0;
-}
diff --git a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp
deleted file mode 100644
index 7945e4f1..00000000
--- a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp
+++ /dev/null
@@ -1,99 +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/>.
- */
-
-#include <gudhi/Persistence_weighted_gaussian.h>
-
-#include <iostream>
-#include <vector>
-#include <utility>
-
-using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram;
-using PWG = Gudhi::Persistence_representations::Persistence_weighted_gaussian;
-
-int main(int argc, char** argv) {
-
-  Persistence_diagram persistence1, persistence2;
-
-  persistence1.push_back(std::make_pair(1, 2));
-  persistence1.push_back(std::make_pair(6, 8));
-  persistence1.push_back(std::make_pair(0, 4));
-  persistence1.push_back(std::make_pair(3, 8));
-
-  persistence2.push_back(std::make_pair(2, 9));
-  persistence2.push_back(std::make_pair(1, 6));
-  persistence2.push_back(std::make_pair(3, 5));
-  persistence2.push_back(std::make_pair(6, 10));
-
-  double sigma = 1;
-  double tau = 1;
-  int m = 10000;
-
-  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, Gudhi::Persistence_representations::arctan_weight(1,1));
-  PWG PWGex2(persistence2, sigma, -1, Gudhi::Persistence_representations::arctan_weight(1,1));
-
-
-  // Linear PWG
-
-  std::cout << "Approx PWG kernel: " << PWG1.compute_scalar_product  (PWG2) << std::endl;
-  std::cout << "Exact  PWG kernel: " << PWGex1.compute_scalar_product  (PWGex2) << std::endl;
-
-  std::cout << "Distance induced by approx PWG kernel: " << PWG1.distance  (PWG2) << std::endl;
-  std::cout << "Distance induced by exact  PWG kernel: " << PWGex1.distance  (PWGex2) << std::endl;
-
-
-
-
-
-
-
-  // Gaussian PWG
-
-  std::cout << "Approx Gaussian PWG kernel: " << std::exp( -PWG1.distance (PWG2) ) / (2*tau*tau) << std::endl;
-  std::cout << "Exact  Gaussian PWG kernel: " << std::exp( -PWGex1.distance (PWGex2) ) / (2*tau*tau) << std::endl;
-
-
-
-
-
-
-
-  // PSS
-
-  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, 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, 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;
-
-
-
-  return 0;
-}
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;