From ef407f0e40099a832f13371401b44ac565325aff Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 15 Mar 2019 15:09:32 +0100 Subject: Remove comments and identify both contributions. Include what is used and add some std::. Use std::sqrt instead of std::pow(x,0.5) --- .../example/persistence_heat_maps.cpp | 26 +- .../include/gudhi/Persistence_heat_maps.h | 279 ++++++++------- .../include/gudhi/Persistence_landscape_on_grid.h | 6 +- .../include/gudhi/Sliced_Wasserstein.h | 373 +++++++++++++-------- 4 files changed, 381 insertions(+), 303 deletions(-) (limited to 'src/Persistence_representations') diff --git a/src/Persistence_representations/example/persistence_heat_maps.cpp b/src/Persistence_representations/example/persistence_heat_maps.cpp index 0beb1052..45208b68 100644 --- a/src/Persistence_representations/example/persistence_heat_maps.cpp +++ b/src/Persistence_representations/example/persistence_heat_maps.cpp @@ -2,9 +2,12 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Pawel Dlotko + * Author(s): Pawel Dlotko and Mathieu Carriere * - * Copyright (C) 2016 Inria + * Copyright (C) 2019 Inria + * + * Modifications: + * - 2018/04 MC: Add persistence heat maps computation * * 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 @@ -26,10 +29,14 @@ #include #include #include - - -std::function, std::pair)> Gaussian_function(double sigma){ - return [=](std::pair p, std::pair q){return std::exp( -( (p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second) )/(sigma) );}; +#include +#include + +std::function, std::pair)> Gaussian_function(double sigma) { + return [=](std::pair p, std::pair q) { + return std::exp(-((p.first - q.first) * (p.first - q.first) + (p.second - q.second) * (p.second - q.second)) / + (sigma)); + }; } using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; @@ -82,14 +89,13 @@ 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; - // Mathieu's code ************************************************************************************************************ Persistence_heat_maps hm1k(persistence1, Gaussian_function(1.0)); Persistence_heat_maps hm2k(persistence2, Gaussian_function(1.0)); Persistence_heat_maps hm1i(persistence1, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11); Persistence_heat_maps hm2i(persistence2, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11); - std::cout << "Scalar product computed with exact 2D kernel on grid is : " << hm1i.compute_scalar_product(hm2i) << std::endl; - std::cout << "Scalar product computed with exact 2D kernel is : " << hm1k.compute_scalar_product(hm2k) << std::endl; - // *************************************************************************************************************************** + std::cout << "Scalar product computed with exact 2D kernel on grid is : " << hm1i.compute_scalar_product(hm2i) + << std::endl; + std::cout << "Scalar product computed with exact 2D kernel is : " << hm1k.compute_scalar_product(hm2k) << std::endl; return 0; } diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h index db51cc14..a8458bda 100644 --- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h +++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h @@ -2,9 +2,12 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Pawel Dlotko + * Author(s): Pawel Dlotko and Mathieu Carriere * - * Copyright (C) 2016 Inria + * Modifications: + * - 2018/04 MC: Add discrete/non-discrete mechanism and non-discrete version + * + * Copyright (C) 2019 Inria * * 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 @@ -44,7 +47,7 @@ namespace Persistence_representations { /** * This is a simple procedure to create n by n (or 2*pixel_radius times 2*pixel_radius cubical approximation of a *Gaussian kernel. -**/ + **/ std::vector > create_Gaussian_filter(size_t pixel_radius, double sigma) { bool dbg = false; // we are computing the kernel mask to 2 standard deviations away from the center. We discretize it in a grid of a @@ -74,7 +77,7 @@ std::vector > create_Gaussian_filter(size_t pixel_radius, do for (int y = -pixel_radius; y <= static_cast(pixel_radius); y++) { double real_x = 2 * sigma * x / pixel_radius; double real_y = 2 * sigma * y / pixel_radius; - r = sqrt(real_x * real_x + real_y * real_y); + r = std::sqrt(real_x * real_x + real_y * real_y); kernel[x + pixel_radius][y + pixel_radius] = (exp(-(r * r) / sigma_sqr)) / (3.141592 * sigma_sqr); sum += kernel[x + pixel_radius][y + pixel_radius]; } @@ -100,18 +103,18 @@ std::vector > create_Gaussian_filter(size_t pixel_radius, do } /* -* There are various options to scale the points depending on their location. One can for instance: -* (1) do nothing (scale all of them with the weight 1), as in the function constant_function -* (2) Scale them by the distance to the diagonal. This is implemented in function -* (3) Scale them with the square of their distance to diagonal. This is implemented in function -* (4) Scale them with -*/ + * There are various options to scale the points depending on their location. One can for instance: + * (1) do nothing (scale all of them with the weight 1), as in the function constant_function + * (2) Scale them by the distance to the diagonal. This is implemented in function + * (3) Scale them with the square of their distance to diagonal. This is implemented in function + * (4) Scale them with + */ /** * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * This particular functionality is a function which always assign value 1 to a point in the diagram. -**/ + **/ class constant_scaling_function { public: double operator()(const std::pair& point_in_diagram) { return 1; } @@ -121,13 +124,13 @@ class constant_scaling_function { * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * The scaling given by this function to a point (b,d) is Euclidean distance of (b,d) from diagonal. -**/ + **/ class distance_from_diagonal_scaling { public: double operator()(const std::pair& point_in_diagram) { // (point_in_diagram.first+point_in_diagram.second)/2.0 - return sqrt(pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + - pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2)); + return std::sqrt(std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + + std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2)); } }; @@ -135,12 +138,12 @@ class distance_from_diagonal_scaling { * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * The scaling given by this function to a point (b,d) is a square of Euclidean distance of (b,d) from diagonal. -**/ + **/ class squared_distance_from_diagonal_scaling { public: double operator()(const std::pair& point_in_diagram) { - return pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + - pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2); + return std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + + std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2); } }; @@ -148,7 +151,7 @@ class squared_distance_from_diagonal_scaling { * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * The scaling given by this function to a point (b,d) is an arctan of a persistence of a point (i.e. arctan( b-d ). -**/ + **/ class arc_tan_of_persistence_of_point { public: double operator()(const std::pair& point_in_diagram) { @@ -162,32 +165,24 @@ class arc_tan_of_persistence_of_point { * This scaling function do not only depend on a point (p,d) in the diagram, but it depends on the whole diagram. * The longest persistence pair get a scaling 1. Any other pair get a scaling belong to [0,1], which is proportional * to the persistence of that pair. -**/ + **/ class weight_by_setting_maximal_interval_to_have_length_one { public: - weight_by_setting_maximal_interval_to_have_length_one(double len) : letngth_of_maximal_interval(len) {} + weight_by_setting_maximal_interval_to_have_length_one(double len) : length_of_maximal_interval(len) {} double operator()(const std::pair& point_in_diagram) { - return (point_in_diagram.second - point_in_diagram.first) / this->letngth_of_maximal_interval; + return (point_in_diagram.second - point_in_diagram.first) / this->length_of_maximal_interval; } private: - double letngth_of_maximal_interval; + double length_of_maximal_interval; }; -// Mathieu's code *********************************************************************************************************************** -//class Gaussian_kernel_2D { -// public: -// double operator()(const std::pair& p, const std::pair& q) const { -// return (1.0/(std::sqrt(2*pi))) * std::exp(-((p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second)) / 2); } -//}; -// ************************************************************************************************************************************** - /** * \class Persistence_heat_maps Persistence_heat_maps.h gudhi/Persistence_heat_maps.h * \brief A class implementing persistence heat maps. * * \ingroup Persistence_representations -**/ + **/ // This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, // Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product @@ -197,7 +192,7 @@ class Persistence_heat_maps { /** * The default constructor. A scaling function from the diagonal is set up to a constant function. The image is not *erased below the diagonal. The Gaussian have diameter 5. - **/ + **/ Persistence_heat_maps() { Scalling_of_kernels f; this->f = f; @@ -218,7 +213,7 @@ class Persistence_heat_maps { *std::numeric_limits::max(), in which case the program compute the values based on the data, * (6) a max x and y value of points that are to be taken into account. By default it is set to *std::numeric_limits::max(), in which case the program compute the values based on the data. - **/ + **/ Persistence_heat_maps(const std::vector >& interval, std::vector > filter = create_Gaussian_filter(5, 1), bool erase_below_diagonal = false, size_t number_of_pixels = 1000, @@ -226,12 +221,12 @@ class Persistence_heat_maps { double max_ = std::numeric_limits::max()); /** - * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by - *default), a scaling function (constant by default), a boolean value which determines if the area of image below - *diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each - *direction (set to 1000 by default) and min and max values of images (both set to std::numeric_limits::max() - *by default. If this is the case, the program will pick the right values based on the data). - **/ + * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by + *default), a scaling function (constant by default), a boolean value which determines if the area of image below + *diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each + *direction (set to 1000 by default) and min and max values of images (both set to std::numeric_limits::max() + *by default. If this is the case, the program will pick the right values based on the data). + **/ /** * Construction that takes at the input the following parameters: * (1) A name of a file with persistence intervals. The file should be readable by the function @@ -245,48 +240,45 @@ class Persistence_heat_maps { *std::numeric_limits::max(), in which case the program compute the values based on the data, * (6) a max x and y value of points that are to be taken into account. By default it is set to *std::numeric_limits::max(), in which case the program compute the values based on the data. - **/ + **/ Persistence_heat_maps(const char* filename, std::vector > filter = create_Gaussian_filter(5, 1), bool erase_below_diagonal = false, size_t number_of_pixels = 1000, double min_ = std::numeric_limits::max(), double max_ = std::numeric_limits::max(), unsigned dimension = std::numeric_limits::max()); -// Mathieu's code *********************************************************************************************************************** /** - * Construction that takes as inputs (1) the diagram, and (2) the grid parameters (min, max and number of samples for x and y axes) - **/ + * Construction that takes as inputs (1) the diagram, and (2) the grid parameters (min, max and number of samples for + * x and y axes) + **/ Persistence_heat_maps(const std::vector >& interval, const std::function, std::pair)>& 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); + 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); /** * Construction that takes only the diagram as input (weight and 2D kernel are template parameters) - **/ + **/ Persistence_heat_maps(const std::vector >& interval, const std::function, std::pair)>& kernel); -// ************************************************************************************************************************************** - - /** * 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. - **/ + **/ void compute_mean(const std::vector& maps); /** * Compute a median 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. - **/ + **/ void compute_median(const std::vector& maps); /** * Compute a percentage of active (i.e) values above the cutoff of a collection of heat maps. - **/ + **/ void compute_percentage_of_active(const std::vector& maps, size_t cutoff = 1); // put to file subroutine @@ -294,18 +286,18 @@ class Persistence_heat_maps { * The function outputs the persistence image to a text file. The format as follow: * In the first line, the values min and max of the image are stored * In the next lines, we have the persistence images in a form of a bitmap image. - **/ + **/ void print_to_file(const char* filename) const; /** * A function that load a heat map from file to the current object (and erase whatever was stored in the current *object before). - **/ + **/ void load_from_file(const char* filename); /** * The procedure checks if min_, max_ and this->heat_maps sizes are the same. - **/ + **/ inline bool check_if_the_same(const Persistence_heat_maps& second) const { bool dbg = false; if (this->heat_map.size() != second.heat_map.size()) { @@ -328,17 +320,17 @@ class Persistence_heat_maps { /** * Return minimal range value of persistent image. - **/ + **/ inline double get_min() const { return this->min_; } /** * Return maximal range value of persistent image. - **/ + **/ inline double get_max() const { return this->max_; } /** * Operator == to check if to persistence heat maps are the same. - **/ + **/ bool operator==(const Persistence_heat_maps& rhs) const { bool dbg = false; if (!this->check_if_the_same(rhs)) { @@ -361,12 +353,12 @@ class Persistence_heat_maps { /** * Operator != to check if to persistence heat maps are different. - **/ + **/ bool operator!=(const Persistence_heat_maps& rhs) const { return !((*this) == rhs); } /** * A function to generate a gnuplot script to visualize the persistent image. - **/ + **/ void plot(const char* filename) const; template @@ -396,7 +388,7 @@ class Persistence_heat_maps { /** * Multiplication of Persistence_heat_maps by scalar (so that all values of the heat map gets multiplied by that *scalar). - **/ + **/ Persistence_heat_maps multiply_by_scalar(double scalar) const { Persistence_heat_maps result; result.min_ = this->min_; @@ -415,56 +407,56 @@ class Persistence_heat_maps { /** * This function computes a sum of two objects of a type Persistence_heat_maps. - **/ + **/ friend Persistence_heat_maps operator+(const Persistence_heat_maps& first, const Persistence_heat_maps& second) { return operation_on_pair_of_heat_maps(first, second, std::plus()); } /** -* This function computes a difference of two objects of a type Persistence_heat_maps. -**/ + * This function computes a difference of two objects of a type Persistence_heat_maps. + **/ friend Persistence_heat_maps operator-(const Persistence_heat_maps& first, const Persistence_heat_maps& second) { return operation_on_pair_of_heat_maps(first, second, std::minus()); } /** -* This function computes a product of an object of a type Persistence_heat_maps with real number. -**/ + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ friend Persistence_heat_maps operator*(double scalar, const Persistence_heat_maps& A) { return A.multiply_by_scalar(scalar); } /** -* This function computes a product of an object of a type Persistence_heat_maps with real number. -**/ + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ friend Persistence_heat_maps operator*(const Persistence_heat_maps& A, double scalar) { return A.multiply_by_scalar(scalar); } /** -* This function computes a product of an object of a type Persistence_heat_maps with real number. -**/ + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ Persistence_heat_maps operator*(double scalar) { return this->multiply_by_scalar(scalar); } /** * += operator for Persistence_heat_maps. - **/ + **/ Persistence_heat_maps operator+=(const Persistence_heat_maps& rhs) { *this = *this + rhs; return *this; } /** - * -= operator for Persistence_heat_maps. - **/ + * -= operator for Persistence_heat_maps. + **/ Persistence_heat_maps operator-=(const Persistence_heat_maps& rhs) { *this = *this - rhs; return *this; } /** - * *= operator for Persistence_heat_maps. - **/ + * *= operator for Persistence_heat_maps. + **/ Persistence_heat_maps operator*=(double x) { *this = *this * x; return *this; } /** - * /= operator for Persistence_heat_maps. - **/ + * /= operator for Persistence_heat_maps. + **/ Persistence_heat_maps operator/=(double x) { if (x == 0) throw("In operator /=, division by 0. Program terminated."); *this = *this * (1 / x); @@ -474,14 +466,14 @@ class Persistence_heat_maps { // Implementations of functions for various concepts. /** - * This function produce a vector of doubles based on a persistence heat map. It is required in a concept + * This function produce a vector of doubles based on a persistence heat map. It is required in a concept * Vectorized_topological_data - */ + */ std::vector vectorize(int number_of_function) const; /** - * This function return the number of functions that allows vectorization of persistence heat map. It is required - *in a concept Vectorized_topological_data. - **/ + * This function return the number of functions that allows vectorization of persistence heat map. It is required + *in a concept Vectorized_topological_data. + **/ size_t number_of_vectorize_functions() const { return this->number_of_functions_for_vectorization; } /** @@ -490,45 +482,45 @@ class Persistence_heat_maps { * At the moment this function is not tested, since it is quite likely to be changed in the future. Given this, when *using it, keep in mind that it * will be most likely changed in the next versions. - **/ + **/ double project_to_R(int number_of_function) const; /** * The function gives the number of possible projections to R. This function is required by the *Real_valued_topological_data concept. - **/ + **/ size_t number_of_projections_to_R() const { return this->number_of_functions_for_projections_to_reals; } /** - * A function to compute distance between persistence heat maps. - * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. - * This function is required in Topological_data_with_distances concept. -* For max norm distance, set power to std::numeric_limits::max() - **/ + * A function to compute distance between persistence heat maps. + * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. + * This function is required in Topological_data_with_distances concept. + * For max norm distance, set power to std::numeric_limits::max() + **/ double distance(const Persistence_heat_maps& second_, double power = 1) const; /** * A function to compute averaged persistence heat map, based on vector of persistence heat maps. * This function is required by Topological_data_with_averages concept. - **/ + **/ void compute_average(const std::vector& to_average); /** - * A function to compute scalar product of persistence heat maps. - * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. - * This function is required in Topological_data_with_scalar_product concept. - **/ + * A function to compute scalar product of persistence heat maps. + * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. + * This function is required in Topological_data_with_scalar_product concept. + **/ double compute_scalar_product(const Persistence_heat_maps& second_) const; // end of implementation of functions needed for concepts. /** * The x-range of the persistence heat map. - **/ + **/ std::pair get_x_range() const { return std::make_pair(this->min_, this->max_); } /** * The y-range of the persistence heat map. - **/ + **/ std::pair get_y_range() const { return this->get_x_range(); } protected: @@ -546,12 +538,11 @@ class Persistence_heat_maps { this->number_of_functions_for_projections_to_reals = 1; } -// Mathieu's code *********************************************************************************************************************** - bool discrete = true; // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false) - std::function, std::pair)> kernel; + // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false) + bool discrete = true; + std::function, std::pair)> kernel; std::vector > interval; std::vector weights; -// ************************************************************************************************************************************** // data Scalling_of_kernels f; @@ -561,27 +552,27 @@ class Persistence_heat_maps { std::vector > heat_map; }; - - -// Mathieu's code *********************************************************************************************************************** template -Persistence_heat_maps::Persistence_heat_maps(const std::vector >& interval, - const std::function, std::pair)>& kernel, - size_t number_of_x_pixels, size_t number_of_y_pixels, - double min_x, double max_x, - double min_y, double max_y) { - - this->discrete = true; this->min_ = min_x; this->max_ = max_x; +Persistence_heat_maps::Persistence_heat_maps( + const std::vector >& interval, + const std::function, std::pair)>& kernel, + size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x, double max_x, double min_y, double max_y) { + this->discrete = true; + this->min_ = min_x; + this->max_ = max_x; 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); + 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 = interval.size(); - for(size_t i = 0; i < number_of_y_pixels; i++){ - 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 grid_point(x,y); double pixel_value = 0; - for(int k = 0; k < num_pts; k++) pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point); + for (size_t i = 0; i < number_of_y_pixels; i++) { + 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 grid_point(x, y); + double pixel_value = 0; + for (int k = 0; k < num_pts; k++) pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point); this->heat_map[i].push_back(pixel_value); } } @@ -589,16 +580,16 @@ Persistence_heat_maps::Persistence_heat_maps(const std::vec } template -Persistence_heat_maps::Persistence_heat_maps(const std::vector >& interval, - const std::function, std::pair)>& kernel) { - this->discrete = false; this->interval = interval; this->kernel = kernel; +Persistence_heat_maps::Persistence_heat_maps( + const std::vector >& interval, + const std::function, std::pair)>& kernel) { + this->discrete = false; + this->interval = interval; + this->kernel = kernel; int num_pts = this->interval.size(); - for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->interval[i])); + for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->interval[i])); 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 @@ -897,9 +888,11 @@ void Persistence_heat_maps::load_from_file(const char* file // Concretizations of virtual methods: template std::vector Persistence_heat_maps::vectorize(int number_of_function) const { - std::vector result; - if(!discrete){std::cout << "No vectorize method in case of infinite dimensional vectorization" << std::endl; return 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; @@ -920,7 +913,7 @@ std::vector Persistence_heat_maps::vectorize(int nu template double Persistence_heat_maps::distance(const Persistence_heat_maps& second, double power) const { - if(this->discrete){ + 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 " @@ -928,30 +921,30 @@ double Persistence_heat_maps::distance(const Persistence_he 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::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); + distance += std::pow(std::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]); + if (distance < std::fabs(this->heat_map[i][j] - second.heat_map[i][j])) { + distance = std::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 std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) - + 2 * this->compute_scalar_product(second)); } } @@ -974,8 +967,7 @@ void Persistence_heat_maps::compute_average( template double Persistence_heat_maps::compute_scalar_product(const Persistence_heat_maps& second) const { - - if(discrete){ + 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 " @@ -992,22 +984,19 @@ double Persistence_heat_maps::compute_scalar_product(const } } return scalar_prod; - } - -// Mathieu's code *********************************************************************************************************************** - else{ - int num_pts1 = this->interval.size(); int num_pts2 = second.interval.size(); double kernel_val = 0; - for(int i = 0; i < num_pts1; i++){ + } else { + int num_pts1 = this->interval.size(); + int num_pts2 = second.interval.size(); + double kernel_val = 0; + for (int i = 0; i < num_pts1; i++) { std::pair pi = this->interval[i]; - for(int j = 0; j < num_pts2; j++){ + for (int j = 0; j < num_pts2; j++) { std::pair pj = second.interval[j]; kernel_val += this->weights[i] * second.weights[j] * this->kernel(pi, pj); } } return kernel_val; } -// ************************************************************************************************************************************** - } } // namespace Persistence_representations 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 db0e362a..fd8a181c 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 * (p[int_no].first + p[int_no].second) - grid_min + 1); + size_t grid_interval_midpoint = (size_t)(0.5 * (grid_interval_begin + grid_interval_end)); 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 = grid_min + dx * (grid_interval_begin + 1) - p[int_no].first; + double landscape_value = dx; 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,8 +1030,6 @@ 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::max()) { diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h index 5d0f4a5d..18165c5f 100644 --- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carriere * - * Copyright (C) 2018 INRIA (France) + * Copyright (C) 2018 Inria * * 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 @@ -28,6 +28,13 @@ #include #include +#include // for std::vector<> +#include // for std::pair<>, std::move +#include // for std::sort, std::max, std::merge +#include // for std::abs, std::sqrt +#include // for std::invalid_argument +#include // for std::random_device + namespace Gudhi { namespace Persistence_representations { @@ -39,31 +46,33 @@ namespace Persistence_representations { * * \details * In this class, we compute infinite-dimensional representations of persistence diagrams by using the - * Sliced Wasserstein kernel (see \ref sec_persistence_kernels for more details on kernels). We recall that infinite-dimensional - * representations are defined implicitly, so only scalar products and distances are available for the representations defined in this class. - * The Sliced Wasserstein kernel is defined as a Gaussian-like kernel between persistence diagrams, where the distance used for - * comparison is the Sliced Wasserstein distance \f$SW\f$ between persistence diagrams, defined as the integral of the 1-norm - * between the sorted projections of the diagrams onto all lines passing through the origin: + * Sliced Wasserstein kernel (see \ref sec_persistence_kernels for more details on kernels). We recall that + * infinite-dimensional representations are defined implicitly, so only scalar products and distances are available for + * the representations defined in this class. + * The Sliced Wasserstein kernel is defined as a Gaussian-like kernel between persistence diagrams, where the distance + * used for comparison is the Sliced Wasserstein distance \f$SW\f$ between persistence diagrams, defined as the + * integral of the 1-norm between the sorted projections of the diagrams onto all lines passing through the origin: * - * \f$ SW(D_1,D_2)=\int_{\theta\in\mathbb{S}}\,\|\pi_\theta(D_1\cup\pi_\Delta(D_2))-\pi_\theta(D_2\cup\pi_\Delta(D_1))\|_1{\rm d}\theta\f$, + * \f$ SW(D_1,D_2)=\int_{\theta\in\mathbb{S}}\,\|\pi_\theta(D_1\cup\pi_\Delta(D_2))-\pi_\theta(D_2\cup\pi_\Delta(D_1))\ + * |_1{\rm d}\theta\f$, * - * where \f$\pi_\theta\f$ is the projection onto the line defined with angle \f$\theta\f$ in the unit circle \f$\mathbb{S}\f$, - * and \f$\pi_\Delta\f$ is the projection onto the diagonal. - * Assuming that the diagrams are in general position (i.e. there is no collinear triple), the integral can be computed exactly in \f$O(n^2{\rm log}(n))\f$ time, where \f$n\f$ is the number of points - * in the diagrams. We provide two approximations of the integral: one in which we slightly perturb the diagram points so that they are in general position, - * and another in which we approximate the integral 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: + * where \f$\pi_\theta\f$ is the projection onto the line defined with angle \f$\theta\f$ in the unit circle + * \f$\mathbb{S}\f$, and \f$\pi_\Delta\f$ is the projection onto the diagonal. + * Assuming that the diagrams are in general position (i.e. there is no collinear triple), the integral can be computed + * exactly in \f$O(n^2{\rm log}(n))\f$ time, where \f$n\f$ is the number of points in the diagrams. We provide two + * approximations of the integral: one in which we slightly perturb the diagram points so that they are in general + * position, and another in which we approximate the integral 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$ * * The first method is usually much more accurate but also * much slower. For more details, please see \cite pmlr-v70-carriere17a . * -**/ + **/ class Sliced_Wasserstein { - protected: - Persistence_diagram diagram; int approx; double sigma; @@ -73,27 +82,26 @@ class Sliced_Wasserstein { // Utils. // ********************************** - void build_rep(){ - - if(approx > 0){ - - double step = pi/this->approx; + void build_rep() { + if (approx > 0) { + double step = pi / this->approx; int n = diagram.size(); - for (int i = 0; i < this->approx; i++){ - std::vector l,l_diag; - for (int j = 0; j < n; j++){ - - double px = diagram[j].first; double py = diagram[j].second; - double proj_diag = (px+py)/2; + for (int i = 0; i < this->approx; i++) { + std::vector l, l_diag; + for (int j = 0; j < n; j++) { + double px = diagram[j].first; + double py = diagram[j].second; + double proj_diag = (px + py) / 2; - l.push_back ( px * cos(-pi/2+i*step) + py * sin(-pi/2+i*step) ); - l_diag.push_back ( proj_diag * cos(-pi/2+i*step) + proj_diag * sin(-pi/2+i*step) ); + l.push_back(px * cos(-pi / 2 + i * step) + py * sin(-pi / 2 + i * step)); + l_diag.push_back(proj_diag * cos(-pi / 2 + i * step) + proj_diag * sin(-pi / 2 + i * step)); } - std::sort(l.begin(), l.end()); std::sort(l_diag.begin(), l_diag.end()); - projections.push_back(std::move(l)); projections_diagonal.push_back(std::move(l_diag)); - + std::sort(l.begin(), l.end()); + std::sort(l_diag.begin(), l_diag.end()); + projections.push_back(std::move(l)); + projections_diagonal.push_back(std::move(l_diag)); } diagram.clear(); @@ -101,179 +109,254 @@ class Sliced_Wasserstein { } // Compute the angle formed by two points of a PD - double compute_angle(const Persistence_diagram & diag, int i, int j) const { - if(diag[i].second == diag[j].second) return pi/2; else return atan((diag[j].first-diag[i].first)/(diag[i].second-diag[j].second)); + double compute_angle(const Persistence_diagram& diag, int i, int j) const { + if (diag[i].second == diag[j].second) + return pi / 2; + else + return atan((diag[j].first - diag[i].first) / (diag[i].second - diag[j].second)); } - // Compute the integral of |cos()| between alpha and beta, valid only if alpha is in [-pi,pi] and beta-alpha is in [0,pi] + // Compute the integral of |cos()| between alpha and beta, valid only if alpha is in [-pi,pi] and beta-alpha is in + // [0,pi] double compute_int_cos(double alpha, double beta) const { double res = 0; - if (alpha >= 0 && alpha <= pi){ - if (cos(alpha) >= 0){ - if(pi/2 <= beta){res = 2-sin(alpha)-sin(beta);} - else{res = sin(beta)-sin(alpha);} - } - else{ - if(1.5*pi <= beta){res = 2+sin(alpha)+sin(beta);} - else{res = sin(alpha)-sin(beta);} + if (alpha >= 0 && alpha <= pi) { + if (cos(alpha) >= 0) { + if (pi / 2 <= beta) { + res = 2 - sin(alpha) - sin(beta); + } else { + res = sin(beta) - sin(alpha); + } + } else { + if (1.5 * pi <= beta) { + res = 2 + sin(alpha) + sin(beta); + } else { + res = sin(alpha) - sin(beta); + } } } - if (alpha >= -pi && alpha <= 0){ - if (cos(alpha) <= 0){ - if(-pi/2 <= beta){res = 2+sin(alpha)+sin(beta);} - else{res = sin(alpha)-sin(beta);} - } - else{ - if(pi/2 <= beta){res = 2-sin(alpha)-sin(beta);} - else{res = sin(beta)-sin(alpha);} + if (alpha >= -pi && alpha <= 0) { + if (cos(alpha) <= 0) { + if (-pi / 2 <= beta) { + res = 2 + sin(alpha) + sin(beta); + } else { + res = sin(alpha) - sin(beta); + } + } else { + if (pi / 2 <= beta) { + res = 2 - sin(alpha) - sin(beta); + } else { + res = sin(beta) - sin(alpha); + } } } return res; } - double compute_int(double theta1, double theta2, int p, int q, const Persistence_diagram & diag1, const Persistence_diagram & diag2) const { - double norm = std::sqrt( (diag1[p].first-diag2[q].first)*(diag1[p].first-diag2[q].first) + (diag1[p].second-diag2[q].second)*(diag1[p].second-diag2[q].second) ); + double compute_int(double theta1, double theta2, int p, int q, const Persistence_diagram& diag1, + const Persistence_diagram& diag2) const { + double norm = std::sqrt((diag1[p].first - diag2[q].first) * (diag1[p].first - diag2[q].first) + + (diag1[p].second - diag2[q].second) * (diag1[p].second - diag2[q].second)); double angle1; - if (diag1[p].first == diag2[q].first) angle1 = theta1 - pi/2; else angle1 = theta1 - atan((diag1[p].second-diag2[q].second)/(diag1[p].first-diag2[q].first)); + if (diag1[p].first == diag2[q].first) + angle1 = theta1 - pi / 2; + else + angle1 = theta1 - atan((diag1[p].second - diag2[q].second) / (diag1[p].first - diag2[q].first)); double angle2 = angle1 + theta2 - theta1; - double integral = compute_int_cos(angle1,angle2); - return norm*integral; + double integral = compute_int_cos(angle1, angle2); + return norm * integral; } // 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")); + 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")); - Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram; double sw = 0; - - if(this->approx == -1){ + Persistence_diagram diagram1 = this->diagram; + Persistence_diagram diagram2 = second.diagram; + double sw = 0; + if (this->approx == -1) { // Add projections onto diagonal. - int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); - double min_ordinate = std::numeric_limits::max(); double min_abscissa = std::numeric_limits::max(); - double max_ordinate = std::numeric_limits::lowest(); double max_abscissa = std::numeric_limits::lowest(); - for (int i = 0; i < n2; i++){ - min_ordinate = std::min(min_ordinate, diagram2[i].second); min_abscissa = std::min(min_abscissa, diagram2[i].first); - max_ordinate = std::max(max_ordinate, diagram2[i].second); max_abscissa = std::max(max_abscissa, diagram2[i].first); - diagram1.emplace_back( (diagram2[i].first+diagram2[i].second)/2, (diagram2[i].first+diagram2[i].second)/2 ); + int n1, n2; + n1 = diagram1.size(); + n2 = diagram2.size(); + double min_ordinate = std::numeric_limits::max(); + double min_abscissa = std::numeric_limits::max(); + double max_ordinate = std::numeric_limits::lowest(); + double max_abscissa = std::numeric_limits::lowest(); + for (int i = 0; i < n2; i++) { + min_ordinate = std::min(min_ordinate, diagram2[i].second); + min_abscissa = std::min(min_abscissa, diagram2[i].first); + max_ordinate = std::max(max_ordinate, diagram2[i].second); + max_abscissa = std::max(max_abscissa, diagram2[i].first); + diagram1.emplace_back((diagram2[i].first + diagram2[i].second) / 2, + (diagram2[i].first + diagram2[i].second) / 2); } - for (int i = 0; i < n1; i++){ - min_ordinate = std::min(min_ordinate, diagram1[i].second); min_abscissa = std::min(min_abscissa, diagram1[i].first); - max_ordinate = std::max(max_ordinate, diagram1[i].second); max_abscissa = std::max(max_abscissa, diagram1[i].first); - diagram2.emplace_back( (diagram1[i].first+diagram1[i].second)/2, (diagram1[i].first+diagram1[i].second)/2 ); + for (int i = 0; i < n1; i++) { + min_ordinate = std::min(min_ordinate, diagram1[i].second); + min_abscissa = std::min(min_abscissa, diagram1[i].first); + max_ordinate = std::max(max_ordinate, diagram1[i].second); + max_abscissa = std::max(max_abscissa, diagram1[i].first); + diagram2.emplace_back((diagram1[i].first + diagram1[i].second) / 2, + (diagram1[i].first + diagram1[i].second) / 2); } int num_pts_dgm = diagram1.size(); // Slightly perturb the points so that the PDs are in generic positions. double epsilon = 0.0001; - double thresh_y = (max_ordinate-min_ordinate) * epsilon; double thresh_x = (max_abscissa-min_abscissa) * epsilon; - std::random_device rd; std::default_random_engine re(rd()); std::uniform_real_distribution uni(-1,1); - for (int i = 0; i < num_pts_dgm; i++){ + double thresh_y = (max_ordinate - min_ordinate) * epsilon; + double thresh_x = (max_abscissa - min_abscissa) * epsilon; + std::random_device rd; + std::default_random_engine re(rd()); + std::uniform_real_distribution uni(-1, 1); + for (int i = 0; i < num_pts_dgm; i++) { double u = uni(re); - diagram1[i].first += u*thresh_x; diagram1[i].second += u*thresh_y; - diagram2[i].first += u*thresh_x; diagram2[i].second += u*thresh_y; + diagram1[i].first += u * thresh_x; + diagram1[i].second += u * thresh_y; + diagram2[i].first += u * thresh_x; + diagram2[i].second += u * thresh_y; } // Compute all angles in both PDs. - std::vector > > angles1, angles2; - for (int i = 0; i < num_pts_dgm; i++){ - for (int j = i+1; j < num_pts_dgm; j++){ - double theta1 = compute_angle(diagram1,i,j); double theta2 = compute_angle(diagram2,i,j); - angles1.emplace_back(theta1, std::pair(i,j)); - angles2.emplace_back(theta2, std::pair(i,j)); + std::vector > > angles1, angles2; + for (int i = 0; i < num_pts_dgm; i++) { + for (int j = i + 1; j < num_pts_dgm; j++) { + double theta1 = compute_angle(diagram1, i, j); + double theta2 = compute_angle(diagram2, i, j); + angles1.emplace_back(theta1, std::pair(i, j)); + angles2.emplace_back(theta2, std::pair(i, j)); } } // Sort angles. - std::sort(angles1.begin(), angles1.end(), [](const std::pair >& p1, const std::pair >& p2){return (p1.first < p2.first);}); - std::sort(angles2.begin(), angles2.end(), [](const std::pair >& p1, const std::pair >& p2){return (p1.first < p2.first);}); + std::sort(angles1.begin(), angles1.end(), + [](const std::pair >& p1, + const std::pair >& p2) { return (p1.first < p2.first); }); + std::sort(angles2.begin(), angles2.end(), + [](const std::pair >& p1, + const std::pair >& p2) { return (p1.first < p2.first); }); // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2). std::vector orderp1, orderp2; - for (int i = 0; i < num_pts_dgm; i++){ orderp1.push_back(i); orderp2.push_back(i); } - std::sort( orderp1.begin(), orderp1.end(), [&](int i, int j){ if(diagram1[i].second != diagram1[j].second) return (diagram1[i].second < diagram1[j].second); else return (diagram1[i].first > diagram1[j].first); } ); - std::sort( orderp2.begin(), orderp2.end(), [&](int i, int j){ if(diagram2[i].second != diagram2[j].second) return (diagram2[i].second < diagram2[j].second); else return (diagram2[i].first > diagram2[j].first); } ); + for (int i = 0; i < num_pts_dgm; i++) { + orderp1.push_back(i); + orderp2.push_back(i); + } + std::sort(orderp1.begin(), orderp1.end(), [&](int i, int j) { + if (diagram1[i].second != diagram1[j].second) + return (diagram1[i].second < diagram1[j].second); + else + return (diagram1[i].first > diagram1[j].first); + }); + std::sort(orderp2.begin(), orderp2.end(), [&](int i, int j) { + if (diagram2[i].second != diagram2[j].second) + return (diagram2[i].second < diagram2[j].second); + else + return (diagram2[i].first > diagram2[j].first); + }); // Find the inverses of the orders. - std::vector order1(num_pts_dgm); std::vector order2(num_pts_dgm); - for(int i = 0; i < num_pts_dgm; i++){ order1[orderp1[i]] = i; order2[orderp2[i]] = i; } + std::vector order1(num_pts_dgm); + std::vector order2(num_pts_dgm); + for (int i = 0; i < num_pts_dgm; i++) { + order1[orderp1[i]] = i; + order2[orderp2[i]] = i; + } // Record all inversions of points in the orders as theta varies along the positive half-disk. - std::vector > > anglePerm1(num_pts_dgm); - std::vector > > anglePerm2(num_pts_dgm); + std::vector > > anglePerm1(num_pts_dgm); + std::vector > > anglePerm2(num_pts_dgm); int m1 = angles1.size(); - for (int i = 0; i < m1; i++){ - double theta = angles1[i].first; int p = angles1[i].second.first; int q = angles1[i].second.second; - anglePerm1[order1[p]].emplace_back(p,theta); - anglePerm1[order1[q]].emplace_back(q,theta); - int a = order1[p]; int b = order1[q]; order1[p] = b; order1[q] = a; + for (int i = 0; i < m1; i++) { + double theta = angles1[i].first; + int p = angles1[i].second.first; + int q = angles1[i].second.second; + anglePerm1[order1[p]].emplace_back(p, theta); + anglePerm1[order1[q]].emplace_back(q, theta); + int a = order1[p]; + int b = order1[q]; + order1[p] = b; + order1[q] = a; } int m2 = angles2.size(); - for (int i = 0; i < m2; i++){ - double theta = angles2[i].first; int p = angles2[i].second.first; int q = angles2[i].second.second; - anglePerm2[order2[p]].emplace_back(p,theta); - anglePerm2[order2[q]].emplace_back(q,theta); - int a = order2[p]; int b = order2[q]; order2[p] = b; order2[q] = a; + for (int i = 0; i < m2; i++) { + double theta = angles2[i].first; + int p = angles2[i].second.first; + int q = angles2[i].second.second; + anglePerm2[order2[p]].emplace_back(p, theta); + anglePerm2[order2[q]].emplace_back(q, theta); + int a = order2[p]; + int b = order2[q]; + order2[p] = b; + order2[q] = a; } - for (int i = 0; i < num_pts_dgm; i++){ - anglePerm1[order1[i]].emplace_back(i,pi/2); - anglePerm2[order2[i]].emplace_back(i,pi/2); + for (int i = 0; i < num_pts_dgm; i++) { + anglePerm1[order1[i]].emplace_back(i, pi / 2); + anglePerm2[order2[i]].emplace_back(i, pi / 2); } // Compute the SW distance with the list of inversions. - for (int i = 0; i < num_pts_dgm; i++){ - std::vector > u,v; u = anglePerm1[i]; v = anglePerm2[i]; - double theta1, theta2; theta1 = -pi/2; - unsigned int ku, kv; ku = 0; kv = 0; theta2 = std::min(u[ku].second,v[kv].second); - while(theta1 != pi/2){ - if(diagram1[u[ku].first].first != diagram2[v[kv].first].first || diagram1[u[ku].first].second != diagram2[v[kv].first].second) - if(theta1 != theta2) - sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2); + for (int i = 0; i < num_pts_dgm; i++) { + std::vector > u, v; + u = anglePerm1[i]; + v = anglePerm2[i]; + double theta1, theta2; + theta1 = -pi / 2; + unsigned int ku, kv; + ku = 0; + kv = 0; + theta2 = std::min(u[ku].second, v[kv].second); + while (theta1 != pi / 2) { + if (diagram1[u[ku].first].first != diagram2[v[kv].first].first || + diagram1[u[ku].first].second != diagram2[v[kv].first].second) + if (theta1 != theta2) sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2); theta1 = theta2; - if ( (theta2 == u[ku].second) && ku < u.size()-1 ) ku++; - if ( (theta2 == v[kv].second) && kv < v.size()-1 ) kv++; + if ((theta2 == u[ku].second) && ku < u.size() - 1) ku++; + if ((theta2 == v[kv].second) && kv < v.size() - 1) kv++; theta2 = std::min(u[ku].second, v[kv].second); } } - } - - - else{ - - double step = pi/this->approx; std::vector v1, v2; - for (int i = 0; i < this->approx; i++){ - - 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; - + } else { + double step = pi / this->approx; + std::vector v1, v2; + for (int i = 0; i < this->approx; i++) { + 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; } } - return sw/pi; + return sw / pi; } public: - /** \brief Sliced Wasserstein kernel constructor. * \implements Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_scalar_product * \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 random perturbation. If positive, then projections of the diagram - * points on all directions are stored in memory to reduce computation time. + * @param[in] _approx number of directions used to approximate the integral in the Sliced Wasserstein distance, set + * to -1 for random perturbation. If positive, then projections of the diagram points on all + * directions are stored in memory to reduce computation time. * */ - Sliced_Wasserstein(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 10):diagram(_diagram), approx(_approx), sigma(_sigma) {build_rep();} + 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 @@ -282,9 +365,10 @@ class Sliced_Wasserstein { * @param[in] second other instance of class Sliced_Wasserstein. * */ - double compute_scalar_product(const Sliced_Wasserstein & second) const { - GUDHI_CHECK(this->sigma == second.sigma, std::invalid_argument("Error: different sigma values for representations")); - return std::exp(-compute_sliced_wasserstein_distance(second)/(2*this->sigma*this->sigma)); + double compute_scalar_product(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->sigma == second.sigma, + std::invalid_argument("Error: different sigma values for representations")); + return std::exp(-compute_sliced_wasserstein_distance(second) / (2 * this->sigma * this->sigma)); } /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. @@ -294,13 +378,14 @@ class Sliced_Wasserstein { * @param[in] second other instance of class Sliced_Wasserstein. * */ - double distance(const Sliced_Wasserstein & second) const { - GUDHI_CHECK(this->sigma == second.sigma, std::invalid_argument("Error: different sigma values for representations")); - return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), 0.5); + double distance(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->sigma == second.sigma, + std::invalid_argument("Error: different sigma values for representations")); + return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) - + 2 * this->compute_scalar_product(second)); } - -}; // class Sliced_Wasserstein +}; // class Sliced_Wasserstein } // namespace Persistence_representations } // namespace Gudhi -- cgit v1.2.3