Persistence_heat_maps.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): Pawel Dlotko and Mathieu Carriere
4  *
5  * Modifications:
6  * - 2018/04 MC: Add discrete/non-discrete mechanism and non-discrete version
7  *
8  * Copyright (C) 2019 Inria
9  *
10  * Modification(s):
11  * - YYYY/MM Author: Description of the modification
12  */
13 
14 #ifndef PERSISTENCE_HEAT_MAPS_H_
15 #define PERSISTENCE_HEAT_MAPS_H_
16 
17 // gudhi include
18 #include <gudhi/read_persistence_from_file.h>
19 #include <gudhi/common_persistence_representations.h>
20 
21 // standard include
22 #include <vector>
23 #include <sstream>
24 #include <iostream>
25 #include <cmath>
26 #include <limits>
27 #include <algorithm>
28 #include <utility>
29 #include <string>
30 #include <functional>
31 
32 namespace Gudhi {
33 namespace Persistence_representations {
34 
39 std::vector<std::vector<double> > create_Gaussian_filter(size_t pixel_radius, double sigma) {
40  bool dbg = false;
41  // we are computing the kernel mask to 2 standard deviations away from the center. We discretize it in a grid of a
42  // size 2*pixel_radius times 2*pixel_radius.
43 
44  double r = 0;
45  double sigma_sqr = sigma * sigma;
46 
47  // sum is for normalization
48  double sum = 0;
49 
50  // initialization of a kernel:
51  std::vector<std::vector<double> > kernel(2 * pixel_radius + 1);
52  for (size_t i = 0; i != kernel.size(); ++i) {
53  std::vector<double> v(2 * pixel_radius + 1, 0);
54  kernel[i] = v;
55  }
56 
57  if (dbg) {
58  std::clog << "Kernel initialize \n";
59  std::clog << "pixel_radius : " << pixel_radius << std::endl;
60  std::clog << "kernel.size() : " << kernel.size() << std::endl;
61  getchar();
62  }
63 
64  for (int x = -pixel_radius; x <= static_cast<int>(pixel_radius); x++) {
65  for (int y = -pixel_radius; y <= static_cast<int>(pixel_radius); y++) {
66  double real_x = 2 * sigma * x / pixel_radius;
67  double real_y = 2 * sigma * y / pixel_radius;
68  r = std::sqrt(real_x * real_x + real_y * real_y);
69  kernel[x + pixel_radius][y + pixel_radius] = (exp(-(r * r) / sigma_sqr)) / (3.141592 * sigma_sqr);
70  sum += kernel[x + pixel_radius][y + pixel_radius];
71  }
72  }
73 
74  // normalize the kernel
75  for (size_t i = 0; i != kernel.size(); ++i) {
76  for (size_t j = 0; j != kernel[i].size(); ++j) {
77  kernel[i][j] /= sum;
78  }
79  }
80 
81  if (dbg) {
82  std::clog << "Here is the kernel : \n";
83  for (size_t i = 0; i != kernel.size(); ++i) {
84  for (size_t j = 0; j != kernel[i].size(); ++j) {
85  std::clog << kernel[i][j] << " ";
86  }
87  std::clog << std::endl;
88  }
89  }
90  return kernel;
91 }
92 
93 /*
94  * There are various options to scale the points depending on their location. One can for instance:
95  * (1) do nothing (scale all of them with the weight 1), as in the function constant_function
96  * (2) Scale them by the distance to the diagonal. This is implemented in function
97  * (3) Scale them with the square of their distance to diagonal. This is implemented in function
98  * (4) Scale them with
99  */
100 
107  public:
108  double operator()(const std::pair<double, double>& point_in_diagram) { return 1; }
109 };
110 
117  public:
118  double operator()(const std::pair<double, double>& point_in_diagram) {
119  // (point_in_diagram.first+point_in_diagram.second)/2.0
120  return std::sqrt(std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) +
121  std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2));
122  }
123 };
124 
131  public:
132  double operator()(const std::pair<double, double>& point_in_diagram) {
133  return std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) +
134  std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2);
135  }
136 };
137 
144  public:
145  double operator()(const std::pair<double, double>& point_in_diagram) {
146  return atan(point_in_diagram.second - point_in_diagram.first);
147  }
148 };
149 
158  public:
159  weight_by_setting_maximal_interval_to_have_length_one(double len) : length_of_maximal_interval(len) {}
160  double operator()(const std::pair<double, double>& point_in_diagram) {
161  return (point_in_diagram.second - point_in_diagram.first) / this->length_of_maximal_interval;
162  }
163 
164  private:
165  double length_of_maximal_interval;
166 };
167 
175 // This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances,
176 // Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product
177 template <typename Scalling_of_kernels = constant_scaling_function>
179  public:
185  Scalling_of_kernels f;
186  this->f = f;
187  this->erase_below_diagonal = false;
188  this->min_ = this->max_ = 0;
189  this->set_up_parameters_for_basic_classes();
190  }
191 
205  Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
206  std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1),
207  bool erase_below_diagonal = false, size_t number_of_pixels = 1000,
208  double min_ = std::numeric_limits<double>::max(),
209  double max_ = std::numeric_limits<double>::max());
210 
232  Persistence_heat_maps(const char* filename, std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1),
233  bool erase_below_diagonal = false, size_t number_of_pixels = 1000,
234  double min_ = std::numeric_limits<double>::max(),
235  double max_ = std::numeric_limits<double>::max(),
236  unsigned dimension = std::numeric_limits<unsigned>::max());
237 
242  Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
243  const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel,
244  size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x = 0, double max_x = 1,
245  double min_y = 0, double max_y = 1);
246 
250  Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval,
251  const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel);
252 
258  void compute_mean(const std::vector<Persistence_heat_maps*>& maps);
259 
265  void compute_median(const std::vector<Persistence_heat_maps*>& maps);
266 
270  void compute_percentage_of_active(const std::vector<Persistence_heat_maps*>& maps, size_t cutoff = 1);
271 
272  // put to file subroutine
278  void print_to_file(const char* filename) const;
279 
284  void load_from_file(const char* filename);
285 
289  inline bool check_if_the_same(const Persistence_heat_maps& second) const {
290  bool dbg = false;
291  if (this->heat_map.size() != second.heat_map.size()) {
292  if (dbg)
293  std::clog << "this->heat_map.size() : " << this->heat_map.size()
294  << " \n second.heat_map.size() : " << second.heat_map.size() << std::endl;
295  return false;
296  }
297  if (this->min_ != second.min_) {
298  if (dbg) std::clog << "this->min_ : " << this->min_ << ", second.min_ : " << second.min_ << std::endl;
299  return false;
300  }
301  if (this->max_ != second.max_) {
302  if (dbg) std::clog << "this->max_ : " << this->max_ << ", second.max_ : " << second.max_ << std::endl;
303  return false;
304  }
305  // in the other case we may assume that the persistence images are defined on the same domain.
306  return true;
307  }
308 
312  inline double get_min() const { return this->min_; }
313 
317  inline double get_max() const { return this->max_; }
318 
322  bool operator==(const Persistence_heat_maps& rhs) const {
323  bool dbg = false;
324  if (!this->check_if_the_same(rhs)) {
325  if (dbg) std::clog << "The domains are not the same \n";
326  return false; // in this case, the domains are not the same, so the maps cannot be the same.
327  }
328  for (size_t i = 0; i != this->heat_map.size(); ++i) {
329  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
330  if (!almost_equal(this->heat_map[i][j], rhs.heat_map[i][j])) {
331  if (dbg) {
332  std::clog << "this->heat_map[" << i << "][" << j << "] = " << this->heat_map[i][j] << std::endl;
333  std::clog << "rhs.heat_map[" << i << "][" << j << "] = " << rhs.heat_map[i][j] << std::endl;
334  }
335  return false;
336  }
337  }
338  }
339  return true;
340  }
341 
345  bool operator!=(const Persistence_heat_maps& rhs) const { return !((*this) == rhs); }
346 
350  void plot(const char* filename) const;
351 
352  template <typename Operation_type>
353  friend Persistence_heat_maps operation_on_pair_of_heat_maps(const Persistence_heat_maps& first,
354  const Persistence_heat_maps& second,
355  Operation_type operation) {
356  // first check if the heat maps are compatible
357  if (!first.check_if_the_same(second)) {
358  std::cerr << "Sizes of the heat maps are not compatible. The program will now terminate \n";
359  throw "Sizes of the heat maps are not compatible. The program will now terminate \n";
360  }
361  Persistence_heat_maps result;
362  result.min_ = first.min_;
363  result.max_ = first.max_;
364  result.heat_map.reserve(first.heat_map.size());
365  for (size_t i = 0; i != first.heat_map.size(); ++i) {
366  std::vector<double> v;
367  v.reserve(first.heat_map[i].size());
368  for (size_t j = 0; j != first.heat_map[i].size(); ++j) {
369  v.push_back(operation(first.heat_map[i][j], second.heat_map[i][j]));
370  }
371  result.heat_map.push_back(v);
372  }
373  return result;
374  } // operation_on_pair_of_heat_maps
375 
381  Persistence_heat_maps result;
382  result.min_ = this->min_;
383  result.max_ = this->max_;
384  result.heat_map.reserve(this->heat_map.size());
385  for (size_t i = 0; i != this->heat_map.size(); ++i) {
386  std::vector<double> v;
387  v.reserve(this->heat_map[i].size());
388  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
389  v.push_back(this->heat_map[i][j] * scalar);
390  }
391  result.heat_map.push_back(v);
392  }
393  return result;
394  }
395 
400  return operation_on_pair_of_heat_maps(first, second, std::plus<double>());
401  }
406  return operation_on_pair_of_heat_maps(first, second, std::minus<double>());
407  }
411  friend Persistence_heat_maps operator*(double scalar, const Persistence_heat_maps& A) {
412  return A.multiply_by_scalar(scalar);
413  }
417  friend Persistence_heat_maps operator*(const Persistence_heat_maps& A, double scalar) {
418  return A.multiply_by_scalar(scalar);
419  }
423  Persistence_heat_maps operator*(double scalar) { return this->multiply_by_scalar(scalar); }
428  *this = *this + rhs;
429  return *this;
430  }
435  *this = *this - rhs;
436  return *this;
437  }
442  *this = *this * x;
443  return *this;
444  }
449  if (x == 0) throw("In operator /=, division by 0. Program terminated.");
450  *this = *this * (1 / x);
451  return *this;
452  }
453 
454  // Implementations of functions for various concepts.
455 
460  std::vector<double> vectorize(int number_of_function) const;
465  size_t number_of_vectorize_functions() const { return this->number_of_functions_for_vectorization; }
466 
474  double project_to_R(int number_of_function) const;
479  size_t number_of_projections_to_R() const { return this->number_of_functions_for_projections_to_reals; }
480 
487  double distance(const Persistence_heat_maps& second_, double power = 1) const;
488 
493  void compute_average(const std::vector<Persistence_heat_maps*>& to_average);
494 
500  double compute_scalar_product(const Persistence_heat_maps& second_) const;
501 
502  // end of implementation of functions needed for concepts.
503 
507  std::pair<double, double> get_x_range() const { return std::make_pair(this->min_, this->max_); }
508 
512  std::pair<double, double> get_y_range() const { return this->get_x_range(); }
513 
514  protected:
515  // private methods
516  std::vector<std::vector<double> > check_and_initialize_maps(const std::vector<Persistence_heat_maps*>& maps);
517  size_t number_of_functions_for_vectorization;
518  size_t number_of_functions_for_projections_to_reals;
519  void construct(const std::vector<std::pair<double, double> >& intervals_,
520  std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1),
521  bool erase_below_diagonal = false, size_t number_of_pixels = 1000,
522  double min_ = std::numeric_limits<double>::max(), double max_ = std::numeric_limits<double>::max());
523 
524  void set_up_parameters_for_basic_classes() {
525  this->number_of_functions_for_vectorization = 1;
526  this->number_of_functions_for_projections_to_reals = 1;
527  }
528 
529  // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false)
530  bool discrete = true;
531  std::function<double(std::pair<double, double>, std::pair<double, double>)> kernel;
532  std::vector<std::pair<double, double> > interval;
533  std::vector<double> weights;
534 
535  // data
536  Scalling_of_kernels f;
537  bool erase_below_diagonal;
538  double min_;
539  double max_;
540  std::vector<std::vector<double> > heat_map;
541 };
542 
543 template <typename Scalling_of_kernels>
545  const std::vector<std::pair<double, double> >& interval,
546  const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel,
547  size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x, double max_x, double min_y, double max_y) {
548  this->discrete = true;
549  this->min_ = min_x;
550  this->max_ = max_x;
551  this->heat_map.resize(number_of_y_pixels);
552  double step_x = (max_x - min_x) / (number_of_x_pixels - 1);
553  double step_y = (max_y - min_y) / (number_of_y_pixels - 1);
554 
555  int num_pts = interval.size();
556  for (size_t i = 0; i < number_of_y_pixels; i++) {
557  double y = min_y + i * step_y;
558  this->heat_map[i].reserve(number_of_x_pixels);
559  for (size_t j = 0; j < number_of_x_pixels; j++) {
560  double x = min_x + j * step_x;
561  std::pair<double, double> grid_point(x, y);
562  double pixel_value = 0;
563  for (int k = 0; k < num_pts; k++) pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point);
564  this->heat_map[i].push_back(pixel_value);
565  }
566  }
567  this->set_up_parameters_for_basic_classes();
568 }
569 
570 template <typename Scalling_of_kernels>
572  const std::vector<std::pair<double, double> >& interval,
573  const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel) {
574  this->discrete = false;
575  this->interval = interval;
576  this->kernel = kernel;
577  int num_pts = this->interval.size();
578  for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->interval[i]));
579  this->set_up_parameters_for_basic_classes();
580 }
581 
582 // if min_ == max_, then the program is requested to set up the values itself based on persistence intervals
583 template <typename Scalling_of_kernels>
584 void Persistence_heat_maps<Scalling_of_kernels>::construct(const std::vector<std::pair<double, double> >& intervals_,
585  std::vector<std::vector<double> > filter,
586  bool erase_below_diagonal, size_t number_of_pixels,
587  double min_, double max_) {
588  bool dbg = false;
589  if (dbg) std::clog << "Entering construct procedure \n";
590  Scalling_of_kernels f;
591  this->f = f;
592 
593  if (dbg) std::clog << "min and max passed to construct() procedure: " << min_ << " " << max_ << std::endl;
594 
595  if (min_ == max_) {
596  if (dbg) std::clog << "min and max parameters will be determined based on intervals \n";
597  // in this case, we want the program to set up the min_ and max_ values by itself.
598  min_ = std::numeric_limits<int>::max();
599  max_ = -std::numeric_limits<int>::max();
600 
601  for (size_t i = 0; i != intervals_.size(); ++i) {
602  if (intervals_[i].first < min_) min_ = intervals_[i].first;
603  if (intervals_[i].second > max_) max_ = intervals_[i].second;
604  }
605  // now we have the structure filled in, and moreover we know min_ and max_ values of the interval, so we know the
606  // range.
607 
608  // add some more space:
609  min_ -= fabs(max_ - min_) / 100;
610  max_ += fabs(max_ - min_) / 100;
611  }
612 
613  if (dbg) {
614  std::clog << "min_ : " << min_ << std::endl;
615  std::clog << "max_ : " << max_ << std::endl;
616  std::clog << "number_of_pixels : " << number_of_pixels << std::endl;
617  getchar();
618  }
619 
620  this->min_ = min_;
621  this->max_ = max_;
622 
623  // initialization of the structure heat_map
624  std::vector<std::vector<double> > heat_map_;
625  for (size_t i = 0; i != number_of_pixels; ++i) {
626  std::vector<double> v(number_of_pixels, 0);
627  heat_map_.push_back(v);
628  }
629  this->heat_map = heat_map_;
630 
631  if (dbg) std::clog << "Done creating of the heat map, now we will fill in the structure \n";
632 
633  for (size_t pt_nr = 0; pt_nr != intervals_.size(); ++pt_nr) {
634  // compute the value of intervals_[pt_nr] in the grid:
635  int x_grid =
636  static_cast<int>((intervals_[pt_nr].first - this->min_) / (this->max_ - this->min_) * number_of_pixels);
637  int y_grid =
638  static_cast<int>((intervals_[pt_nr].second - this->min_) / (this->max_ - this->min_) * number_of_pixels);
639 
640  if (dbg) {
641  std::clog << "point : " << intervals_[pt_nr].first << " , " << intervals_[pt_nr].second << std::endl;
642  std::clog << "x_grid : " << x_grid << std::endl;
643  std::clog << "y_grid : " << y_grid << std::endl;
644  }
645 
646  // x_grid and y_grid gives a center of the kernel. We want to have its lower left corner. To get this, we need to
647  // shift x_grid and y_grid by a grid diameter.
648  x_grid -= filter.size() / 2;
649  y_grid -= filter.size() / 2;
650  // note that the numbers x_grid and y_grid may be negative.
651 
652  if (dbg) {
653  std::clog << "After shift : \n";
654  std::clog << "x_grid : " << x_grid << std::endl;
655  std::clog << "y_grid : " << y_grid << std::endl;
656  }
657 
658  double scaling_value = this->f(intervals_[pt_nr]);
659 
660  for (size_t i = 0; i != filter.size(); ++i) {
661  for (size_t j = 0; j != filter.size(); ++j) {
662  // if the point (x_grid+i,y_grid+j) is the correct point in the grid.
663  if (((x_grid + i) >= 0) && (x_grid + i < this->heat_map.size()) && ((y_grid + j) >= 0) &&
664  (y_grid + j < this->heat_map.size())) {
665  if (dbg) {
666  std::clog << y_grid + j << " " << x_grid + i << std::endl;
667  }
668  this->heat_map[y_grid + j][x_grid + i] += scaling_value * filter[i][j];
669  if (dbg) {
670  std::clog << "Position : (" << x_grid + i << "," << y_grid + j
671  << ") got increased by the value : " << filter[i][j] << std::endl;
672  }
673  }
674  }
675  }
676  }
677 
678  // now it remains to cut everything below diagonal if the user wants us to.
679  if (erase_below_diagonal) {
680  for (size_t i = 0; i != this->heat_map.size(); ++i) {
681  for (size_t j = i; j != this->heat_map.size(); ++j) {
682  this->heat_map[i][j] = 0;
683  }
684  }
685  }
686 } // construct
687 
688 template <typename Scalling_of_kernels>
690  const std::vector<std::pair<double, double> >& interval, std::vector<std::vector<double> > filter,
691  bool erase_below_diagonal, size_t number_of_pixels, double min_, double max_) {
692  this->construct(interval, filter, erase_below_diagonal, number_of_pixels, min_, max_);
693  this->set_up_parameters_for_basic_classes();
694 }
695 
696 template <typename Scalling_of_kernels>
698  std::vector<std::vector<double> > filter,
699  bool erase_below_diagonal, size_t number_of_pixels,
700  double min_, double max_, unsigned dimension) {
701  std::vector<std::pair<double, double> > intervals_;
702  if (dimension == std::numeric_limits<unsigned>::max()) {
703  intervals_ = read_persistence_intervals_in_one_dimension_from_file(filename);
704  } else {
705  intervals_ = read_persistence_intervals_in_one_dimension_from_file(filename, dimension);
706  }
707  this->construct(intervals_, filter, erase_below_diagonal, number_of_pixels, min_, max_);
708  this->set_up_parameters_for_basic_classes();
709 }
710 
711 template <typename Scalling_of_kernels>
713  const std::vector<Persistence_heat_maps*>& maps) {
714  // checking if all the heat maps are of the same size:
715  for (size_t i = 0; i != maps.size(); ++i) {
716  if (maps[i]->heat_map.size() != maps[0]->heat_map.size()) {
717  std::cerr << "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";
718  throw "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";
719  }
720  if (maps[i]->heat_map[0].size() != maps[0]->heat_map[0].size()) {
721  std::cerr << "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";
722  throw "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";
723  }
724  }
725  std::vector<std::vector<double> > heat_maps(maps[0]->heat_map.size());
726  for (size_t i = 0; i != maps[0]->heat_map.size(); ++i) {
727  std::vector<double> v(maps[0]->heat_map[0].size(), 0);
728  heat_maps[i] = v;
729  }
730  return heat_maps;
731 }
732 
733 template <typename Scalling_of_kernels>
734 void Persistence_heat_maps<Scalling_of_kernels>::compute_median(const std::vector<Persistence_heat_maps*>& maps) {
735  std::vector<std::vector<double> > heat_maps = this->check_and_initialize_maps(maps);
736 
737  std::vector<double> to_compute_median(maps.size());
738  for (size_t i = 0; i != heat_maps.size(); ++i) {
739  for (size_t j = 0; j != heat_maps[i].size(); ++j) {
740  for (size_t map_no = 0; map_no != maps.size(); ++map_no) {
741  to_compute_median[map_no] = maps[map_no]->heat_map[i][j];
742  }
743  std::nth_element(to_compute_median.begin(), to_compute_median.begin() + to_compute_median.size() / 2,
744  to_compute_median.end());
745  heat_maps[i][j] = to_compute_median[to_compute_median.size() / 2];
746  }
747  }
748  this->heat_map = heat_maps;
749  this->min_ = maps[0]->min_;
750  this->max_ = maps[0]->max_;
751 }
752 
753 template <typename Scalling_of_kernels>
754 void Persistence_heat_maps<Scalling_of_kernels>::compute_mean(const std::vector<Persistence_heat_maps*>& maps) {
755  std::vector<std::vector<double> > heat_maps = this->check_and_initialize_maps(maps);
756  for (size_t i = 0; i != heat_maps.size(); ++i) {
757  for (size_t j = 0; j != heat_maps[i].size(); ++j) {
758  double mean = 0;
759  for (size_t map_no = 0; map_no != maps.size(); ++map_no) {
760  mean += maps[map_no]->heat_map[i][j];
761  }
762  heat_maps[i][j] = mean / static_cast<double>(maps.size());
763  }
764  }
765  this->heat_map = heat_maps;
766  this->min_ = maps[0]->min_;
767  this->max_ = maps[0]->max_;
768 }
769 
770 template <typename Scalling_of_kernels>
772  const std::vector<Persistence_heat_maps*>& maps, size_t cutoff) {
773  std::vector<std::vector<double> > heat_maps = this->check_and_initialize_maps(maps);
774 
775  for (size_t i = 0; i != heat_maps.size(); ++i) {
776  for (size_t j = 0; j != heat_maps[i].size(); ++j) {
777  size_t number_of_active_levels = 0;
778  for (size_t map_no = 0; map_no != maps.size(); ++map_no) {
779  if (maps[map_no]->heat_map[i][j]) number_of_active_levels++;
780  }
781  if (number_of_active_levels > cutoff) {
782  heat_maps[i][j] = number_of_active_levels;
783  } else {
784  heat_maps[i][j] = 0;
785  }
786  }
787  }
788  this->heat_map = heat_maps;
789  this->min_ = maps[0]->min_;
790  this->max_ = maps[0]->max_;
791 }
792 
793 template <typename Scalling_of_kernels>
794 void Persistence_heat_maps<Scalling_of_kernels>::plot(const char* filename) const {
795  std::ofstream out;
796  std::stringstream gnuplot_script;
797  gnuplot_script << filename << "_GnuplotScript";
798 
799  out.open(gnuplot_script.str().c_str());
800  out << "plot '-' matrix with image" << std::endl;
801  for (size_t i = 0; i != this->heat_map.size(); ++i) {
802  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
803  out << this->heat_map[i][j] << " ";
804  }
805  out << std::endl;
806  }
807  out.close();
808  std::clog << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'"
809  << gnuplot_script.str().c_str() << "\'\"" << std::endl;
810 }
811 
812 template <typename Scalling_of_kernels>
814  std::ofstream out;
815  out.open(filename);
816 
817  // First we store this->min_ and this->max_ values:
818  out << this->min_ << " " << this->max_ << std::endl;
819  for (size_t i = 0; i != this->heat_map.size(); ++i) {
820  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
821  out << this->heat_map[i][j] << " ";
822  }
823  out << std::endl;
824  }
825  out.close();
826 }
827 
828 template <typename Scalling_of_kernels>
830  bool dbg = false;
831 
832  std::ifstream in;
833  in.open(filename);
834 
835  // checking if the file exist / if it was open.
836  if (!in.good()) {
837  std::cerr << "The file : " << filename << " do not exist. The program will now terminate \n";
838  throw "The persistence landscape file do not exist. The program will now terminate \n";
839  }
840 
841  // now we read the file one by one.
842 
843  in >> this->min_ >> this->max_;
844  if (dbg) {
845  std::clog << "Reading the following values of min and max : " << this->min_ << " , " << this->max_ << std::endl;
846  }
847 
848  std::string temp;
849  std::getline(in, temp);
850  while (in.good()) {
851  std::getline(in, temp);
852  std::stringstream lineSS;
853  lineSS << temp;
854 
855  std::vector<double> line_of_heat_map;
856  while (lineSS.good()) {
857  double point;
858 
859  lineSS >> point;
860  line_of_heat_map.push_back(point);
861  if (dbg) {
862  std::clog << point << " ";
863  }
864  }
865  if (dbg) {
866  std::clog << std::endl;
867  getchar();
868  }
869 
870  if (in.good()) this->heat_map.push_back(line_of_heat_map);
871  }
872  in.close();
873  if (dbg) std::clog << "Done \n";
874 }
875 
876 // Concretizations of virtual methods:
877 template <typename Scalling_of_kernels>
878 std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize(int number_of_function) const {
879  std::vector<double> result;
880  if (!discrete) {
881  std::cerr << "No vectorize method in case of infinite dimensional vectorization" << std::endl;
882  return result;
883  }
884 
885  // convert this->heat_map into one large vector:
886  size_t size_of_result = 0;
887  for (size_t i = 0; i != this->heat_map.size(); ++i) {
888  size_of_result += this->heat_map[i].size();
889  }
890 
891  result.reserve(size_of_result);
892 
893  for (size_t i = 0; i != this->heat_map.size(); ++i) {
894  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
895  result.push_back(this->heat_map[i][j]);
896  }
897  }
898 
899  return result;
900 }
901 
902 template <typename Scalling_of_kernels>
904  if (this->discrete) {
905  // first we need to check if (*this) and second are defined on the same domain and have the same dimensions:
906  if (!this->check_if_the_same(second)) {
907  std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between "
908  "them. The program will now terminate";
909  throw "The persistence images are of non compatible sizes. The program will now terminate";
910  }
911 
912  // if we are here, we know that the two persistence images are defined on the same domain, so we can start
913  // computing their distances:
914 
915  double distance = 0;
916  if (power < std::numeric_limits<double>::max()) {
917  for (size_t i = 0; i != this->heat_map.size(); ++i) {
918  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
919  distance += std::pow(std::fabs(this->heat_map[i][j] - second.heat_map[i][j]), power);
920  }
921  }
922  } else {
923  // in this case, we compute max norm distance
924  for (size_t i = 0; i != this->heat_map.size(); ++i) {
925  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
926  if (distance < std::fabs(this->heat_map[i][j] - second.heat_map[i][j])) {
927  distance = std::fabs(this->heat_map[i][j] - second.heat_map[i][j]);
928  }
929  }
930  }
931  }
932  return distance;
933  } else {
934  return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) -
935  2 * this->compute_scalar_product(second));
936  }
937 }
938 
939 template <typename Scalling_of_kernels>
940 double Persistence_heat_maps<Scalling_of_kernels>::project_to_R(int number_of_function) const {
941  double result = 0;
942  for (size_t i = 0; i != this->heat_map.size(); ++i) {
943  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
944  result += this->heat_map[i][j];
945  }
946  }
947  return result;
948 }
949 
950 template <typename Scalling_of_kernels>
952  const std::vector<Persistence_heat_maps*>& to_average) {
953  this->compute_mean(to_average);
954 }
955 
956 template <typename Scalling_of_kernels>
958  if (discrete) {
959  // first we need to check if (*this) and second are defined on the same domain and have the same dimensions:
960  if (!this->check_if_the_same(second)) {
961  std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between "
962  "them. The program will now terminate";
963  throw "The persistence images are of non compatible sizes. The program will now terminate";
964  }
965 
966  // if we are here, we know that the two persistence images are defined on the same domain, so we can start computing
967  // their scalar product:
968  double scalar_prod = 0;
969  for (size_t i = 0; i != this->heat_map.size(); ++i) {
970  for (size_t j = 0; j != this->heat_map[i].size(); ++j) {
971  scalar_prod += this->heat_map[i][j] * second.heat_map[i][j];
972  }
973  }
974  return scalar_prod;
975  } else {
976  int num_pts1 = this->interval.size();
977  int num_pts2 = second.interval.size();
978  double kernel_val = 0;
979  for (int i = 0; i < num_pts1; i++) {
980  std::pair<double, double> pi = this->interval[i];
981  for (int j = 0; j < num_pts2; j++) {
982  std::pair<double, double> pj = second.interval[j];
983  kernel_val += this->weights[i] * second.weights[j] * this->kernel(pi, pj);
984  }
985  }
986  return kernel_val;
987  }
988 }
989 
990 } // namespace Persistence_representations
991 } // namespace Gudhi
992 
993 #endif // PERSISTENCE_HEAT_MAPS_H_
A class implementing persistence heat maps.
Definition: Persistence_heat_maps.h:178
size_t number_of_vectorize_functions() const
Definition: Persistence_heat_maps.h:465
std::pair< double, double > get_x_range() const
Definition: Persistence_heat_maps.h:507
void compute_mean(const std::vector< Persistence_heat_maps * > &maps)
Definition: Persistence_heat_maps.h:754
bool check_if_the_same(const Persistence_heat_maps &second) const
Definition: Persistence_heat_maps.h:289
std::vector< double > vectorize(int number_of_function) const
Definition: Persistence_heat_maps.h:878
friend Persistence_heat_maps operator*(const Persistence_heat_maps &A, double scalar)
Definition: Persistence_heat_maps.h:417
friend Persistence_heat_maps operator*(double scalar, const Persistence_heat_maps &A)
Definition: Persistence_heat_maps.h:411
Persistence_heat_maps operator/=(double x)
Definition: Persistence_heat_maps.h:448
Persistence_heat_maps operator*(double scalar)
Definition: Persistence_heat_maps.h:423
bool operator!=(const Persistence_heat_maps &rhs) const
Definition: Persistence_heat_maps.h:345
Persistence_heat_maps operator-=(const Persistence_heat_maps &rhs)
Definition: Persistence_heat_maps.h:434
void print_to_file(const char *filename) const
Definition: Persistence_heat_maps.h:813
Persistence_heat_maps(const char *filename, 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(), unsigned dimension=std::numeric_limits< unsigned >::max())
Definition: Persistence_heat_maps.h:697
Persistence_heat_maps multiply_by_scalar(double scalar) const
Definition: Persistence_heat_maps.h:380
Persistence_heat_maps operator*=(double x)
Definition: Persistence_heat_maps.h:441
Persistence_heat_maps(const std::vector< std::pair< double, double > > &interval, 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())
Definition: Persistence_heat_maps.h:689
void compute_percentage_of_active(const std::vector< Persistence_heat_maps * > &maps, size_t cutoff=1)
Definition: Persistence_heat_maps.h:771
size_t number_of_projections_to_R() const
Definition: Persistence_heat_maps.h:479
bool operator==(const Persistence_heat_maps &rhs) const
Definition: Persistence_heat_maps.h:322
double project_to_R(int number_of_function) const
Definition: Persistence_heat_maps.h:940
friend Persistence_heat_maps operator-(const Persistence_heat_maps &first, const Persistence_heat_maps &second)
Definition: Persistence_heat_maps.h:405
Persistence_heat_maps()
Definition: Persistence_heat_maps.h:184
void compute_average(const std::vector< Persistence_heat_maps * > &to_average)
Definition: Persistence_heat_maps.h:951
Persistence_heat_maps(const std::vector< std::pair< double, double > > &interval, const std::function< double(std::pair< double, double >, std::pair< double, double >)> &kernel, 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)
Definition: Persistence_heat_maps.h:544
double get_min() const
Definition: Persistence_heat_maps.h:312
double get_max() const
Definition: Persistence_heat_maps.h:317
double distance(const Persistence_heat_maps &second_, double power=1) const
Definition: Persistence_heat_maps.h:903
void load_from_file(const char *filename)
Definition: Persistence_heat_maps.h:829
Persistence_heat_maps(const std::vector< std::pair< double, double > > &interval, const std::function< double(std::pair< double, double >, std::pair< double, double >)> &kernel)
Definition: Persistence_heat_maps.h:571
void compute_median(const std::vector< Persistence_heat_maps * > &maps)
Definition: Persistence_heat_maps.h:734
Persistence_heat_maps operator+=(const Persistence_heat_maps &rhs)
Definition: Persistence_heat_maps.h:427
void plot(const char *filename) const
Definition: Persistence_heat_maps.h:794
std::pair< double, double > get_y_range() const
Definition: Persistence_heat_maps.h:512
double compute_scalar_product(const Persistence_heat_maps &second_) const
Definition: Persistence_heat_maps.h:957
friend Persistence_heat_maps operator+(const Persistence_heat_maps &first, const Persistence_heat_maps &second)
Definition: Persistence_heat_maps.h:399
Gudhi namespace.
Definition: SimplicialComplexForAlpha.h:14