1 files changed, 138 insertions, 36 deletions

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