From 82dca01519298502af4d478f47cd8e42d0d9478f Mon Sep 17 00:00:00 2001 From: mcarrier Date: Wed, 3 Jan 2018 15:30:49 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3110 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 95fd5b9e4cd4738ecdeaeb27f20d267816b7bead --- biblio/bibliography.bib | 22 +++++++ data/persistence_diagram/PD1 | 3 + data/persistence_diagram/PD2 | 2 + src/Kernels/doc/Intro_kernels.h | 53 ++++++++++++++++- src/Kernels/example/kernel.txt | 8 +++ src/Kernels/example/kernel_basic_example.cpp | 48 ++++++++++------ src/Kernels/include/gudhi/kernel.h | 86 +++++++++++++++++----------- src/Kernels/test/CMakeLists.txt | 14 +++++ src/Kernels/test/test_kernel.cpp | 56 ++++++++++++++++++ 9 files changed, 239 insertions(+), 53 deletions(-) create mode 100644 data/persistence_diagram/PD1 create mode 100644 data/persistence_diagram/PD2 create mode 100644 src/Kernels/example/kernel.txt create mode 100644 src/Kernels/test/CMakeLists.txt create mode 100644 src/Kernels/test/test_kernel.cpp diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index b101cb76..e56734e4 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -1072,3 +1072,25 @@ language={English} + +@InProceedings{pmlr-v70-carriere17a, + title = {Sliced {W}asserstein Kernel for Persistence Diagrams}, + author = {Mathieu Carri{\`e}re and Marco Cuturi and Steve Oudot}, + booktitle = {Proceedings of the 34th International Conference on Machine Learning}, + pages = {664--673}, + year = {2017}, + editor = {Doina Precup and Yee Whye Teh}, + volume = {70}, + series = {Proceedings of Machine Learning Research}, + address = {International Convention Centre, Sydney, Australia}, + month = {06--11 Aug}, + publisher = {PMLR}, +} + +@INPROCEEDINGS{Rahimi07randomfeatures, + author = {Ali Rahimi and Ben Recht}, + title = {Random features for large-scale kernel machines}, + booktitle = {In Neural Information Processing Systems}, + year = {2007} +} + diff --git a/data/persistence_diagram/PD1 b/data/persistence_diagram/PD1 new file mode 100644 index 00000000..404199b4 --- /dev/null +++ b/data/persistence_diagram/PD1 @@ -0,0 +1,3 @@ +2.7 3.7 +9.6 14 +34.2 34.974 \ No newline at end of file diff --git a/data/persistence_diagram/PD2 b/data/persistence_diagram/PD2 new file mode 100644 index 00000000..125d8e4b --- /dev/null +++ b/data/persistence_diagram/PD2 @@ -0,0 +1,2 @@ +2.8 4.45 +9.5 14.1 \ No newline at end of file diff --git a/src/Kernels/doc/Intro_kernels.h b/src/Kernels/doc/Intro_kernels.h index be97a6cf..163690b1 100644 --- a/src/Kernels/doc/Intro_kernels.h +++ b/src/Kernels/doc/Intro_kernels.h @@ -37,17 +37,64 @@ namespace kernel { * to the scalar products of the images of the diagrams under some feature map into a (generally unknown and infinite dimensional) * Hilbert space. Kernels are * very useful to handle any type of data for algorithms that require at least a Hilbert structure, such as Principal Component Analysis - * or Support Vector Machines. In this package, we implement three kernels for persistence diagrams: the Persistence Scale Space kernel, - * the Persistence Weighted Gaussian kernel and the Sliced Wasserstein kernel. + * or Support Vector Machines. In this package, we implement three kernels for persistence diagrams: + * the Persistence Scale Space Kernel (PSSK)---see \cite Reininghaus_Huber_ALL_PSSK, + * the Persistence Weighted Gaussian Kernel (PWGK)---see \cite Kusano_Fukumizu_Hiraoka_PWGK, + * and the Sliced Wasserstein Kernel (SWK)---see \cite pmlr-v70-carriere17a. * + * \section pwg Persistence Weighted Gaussian Kernel and Persistence Scale Space Kernel + * + * The PWGK 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, either their scalar product in this space is + * computed (Linear Persistence Weighted Gaussian Kernel): + * + * \f$ LPWGK(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$, + * + * or a second Gaussian kernel with bandwidth parameter \f$\tau >0\f$ is applied to their distance in this space + * (Gaussian Persistence Weighted Gaussian Kernel): + * + * \f$ GPWGK(D_1,D_2)={\rm exp}\left(-\frac{\|\Phi(D_1)-\Phi(D_2)\|^2}{2\tau^2} \right)\f$, + * where \f$\|\Phi(D_1)-\Phi(D_2)\|^2 = \langle\Phi(D_1)-\Phi(D_2),\Phi(D_1)-\Phi(D_2)\rangle\f$. + * + * 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 PSSK is a Linear 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. + * + * \section sw Sliced Wasserstein Kernel + * + * 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$m\f$ lines in the circle in \f$O(mn{\rm log}(n))\f$ time. The SWK is then computed as: + * + * \f$ SWK(D_1,D_2) = {\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right).\f$ * * When launching: * - * \code $> ./BasicEx + * \code $> ./BasicEx ../../../../data/persistence_diagram/PD1 ../../../../data/persistence_diagram/PD2 * \endcode * * the program output is: * + * \include Kernels/kernel.txt + * * * \copyright GNU General Public License v3. * \verbatim Contact: gudhi-users@lists.gforge.inria.fr \endverbatim diff --git a/src/Kernels/example/kernel.txt b/src/Kernels/example/kernel.txt new file mode 100644 index 00000000..5fb8b504 --- /dev/null +++ b/src/Kernels/example/kernel.txt @@ -0,0 +1,8 @@ +SWK exact = 0.875446 +SWK approx = 0.875204 +PSSK exact = 0.0218669 +PSSK approx = 0.0213766 +LPWGK exact = 2.57351 +LPWGK approx = 2.49102 +GPWGK exact = 0.98783 +GPWGK approx = 0.987591 \ No newline at end of file diff --git a/src/Kernels/example/kernel_basic_example.cpp b/src/Kernels/example/kernel_basic_example.cpp index 8e9925c5..46e42c9d 100644 --- a/src/Kernels/example/kernel_basic_example.cpp +++ b/src/Kernels/example/kernel_basic_example.cpp @@ -20,29 +20,41 @@ * along with this program. If not, see . */ -#define NUMPI 3.14159265359 #include -int main() { +void usage(int nbArgs, char *const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " PD1 PD2 \n"; + std::cerr << " i.e.: " << progName << " ../../../../data/persistence_diagram/PD1 ../../../../data/persistence_diagram/PD2 \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { - std::vector< std::pair > v1, v2; + if (argc != 3) usage(argc, argv[0]); double sigma = 2; double tau = 5; - v1.emplace_back(2.7, 3.7); - v1.emplace_back(9.6, 14.); - v1.emplace_back(34.2, 34.974); - - v2.emplace_back(2.8, 4.45); - v2.emplace_back(9.5, 14.1); - - std::cout << "SW exact = " << Gudhi::kernel::sw (v1, v2) << std::endl; - std::cout << "SW approx = " << Gudhi::kernel::approx_sw (v1, v2) << std::endl; - std::cout << "PSS exact = " << Gudhi::kernel::pss (v1,v2,sigma) << std::endl; - std::cout << "PSS approx = " << Gudhi::kernel::approx_pss (v1,v2,sigma) << std::endl; - std::cout << "PWG exact = " << Gudhi::kernel::lpwg (v1,v2,sigma) << std::endl; - std::cout << "PWG approx = " << Gudhi::kernel::approx_lpwg (v1,v2,sigma) << std::endl; - std::cout << "GPWG exact = " << Gudhi::kernel::gpwg (v1,v2,sigma,tau) << std::endl; - std::cout << "GPWG approx = " << Gudhi::kernel::approx_gpwg (v1,v2,sigma,tau) << std::endl; + std::string PDname1(argv[1]); std::string PDname2(argv[2]); + std::vector< std::pair > v1, v2; std::string line; double b,d; + + std::ifstream input1(PDname1); + while(std::getline(input1,line)){ + std::stringstream stream(line); stream >> b; stream >> d; v1.push_back(std::pair(b,d)); + } + + std::ifstream input2(PDname2); + while(std::getline(input2,line)){ + std::stringstream stream(line); stream >> b; stream >> d; v2.push_back(std::pair(b,d)); + } + + std::cout << "SWK exact = " << Gudhi::kernel::swk (v1,v2,sigma) << std::endl; + std::cout << "SWK approx = " << Gudhi::kernel::approx_swk (v1,v2,sigma) << std::endl; + std::cout << "PSSK exact = " << Gudhi::kernel::pssk (v1,v2,sigma) << std::endl; + std::cout << "PSSK approx = " << Gudhi::kernel::approx_pssk (v1,v2,sigma) << std::endl; + std::cout << "LPWGK exact = " << Gudhi::kernel::lpwgk (v1,v2,sigma) << std::endl; + std::cout << "LPWGK approx = " << Gudhi::kernel::approx_lpwgk (v1,v2,sigma) << std::endl; + std::cout << "GPWGK exact = " << Gudhi::kernel::gpwgk (v1,v2,sigma,tau) << std::endl; + std::cout << "GPWGK approx = " << Gudhi::kernel::approx_gpwgk (v1,v2,sigma,tau) << std::endl; } diff --git a/src/Kernels/include/gudhi/kernel.h b/src/Kernels/include/gudhi/kernel.h index 6429efed..44d984bd 100644 --- a/src/Kernels/include/gudhi/kernel.h +++ b/src/Kernels/include/gudhi/kernel.h @@ -23,8 +23,6 @@ #ifndef KERNEL_H_ #define KERNEL_H_ -#define NUMPI 3.14159265359 - #include #include #include @@ -80,7 +78,7 @@ double pss_weight(std::pair P){ * @param[in] weight weight function for the points in the diagrams. * */ -double lpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ +double lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ int num_pts1 = PD1.size(); int num_pts2 = PD2.size(); double k = 0; for(int i = 0; i < num_pts1; i++) for(int j = 0; j < num_pts2; j++) @@ -96,10 +94,10 @@ double lpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair(PD1[i].second,PD1[i].first)); PD pd2 = PD2; numpts = PD2.size(); for(int i = 0; i < numpts; i++) pd2.push_back(std::pair(PD2[i].second,PD2[i].first)); - return lpwg(pd1, pd2, 2*sqrt(sigma), &pss_weight) / (2*8*NUMPI*sigma); + return lpwgk(pd1, pd2, 2*sqrt(sigma), &pss_weight) / (2*8*3.14159265359*sigma); } /** \brief Computes the Gaussian Persistence Weighted Gaussian Kernel between two persistence diagrams. @@ -112,10 +110,10 @@ double pss(PD PD1, PD PD2, double sigma){ * @param[in] weight weight function for the points in the diagrams. * */ -double gpwg(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ - double k1 = lpwg(PD1,PD1,sigma,weight); - double k2 = lpwg(PD2,PD2,sigma,weight); - double k3 = lpwg(PD1,PD2,sigma,weight); +double gpwgk(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ + double k1 = lpwgk(PD1,PD1,sigma,weight); + double k2 = lpwgk(PD2,PD2,sigma,weight); + double k3 = lpwgk(PD1,PD2,sigma,weight); return exp( - (k1+k2-2*k3) / (2*pow(tau,2)) ); } @@ -129,9 +127,9 @@ double gpwg(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pair * */ double dpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ - double k1 = lpwg(PD1,PD1,sigma,weight); - double k2 = lpwg(PD2,PD2,sigma,weight); - double k3 = lpwg(PD1,PD2,sigma,weight); + double k1 = lpwgk(PD1,PD1,sigma,weight); + double k2 = lpwgk(PD2,PD2,sigma,weight); + double k3 = lpwgk(PD1,PD2,sigma,weight); return std::sqrt(k1+k2-2*k3); } @@ -160,24 +158,24 @@ double compute_angle(const PD & PersDiag, const int & i, const int & j){ // 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 res = 0; - assert((alpha >= 0 && alpha <= NUMPI) || (alpha >= -NUMPI && alpha <= 0)); - if (alpha >= 0 && alpha <= NUMPI){ + assert((alpha >= 0 && alpha <= 3.14159265359) || (alpha >= -3.14159265359 && alpha <= 0)); + if (alpha >= 0 && alpha <= 3.14159265359){ if (cos(alpha) >= 0){ - if(NUMPI/2 <= beta){res = 2-sin(alpha)-sin(beta);} + if(3.14159265359/2 <= beta){res = 2-sin(alpha)-sin(beta);} else{res = sin(beta)-sin(alpha);} } else{ - if(1.5*NUMPI <= beta){res = 2+sin(alpha)+sin(beta);} + if(1.5*3.14159265359 <= beta){res = 2+sin(alpha)+sin(beta);} else{res = sin(alpha)-sin(beta);} } } - if (alpha >= -NUMPI && alpha <= 0){ + if (alpha >= -3.14159265359 && alpha <= 0){ if (cos(alpha) <= 0){ - if(-NUMPI/2 <= beta){res = 2+sin(alpha)+sin(beta);} + if(-3.14159265359/2 <= beta){res = 2+sin(alpha)+sin(beta);} else{res = sin(alpha)-sin(beta);} } else{ - if(NUMPI/2 <= beta){res = 2-sin(alpha)-sin(beta);} + if(3.14159265359/2 <= beta){res = 2-sin(alpha)-sin(beta);} else{res = sin(beta)-sin(alpha);} } } @@ -202,9 +200,9 @@ double compute_sw(const std::vector > > & V1, int N = V1.size(); double sw = 0; for (int i = 0; i < N; i++){ std::vector > U,V; U = V1[i]; V = V2[i]; - double theta1, theta2; theta1 = -NUMPI/2; + double theta1, theta2; theta1 = -3.14159265359/2; unsigned int ku, kv; ku = 0; kv = 0; theta2 = std::min(U[ku].second,V[kv].second); - while(theta1 != NUMPI/2){ + while(theta1 != 3.14159265359/2){ if(PD1[U[ku].first].first != PD2[V[kv].first].first || PD1[U[ku].first].second != PD2[V[kv].first].second) if(theta1 != theta2) sw += compute_int(theta1, theta2, U[ku].first, V[kv].first, PD1, PD2); @@ -214,7 +212,7 @@ double compute_sw(const std::vector > > & V1, theta2 = std::min(U[ku].second, V[kv].second); } } - return sw/NUMPI; + return sw/3.14159265359; } /** \brief Computes the Sliced Wasserstein distance between two persistence diagrams. @@ -292,8 +290,8 @@ double compute_sw(const std::vector > > & V1, } for (int i = 0; i < N; i++){ - anglePerm1[order1[i]].push_back(std::pair(i,NUMPI/2)); - anglePerm2[order2[i]].push_back(std::pair(i,NUMPI/2)); + anglePerm1[order1[i]].push_back(std::pair(i,3.14159265359/2)); + anglePerm2[order2[i]].push_back(std::pair(i,3.14159265359/2)); } // Compute the SW distance with the list of inversions. @@ -301,6 +299,17 @@ double compute_sw(const std::vector > > & V1, } + /** \brief Computes the Sliced Wasserstein Kernel between two persistence diagrams. + * \ingroup kernel + * + * @param[in] PD1 first persistence diagram. + * @param[in] PD2 second persistence diagram. + * @param[in] sigma bandwidth parameter. + * + */ + double swk(PD PD1, PD PD2, double sigma){ + return exp( - sw(PD1,PD2) / (2*pow(sigma, 2)) ); + } @@ -370,7 +379,7 @@ std::vector > random_Fourier(double sigma, int M = 1000 * @param[in] M number of Fourier features. * */ -double approx_lpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}, int M = 1000){ +double approx_lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}, int M = 1000){ std::vector > Z = random_Fourier(sigma, M); std::vector > B1 = Fourier_feat(PD1,Z,weight); std::vector > B2 = Fourier_feat(PD2,Z,weight); @@ -386,10 +395,10 @@ double approx_lpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair(PD1[i].second,PD1[i].first)); PD pd2 = PD2; numpts = PD2.size(); for(int i = 0; i < numpts; i++) pd2.push_back(std::pair(PD2[i].second,PD2[i].first)); - return approx_lpwg(pd1, pd2, 2*sqrt(sigma), &pss_weight, M) / (2*8*NUMPI*sigma); + return approx_lpwgk(pd1, pd2, 2*sqrt(sigma), &pss_weight, M) / (2*8*3.14159265359*sigma); } @@ -404,7 +413,7 @@ double approx_pss(PD PD1, PD PD2, double sigma, int M = 1000){ * @param[in] M number of Fourier features. * */ -double approx_gpwg(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}, int M = 1000){ +double approx_gpwgk(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}, int M = 1000){ std::vector > Z = random_Fourier(sigma, M); std::vector > B1 = Fourier_feat(PD1,Z,weight); std::vector > B2 = Fourier_feat(PD2,Z,weight); @@ -422,7 +431,7 @@ double approx_gpwg(PD PD1, PD PD2, double sigma, double tau, double (*weight)(st */ double approx_sw(PD PD1, PD PD2, int N = 100){ - double step = NUMPI/N; double sw = 0; + double step = 3.14159265359/N; double sw = 0; // Add projections onto diagonal. int n1, n2; n1 = PD1.size(); n2 = PD2.size(); @@ -437,14 +446,27 @@ double approx_sw(PD PD1, PD PD2, int N = 100){ for (int i = 0; i < N; i++){ std::vector > L1, L2; for (int j = 0; j < n; j++){ - L1.push_back( std::pair(j, PD1[j].first*cos(-NUMPI/2+i*step) + PD1[j].second*sin(-NUMPI/2+i*step)) ); - L2.push_back( std::pair(j, PD2[j].first*cos(-NUMPI/2+i*step) + PD2[j].second*sin(-NUMPI/2+i*step)) ); + L1.push_back( std::pair(j, PD1[j].first*cos(-3.14159265359/2+i*step) + PD1[j].second*sin(-3.14159265359/2+i*step)) ); + L2.push_back( std::pair(j, PD2[j].first*cos(-3.14159265359/2+i*step) + PD2[j].second*sin(-3.14159265359/2+i*step)) ); } std::sort(L1.begin(),L1.end(), myComp); std::sort(L2.begin(),L2.end(), myComp); double f = 0; for (int j = 0; j < n; j++) f += std::abs(L1[j].second - L2[j].second); sw += f*step; } - return sw/NUMPI; + return sw/3.14159265359; +} + +/** \brief Computes an approximation of the Sliced Wasserstein Kernel between two persistence diagrams. + * \ingroup kernel + * + * @param[in] PD1 first persistence diagram. + * @param[in] PD2 second persistence diagram. + * @param[in] sigma bandwidth parameter. + * @param[in] N number of points sampled on the circle. + * + */ +double approx_swk(PD PD1, PD PD2, double sigma, int N = 100){ + return exp( - approx_sw(PD1,PD2,N) / (2*pow(sigma,2))); } diff --git a/src/Kernels/test/CMakeLists.txt b/src/Kernels/test/CMakeLists.txt new file mode 100644 index 00000000..9dbb9ed4 --- /dev/null +++ b/src/Kernels/test/CMakeLists.txt @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 2.6) +project(kernel_tests) + +include(GUDHI_test_coverage) + +add_executable ( kernel_test_unit test_kernel.cpp ) +target_link_libraries(kernel_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(kernel_test_unit ${TBB_LIBRARIES}) +endif() + +file(COPY data DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +gudhi_add_coverage_test(kernel_test_unit) diff --git a/src/Kernels/test/test_kernel.cpp b/src/Kernels/test/test_kernel.cpp new file mode 100644 index 00000000..db05fd28 --- /dev/null +++ b/src/Kernels/test/test_kernel.cpp @@ -0,0 +1,56 @@ +/* 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) 2017 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 . + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "kernel" + +#include +#include // float comparison +#include +#include +#include +#include // std::max +#include +#include +#include + +BOOST_AUTO_TEST_CASE(check_PSS) { + std::vector< std::pair > v1, v2; + v1.emplace_back(std::pair(0,1)); + v2.emplace_back(std::pair(0,2)); + BOOST_CHECK(std::abs(Gudhi::kernel::pssk(v1,v2,1) - Gudhi::kernel::approx_pssk(v1,v2,1)) <= 1e-1); +} + +BOOST_AUTO_TEST_CASE(check_PWG) { + std::vector< std::pair > v1, v2; + v1.emplace_back(std::pair(0,1)); + v2.emplace_back(std::pair(0,2)); + BOOST_CHECK(std::abs(Gudhi::kernel::lpwgk(v1,v2,1) - Gudhi::kernel::approx_lpwgk(v1,v2,1)) <= 1e-1); + BOOST_CHECK(std::abs(Gudhi::kernel::gpwgk(v1,v2,1,1) - Gudhi::kernel::approx_gpwgk(v1,v2,1,1)) <= 1e-1); +} + +BOOST_AUTO_TEST_CASE(check_SW) { + std::vector< std::pair > v1, v2; + v2.emplace_back(std::pair(0,2)); + BOOST_CHECK(std::abs(Gudhi::kernel::sw(v1,v2) - Gudhi::kernel::approx_sw(v1,v2)) <= 1e-3); + BOOST_CHECK(std::abs(Gudhi::kernel::sw(v1,v2) - 2*std::sqrt(2)/3.1415) <= 1e-3); +} -- cgit v1.2.3