git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3110 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 95fd5b9e4cd4738ecdeaeb27f20d267816b7bead

9 files changed, 239 insertions, 53 deletions

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 <http://www.gnu.org/licenses/>.
  */
 
-#define NUMPI 3.14159265359
 #include <gudhi/kernel.h>
 
-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<double, double> > 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<double, double> > 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<double,double>(b,d));
+  }
+
+  std::ifstream input2(PDname2);
+  while(std::getline(input2,line)){
+    std::stringstream stream(line); stream >> b; stream >> d; v2.push_back(std::pair<double,double>(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 <stdlib.h>
 #include <stdio.h>
 #include <cstdlib>
@@ -80,7 +78,7 @@ double pss_weight(std::pair<double,double> P){
  * @param[in] weight   weight function for the points in the diagrams.
  *
  */
-double lpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair<double,double>) = [](std::pair<double,double> P){return atan(P.second - P.first);}){
+double lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair<double,double>) = [](std::pair<double,double> 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<double,doub
  * @param[in] sigma    bandwidth parameter of the Gaussian Kernel used for the Kernel Mean Embedding of the diagrams.
  *
  */
-double pss(PD PD1, PD PD2, double sigma){
+double pssk(PD PD1, PD PD2, double sigma){
   PD pd1 = PD1; int numpts = PD1.size();    for(int i = 0; i < numpts; i++)  pd1.push_back(std::pair<double,double>(PD1[i].second,PD1[i].first));
   PD pd2 = PD2;     numpts = PD2.size();    for(int i = 0; i < numpts; i++)  pd2.push_back(std::pair<double,double>(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<double,double>) = [](std::pair<double,double> 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<double,double>) = [](std::pair<double,double> 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<double,double>) = [](std::pair<double,double> 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<std::vector<std::pair<int,double> > > & V1,
   int N = V1.size(); double sw = 0;
   for (int i = 0; i < N; i++){
     std::vector<std::pair<int,double> > 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<std::vector<std::pair<int,double> > > & 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<std::vector<std::pair<int,double> > > & V1,
   }
 
   for (int i = 0; i < N; i++){
-    anglePerm1[order1[i]].push_back(std::pair<int, double>(i,NUMPI/2));
-    anglePerm2[order2[i]].push_back(std::pair<int, double>(i,NUMPI/2));
+    anglePerm1[order1[i]].push_back(std::pair<int, double>(i,3.14159265359/2));
+    anglePerm2[order2[i]].push_back(std::pair<int, double>(i,3.14159265359/2));
   }
 
   // Compute the SW distance with the list of inversions.
@@ -301,6 +299,17 @@ double compute_sw(const std::vector<std::vector<std::pair<int,double> > > & 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<std::pair<double,double> > 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<double,double>) = [](std::pair<double,double> P){return atan(P.second - P.first);}, int M = 1000){
+double approx_lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair<double,double>) = [](std::pair<double,double> P){return atan(P.second - P.first);}, int M = 1000){
   std::vector<std::pair<double,double> > Z =  random_Fourier(sigma, M);
   std::vector<std::pair<double,double> > B1 = Fourier_feat(PD1,Z,weight);
   std::vector<std::pair<double,double> > B2 = Fourier_feat(PD2,Z,weight);
@@ -386,10 +395,10 @@ double approx_lpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair<doub
  * @param[in] M         number of Fourier features.
  *
  */
-double approx_pss(PD PD1, PD PD2, double sigma, int M = 1000){
+double approx_pssk(PD PD1, PD PD2, double sigma, int M = 1000){
   PD pd1 = PD1; int numpts = PD1.size();    for(int i = 0; i < numpts; i++)  pd1.push_back(std::pair<double,double>(PD1[i].second,PD1[i].first));
   PD pd2 = PD2;     numpts = PD2.size();    for(int i = 0; i < numpts; i++)  pd2.push_back(std::pair<double,double>(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<double,double>) = [](std::pair<double,double> 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<double,double>) = [](std::pair<double,double> P){return atan(P.second - P.first);}, int M = 1000){
   std::vector<std::pair<double,double> > Z =  random_Fourier(sigma, M);
   std::vector<std::pair<double,double> > B1 = Fourier_feat(PD1,Z,weight);
   std::vector<std::pair<double,double> > 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<std::pair<int,double> > L1, L2;
     for (int j = 0; j < n; j++){
-      L1.push_back(   std::pair<int,double>(j, PD1[j].first*cos(-NUMPI/2+i*step) + PD1[j].second*sin(-NUMPI/2+i*step))   );
-      L2.push_back(   std::pair<int,double>(j, PD2[j].first*cos(-NUMPI/2+i*step) + PD2[j].second*sin(-NUMPI/2+i*step))   );
+      L1.push_back(   std::pair<int,double>(j, PD1[j].first*cos(-3.14159265359/2+i*step) + PD1[j].second*sin(-3.14159265359/2+i*step))   );
+      L2.push_back(   std::pair<int,double>(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 <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/kernel.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/reader_utils.h>
+
+BOOST_AUTO_TEST_CASE(check_PSS) {
+  std::vector< std::pair<double, double> > v1, v2;
+  v1.emplace_back(std::pair<double,double>(0,1));
+  v2.emplace_back(std::pair<double,double>(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<double, double> > v1, v2;
+  v1.emplace_back(std::pair<double,double>(0,1));
+  v2.emplace_back(std::pair<double,double>(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<double, double> > v1, v2;
+  v2.emplace_back(std::pair<double,double>(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);
+}