diff options
Diffstat (limited to 'src/Persistence_representations')
14 files changed, 1305 insertions, 46 deletions
diff --git a/src/Persistence_representations/doc/Persistence_representations_doc.h b/src/Persistence_representations/doc/Persistence_representations_doc.h index 4d850a02..d0b02739 100644 --- a/src/Persistence_representations/doc/Persistence_representations_doc.h +++ b/src/Persistence_representations/doc/Persistence_representations_doc.h @@ -24,7 +24,6 @@ #define DOC_GUDHI_STAT_H_ namespace Gudhi { - namespace Persistence_representations { /** \defgroup Persistence_representations Persistence representations @@ -129,33 +128,35 @@ namespace Persistence_representations { function \f$L : \mathbb{N} \times \mathbb{R} \to [0,\infty)\f$ of two variables, if we define \f$L(k,t) = \lambda_k(t)\f$. - The detailed description of algorithms used to compute persistence landscapes can be found in - \cite bubenik_dlotko_landscapes_2016. - Note that this implementation provides exact representation of landscapes. That have many advantages, but also a few - drawbacks. For instance, as discussed - in \cite bubenik_dlotko_landscapes_2016, the exact representation of landscape may be of quadratic size with respect - to the input persistence diagram. It may therefore happen - that, for very large diagrams, using this representation may be memory--prohibitive. In such a case, there are two - possible ways to proceed: + The detailed description of algorithms used to compute persistence landscapes can be found in \cite bubenik_dlotko_landscapes_2016. + Note that this implementation provides exact representation of landscapes. That have many advantages, but also a few drawbacks. + For instance, as discussed in \cite bubenik_dlotko_landscapes_2016, the exact representation of landscape may be of quadratic size with respect + to the input persistence diagram. It may therefore happen that, for very large diagrams, using this representation may be memory--prohibitive. + In such a case, there are two possible ways to proceed: - \li Use non exact representation on a grid described in the Section \ref sec_landscapes_on_grid. + \li Use representation on a grid---see section \ref sec_landscapes_on_grid. \li Compute just a number of initial nonzero landscapes. This option is available from C++ level as a last parameter of the constructor of persistence landscape (set by default to std::numeric_limits<size_t>::max()). \section sec_landscapes_on_grid Persistence Landscapes on a grid + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_landscape_on_grid <br> - This is an alternative, not--exact, representation of persistence landscapes defined in the Section \ref - sec_persistence_landscapes. Unlike in the Section \ref sec_persistence_landscapes we build a - representation of persistence landscape by sampling its values on a finite, equally distributed grid of points. - Since, the persistence landscapes that originate from persistence diagrams have slope \f$1\f$ or \f$-1\f$, we have an - estimate of a region between the grid points where the landscape cab be located. - That allows to estimate an error make when performing various operations on landscape. Note that for average - landscapes the slope is in range \f$[-1,1]\f$ and similar estimate can be used. + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_landscape_on_grid_exact <br> + + Here, we provide alternative, not exact, representations of persistence landscapes defined in Section \ref sec_persistence_landscapes. + Unlike Section \ref sec_persistence_landscapes, we build representations of persistence landscapes by evaluating the landscape functions on a finite, equally distributed grid of points. + We propose two different representations depending on whether the persistence intervals are also mapped on the grid (Persistence_landscape_on_grid) or not (Persistence_landscape_on_grid_exact). + This makes a big difference since mapping the intervals on the grid makes the computation time smaller but only provides an approximation of the landscape values. - Due to a lack of rigorous description of the algorithms to deal with this non--rigorous representation of persistence - landscapes in the literature, we are providing a short discussion of them in below. + Since persistence landscapes originating from persistence diagrams have slope \f$1\f$ or \f$-1\f$, we have an + estimate of a region between the grid points where the landscapes can be located. + That allows to estimate an error made when performing various operations on landscapes. Note that for average + landscapes the slope is in range \f$[-1,1]\f$ and similar estimates can be used. + + Due to the lack of rigorous description of the algorithms for these non rigorous representations of persistence + landscapes in the literature, we provide a short discussion below. Let us assume that we want to compute persistence landscape on a interval \f$[x,y]\f$. Let us assume that we want to use \f$N\f$ grid points for that purpose. @@ -167,11 +168,11 @@ namespace Persistence_representations { functions) on the i-th point of a grid, i.e. \f$x + i \frac{y-x}{N}\f$. When averaging two persistence landscapes represented by a grid we need to make sure that they are defined in a - compatible grids. I.e. the intervals \f$[x,y]\f$ on which they are defined are + compatible grids, i.e. the intervals \f$[x,y]\f$ on which they are defined are the same, and the numbers of grid points \f$N\f$ are the same in both cases. If this is the case, we simply compute - point-wise averages of the entries of corresponding - vectors (In this whole section we assume that if one vector of numbers is shorter than another, we extend the shorter - one with zeros so that they have the same length.) + point-wise averages of the entries of the corresponding + vectors (in this whole section we assume that if one vector of numbers is shorter than the other, we extend the shortest + one with zeros so that they have the same length). Computations of distances between two persistence landscapes on a grid is not much different than in the rigorous case. In this case, we sum up the distances between the same levels of @@ -180,11 +181,11 @@ namespace Persistence_representations { Similarly as in case of distance, when computing the scalar product of two persistence landscapes on a grid, we sum up the scalar products of corresponding levels of landscapes. For each level, - we assume that the persistence landscape on a grid between two grid points is approximated by linear function. - Therefore to compute scalar product of two corresponding levels of landscapes, + we assume that the persistence landscape on a grid between two grid points is approximated by a linear function. + Therefore to compute the scalar product of two corresponding levels of landscapes, we sum up the integrals of products of line segments for every pair of constitutive grid points. - Note that for this representation we need to specify a few parameters: + Note that for these representations we need to specify a few parameters: \li Begin and end point of a grid -- the interval \f$[x,y]\f$ (real numbers). \li Number of points in a grid (positive integer \f$N\f$). @@ -193,29 +194,33 @@ namespace Persistence_representations { Note that the same representation is used in TDA R-package \cite Fasy_Kim_Lecci_Maria_tda. \section sec_persistence_heat_maps Persistence heat maps + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_heat_maps <br> - This is a general class of discrete structures which are based on idea of placing a kernel in the points of - persistence diagrams. + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_heat_maps_exact <br> + + This is a general class of discrete structures which are based on idea of placing a kernel in the points of persistence diagrams. This idea appeared in work by many authors over the last 15 years. As far as we know this idea was firstly described in the work of Bologna group in \cite Ferri_Frosini_comparision_sheme_1 and \cite Ferri_Frosini_comparision_sheme_2. Later it has been described by Colorado State University group in \cite Persistence_Images_2017. The presented paper - in the first time provide a discussion of stability of the representation. - Also, the same ideas are used in construction of two recent kernels used for machine learning: - \cite Kusano_Fukumizu_Hiraoka_PWGK and \cite Reininghaus_Huber_ALL_PSSK. Both the kernel's construction uses - interesting ideas to ensure stability of the representation with respect to Wasserstein metric. In the kernel + in the first time provided a discussion of stability of this representation. + Also, the same ideas are used in the construction of two recent kernels used for machine learning: + \cite Kusano_Fukumizu_Hiraoka_PWGK and \cite Reininghaus_Huber_ALL_PSSK. Both the kernels use + interesting ideas to ensure stability of the representations with respect to the 1-Wasserstein metric. In the kernel presented in \cite Kusano_Fukumizu_Hiraoka_PWGK, a scaling function is used to multiply the Gaussian kernel in the - way that the points close to diagonal got low weight and consequently do not have a big influence on the resulting + way that the points close to diagonal have low weights and consequently do not have a big influence on the resulting distribution. In \cite Reininghaus_Huber_ALL_PSSK for every point \f$(b,d)\f$ two Gaussian kernels are added: first, with a weight 1 in a point \f$(b,d)\f$, and the second, with the weight -1 for a point \f$(b,d)\f$. In both cases, the representations are stable with respect to 1-Wasserstein distance. - In Persistence\_representations package we currently implement a discretization of the distributions described above. - The base of this implementation is 2-dimensional array of pixels. Each pixel have assigned a real value which - is a sum of values of distributions induced by each point of the persistence diagram. At the moment we compute the - sum of values on a center of a pixels. It can be easily extended to any other function - (like for instance sum of integrals of the intermediate distribution on a pixel). + In Persistence_representations package, we currently implement a discretization of the distributions described above. + The base of this implementation is a 2-dimensional array of pixels. To each pixel is assigned a real value which + is the sum of the distribution values induced by each point of the persistence diagram. + As for Persistence_landscapes, we propose two different representations depending on whether the persistence intervals are also mapped on the pixels + (Persistence_heat_maps) or not (Persistence_heat_maps_exact). + At the moment we compute the sum over the evaluations of the distributions on the pixel centers. It can be easily extended to any other function + (like for instance the sum of the integrals of the distributions over the pixels). - The parameters that determine the structure are the following: + Concerning Persistence_heat_maps, the parameters that determine the structure are the following: \li A positive integer k determining the size of the kernel we used (we always assume that the kernels are square). \li A filter: in practice a square matrix of a size \f$2k+1 \times 2k+1\f$. By default, this is a discretization of @@ -227,6 +232,7 @@ namespace Persistence_representations { to diagonal are given then sometimes the kernel have support that reaches the region below the diagonal. If the value of this parameter is true, then the values below diagonal can be erased. + Concerning Persistence_heat_maps_exact, only Gaussian kernels are implemented, so the parameters are the array of pixels, the weight functions for the Gaussians and the bandwidth of the Gaussians. \section sec_persistence_vectors Persistence vectors <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Vector_distances_in_diagram <br> @@ -250,6 +256,68 @@ namespace Persistence_representations { absolute value of differences between coordinates. A scalar product is a sum of products of values at the corresponding positions of two vectors. + + + + +\section sec_persistence_kernels Kernels on persistence diagrams + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Sliced_Wasserstein <br> + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_weighted_gaussian <br> + + Kernels for persistence diagrams can be regarded as infinite-dimensional vectorizations. More specifically, + they are similarity functions whose evaluations on pairs of persistence diagrams equals the scalar products + between images of these pairs under a map \f$\Phi\f$ taking values in a specific (possibly non Euclidean) Hilbert space \f$k(D_i, D_j) = \langle \Phi(D_i),\Phi(D_j)\rangle\f$. + Reciprocally, classical results of learning theory ensure that such a \f$\Phi\f$ exists for a given similarity function \f$k\f$ if and only if \f$k\f$ is <i>positive semi-definite</i>. + Kernels are designed for algorithms that can be <i>kernelized</i>, i.e., algorithms that only require to know scalar products between instances in order to run. + Examples of such algorithms include Support Vector Machines, Principal Component Analysis and Ridge Regression. + + There have been several attempts at defining kernels, i.e., positive semi-definite functions, between persistence diagrams within the last few years. We provide implementation + for three of them: + + \li the <i>Persistence Scale Space Kernel</i>---see \cite Reininghaus_Huber_ALL_PSSK, which is the classical scalar product between \f$L^2\f$ functions, where persistence diagrams + are turned into functions by centering and summing Gaussian functions over the diagram points and their symmetric counterparts w.r.t. the diagonal: \f$k(D_1,D_2)=\int \Phi(D_1)\Phi(D_2)\f$, + where \f$\Phi(D)=\sum_{p\in D} {\rm exp}\left(-\frac{\|p-\cdot\|_2^2}{2\sigma^2}\right)\f$. + + \li the <i>Persistence Weighted Gaussian Kernel</i>---see \cite Kusano_Fukumizu_Hiraoka_PWGK, which is a slight generalization of the previous kernel, is the scalar product between + weighted Kernel Mean Embeddings of persistence diagrams w.r.t. the Gaussian Kernel \f$k_G\f$ (with corresponding map \f$\Phi_G\f$) in \f$\mathbb{R}^2\f$: + \f$k(D_1,D_2)=\langle\sum_{p\in D_1} w(p)\Phi_G(p), \sum_{q\in D_2} w(q)\Phi_G(q)\rangle\f$ + + \li the <i>Sliced Wasserstein Kernel</i>---see \cite pmlr-v70-carriere17a, which takes the form of a Gaussian kernel with a specific distance between persistence diagrams + called the <i>Sliced Wasserstein Distance</i>: \f$k(D_1,D_2)={\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right)\f$ + + When launching: + + \code $> ./Sliced_Wasserstein + \endcode + + the program output is: + + \code $> Approx SW distance: 5.33648 + $> Exact SW distance: 5.33798 + $> Approx SW kernel: 0.0693743 + $> Exact SW kernel: 0.0693224 + $> Distance induced by approx SW kernel: 1.36428 + $> Distance induced by exact SW kernel: 1.3643 + \endcode + + + and when launching: + + \code $> ./Persistence_weighted_gaussian + \endcode + + the program output is: + + \code $> Approx PWG kernel: 1.21509 + $> Exact PWG kernel: 1.13628 + $> Distance induced by approx PWG kernel: 3.23354 + $> Distance induced by exact PWG kernel: 3.25697 + $> Approx Gaussian PWG kernel: 0.0194222 + $> Exact Gaussian PWG kernel: 0.0192524 + $> Approx PSS kernel: 0.134413 + $> Exact PSS kernel: 0.133394 + \endcode + */ /** @} */ // end defgroup Persistence_representations diff --git a/src/Persistence_representations/example/CMakeLists.txt b/src/Persistence_representations/example/CMakeLists.txt index 54d719ac..3142f19b 100644 --- a/src/Persistence_representations/example/CMakeLists.txt +++ b/src/Persistence_representations/example/CMakeLists.txt @@ -27,3 +27,22 @@ add_test(NAME Persistence_representations_example_heat_maps COMMAND $<TARGET_FILE:Persistence_representations_example_heat_maps>) install(TARGETS Persistence_representations_example_heat_maps DESTINATION bin) +add_executable ( Sliced_Wasserstein sliced_wasserstein.cpp ) +add_test(NAME Sliced_Wasserstein + COMMAND $<TARGET_FILE:Sliced_Wasserstein>) +install(TARGETS Sliced_Wasserstein DESTINATION bin) + +add_executable ( Persistence_weighted_gaussian persistence_weighted_gaussian.cpp ) +add_test(NAME Persistence_weighted_gaussian + COMMAND $<TARGET_FILE:Persistence_weighted_gaussian>) +install(TARGETS Persistence_weighted_gaussian DESTINATION bin) + +add_executable ( Persistence_heat_maps_exact persistence_heat_maps_exact.cpp ) +add_test(NAME Persistence_heat_maps_exact + COMMAND $<TARGET_FILE:Persistence_heat_maps_exact>) +install(TARGETS Persistence_heat_maps_exact DESTINATION bin) + +add_executable ( Persistence_landscape_on_grid_exact persistence_landscape_on_grid_exact.cpp ) +add_test(NAME Persistence_landscape_on_grid_exact + COMMAND $<TARGET_FILE:Persistence_landscape_on_grid_exact>) +install(TARGETS Persistence_landscape_on_grid_exact DESTINATION bin) diff --git a/src/Persistence_representations/example/persistence_heat_maps_exact.cpp b/src/Persistence_representations/example/persistence_heat_maps_exact.cpp new file mode 100644 index 00000000..f15b710d --- /dev/null +++ b/src/Persistence_representations/example/persistence_heat_maps_exact.cpp @@ -0,0 +1,55 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Persistence_heat_maps_exact.h> +#include <gudhi/Weight_functions.h> + +#include <iostream> +#include <vector> +#include <utility> +#include <string> + +using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram; +using PI = Gudhi::Persistence_representations::Persistence_heat_maps_exact; +using Weight = std::function<double (std::pair<double,double>) >; + +int main(int argc, char** argv) { + + Persistence_diagram persistence; + + persistence.push_back(std::make_pair(1, 2)); + persistence.push_back(std::make_pair(6, 8)); + persistence.push_back(std::make_pair(0, 4)); + persistence.push_back(std::make_pair(3, 8)); + + double min_x = 0.0; double max_x = 10.0; int res_x = 100; double min_y = 0.0; double max_y = 10.0; int res_y = 100; double sigma = 1.0; Weight weight = Gudhi::Persistence_representations::linear_weight; + + PI pim(persistence, min_x, max_x, res_x, min_y, max_y, res_y, weight, sigma); + std::vector<std::vector<double> > P = pim.vectorize(); + + for(int i = 0; i < res_y; i++){ + for(int j = 0; j < res_x; j++) std::cout << P[i][j] << " "; + std::cout << std::endl; + } + + return 0; +} diff --git a/src/Persistence_representations/example/persistence_landscape_on_grid_exact.cpp b/src/Persistence_representations/example/persistence_landscape_on_grid_exact.cpp new file mode 100644 index 00000000..da27bc5a --- /dev/null +++ b/src/Persistence_representations/example/persistence_landscape_on_grid_exact.cpp @@ -0,0 +1,52 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Persistence_landscape_on_grid_exact.h> + +#include <iostream> +#include <vector> +#include <utility> + +using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram; +using LS = Gudhi::Persistence_representations::Persistence_landscape_on_grid_exact; + +int main(int argc, char** argv) { + + Persistence_diagram persistence; + + persistence.push_back(std::make_pair(1, 2)); + persistence.push_back(std::make_pair(6, 8)); + persistence.push_back(std::make_pair(0, 4)); + persistence.push_back(std::make_pair(3, 8)); + + int nb_ls = 3; double min_x = 0.0; double max_x = 10.0; int res_x = 100; + + LS ls(persistence, nb_ls, min_x, max_x, res_x); + std::vector<std::vector<double> > L = ls.vectorize(); + + for(int i = 0; i < nb_ls; i++){ + for(int j = 0; j < res_x; j++) std::cout << L[i][j] << " "; + std::cout << std::endl; + } + + return 0; +} diff --git a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp new file mode 100644 index 00000000..7945e4f1 --- /dev/null +++ b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp @@ -0,0 +1,99 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Persistence_weighted_gaussian.h> + +#include <iostream> +#include <vector> +#include <utility> + +using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram; +using PWG = Gudhi::Persistence_representations::Persistence_weighted_gaussian; + +int main(int argc, char** argv) { + + Persistence_diagram persistence1, persistence2; + + persistence1.push_back(std::make_pair(1, 2)); + persistence1.push_back(std::make_pair(6, 8)); + persistence1.push_back(std::make_pair(0, 4)); + persistence1.push_back(std::make_pair(3, 8)); + + persistence2.push_back(std::make_pair(2, 9)); + persistence2.push_back(std::make_pair(1, 6)); + persistence2.push_back(std::make_pair(3, 5)); + persistence2.push_back(std::make_pair(6, 10)); + + double sigma = 1; + double tau = 1; + int m = 10000; + + PWG PWG1(persistence1, sigma, m, Gudhi::Persistence_representations::arctan_weight(1,1)); + PWG PWG2(persistence2, sigma, m, Gudhi::Persistence_representations::arctan_weight(1,1)); + + PWG PWGex1(persistence1, sigma, -1, Gudhi::Persistence_representations::arctan_weight(1,1)); + PWG PWGex2(persistence2, sigma, -1, Gudhi::Persistence_representations::arctan_weight(1,1)); + + + // Linear PWG + + std::cout << "Approx PWG kernel: " << PWG1.compute_scalar_product (PWG2) << std::endl; + std::cout << "Exact PWG kernel: " << PWGex1.compute_scalar_product (PWGex2) << std::endl; + + std::cout << "Distance induced by approx PWG kernel: " << PWG1.distance (PWG2) << std::endl; + std::cout << "Distance induced by exact PWG kernel: " << PWGex1.distance (PWGex2) << std::endl; + + + + + + + + // Gaussian PWG + + std::cout << "Approx Gaussian PWG kernel: " << std::exp( -PWG1.distance (PWG2) ) / (2*tau*tau) << std::endl; + std::cout << "Exact Gaussian PWG kernel: " << std::exp( -PWGex1.distance (PWGex2) ) / (2*tau*tau) << std::endl; + + + + + + + + // PSS + + Persistence_diagram pd1 = persistence1; int numpts = persistence1.size(); for(int i = 0; i < numpts; i++) pd1.emplace_back(persistence1[i].second,persistence1[i].first); + Persistence_diagram pd2 = persistence2; numpts = persistence2.size(); for(int i = 0; i < numpts; i++) pd2.emplace_back(persistence2[i].second,persistence2[i].first); + + PWG pwg1(pd1, 2*std::sqrt(sigma), m, Gudhi::Persistence_representations::pss_weight); + PWG pwg2(pd2, 2*std::sqrt(sigma), m, Gudhi::Persistence_representations::pss_weight); + + PWG pwgex1(pd1, 2*std::sqrt(sigma), -1, Gudhi::Persistence_representations::pss_weight); + PWG pwgex2(pd2, 2*std::sqrt(sigma), -1, Gudhi::Persistence_representations::pss_weight); + + std::cout << "Approx PSS kernel: " << pwg1.compute_scalar_product (pwg2) / (16*Gudhi::Persistence_representations::pi*sigma) << std::endl; + std::cout << "Exact PSS kernel: " << pwgex1.compute_scalar_product (pwgex2) / (16*Gudhi::Persistence_representations::pi*sigma) << std::endl; + + + + return 0; +} diff --git a/src/Persistence_representations/example/sliced_wasserstein.cpp b/src/Persistence_representations/example/sliced_wasserstein.cpp new file mode 100644 index 00000000..2104e2b2 --- /dev/null +++ b/src/Persistence_representations/example/sliced_wasserstein.cpp @@ -0,0 +1,61 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Sliced_Wasserstein.h> + +#include <iostream> +#include <vector> +#include <utility> + +using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram; +using SW = Gudhi::Persistence_representations::Sliced_Wasserstein; + +int main(int argc, char** argv) { + + Persistence_diagram persistence1, persistence2; + + persistence1.push_back(std::make_pair(1, 2)); + persistence1.push_back(std::make_pair(6, 8)); + persistence1.push_back(std::make_pair(0, 4)); + persistence1.push_back(std::make_pair(3, 8)); + + persistence2.push_back(std::make_pair(2, 9)); + persistence2.push_back(std::make_pair(1, 6)); + persistence2.push_back(std::make_pair(3, 5)); + persistence2.push_back(std::make_pair(6, 10)); + + + SW sw1(persistence1, 1, 100); + SW sw2(persistence2, 1, 100); + + SW swex1(persistence1, 1, -1); + SW swex2(persistence2, 1, -1); + + std::cout << "Approx SW distance: " << sw1.compute_sliced_wasserstein_distance(sw2) << std::endl; + std::cout << "Exact SW distance: " << swex1.compute_sliced_wasserstein_distance(swex2) << std::endl; + std::cout << "Approx SW kernel: " << sw1.compute_scalar_product(sw2) << std::endl; + std::cout << "Exact SW kernel: " << swex1.compute_scalar_product(swex2) << std::endl; + std::cout << "Distance induced by approx SW kernel: " << sw1.distance(sw2) << std::endl; + std::cout << "Distance induced by exact SW kernel: " << swex1.distance(swex2) << std::endl; + + return 0; +} diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps_exact.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps_exact.h new file mode 100644 index 00000000..7c5b2fdc --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps_exact.h @@ -0,0 +1,125 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef PERSISTENCE_HEAT_MAPS_EXACT_H_ +#define PERSISTENCE_HEAT_MAPS_EXACT_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> +#include <gudhi/Weight_functions.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 { + +/** + * \class Persistence_heat_maps_exact gudhi/Persistence_heat_maps_exact.h + * \brief A class implementing exact persistence heat maps. + * + * \ingroup Persistence_representations + * + * \details + * + * In this class, we propose a way to approximate persistence heat maps, or persistence surfaces, by centering weighted Gaussians on each point of the persistence diagram, and evaluating these (exact) weighted Gaussian functions + * on the pixels of a 2D grid. Note that this scheme is different from the one proposed in Persistence_heat_maps, which first maps the points of the diagram to a 2D grid, and then evaluates the (approximate) weighted Gaussian functions. + * Hence, the difference is that we do not modify the diagram in this implementation, but the code can be slower to run. +**/ + +class Persistence_heat_maps_exact { + + protected: + Persistence_diagram diagram; + int res_x, res_y; + double min_x, max_x, min_y, max_y; + Weight weight; + double sigma; + + public: + + /** \brief Persistence_heat_maps_exact constructor. + * \ingroup Persistence_heat_maps_exact + * + * @param[in] _diagram persistence diagram. + * @param[in] _min_x minimum value of pixel abscissa. + * @param[in] _max_x maximum value of pixel abscissa. + * @param[in] _res_x number of pixels for the x-direction. + * @param[in] _min_y minimum value of pixel ordinate. + * @param[in] _max_y maximum value of pixel ordinate. + * @param[in] _res_y number of pixels for the y-direction. + * @param[in] _weight weight function for the Gaussians. + * @param[in] _sigma bandwidth parameter for the Gaussians. + * + */ + Persistence_heat_maps_exact(const Persistence_diagram & _diagram, double _min_x = 0.0, double _max_x = 1.0, int _res_x = 10, double _min_y = 0.0, double _max_y = 1.0, int _res_y = 10, const Weight & _weight = arctan_weight(1,1), double _sigma = 1.0){ + diagram = _diagram; min_x = _min_x; max_x = _max_x; res_x = _res_x; min_y = _min_y; max_y = _max_y; res_y = _res_y, weight = _weight; sigma = _sigma; + } + + /** \brief Computes the persistence image of a diagram. + * \ingroup Persistence_heat_maps_exact + * + */ + std::vector<std::vector<double> > vectorize() const { + std::vector<std::vector<double> > im; for(int i = 0; i < res_y; i++) im.emplace_back(); + double step_x = (max_x - min_x)/(res_x - 1); double step_y = (max_y - min_y)/(res_y - 1); + + int num_pts = diagram.size(); + + for(int i = 0; i < res_y; i++){ + double y = min_y + i*step_y; + for(int j = 0; j < res_x; j++){ + double x = min_x + j*step_x; + + double pixel_value = 0; + for(int k = 0; k < num_pts; k++){ + double px = diagram[k].first; double py = diagram[k].second; + pixel_value += weight(std::pair<double,double>(px,py)) * std::exp( -((x-px)*(x-px) + (y-(py-px))*(y-(py-px))) / (2*sigma*sigma) ) / (sigma*std::sqrt(2*pi)); + } + im[i].push_back(pixel_value); + + } + } + + return im; + + } + + + + +}; // class Persistence_heat_maps_exact +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_HEAT_MAPS_EXACT_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid_exact.h b/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid_exact.h new file mode 100644 index 00000000..25f71e27 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid_exact.h @@ -0,0 +1,107 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef LANDSCAPE_H_ +#define LANDSCAPE_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#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 { + +/** + * \class Persistence_landscape_on_grid_exact gudhi/Persistence_landscape_on_grid_exact.h + * \brief A class implementing exact persistence landscapes by approximating them on a collection of grid points + * + * \ingroup Persistence_representations + * + * \details + * In this class, we propose a way to approximate landscapes by sampling the x-axis of the persistence diagram and evaluating the (exact) landscape functions on the sample projections onto the diagonal. Note that this is a different approximation scheme + * from the one proposed in Persistence_landscape_on_grid, which puts a grid on the diagonal, maps the persistence intervals on this grid and computes the (approximate) landscape functions on the samples. + * Hence, the difference is that we do not modify the diagram in this implementation, but the code can be slower to run. +**/ + +class Persistence_landscape_on_grid_exact { + + protected: + Persistence_diagram diagram; + int res_x, nb_ls; + double min_x, max_x; + + public: + + /** \brief Persistence_landscape_on_grid_exact constructor. + * \ingroup Persistence_landscape_on_grid_exact + * + * @param[in] _diagram persistence diagram. + * @param[in] _nb_ls number of landscape functions. + * @param[in] _min_x minimum value of samples. + * @param[in] _max_x maximum value of samples. + * @param[in] _res_x number of samples. + * + */ + Persistence_landscape_on_grid_exact(const Persistence_diagram & _diagram, int _nb_ls = 5, double _min_x = 0.0, double _max_x = 1.0, int _res_x = 10){diagram = _diagram; nb_ls = _nb_ls; min_x = _min_x; max_x = _max_x; res_x = _res_x;} + + /** \brief Computes the landscape approximation of a diagram. + * \ingroup Persistence_landscape_on_grid_exact + * + */ + std::vector<std::vector<double> > vectorize() const { + std::vector<std::vector<double> > ls; for(int i = 0; i < nb_ls; i++) ls.emplace_back(); + int num_pts = diagram.size(); double step = (max_x - min_x)/res_x; + + for(int i = 0; i < res_x; i++){ + double x = min_x + i*step; double t = x / std::sqrt(2); std::vector<double> events; + for(int j = 0; j < num_pts; j++){ + double px = diagram[j].first; double py = diagram[j].second; + if(t >= px && t <= py){ if(t >= (px+py)/2) events.push_back(std::sqrt(2)*(py-t)); else events.push_back(std::sqrt(2)*(t-px)); } + } + + std::sort(events.begin(), events.end(), [](const double & a, const double & b){return a > b;}); int nb_events = events.size(); + for (int j = 0; j < nb_ls; j++){ if(j < nb_events) ls[j].push_back(events[j]); else ls[j].push_back(0); } + } + return ls; + } + + + + +}; // class Persistence_landscape_on_grid_exact +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // LANDSCAPE_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h new file mode 100644 index 00000000..76c43e65 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h @@ -0,0 +1,181 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef PERSISTENCE_WEIGHTED_GAUSSIAN_H_ +#define PERSISTENCE_WEIGHTED_GAUSSIAN_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> +#include <gudhi/Weight_functions.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 { +/** + * \class Persistence_weighted_gaussian gudhi/Persistence_weighted_gaussian.h + * \brief A class implementing the Persistence Weighted Gaussian kernel and a specific case thereof called the Persistence Scale Space kernel. + * + * \ingroup Persistence_representations + * + * \details + * The Persistence Weighted Gaussian kernel is built with Gaussian Kernel Mean Embedding, meaning that each persistence diagram is first + * sent to the Hilbert space of a Gaussian kernel with bandwidth parameter \f$\sigma >0\f$ using a weighted mean embedding \f$\Phi\f$: + * + * \f$ \Phi\,:\,D\,\rightarrow\,\sum_{p\in D}\,w(p)\,{\rm exp}\left(-\frac{\|p-\cdot\|_2^2}{2\sigma^2}\right) \f$, + * + * Usually, the weight function is chosen to be an arctan function of the distance of the point to the diagonal: + * \f$w(p) = {\rm arctan}(C\,|y-x|^\alpha)\f$, for some parameters \f$C,\alpha >0\f$. + * Then, their scalar product in this space is computed: + * + * \f$ k(D_1,D_2)=\langle\Phi(D_1),\Phi(D_2)\rangle + * \,=\,\sum_{p\in D_1}\,\sum_{q\in D_2}\,w(p)\,w(q)\,{\rm exp}\left(-\frac{\|p-q\|_2^2}{2\sigma^2}\right).\f$ + * + * Note that one may apply a second Gaussian kernel to their distance in this space and still get a kernel. + * + * It follows that the computation time is \f$O(n^2)\f$ where \f$n\f$ is the number of points + * in the diagrams. This time can be improved by computing approximations of the kernel + * with \f$m\f$ Fourier features \cite Rahimi07randomfeatures. In that case, the computation time becomes \f$O(mn)\f$. + * + * The Persistence Scale Space kernel is a Persistence Weighted Gaussian kernel between modified diagrams: + * the symmetric of each point with respect to the diagonal is first added in each diagram, and then the weight function + * is set to be +1 if the point is above the diagonal and -1 otherwise. + * + * For more details, please see \cite Kusano_Fukumizu_Hiraoka_PWGK + * and \cite Reininghaus_Huber_ALL_PSSK . + * +**/ +class Persistence_weighted_gaussian{ + + protected: + Persistence_diagram diagram; + Weight weight; + double sigma; + int approx; + + public: + + /** \brief Persistence Weighted Gaussian kernel constructor. + * \ingroup Persistence_weighted_gaussian + * + * @param[in] _diagram persistence diagram. + * @param[in] _sigma bandwidth parameter of the Gaussian kernel used for the Kernel Mean Embedding of the diagrams. + * @param[in] _approx number of random Fourier features in case of approximate computation, set to -1 for exact computation. + * @param[in] _weight weight function for the points in the diagrams. + * + */ + Persistence_weighted_gaussian(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 1000, const Weight & _weight = arctan_weight(1,1)){diagram = _diagram; sigma = _sigma; approx = _approx; weight = _weight;} + + + // ********************************** + // Utils. + // ********************************** + + std::vector<std::pair<double,double> > Fourier_feat(const Persistence_diagram & diag, const std::vector<std::pair<double,double> > & z, const Weight & weight = arctan_weight(1,1)) const { + int md = diag.size(); std::vector<std::pair<double,double> > b; int mz = z.size(); + for(int i = 0; i < mz; i++){ + double d1 = 0; double d2 = 0; double zx = z[i].first; double zy = z[i].second; + for(int j = 0; j < md; j++){ + double x = diag[j].first; double y = diag[j].second; + d1 += weight(diag[j])*cos(x*zx + y*zy); + d2 += weight(diag[j])*sin(x*zx + y*zy); + } + b.emplace_back(d1,d2); + } + return b; + } + + std::vector<std::pair<double,double> > random_Fourier(double sigma, int m = 1000) const { + std::normal_distribution<double> distrib(0,1); std::vector<std::pair<double,double> > z; std::random_device rd; + for(int i = 0; i < m; i++){ + std::mt19937 e1(rd()); std::mt19937 e2(rd()); + double zx = distrib(e1); double zy = distrib(e2); + z.emplace_back(zx/sigma,zy/sigma); + } + return z; + } + + + + // ********************************** + // Scalar product + distance. + // ********************************** + + /** \brief Evaluation of the kernel on a pair of diagrams. + * \ingroup Persistence_weighted_gaussian + * + * @pre sigma, approx and weight attributes need to be the same for both instances. + * @param[in] second other instance of class Persistence_weighted_gaussian. + * + */ + double compute_scalar_product(const Persistence_weighted_gaussian & second) const { + + GUDHI_CHECK(this->sigma != second.sigma || this->approx != second.approx || this->weight != second.weight, std::invalid_argument("Error: different values for representations")); + Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram; + + if(this->approx == -1){ + int num_pts1 = diagram1.size(); int num_pts2 = diagram2.size(); double k = 0; + for(int i = 0; i < num_pts1; i++) + for(int j = 0; j < num_pts2; j++) + k += this->weight(diagram1[i])*this->weight(diagram2[j])*exp(-((diagram1[i].first - diagram2[j].first) * (diagram1[i].first - diagram2[j].first) + + (diagram1[i].second - diagram2[j].second) * (diagram1[i].second - diagram2[j].second)) + /(2*this->sigma*this->sigma)); + return k; + } + else{ + std::vector<std::pair<double,double> > z = random_Fourier(this->sigma, this->approx); + std::vector<std::pair<double,double> > b1 = Fourier_feat(diagram1,z,this->weight); + std::vector<std::pair<double,double> > b2 = Fourier_feat(diagram2,z,this->weight); + double d = 0; for(int i = 0; i < this->approx; i++) d += b1[i].first*b2[i].first + b1[i].second*b2[i].second; + return d/this->approx; + } + } + + /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. + * \ingroup Persistence_weighted_gaussian + * + * @pre sigma, approx and weight attributes need to be the same for both instances. + * @param[in] second other instance of class Persistence_weighted_gaussian. + * + */ + double distance(const Persistence_weighted_gaussian & second) const { + GUDHI_CHECK(this->sigma != second.sigma || this->approx != second.approx || this->weight != second.weight, std::invalid_argument("Error: different values for representations")); + return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), 0.5); + } + + +}; // class Persistence_weighted_gaussian +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_WEIGHTED_GAUSSIAN_H_ diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h new file mode 100644 index 00000000..d8ed0d98 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -0,0 +1,340 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef SLICED_WASSERSTEIN_H_ +#define SLICED_WASSERSTEIN_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#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 { + +/** + * \class Sliced_Wasserstein gudhi/Sliced_Wasserstein.h + * \brief A class implementing the Sliced Wasserstein kernel. + * + * \ingroup Persistence_representations + * + * \details + * 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$, + * + * 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. + * The integral can be either computed exactly in \f$O(n^2{\rm log}(n))\f$ time, where \f$n\f$ is the number of points + * in the diagrams, or approximated 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$ + * + * For more details, please see \cite pmlr-v70-carriere17a . + * +**/ + +class Sliced_Wasserstein { + + protected: + Persistence_diagram diagram; + int approx; + double sigma; + std::vector<std::vector<double> > projections, projections_diagonal; + + public: + + 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<double> 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) ); + } + + std::sort(l.begin(), l.end()); std::sort(l_diag.begin(), l_diag.end()); + projections.push_back(l); projections_diagonal.push_back(l_diag); + + } + + } + + } + + /** \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; + x1 = diag[i].first; y1 = diag[i].second; + x2 = diag[j].first; y2 = diag[j].second; + if (y1 - y2 > 0){ + vect.first = y1 - y2; + vect.second = x2 - x1;} + else{ + if(y1 - y2 < 0){ + vect.first = y2 - y1; + vect.second = x1 - x2; + } + else{ + vect.first = 0; + vect.second = abs(x1 - x2);} + } + double norm = std::sqrt(vect.first*vect.first + vect.second*vect.second); + return asin(vect.second/norm); + } + + // 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 >= -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 angle1; + if (diag1[p].first > diag2[q].first) + angle1 = theta1 - asin( (diag1[p].second-diag2[q].second)/norm ); + else + angle1 = theta1 - asin( (diag2[q].second-diag1[p].second)/norm ); + double angle2 = angle1 + theta2 - theta1; + double integral = compute_int_cos(angle1,angle2); + 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. + * + * + */ + 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){ + + // Add projections onto diagonal. + int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double max_ordinate = std::numeric_limits<double>::lowest(); + for (int i = 0; i < n2; i++){ + max_ordinate = std::max(max_ordinate, diagram2[i].second); + diagram1.emplace_back( (diagram2[i].first+diagram2[i].second)/2, (diagram2[i].first+diagram2[i].second)/2 ); + } + for (int i = 0; i < n1; i++){ + max_ordinate = std::max(max_ordinate, diagram1[i].second); + 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. + int mag = 0; while(max_ordinate > 10){mag++; max_ordinate/=10;} + double thresh = pow(10,-5+mag); + srand(time(NULL)); + for (int i = 0; i < num_pts_dgm; i++){ + diagram1[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); diagram1[i].second += thresh*(1.0-2.0*rand()/RAND_MAX); + diagram2[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); diagram2[i].second += thresh*(1.0-2.0*rand()/RAND_MAX); + } + + // Compute all angles in both PDs. + std::vector<std::pair<double, std::pair<int,int> > > 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<int,int>(i,j)); + angles2.emplace_back(theta2, std::pair<int,int>(i,j)); + } + } + + // 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);}); + + // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2). + std::vector<int> 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); } ); + + // Find the inverses of the orders. + std::vector<int> order1(num_pts_dgm); std::vector<int> order2(num_pts_dgm); + for(int i = 0; i < num_pts_dgm; i++) for (int j = 0; j < num_pts_dgm; j++) if(orderp1[j] == i){ order1[i] = j; break; } + for(int i = 0; i < num_pts_dgm; i++) for (int j = 0; j < num_pts_dgm; j++) if(orderp2[j] == i){ order2[i] = j; break; } + + // Record all inversions of points in the orders as theta varies along the positive half-disk. + std::vector<std::vector<std::pair<int,double> > > anglePerm1(num_pts_dgm); + std::vector<std::vector<std::pair<int,double> > > 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; + } + + 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 < 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<std::pair<int,double> > 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++; + theta2 = std::min(u[ku].second, v[kv].second); + } + } + } + + + else{ + + double step = pi/this->approx; + 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)); + 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; + } + + /** \brief Evaluation of the kernel on a pair of diagrams. + * \ingroup Sliced_Wasserstein + * + * @pre approx and sigma attributes need to be the same for both instances. + * @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)); + } + + /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. + * \ingroup Sliced_Wasserstein + * + * @pre approx and sigma attributes need to be the same for both instances. + * @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); + } + + + + +}; // class Sliced_Wasserstein +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // SLICED_WASSERSTEIN_H_ diff --git a/src/Persistence_representations/include/gudhi/Weight_functions.h b/src/Persistence_representations/include/gudhi/Weight_functions.h new file mode 100644 index 00000000..78de406d --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Weight_functions.h @@ -0,0 +1,81 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef WEIGHT_FUNCTIONS_H_ +#define WEIGHT_FUNCTIONS_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.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 { + +/** \fn static double pss_weight(std::pair<double,double> p) + * \brief Persistence Scale Space kernel weight function. + * \ingroup Persistence_representations + * + * @param[in] p point in 2D. + */ +static double pss_weight(std::pair<double,double> p) {if(p.second > p.first) return 1; else return -1;} + +/** \fn static double linear_weight(std::pair<double,double> p) + * \brief Linear weight function. + * \ingroup Persistence_representations + * + * @param[in] p point in 2D. + */ +static double linear_weight(std::pair<double,double> p) {return std::abs(p.second - p.first);} + +/** \fn static double const_weight(std::pair<double,double> p) + * \brief Constant weight function. + * \ingroup Persistence_representations + * + * @param[in] p point in 2D. + */ +static double const_weight(std::pair<double,double> p) {return 1;} + +/** \fn static std::function<double (std::pair<double,double>) > arctan_weight(double C, double alpha) + * \brief Returns the arctan weight function with parameters C and alpha. + * \ingroup Persistence_representations + * + * @param[in] C positive constant. + * @param[in] alpha positive power. + */ +static std::function<double (std::pair<double,double>) > arctan_weight(double C, double alpha) {return [=](std::pair<double,double> p){return C * atan(std::pow(std::abs(p.second - p.first), alpha));};} + +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // WEIGHT_FUNCTIONS_H_ diff --git a/src/Persistence_representations/include/gudhi/common_persistence_representations.h b/src/Persistence_representations/include/gudhi/common_persistence_representations.h index 3d03f1f6..539eee60 100644 --- a/src/Persistence_representations/include/gudhi/common_persistence_representations.h +++ b/src/Persistence_representations/include/gudhi/common_persistence_representations.h @@ -26,17 +26,32 @@ #include <utility> #include <string> #include <cmath> +#include <boost/math/constants/constants.hpp> + + namespace Gudhi { namespace Persistence_representations { // this file contain an implementation of some common procedures used in Persistence_representations. +static constexpr double pi = boost::math::constants::pi<double>(); + + +/** + * In this module, we use the name Persistence_diagram for the representation of a diagram in a vector of pairs of two 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. + */ +using Weight = std::function<double (std::pair<double,double>) >; // double epsi = std::numeric_limits<double>::epsilon(); double epsi = 0.000005; /** * A procedure used to compare doubles. Typically given two doubles A and B, comparing A == B is not good idea. In this - *case, we use the procedure almostEqual with the epsi defined at + * case, we use the procedure almostEqual with the epsi defined at * the top of the file. Setting up the epsi gives the user a tolerance on what should be consider equal. **/ inline bool almost_equal(double a, double b) { @@ -53,8 +68,7 @@ double birth_plus_deaths(std::pair<double, double> a) { return a.first + a.secon // landscapes /** - * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two - *points. + * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two points. **/ std::pair<double, double> compute_parameters_of_a_line(std::pair<double, double> p1, std::pair<double, double> p2) { double a = (p2.second - p1.second) / (p2.first - p1.first); @@ -64,8 +78,7 @@ std::pair<double, double> compute_parameters_of_a_line(std::pair<double, double> // landscapes /** - * This procedure given two points which lies on the opposite sides of x axis, compute x for which the line connecting - *those two points crosses x axis. + * This procedure given two points which lies on the opposite sides of x axis, compute x for which the line connecting those two points crosses x axis. **/ double find_zero_of_a_line_segment_between_those_two_points(std::pair<double, double> p1, std::pair<double, double> p2) { @@ -89,8 +102,7 @@ double find_zero_of_a_line_segment_between_those_two_points(std::pair<double, do // landscapes /** * This method provides a comparison of points that is used in construction of persistence landscapes. The ordering is - *lexicographical for the first coordinate, and reverse-lexicographical for the - * second coordinate. + * lexicographical for the first coordinate, and reverse-lexicographical for the second coordinate. **/ bool compare_points_sorting(std::pair<double, double> f, std::pair<double, double> s) { if (f.first < s.first) { diff --git a/src/Persistence_representations/test/CMakeLists.txt b/src/Persistence_representations/test/CMakeLists.txt index 335a71ef..ca1713c3 100644 --- a/src/Persistence_representations/test/CMakeLists.txt +++ b/src/Persistence_representations/test/CMakeLists.txt @@ -35,6 +35,11 @@ target_link_libraries(Read_persistence_from_file_test_unit ${Boost_UNIT_TEST_FRA gudhi_add_coverage_test(Read_persistence_from_file_test_unit) +add_executable ( kernels_unit kernels.cpp ) +target_link_libraries(kernels_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +gudhi_add_coverage_test(kernels_unit) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) add_executable (Persistence_intervals_with_distances_test_unit persistence_intervals_with_distances_test.cpp ) target_link_libraries(Persistence_intervals_with_distances_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) diff --git a/src/Persistence_representations/test/kernels.cpp b/src/Persistence_representations/test/kernels.cpp new file mode 100644 index 00000000..9db19123 --- /dev/null +++ b/src/Persistence_representations/test/kernels.cpp @@ -0,0 +1,54 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carrière + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "kernel" + +#include <boost/test/unit_test.hpp> +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <algorithm> // std::max +#include <gudhi/common_persistence_representations.h> +#include <gudhi/Weight_functions.h> +#include <gudhi/Persistence_weighted_gaussian.h> +#include <gudhi/Sliced_Wasserstein.h> +#include <gudhi/distance_functions.h> +#include <gudhi/reader_utils.h> + +using SW = Gudhi::Persistence_representations::Sliced_Wasserstein; +using PWG = Gudhi::Persistence_representations::Persistence_weighted_gaussian; + +BOOST_AUTO_TEST_CASE(check_PWG) { + Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); + PWG pwg1(v1, 1.0, 1000, Gudhi::Persistence_representations::arctan_weight(1,1)); PWG pwgex1(v1, 1.0, -1, Gudhi::Persistence_representations::arctan_weight(1,1)); + PWG pwg2(v2, 1.0, 1000, Gudhi::Persistence_representations::arctan_weight(1,1)); PWG pwgex2(v2, 1.0, -1, Gudhi::Persistence_representations::arctan_weight(1,1)); + BOOST_CHECK(std::abs(pwg1.compute_scalar_product(pwg2) - pwgex1.compute_scalar_product(pwgex2)) <= 1e-1); +} + +BOOST_AUTO_TEST_CASE(check_SW) { + Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); + SW sw1(v1, 1.0, 100); SW swex1(v1, 1.0, -1); + SW sw2(v2, 1.0, 100); SW swex2(v2, 1.0, -1); + BOOST_CHECK(std::abs(sw1.compute_scalar_product(sw2) - swex1.compute_scalar_product(swex2)) <= 1e-1); +} |