From 541284f6f1bf7d4a76daac8a52850c7162a765cb Mon Sep 17 00:00:00 2001 From: mcarrier Date: Mon, 23 Apr 2018 15:22:13 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3387 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 3fe2ae4af0c7cadf507fc5148c05dcf664c5e151 --- .../doc/Persistence_representations_doc.h | 5 +- .../example/persistence_image.cpp | 2 +- .../example/persistence_weighted_gaussian.cpp | 24 +- .../example/sliced_wasserstein.cpp | 3 +- .../include/gudhi/Landscape.h | 87 ++++---- .../include/gudhi/Persistence_image.h | 41 ++-- .../include/gudhi/Persistence_weighted_gaussian.h | 70 +++--- .../include/gudhi/Sliced_Wasserstein.h | 243 ++++++++++----------- .../gudhi/common_persistence_representations.h | 22 +- 9 files changed, 242 insertions(+), 255 deletions(-) (limited to 'src/Persistence_representations') diff --git a/src/Persistence_representations/doc/Persistence_representations_doc.h b/src/Persistence_representations/doc/Persistence_representations_doc.h index 6d4cc96c..ca283017 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 @@ -254,11 +253,11 @@ namespace Persistence_representations { -\section sec_persistence_kernels Kernels on Persistence Diagrams +\section sec_persistence_kernels Kernels on persistence diagrams Reference manual: \ref Gudhi::Persistence_representations::Sliced_Wasserstein
Reference manual: \ref Gudhi::Persistence_representations::Persistence_weighted_gaussian
- Kernels for Persistence Diagrams can be regarded as infinite-dimensional vectorizations. More specifically, + 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 positive semi-definite. diff --git a/src/Persistence_representations/example/persistence_image.cpp b/src/Persistence_representations/example/persistence_image.cpp index dfa469d4..cdce3bbf 100644 --- a/src/Persistence_representations/example/persistence_image.cpp +++ b/src/Persistence_representations/example/persistence_image.cpp @@ -40,7 +40,7 @@ int main(int argc, char** argv) { 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::Persistence_weighted_gaussian::linear_weight; + 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 > P = pim.vectorize(); diff --git a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp index 234f6323..db60755f 100644 --- a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp +++ b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp @@ -26,13 +26,11 @@ #include #include -using PD = std::vector >; using PWG = Gudhi::Persistence_representations::Persistence_weighted_gaussian; int main(int argc, char** argv) { - std::vector > persistence1; - std::vector > persistence2; + Persistence_diagram persistence1, persistence2; persistence1.push_back(std::make_pair(1, 2)); persistence1.push_back(std::make_pair(6, 8)); @@ -48,11 +46,11 @@ int main(int argc, char** argv) { double tau = 1; int m = 10000; - PWG PWG1(persistence1, sigma, m, PWG::arctan_weight(1,1)); - PWG PWG2(persistence2, sigma, m, PWG::arctan_weight(1,1)); + 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, PWG::arctan_weight(1,1)); - PWG PWGex2(persistence2, sigma, -1, PWG::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 @@ -82,14 +80,14 @@ int main(int argc, char** argv) { // PSS - PD pd1 = persistence1; int numpts = persistence1.size(); for(int i = 0; i < numpts; i++) pd1.emplace_back(persistence1[i].second,persistence1[i].first); - PD pd2 = persistence2; numpts = persistence2.size(); for(int i = 0; i < numpts; i++) pd2.emplace_back(persistence2[i].second,persistence2[i].first); + 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, PWG::pss_weight); - PWG pwg2(pd2, 2*std::sqrt(sigma), m, PWG::pss_weight); + 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, PWG::pss_weight); - PWG pwgex2(pd2, 2*std::sqrt(sigma), -1, PWG::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; diff --git a/src/Persistence_representations/example/sliced_wasserstein.cpp b/src/Persistence_representations/example/sliced_wasserstein.cpp index f1aeea5c..d37cb23c 100644 --- a/src/Persistence_representations/example/sliced_wasserstein.cpp +++ b/src/Persistence_representations/example/sliced_wasserstein.cpp @@ -30,8 +30,7 @@ using SW = Gudhi::Persistence_representations::Sliced_Wasserstein; int main(int argc, char** argv) { - std::vector > persistence1; - std::vector > persistence2; + Persistence_diagram persistence1, persistence2; persistence1.push_back(std::make_pair(1, 2)); persistence1.push_back(std::make_pair(6, 8)); diff --git a/src/Persistence_representations/include/gudhi/Landscape.h b/src/Persistence_representations/include/gudhi/Landscape.h index d6608a57..bbbca36b 100644 --- a/src/Persistence_representations/include/gudhi/Landscape.h +++ b/src/Persistence_representations/include/gudhi/Landscape.h @@ -40,64 +40,69 @@ #include #include -using PD = std::vector >; - namespace Gudhi { namespace Persistence_representations { /** * \class Landscape gudhi/Landscape.h - * \brief A class implementing the Landscapes. + * \brief A class implementing landscapes. * * \ingroup Persistence_representations * * \details * + * The landscape is a way to turn a persistence diagram into \f$L^2\f$ functions. Roughly, the idea is to see the boundaries of the rank functions as scalar functions taking values on the diagonal. + * See \cite bubenik_landscapes_2015 for more details. Here we provide a way to approximate such functions by computing their values on a set of samples. + * **/ class Landscape { - protected: - PD diagram; - int res_x, nb_ls; - double min_x, max_x; + protected: + Persistence_diagram diagram; + int res_x, nb_ls; + double min_x, max_x; public: - /** \brief Landscape constructor. - * \ingroup Landscape - * - */ - Landscape(PD _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 of a diagram. - * \ingroup Landscape - * - */ - std::vector > vectorize() { - std::vector > 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 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; - } - - - - -}; - -} // namespace Landscape + /** \brief Landscape constructor. + * \ingroup Landscape + * + * @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. + * + */ + Landscape(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 of a diagram. + * \ingroup Landscape + * + */ + std::vector > vectorize() const { + std::vector > 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 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 Landscape +} // namespace Persistence_representations } // namespace Gudhi #endif // LANDSCAPE_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_image.h b/src/Persistence_representations/include/gudhi/Persistence_image.h index 6c9f75b7..76b34d8d 100644 --- a/src/Persistence_representations/include/gudhi/Persistence_image.h +++ b/src/Persistence_representations/include/gudhi/Persistence_image.h @@ -26,8 +26,8 @@ // gudhi include #include #include +#include #include -#include // standard include #include @@ -41,39 +41,49 @@ #include #include -using PD = std::vector >; -using Weight = std::function) >; - namespace Gudhi { namespace Persistence_representations { /** * \class Persistence_image gudhi/Persistence_image.h - * \brief A class implementing the Persistence Images. + * \brief A class implementing the persistence images. * * \ingroup Persistence_representations * * \details * + * Persistence images are a way to build images from persistence diagrams. Roughly, the idea is to center Gaussians on each diagram point, with a weight that usually depends on + * the distance to the diagonal, so that the diagram is turned into a function, and then to discretize the plane into pixels, and integrate this function on each pixel. + * See \cite Persistence_Images_2017 for more details. + * **/ class Persistence_image { protected: - PD diagram; - int res_x, res_y; - double min_x, max_x, min_y, max_y; - Weight weight; - double sigma; + Persistence_diagram diagram; + int res_x, res_y; + double min_x, max_x, min_y, max_y; + Weight weight; + double sigma; public: /** \brief Persistence Image constructor. * \ingroup Persistence_image * + * @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_image(PD _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, - Weight _weight = Gudhi::Persistence_representations::Persistence_weighted_gaussian::arctan_weight(1,1), double _sigma = 1.0){ + Persistence_image(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; } @@ -81,7 +91,7 @@ class Persistence_image { * \ingroup Persistence_image * */ - std::vector > vectorize() { + std::vector > vectorize() const { std::vector > im; for(int i = 0; i < res_y; i++) im.emplace_back(); double step_x = (max_x - min_x)/res_x; double step_y = (max_y - min_y)/res_y; @@ -109,9 +119,8 @@ class Persistence_image { -}; - -} // namespace Persistence_image +}; // class Persistence_image +} // namespace Persistence_representations } // namespace Gudhi #endif // PERSISTENCE_IMAGE_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h index 9a63fccd..76c43e65 100644 --- a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h +++ b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h @@ -26,6 +26,7 @@ // gudhi include #include #include +#include // standard include #include @@ -39,19 +40,16 @@ #include #include -using PD = std::vector >; -using Weight = std::function) >; - 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 of it called the Persistence Scale Space Kernel. + * \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 + * 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$, @@ -69,59 +67,41 @@ namespace Persistence_representations { * 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 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 consult Persistence Weighted Kernel for Topological Data Analysis\cite Kusano_Fukumizu_Hiraoka_PWGK - * and A Stable Multi-Scale Kernel for Topological Machine Learning\cite Reininghaus_Huber_ALL_PSSK . - * It implements the following concepts: Topological_data_with_distances, Topological_data_with_scalar_product. + * For more details, please see \cite Kusano_Fukumizu_Hiraoka_PWGK + * and \cite Reininghaus_Huber_ALL_PSSK . * **/ class Persistence_weighted_gaussian{ protected: - PD diagram; + Persistence_diagram diagram; Weight weight; double sigma; int approx; public: - /** \brief Persistence Weighted Gaussian Kernel constructor. + /** \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] _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(PD _diagram, double _sigma = 1.0, int _approx = 1000, Weight _weight = arctan_weight(1,1)){diagram = _diagram; sigma = _sigma; approx = _approx; weight = _weight;} - - PD get_diagram() const {return this->diagram;} - double get_sigma() const {return this->sigma;} - int get_approx() const {return this->approx;} - Weight get_weight() const {return this->weight;} + 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. // ********************************** - /** \brief Specific weight of Persistence Scale Space Kernel. - * \ingroup Persistence_weighted_gaussian - * - * @param[in] p point in 2D. - * - */ - static double pss_weight(std::pair p) {if(p.second > p.first) return 1; else return -1;} - static double linear_weight(std::pair p) {return std::abs(p.second - p.first);} - static double const_weight(std::pair p) {return 1;} - static std::function) > arctan_weight(double C, double power) {return [=](std::pair p){return C * atan(std::pow(std::abs(p.second - p.first), power));};} - - - std::vector > Fourier_feat(PD diag, std::vector > z, Weight weight = arctan_weight(1,1)){ + std::vector > Fourier_feat(const Persistence_diagram & diag, const std::vector > & z, const Weight & weight = arctan_weight(1,1)) const { int md = diag.size(); std::vector > 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; @@ -135,7 +115,7 @@ class Persistence_weighted_gaussian{ return b; } - std::vector > random_Fourier(double sigma, int m = 1000){ + std::vector > random_Fourier(double sigma, int m = 1000) const { std::normal_distribution distrib(0,1); std::vector > z; std::random_device rd; for(int i = 0; i < m; i++){ std::mt19937 e1(rd()); std::mt19937 e2(rd()); @@ -154,12 +134,14 @@ class Persistence_weighted_gaussian{ /** \brief Evaluation of the kernel on a pair of diagrams. * \ingroup Persistence_weighted_gaussian * - * @param[in] second other instance of class Persistence_weighted_gaussian. Warning: sigma, approx and weight parameters need to be the same for both instances!!! + * @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(Persistence_weighted_gaussian second){ + double compute_scalar_product(const Persistence_weighted_gaussian & second) const { - PD diagram1 = this->diagram; PD diagram2 = second.diagram; + 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; @@ -171,7 +153,7 @@ class Persistence_weighted_gaussian{ return k; } else{ - std::vector > z = random_Fourier(this->sigma, this->approx); + std::vector > z = random_Fourier(this->sigma, this->approx); std::vector > b1 = Fourier_feat(diagram1,z,this->weight); std::vector > 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; @@ -182,20 +164,18 @@ class Persistence_weighted_gaussian{ /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. * \ingroup Persistence_weighted_gaussian * - * @param[in] second other instance of class Persistence_weighted_gaussian. Warning: sigma, approx and weight parameters need to be the same for both instances!!! + * @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(Persistence_weighted_gaussian second) { - if(this->sigma != second.get_sigma() || this->approx != second.get_approx()){ - std::cout << "Error: different representations!" << std::endl; return 0; - } - else return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), 0.5); + 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); } -}; - -} // namespace Persistence_weighted_gaussian +}; // 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 index 6a9a607e..235918fe 100644 --- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -40,19 +40,17 @@ #include #include -using PD = std::vector >; - namespace Gudhi { namespace Persistence_representations { /** * \class Sliced_Wasserstein gudhi/Sliced_Wasserstein.h - * \brief A class implementing the Sliced Wasserstein Kernel. + * \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 + * 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: * @@ -65,15 +63,14 @@ namespace Persistence_representations { * * \f$ k(D_1,D_2) = {\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right).\f$ * - * For more details, please consult Sliced Wasserstein Kernel for Persistence Diagrams\cite pmlr-v70-carriere17a . - * It implements the following concepts: Topological_data_with_distances, Topological_data_with_scalar_product. + * For more details, please see \cite pmlr-v70-carriere17a . * **/ class Sliced_Wasserstein { protected: - PD diagram; + Persistence_diagram diagram; int approx; double sigma; std::vector > projections, projections_diagonal; @@ -107,7 +104,7 @@ class Sliced_Wasserstein { } - /** \brief Sliced Wasserstein Kernel constructor. + /** \brief Sliced Wasserstein kernel constructor. * \ingroup Sliced_Wasserstein * * @param[in] _diagram persistence diagram. @@ -115,21 +112,14 @@ class Sliced_Wasserstein { * @param[in] _approx number of directions used to approximate the integral in the Sliced Wasserstein distance, set to -1 for exact computation. * */ - Sliced_Wasserstein(PD _diagram, double _sigma = 1.0, int _approx = 100){diagram = _diagram; approx = _approx; sigma = _sigma; build_rep();} - - PD get_diagram() const {return this->diagram;} - int get_approx() const {return this->approx;} - double get_sigma() const {return this->sigma;} - - - + 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(PD diag, int i, int j){ + double compute_angle(const Persistence_diagram & diag, int i, int j) const { std::pair vect; double x1,y1, x2,y2; x1 = diag[i].first; y1 = diag[i].second; x2 = diag[j].first; y2 = diag[j].second; @@ -150,7 +140,7 @@ class Sliced_Wasserstein { } // 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(const double & alpha, const double & beta){ + double compute_int_cos(double alpha, double beta) const { double res = 0; if (alpha >= 0 && alpha <= pi){ if (cos(alpha) >= 0){ @@ -175,13 +165,13 @@ class Sliced_Wasserstein { return res; } - double compute_int(const double & theta1, const double & theta2, const int & p, const int & q, const PD & PD1, const PD & PD2){ - double norm = std::sqrt( (PD1[p].first-PD2[q].first)*(PD1[p].first-PD2[q].first) + (PD1[p].second-PD2[q].second)*(PD1[p].second-PD2[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 (PD1[p].first > PD2[q].first) - angle1 = theta1 - asin( (PD1[p].second-PD2[q].second)/norm ); + if (diag1[p].first > diag2[q].first) + angle1 = theta1 - asin( (diag1[p].second-diag2[q].second)/norm ); else - angle1 = theta1 - asin( (PD2[q].second-PD1[p].second)/norm ); + 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; @@ -197,134 +187,133 @@ class Sliced_Wasserstein { /** \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. - * For warning in red: - * @warning approx parameter needs to be the same for both instances. + * * */ - double compute_sliced_wasserstein_distance(Sliced_Wasserstein second) { + 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")); + GUDHI_CHECK(this->approx != second.approx, std::invalid_argument("Error: different approx values for representations")); - PD diagram1 = this->diagram; PD diagram2 = second.diagram; double sw = 0; + Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram; double sw = 0; - if(this->approx == -1){ + if(this->approx == -1){ - // Add projections onto diagonal. - int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double max_ordinate = std::numeric_limits::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); - } + // Add projections onto diagonal. + int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double max_ordinate = std::numeric_limits::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 > > 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)); - } + // 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)); } + } - // Sort angles. - std::sort(angles1.begin(), angles1.end(), [=](std::pair >& p1, const std::pair >& p2){return (p1.first < p2.first);}); - std::sort(angles2.begin(), angles2.end(), [=](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); } ); - - // 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++) 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 > > 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; - } + // Sort angles. + std::sort(angles1.begin(), angles1.end(), [=](std::pair >& p1, const std::pair >& p2){return (p1.first < p2.first);}); + std::sort(angles2.begin(), angles2.end(), [=](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); } ); + + // 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++) 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 > > 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; + } - 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; - } + 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); - } + 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); - 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); - } + // 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); + 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++){ + else{ - std::vector v1; std::vector l1 = this->projections[i]; std::vector l1bis = second.projections_diagonal[i]; std::merge(l1.begin(), l1.end(), l1bis.begin(), l1bis.end(), std::back_inserter(v1)); - std::vector v2; std::vector l2 = second.projections[i]; std::vector 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; + double step = pi/this->approx; + for (int i = 0; i < this->approx; i++){ - } + std::vector v1; std::vector l1 = this->projections[i]; std::vector l1bis = second.projections_diagonal[i]; std::merge(l1.begin(), l1.end(), l1bis.begin(), l1bis.end(), std::back_inserter(v1)); + std::vector v2; std::vector l2 = second.projections[i]; std::vector 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; + return sw/pi; } /** \brief Evaluation of the kernel on a pair of diagrams. * \ingroup Sliced_Wasserstein * - * @param[in] second other instance of class Sliced_Wasserstein. Warning: sigma and approx parameters need to be the same for both instances!!! + * @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(Sliced_Wasserstein second){ + 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)); } @@ -332,10 +321,11 @@ class Sliced_Wasserstein { /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. * \ingroup Sliced_Wasserstein * - * @param[in] second other instance of class Sliced_Wasserstein. Warning: sigma and approx parameters need to be the same for both instances!!! + * @pre approx and sigma attributes need to be the same for both instances. + * @param[in] second other instance of class Sliced_Wasserstein. * */ - double distance(Sliced_Wasserstein second) { + 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); } @@ -343,9 +333,8 @@ class Sliced_Wasserstein { -}; - -} // namespace Sliced_Wasserstein +}; // class Sliced_Wasserstein +} // namespace Persistence_representations } // namespace Gudhi #endif // SLICED_WASSERSTEIN_H_ diff --git a/src/Persistence_representations/include/gudhi/common_persistence_representations.h b/src/Persistence_representations/include/gudhi/common_persistence_representations.h index 90f2626d..884fce58 100644 --- a/src/Persistence_representations/include/gudhi/common_persistence_representations.h +++ b/src/Persistence_representations/include/gudhi/common_persistence_representations.h @@ -28,17 +28,28 @@ #include #include +/** + * 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 >; + +/** + * 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) >; + 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 epsi = std::numeric_limits::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) { @@ -55,8 +66,7 @@ double birth_plus_deaths(std::pair 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 compute_parameters_of_a_line(std::pair p1, std::pair p2) { double a = (p2.second - p1.second) / (p2.first - p1.first); @@ -66,8 +76,7 @@ std::pair compute_parameters_of_a_line(std::pair // 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 p1, std::pair p2) { @@ -91,8 +100,7 @@ double find_zero_of_a_line_segment_between_those_two_points(std::pair f, std::pair s) { if (f.first < s.first) { -- cgit v1.2.3