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

Former-commit-id: 189ac5572f69842e1d8d1cec68ca6a4f62e39bd4

6 files changed, 54 insertions, 80 deletions

diff --git a/src/Persistence_representations/doc/Persistence_representations_doc.h b/src/Persistence_representations/doc/Persistence_representations_doc.h
index 73800d0d..3aa99315 100644
--- a/src/Persistence_representations/doc/Persistence_representations_doc.h
+++ b/src/Persistence_representations/doc/Persistence_representations_doc.h
@@ -227,6 +227,9 @@ 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.
 
+ We also provide two methods to perform exact calculations. In both methods, the kernel is no longer provided as a filter (i.e. a square matrix---see parameters above), but rather as
+ a function assigning a real value to a 2D point. In the first method, calculations are done on a grid, leading to a finite-dimensional representation.
+ On the other hand, in the second method, we do not use grids, meaning that diagrams are represented as functions. Thus, only scalar products are available.
 
  \section sec_persistence_vectors Persistence vectors
  <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Vector_distances_in_diagram <br>
@@ -273,10 +276,8 @@ namespace Persistence_representations {
 
  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
+ \code $> Approx SW kernel: 0.0693743
+ $> Exact  SW kernel: 0.0693218
  $> Distance induced by approx SW kernel: 1.36428
  $> Distance induced by exact  SW kernel: 1.3643
  \endcode
diff --git a/src/Persistence_representations/example/persistence_heat_maps.cpp b/src/Persistence_representations/example/persistence_heat_maps.cpp
index f1791e97..5874336e 100644
--- a/src/Persistence_representations/example/persistence_heat_maps.cpp
+++ b/src/Persistence_representations/example/persistence_heat_maps.cpp
@@ -77,7 +77,7 @@ 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);
+  Gudhi::Persistence_representations::Kernel2D k = Gudhi::Persistence_representations::Gaussian_kernel(1.0);
 
   Persistence_heat_maps hm1k(persistence1, k);
   Persistence_heat_maps hm2k(persistence2, k);
diff --git a/src/Persistence_representations/example/sliced_wasserstein.cpp b/src/Persistence_representations/example/sliced_wasserstein.cpp
index 2104e2b2..089172a0 100644
--- a/src/Persistence_representations/example/sliced_wasserstein.cpp
+++ b/src/Persistence_representations/example/sliced_wasserstein.cpp
@@ -50,8 +50,6 @@ int main(int argc, char** argv) {
   SW swex1(persistence1, 1, -1);
   SW swex2(persistence2, 1, -1);
 
-  std::cout << "Approx SW distance: " << sw1.compute_sliced_wasserstein_distance(sw2) << std::endl;
-  std::cout << "Exact SW distance: " << swex1.compute_sliced_wasserstein_distance(swex2) << std::endl;
   std::cout << "Approx SW kernel: " << sw1.compute_scalar_product(sw2) << std::endl;
   std::cout << "Exact  SW kernel: " << swex1.compute_scalar_product(swex2) << std::endl;
   std::cout << "Distance induced by approx SW kernel: " << sw1.distance(sw2) << std::endl;
diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
index 63c6e239..12188526 100644
--- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
+++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
@@ -249,14 +249,14 @@ class Persistence_heat_maps {
    * 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));
+                        double min_x = 0, double max_x = 1, double min_y = 0, double max_y = 1, const Kernel2D & 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));
+  Persistence_heat_maps(const Persistence_diagram & interval, const Kernel2D & 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
@@ -531,8 +531,8 @@ class Persistence_heat_maps {
 
   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));
