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

Former-commit-id: 6fd757a9dfd0abe7ddfbaeca1af667a4c93af34b

3 files changed, 600 insertions, 0 deletions

diff --git a/src/Kernels/include/gudhi/PSS.h b/src/Kernels/include/gudhi/PSS.h
new file mode 100644
index 00000000..70743c47
--- /dev/null
+++ b/src/Kernels/include/gudhi/PSS.h
@@ -0,0 +1,108 @@
+/*    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 (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 PSS_H_
+#define PSS_H_
+
+#define NUMPI 3.14159265359
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <cstdlib>
+#include <cstdio>
+#include <iomanip>
+
+#include <string>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+
+#include <set>
+#include <map>
+#include <vector>
+#include <algorithm>
+#include <limits>
+#include <assert.h>
+
+#include <cmath>
+#include <math.h>
+
+#include <memory>
+#include <stdexcept>
+#include <omp.h>
+
+#include <gmp.h>
+#include <gmpxx.h>
+#include <random>
+#include <chrono>
+#include <ctime>
+
+#include "../../figtree-0.9.3/include/figtree.h"
+#include "../../figtree-0.9.3/external/ann_1.1.1/include/ANN/ANN.h"
+
+using PD = std::vector<std::pair<double,double> >;
+
+namespace Gudhi {
+namespace persistence_scale_space {
+
+double compute_exact_pss(PD PD1, PD PD2, double sigma = 1){
+  double k = 0;
+  for(int i = 0; i < PD1.size(); i++){
+    for(int j = 0; j < PD2.size(); j++){
+      k += exp(  -(    pow(PD1[i].first - PD2[j].first,  2)  +  pow(PD1[i].second - PD2[j].second, 2)     )/(8*sigma))   -\
+           exp(  -(    pow(PD1[i].first - PD2[j].second, 2)  +  pow(PD1[i].second - PD2[j].first,  2)     )/(8*sigma));
+    }
+  }
+  return k/(8*NUMPI*sigma);
+}
+
+double compute_approximate_pss(PD PD1, PD PD2, double sigma = 1, double error = 1e-2){
+
+  double k = 0;
+
+  int d = 2; int N = PD1.size(); int M = PD2.size(); double h = std::sqrt(8*sigma);
+  double* x = new double[2*N]; double* y = new double[2*M]; double* q = new double[N];
+  for(int i = 0; i < N; i++){
+    q[i] = 1.0/(8*NUMPI*sigma);
+    x[2*i] = PD1[i].first; x[2*i+1] = PD1[i].second;
+  }
+  for(int i = 0; i < M; i++){  y[2*i] = PD2[i].first; y[2*i+1] = PD2[i].second;  }
+  double* g_auto = new double[M];
+  memset(g_auto, 0, sizeof(double)*M);
+
+  figtree(d, N, M, 1, x, h, q, y, error, g_auto);
+  for(int i = 0; i < M; i++)  k += g_auto[i];
+
+  for(int i = 0; i < M; i++){  y[2*i] = PD2[i].second; y[2*i+1] = PD2[i].first;  }
+
+  figtree(d, N, M, 1, x, h, q, y, error, g_auto);
+  for(int i = 0; i < M; i++)  k -= g_auto[i];
+
+  delete[] x; delete[] y; delete[] q; delete[] g_auto;
+  return k;
+}
+
+} // namespace persistence_scale_space
+
+} // namespace Gudhi
+
+#endif // PSS_H_
diff --git a/src/Kernels/include/gudhi/PWG.h b/src/Kernels/include/gudhi/PWG.h
new file mode 100644
index 00000000..bc491ae7
--- /dev/null
+++ b/src/Kernels/include/gudhi/PWG.h
@@ -0,0 +1,204 @@
+/*    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 (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 PWG_H_
+#define PWG_H_
+
+#define NUMPI 3.14159265359
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <cstdlib>
+#include <cstdio>
+#include <iomanip>
+
+#include <string>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+
+#include <set>
+#include <map>
+#include <vector>
+#include <algorithm>
+#include <limits>
+#include <assert.h>
+
+#include <cmath>
+#include <math.h>
+
+#include <memory>
+#include <stdexcept>
+#include <omp.h>
+
+#include <gmp.h>
+#include <gmpxx.h>
+#include <random>
+#include <chrono>
+#include <ctime>
+
+using PD = std::vector<std::pair<double,double> >;
+
+namespace Gudhi {
+namespace persistence_weighted_gaussian {
+
+double compute_exact_linear_pwg(PD PD1, PD PD2, double sigma, double C, int p){
+
+  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++){
+      k += atan(C*pow(PD1[i].second-PD1[i].first,p))*atan(C*pow(PD2[j].second-PD2[j].first,p))*\
+            exp(  -(    pow(PD1[i].first-PD2[j].first,2) + pow(PD1[i].second-PD2[j].second,2)     )/(2*pow(sigma,2))  );
+    }
+  }
+
+  return k;
+
+}
+
+double compute_exact_gaussian_pwg(PD PD1, PD PD2, double sigma, double C, int p, double tau){
+
+  int num_pts1 = PD1.size();
+  int num_pts2 = PD2.size();
+
+  double k1 = 0;
+  for(int i = 0; i < num_pts1; i++){
+    for(int j = 0; j < num_pts1; j++){
+      k1 += atan(C*pow(PD1[i].second-PD1[i].first,p))*atan(C*pow(PD1[j].second-PD1[j].first,p))*\
+            exp(  -(    pow(PD1[i].first-PD1[j].first,2) + pow(PD1[i].second-PD1[j].second,2)     )/(2*pow(sigma,2))  );
+    }
+  }
+
+  double k2 = 0;
+  for(int i = 0; i < num_pts2; i++){
+    for(int j = 0; j < num_pts2; j++){
+      k2 += atan(C*pow(PD2[i].second-PD2[i].first,p))*atan(C*pow(PD2[j].second-PD2[j].first,p))*\
+            exp(  -(    pow(PD2[i].first-PD2[j].first,2) + pow(PD2[i].second-PD2[j].second,2)     )/(2*pow(sigma,2))  );
+    }
+  }
+
+  double k3 = compute_exact_linear_pwg(PD1,PD2,sigma,C,p);
+  return exp( - (k1+k2-2*k3) / (2*pow(tau,2))  );
+
+}
+
+double compute_exact_gaussian_RKHSdist(PD PD1, PD PD2, double sigma, double C, int p){
+
+  int num_pts1 = PD1.size();
+  int num_pts2 = PD2.size();
+
+  double k1 = 0;
+  for(int i = 0; i < num_pts1; i++){
+    for(int j = 0; j < num_pts1; j++){
+      k1 += atan(C*pow(PD1[i].second-PD1[i].first,p))*atan(C*pow(PD1[j].second-PD1[j].first,p))*\
+            exp(  -(    pow(PD1[i].first-PD1[j].first,2) + pow(PD1[i].second-PD1[j].second,2)     )/(2*pow(sigma,2))  );
+    }
+  }
+
+  double k2 = 0;
+  for(int i = 0; i < num_pts2; i++){
+    for(int j = 0; j < num_pts2; j++){
+      k2 += atan(C*pow(PD2[i].second-PD2[i].first,p))*atan(C*pow(PD2[j].second-PD2[j].first,p))*\
+            exp(  -(    pow(PD2[i].first-PD2[j].first,2) + pow(PD2[i].second-PD2[j].second,2)     )/(2*pow(sigma,2))  );
+    }
+  }
+
+  double k3 = compute_exact_linear_pwg(PD1,PD2,sigma,C,p);
+  return std::sqrt(k1+k2-2*k3);
+
+}
+
+double compute_approximate_linear_pwg_from_Fourier_features(const std::vector<std::pair<double,double> >& B1, \
+                                                            const std::vector<std::pair<double,double> >& B2){
+  double d = 0; int M  = B1.size();
+  for(int i = 0; i < M; i++) d += B1[i].first*B2[i].first + B1[i].second*B2[i].second;
+  return (1.0/M)*d;
+}
+
+double compute_approximate_gaussian_pwg_from_Fourier_features(const std::vector<std::pair<double,double> >& B1, \
+                                                              const std::vector<std::pair<double,double> >& B2, double tau){
+  int M  = B1.size();
+  double d3 = compute_approximate_linear_pwg_from_Fourier_features(B1, B2);
+  double d1 = 0; double d2 = 0;
+  for(int i = 0; i < M; i++){d1 += pow(B1[i].first,2) + pow(B1[i].second,2); d2 += pow(B2[i].first,2) + pow(B2[i].second,2);}
+  return exp( -((1.0/M)*(d1+d2)-2*d3) / (2*pow(tau,2))  );
+}
+
+double compute_approximate_gaussian_RKHSdist_from_Fourier_features(const std::vector<std::pair<double,double> >& B1, \
+                                                              const std::vector<std::pair<double,double> >& B2){
+  int M  = B1.size();
+  double d3 = compute_approximate_linear_pwg_from_Fourier_features(B1, B2);
+  double d1 = 0; double d2 = 0;
+  for(int i = 0; i < M; i++){d1 += pow(B1[i].first,2) + pow(B1[i].second,2); d2 += pow(B2[i].first,2) + pow(B2[i].second,2);}
+  return std::sqrt((1.0/M)*(d1+d2)-2*d3);
+}
+
+std::vector<std::pair<double,double> > compute_Fourier_features(double C, int p, PD D, std::vector<std::pair<double,double> > Z){
+  int m = D.size(); std::vector<std::pair<double,double> > B; int M = Z.size();
+  for(int i = 0; i < M; i++){
+    double d1 = 0; double d2 = 0; double zx = Z[i].first; double zy = Z[i].second;
+    for(int j = 0; j < m; j++){
+      double x = D[j].first; double y = D[j].second;
+      d1 += atan(C*pow(y-x,p))*cos(x*zx + y*zy);
+      d2 += atan(C*pow(y-x,p))*sin(x*zx + y*zy);
+    }
+    B.push_back(std::pair<double,double>(d1,d2));
+  }
+  return B;
+}
+
+std::vector<std::pair<double,double> > random_Fourier(double sigma, int M = 1000){
+  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++){
+    //unsigned seedx = 2*i; unsigned seedy = 2*i+1;
+    //std::default_random_engine generatorx(seedx); std::default_random_engine generatory(seedy);
+    std::mt19937 e1(rd()); std::mt19937 e2(rd());
+    double zx = distrib(e1/*generatorx*/); double zy = distrib(e2/*generatory*/);
+    Z.push_back(std::pair<double,double>((1/sigma)*zx,(1/sigma)*zy));
+  }
+  return Z;
+}
+
+double compute_approximate_linear_pwg(PD PD1, PD PD2, double sigma, double C, int p, int M = 1000){
+  std::vector<std::pair<double,double> > Z = random_Fourier(sigma, M);
+  std::vector<std::pair<double,double> > B1 = compute_Fourier_features(C,p,PD1,Z);
+  std::vector<std::pair<double,double> > B2 = compute_Fourier_features(C,p,PD2,Z);
+  return compute_approximate_linear_pwg_from_Fourier_features(B1,B2);
+}
+
+double compute_approximate_gaussian_pwg(PD PD1, PD PD2, double sigma, double C, int p, double tau, int M = 1000){
+  std::vector<std::pair<double,double> > Z = random_Fourier(sigma, M);
+  std::vector<std::pair<double,double> > B1 = compute_Fourier_features(C,p,PD1,Z);
+  std::vector<std::pair<double,double> > B2 = compute_Fourier_features(C,p,PD2,Z);
+  return compute_approximate_gaussian_pwg_from_Fourier_features(B1,B2,tau);
+}
+
+
+} // namespace persistence_weighted_gaussian
+
+} // namespace Gudhi
+
+#endif //PWG_H_
diff --git a/src/Kernels/include/gudhi/SW.h b/src/Kernels/include/gudhi/SW.h
new file mode 100644
index 00000000..6871d990
--- /dev/null
+++ b/src/Kernels/include/gudhi/SW.h
@@ -0,0 +1,288 @@
+/*    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 (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 SW_H_
+#define SW_H_
+
+#define NUMPI 3.14159265359
+
+#include <stdlib.h>
+#include <cstdlib>
+#include <string>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include <set>
+#include <map>
+#include <limits>
+#include <cmath>
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+#include <cstdio>
+#include <memory>
+#include <stdexcept>
+#include <omp.h>
+#include <assert.h>
+#include <iomanip>
+#include <gmp.h>
+#include <gmpxx.h>
+#include <random>
+#include <chrono>
+#include <ctime>
+#include <math.h>
+
+using PD = std::vector<std::pair<double,double> >;
+
+std::vector<std::pair<double,double> > PDi, PDj;
+
+bool compOri(const int& p, const int& q){
+  if(PDi[p].second != PDi[q].second)
+    return (PDi[p].second < PDi[q].second);
+  else
+    return (PDi[p].first > PDi[q].first);
+}
+
+bool compOrj(const int& p, const int& q){
+  if(PDj[p].second != PDj[q].second)
+    return (PDj[p].second < PDj[q].second);
+  else
+    return (PDj[p].first > PDj[q].first);
+}
+
+bool sortAngle(const std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){
+  return p1.first < p2.first;
+}
+
+bool myComp(const std::pair<int,double> & P1, const std::pair<int,double> & P2){return P1.second < P2.second;}
+
+namespace Gudhi {
+namespace sliced_wasserstein {
+
+
+double compute_approximate_SW(PD PD1, PD PD2, int N = 100){
+
+  double step = NUMPI/N; double sw = 0;
+
+  // Add projections onto diagonal.
+  // ******************************
+  int n1, n2; n1 = PD1.size(); n2 = PD2.size();
+  for (int i = 0; i < n2; i++)
+    PD1.push_back(std::pair<double,double>(   (PD2[i].first+PD2[i].second)/2,  (PD2[i].first+PD2[i].second)/2)   );
+  for (int i = 0; i < n1; i++)
+    PD2.push_back(std::pair<double,double>(   (PD1[i].first+PD1[i].second)/2,  (PD1[i].first+PD1[i].second)/2)   );
+  int n = PD1.size();
+
+  // Sort and compare all projections.
+  // *********************************
+  //#pragma omp parallel for
+  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))   );
+    }
+    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;
+}
+
+double compute_int_cos(const double& alpha, const double& beta){ // Valid only if alpha is in [-pi,pi] and beta-alpha is in [0,pi]
+  double res;
+  assert((alpha >= 0 && alpha <= NUMPI) || (alpha >= -NUMPI && alpha <= 0));
+  if (alpha >= 0 && alpha <= NUMPI){
+    if (cos(alpha) >= 0){
+      if(NUMPI/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);}
+      else{res = sin(alpha)-sin(beta);}
+    }
+  }
+  if (alpha >= -NUMPI && alpha <= 0){
+    if (cos(alpha) <= 0){
+      if(-NUMPI/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);}
+      else{res = sin(beta)-sin(alpha);}
+    }
+  }
+  return res;
+}
+
+double compute_int(const double& theta1, const double& theta2, const int& p, const int& q){
+  double norm = std::sqrt(pow(PDi[p].first-PDj[q].first,2) + pow(PDi[p].second-PDj[q].second,2));
+  double angle1;
+  if (PDi[p].first > PDj[q].first)
+    angle1 = theta1 - asin( (PDi[p].second-PDj[q].second)/norm  );
+  else
+    angle1 = theta1 - asin( (PDj[q].second-PDi[p].second)/norm  );
+  double angle2 = angle1+theta2-theta1;
+  double integral = compute_int_cos(angle1,angle2);
+  return norm*integral;
+}
+
+double compute_sw(const std::vector<std::vector<std::pair<int,double> > >& V1, \
+                  const std::vector<std::vector<std::pair<int,double> > >& V2){
+  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;
+    int ku, kv; ku = 0; kv = 0; theta2 = std::min(U[ku].second,V[kv].second);
+    while(theta1 != NUMPI/2){
+      if(PDi[U[ku].first].first != PDj[V[kv].first].first || PDi[U[ku].first].second != PDj[V[kv].first].second)
+        if(theta1 != theta2)
+          sw += compute_int(theta1,theta2,U[ku].first,V[kv].first);
+      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);
+    }
+  }
+  return sw/NUMPI;
+}
+
+double compute_angle(const PD& PersDiag, const int& i, const int& j){
+  std::pair<double,double> vect; double x1,y1, x2,y2;
+  x1 = PersDiag[i].first; y1 = PersDiag[i].second;
+  x2 = PersDiag[j].first; y2 = PersDiag[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(pow(vect.first,2) + pow(vect.second,2));
+  return asin(vect.second/norm);
+}
+
+double compute_exact_SW(PD PD1, PD PD2){
+
+  // Add projections onto diagonal.
+  // ******************************
+  int n1, n2; n1 = PD1.size(); n2 = PD2.size(); double max_ordinate = std::numeric_limits<double>::min();
+  for (int i = 0; i < n2; i++){
+    max_ordinate = std::max(max_ordinate, PD2[i].second);
+    PD1.push_back(std::pair<double,double>(  ((PD2[i].first+PD2[i].second)/2), ((PD2[i].first+PD2[i].second)/2))   );
+  }
+  for (int i = 0; i < n1; i++){
+    max_ordinate = std::max(max_ordinate, PD1[i].second);
+    PD2.push_back(std::pair<double,double>(  ((PD1[i].first+PD1[i].second)/2), ((PD1[i].first+PD1[i].second)/2))   );
+  }
+  int N = PD1.size(); assert(N==PD2.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 < N; i++){
+    PD1[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); PD1[i].second += thresh*(1.0-2.0*rand()/RAND_MAX);
+    PD2[i].first += thresh*(1.0-2.0*rand()/RAND_MAX); PD2[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 < N; i++){
+    for (int j = i+1; j < N; j++){
+      double theta1 = compute_angle(PD1,i,j); double theta2 = compute_angle(PD2,i,j);
+      angles1.push_back(std::pair<double, std::pair<int,int> >(theta1, std::pair<int,int>(i,j)));
+      angles2.push_back(std::pair<double, std::pair<int,int> >(theta2, std::pair<int,int>(i,j)));
+    }
+  }
+
+  // Sort angles.
+  // ************
+  std::sort(angles1.begin(), angles1.end(), sortAngle); std::sort(angles2.begin(), angles2.end(), sortAngle);
+
+  // Initialize orders of the points of both PD (given by ordinates when theta = -pi/2).
+  // ***********************************************************************************
+  PDi = PD1; PDj = PD2;
+  std::vector<int> orderp1, orderp2;
+  for (int i = 0; i < N; i++){orderp1.push_back(i); orderp2.push_back(i);}
+  std::sort(orderp1.begin(),orderp1.end(),compOri); std::sort(orderp2.begin(),orderp2.end(),compOrj);
+
+  // Find the inverses of the orders.
+  // ********************************
+  std::vector<int> order1(N); std::vector<int> order2(N);
+  for(int i = 0; i < N; i++){
+    for (int j = 0; j < N; j++)
+      if(orderp1[j] == i)
+        order1[i] = j;
+  }
+  for(int i = 0; i < N; i++){
+    for (int j = 0; j < N; j++)
+      if(orderp2[j] == i)
+        order2[i] = j;
+  }
+
+  // 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(N);
+  std::vector<std::vector<std::pair<int,double> > > anglePerm2(N);
+
+  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]].push_back(std::pair<int, double>(p,theta));
+    anglePerm1[order1[q]].push_back(std::pair<int, double>(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]].push_back(std::pair<int, double>(p,theta));
+    anglePerm2[order2[q]].push_back(std::pair<int, double>(q,theta));
+    int a = order2[p]; int b = order2[q]; order2[p] = b; order2[q] = a;
+  }
+
+  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));
+  }
+
+  // Compute the SW distance with the list of inversions.
+  // ****************************************************
+  return compute_sw(anglePerm1,anglePerm2);
+
+}
+
+}}
+
+#endif
+
+