diff options
Diffstat (limited to 'src/Persistence_representations/include/gudhi')
3 files changed, 48 insertions, 73 deletions
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);}; } |