+                                                   double min_x = 0, double max_x = 1, double min_y = 0, double max_y = 1, const Kernel2D & kernel = Gaussian_kernel(1.0));
+  void construct_kernel_from_exact_universal_kernel(const Persistence_diagram & interval, const Kernel2D & kernel = Gaussian_kernel(1.0));
 
   void set_up_parameters_for_basic_classes() {
     this->number_of_functions_for_vectorization = 1;
@@ -543,7 +543,7 @@ class Persistence_heat_maps {
   bool discrete = true;
 
   // PWGK
-  Kernel k;
+  Kernel2D k;
   Persistence_diagram d;
   std::vector<double> weights;
 
@@ -559,23 +559,22 @@ 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) {
+                                                                                             double min_y, double max_y, const Kernel2D & 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();
+  this->heat_map.resize(number_of_y_pixels);
   double step_x = (max_x - min_x)/(number_of_x_pixels - 1); double step_y = (max_y - min_y)/(number_of_y_pixels - 1);
 
   int num_pts = diagram.size();
 
   for(size_t i = 0; i < number_of_y_pixels; i++){
-    double y = min_y + i*step_y;
+    double y = min_y + i*step_y; this->heat_map[i].reserve(number_of_x_pixels);
     for(size_t j = 0; j < number_of_x_pixels; j++){
       double x = min_x + j*step_x;
 
       std::pair<double, double> grid_point(x,y); double pixel_value = 0;
       for(int k = 0; k < num_pts; k++){
-        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);
+        pixel_value += this->f(diagram[k]) * kernel(diagram[k], grid_point);
       }
       this->heat_map[i].push_back(pixel_value);
 
@@ -589,13 +588,13 @@ 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) {
+                                                                  double min_y, double max_y, const Kernel2D & 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){
+void Persistence_heat_maps<Scalling_of_kernels>::construct_kernel_from_exact_universal_kernel(const Persistence_diagram & diagram, const Kernel2D & 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]));
@@ -603,7 +602,7 @@ void Persistence_heat_maps<Scalling_of_kernels>::construct_kernel_from_exact_uni
 
 
 template <typename Scalling_of_kernels>
-Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(const Persistence_diagram& diagram, const Kernel & kernel) {
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps(const Persistence_diagram& diagram, const Kernel2D & kernel) {
   this->construct_kernel_from_exact_universal_kernel(diagram, kernel);
   this->set_up_parameters_for_basic_classes();
 }
diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
index 8c92ab54..a3c0dc2f 100644
--- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
+++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
@@ -28,18 +28,6 @@
 #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 {
 
@@ -70,12 +58,15 @@ namespace Persistence_representations {
 class Sliced_Wasserstein {
 
  protected:
-    Persistence_diagram diagram;
-    int approx;
-    double sigma;
-    std::vector<std::vector<double> > projections, projections_diagonal;
 
- public:
+  Persistence_diagram diagram;
+  int approx;
+  double sigma;
+  std::vector<std::vector<double> > projections, projections_diagonal;
+
+  // **********************************
+  // Utils.
+  // **********************************
 
   void build_rep(){
 
@@ -96,28 +87,14 @@ class Sliced_Wasserstein {
         }
 
         std::sort(l.begin(), l.end()); std::sort(l_diag.begin(), l_diag.end());
-        projections.push_back(l); projections_diagonal.push_back(l_diag);
+        projections.push_back(std::move(l)); projections_diagonal.push_back(std::move(l_diag));
 
       }
 
+      diagram.clear();
     }
-
   }
 
-  /** \brief Sliced Wasserstein kernel constructor.
-   * \ingroup Sliced_Wasserstein
-   *
-   * @param[in] _diagram  persistence diagram.
-   * @param[in] _sigma  bandwidth parameter.
-   * @param[in] _approx number of directions used to approximate the integral in the Sliced Wasserstein distance, set to -1 for exact computation.
-   *
-   */
-  Sliced_Wasserstein(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 100){diagram = _diagram; approx = _approx; sigma = _sigma; build_rep();}
-
-  // **********************************
-  // Utils.
-  // **********************************
-
   // Compute the angle formed by two points of a PD
   double compute_angle(const Persistence_diagram & diag, int i, int j) const {
     std::pair<double,double> vect; double x1,y1, x2,y2;
@@ -177,21 +154,7 @@ class Sliced_Wasserstein {
     return norm*integral;
   }
 
-
-
-
-  // **********************************
-  // Scalar product + distance.
-  // **********************************
-
-  /** \brief Evaluation of the Sliced Wasserstein Distance between a pair of diagrams.
-   * \ingroup Sliced_Wasserstein
-   *
-   * @pre       approx attribute needs to be the same for both instances.
-   * @param[in] second other instance of class Sliced_Wasserstein.
-   *
-   *
-   */
+  // Evaluation of the Sliced Wasserstein Distance between a pair of diagrams.
   double compute_sliced_wasserstein_distance(const Sliced_Wasserstein & second) const {
 
     GUDHI_CHECK(this->approx != second.approx, std::invalid_argument("Error: different approx values for representations"));
@@ -232,8 +195,8 @@ class Sliced_Wasserstein {
       }
 
       // Sort angles.
-      std::sort(angles1.begin(), angles1.end(), [=](const std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);});
-      std::sort(angles2.begin(), angles2.end(), [=](const std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);});
+      std::sort(angles1.begin(), angles1.end(), [](const std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);});
+      std::sort(angles2.begin(), angles2.end(), [](const std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);});
 
       // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2).
       std::vector<int> orderp1, orderp2;
