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

Former-commit-id: f7ac7ffce4fbfb04b099797a02999d4a93e3f22b

3 files changed, 115 insertions, 87 deletions

diff --git a/src/Kernels/include/gudhi/PSS.h b/src/Kernels/include/gudhi/PSS.h
index 70743c47..5111a09f 100644
--- a/src/Kernels/include/gudhi/PSS.h
+++ b/src/Kernels/include/gudhi/PSS.h
@@ -56,8 +56,8 @@
 #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"
+#include "figtree.h"
+#include "ANN.h"
 
 using PD = std::vector<std::pair<double,double> >;
 
diff --git a/src/Kernels/include/gudhi/SW.h b/src/Kernels/include/gudhi/SW.h
index 6871d990..0b041252 100644
--- a/src/Kernels/include/gudhi/SW.h
+++ b/src/Kernels/include/gudhi/SW.h
@@ -55,47 +55,36 @@
 
 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 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 {
 
 
+// ********************************************************************
+// Approximate computation.
+// ********************************************************************
+
+/** \brief Computes an approximation of the Sliced Wasserstein distance between two persistence diagrams
+ *
+ * @param[in] N number of points sampled on the circle.
+ *
+ */
+
 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)   );
+    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)   );
+    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;
@@ -110,8 +99,55 @@ double compute_approximate_SW(PD PD1, PD PD2, int N = 100){
   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;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+// ********************************************************************
+// Exact computation.
+// ********************************************************************
+
+
+
+// Compute the angle formed by two points of a PD
+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);
+}
+
+// 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 res = 0;
   assert((alpha >= 0 && alpha <= NUMPI) || (alpha >= -NUMPI && alpha <= 0));
   if (alpha >= 0 && alpha <= NUMPI){
     if (cos(alpha) >= 0){
@@ -136,75 +172,76 @@ double compute_int_cos(const double& alpha, const double& beta){ // Valid only i
   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 compute_int(const double & theta1, const double & theta2, const int & p, const int & q, const PD & PD1, const PD & PD2){
+  double norm = std::sqrt(pow(PD1[p].first-PD2[q].first,2) + pow(PD1[p].second-PD2[q].second,2));
   double angle1;
-  if (PDi[p].first > PDj[q].first)
-    angle1 = theta1 - asin( (PDi[p].second-PDj[q].second)/norm  );
+  if (PD1[p].first > PD2[q].first)
+    angle1 = theta1 - asin( (PD1[p].second-PD2[q].second)/norm  );
   else
-    angle1 = theta1 - asin( (PDj[q].second-PDi[p].second)/norm  );
-  double angle2 = angle1+theta2-theta1;
+    angle1 = theta1 - asin( (PD2[q].second-PD1[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){
+
+
+double compute_sw(const std::vector<std::vector<std::pair<int,double> > > & V1, const std::vector<std::vector<std::pair<int,double> > > & V2, const PD & PD1, const PD & PD2){
   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);
+    unsigned 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(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);
+          sw += Gudhi::sliced_wasserstein::compute_int(theta1, theta2, U[ku].first, V[kv].first, PD1, PD2);
       theta1 = theta2;
-      if (  (theta2 == U[ku].second)  &&  ku < U.size()-1  ){ku++;}
-      if (  (theta2 == V[kv].second)  &&  kv < V.size()-1  ){kv++;}
+      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);
-}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/** \brief Computes the Sliced Wasserstein distance between two persistence diagrams
+ *
+ */
 
 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();
+  int n1, n2; n1 = PD1.size(); n2 = PD2.size(); double max_ordinate = std::numeric_limits<double>::lowest();
   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))   );
+    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))   );
+    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));
@@ -214,43 +251,30 @@ double compute_exact_SW(PD PD1, PD PD2){
   }
 
   // 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);
+      double theta1 = Gudhi::sliced_wasserstein::compute_angle(PD1,i,j); double theta2 = Gudhi::sliced_wasserstein::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;
+  // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2).
   std::vector<int> orderp1, orderp2;
-  for (int i = 0; i < N; i++){orderp1.push_back(i); orderp2.push_back(i);}
-  std::sort(orderp1.begin(),orderp1.end(),compOri); std::sort(orderp2.begin(),orderp2.end(),compOrj);
+  for (int i = 0; i < N; i++){ orderp1.push_back(i); orderp2.push_back(i); }
+  std::sort( orderp1.begin(), orderp1.end(), [=](int i, int j){ if(PD1[i].second != PD1[j].second) return (PD1[i].second < PD1[j].second); else return (PD1[i].first > PD1[j].first); } );
+  std::sort( orderp2.begin(), orderp2.end(), [=](int i, int j){ if(PD2[i].second != PD2[j].second) return (PD2[i].second < PD2[j].second); else return (PD2[i].first > PD2[j].first); } );
 
   // 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;
-  }
+  for(int i = 0; i < N; i++)  for (int j = 0; j < N; j++)  if(orderp1[j] == i){  order1[i] = j; break;  }
+  for(int i = 0; i < N; i++)  for (int j = 0; j < N; j++)  if(orderp2[j] == i){  order2[i] = j; break;  }
 
   // Record all inversions of points in the orders as theta varies along the positive half-disk.
-  // *******************************************************************************************
   std::vector<std::vector<std::pair<int,double> > > anglePerm1(N);
   std::vector<std::vector<std::pair<int,double> > > anglePerm2(N);
 
@@ -276,11 +300,15 @@ double compute_exact_SW(PD PD1, PD PD2){
   }
 
   // Compute the SW distance with the list of inversions.
-  // ****************************************************
-  return compute_sw(anglePerm1,anglePerm2);
+  return Gudhi::sliced_wasserstein::compute_sw(anglePerm1, anglePerm2, PD1, PD2);
 
 }
 
+
+
+
+
+
 }}
 
 #endif
diff --git a/src/Kernels/include/gudhi/figtree-0.9.3.zip b/src/Kernels/include/gudhi/figtree-0.9.3.zip
new file mode 100644
index 00000000..a9468274
--- /dev/null
+++ b/src/Kernels/include/gudhi/figtree-0.9.3.zip