@@ -291,11 +254,13 @@ class Sliced_Wasserstein {
 
     else{
 
-      double step = pi/this->approx;
+      double step = pi/this->approx; std::vector<double> v1, v2;
       for (int i = 0; i < this->approx; i++){
 
-        std::vector<double> v1; std::vector<double> l1 = this->projections[i]; std::vector<double> l1bis = second.projections_diagonal[i]; std::merge(l1.begin(), l1.end(), l1bis.begin(), l1bis.end(), std::back_inserter(v1));
-        std::vector<double> v2; std::vector<double> l2 = second.projections[i]; std::vector<double> l2bis = this->projections_diagonal[i]; std::merge(l2.begin(), l2.end(), l2bis.begin(), l2bis.end(), std::back_inserter(v2));
+        v1.clear(); v2.clear();
+        std::merge(this->projections[i].begin(), this->projections[i].end(), second.projections_diagonal[i].begin(), second.projections_diagonal[i].end(), std::back_inserter(v1));
+        std::merge(second.projections[i].begin(), second.projections[i].end(), this->projections_diagonal[i].begin(), this->projections_diagonal[i].end(), std::back_inserter(v2));
+
         int n = v1.size(); double f = 0;
         for (int j = 0; j < n; j++)  f += std::abs(v1[j] - v2[j]);
         sw += f*step;
@@ -306,6 +271,20 @@ class Sliced_Wasserstein {
     return sw/pi;
   }
 
+ public:
+
+  /** \brief Sliced Wasserstein kernel constructor.
+   * \ingroup Sliced_Wasserstein
+   *
+   * @param[in] _diagram  persistence diagram.
+   * @param[in] _sigma    bandwidth parameter.
+   * @param[in] _approx   number of directions used to approximate the integral in the Sliced Wasserstein distance, set to -1 for exact computation. If positive, then projections of the diagram
+   *                      points on all directions are stored in memory to reduce computation time.
+   *
+   */
+  // This class implements the following concepts: Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_scalar_product
+  Sliced_Wasserstein(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 10):diagram(_diagram), approx(_approx), sigma(_sigma) {build_rep();}
+
   /** \brief Evaluation of the kernel on a pair of diagrams.
    * \ingroup Sliced_Wasserstein
    *
@@ -331,8 +310,6 @@ class Sliced_Wasserstein {
   }
 
 
-
-
 }; // class Sliced_Wasserstein
 }  // namespace Persistence_representations
 }  // namespace Gudhi
diff --git a/src/Persistence_representations/include/gudhi/common_persistence_representations.h b/src/Persistence_representations/include/gudhi/common_persistence_representations.h
index 66ed3bf8..f0cd8146 100644
--- a/src/Persistence_representations/include/gudhi/common_persistence_representations.h
+++ b/src/Persistence_representations/include/gudhi/common_persistence_representations.h
@@ -43,16 +43,15 @@ static constexpr double pi = boost::math::constants::pi<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.
+ * In this module, we use the name Kernel for the representation of a function taking a pair of two points in the plane and returning a double.
  */
-using Weight = std::function<double (std::pair<double, double>) >;
-using Kernel = std::function<double (std::pair<double, double>, std::pair<double, double> )>;
+using Kernel2D = std::function<double (std::pair<double, double>, std::pair<double, double> )>;
 
-Kernel Gaussian_kernel(double sigma){
+inline Kernel2D Gaussian_kernel(double sigma){
   return [=](std::pair<double, double> p, std::pair<double, double> q){return (1.0 / (std::sqrt(2*pi)*sigma)) * std::exp(  -((p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second)) / (2*sigma*sigma)   );};
 }
 
-Kernel polynomial_kernel(double c, double d){
+inline Kernel2D 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);};
 }