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

Former-commit-id: a6a297a7f14703e55954706328072acdd447484c

2 files changed, 11 insertions, 33 deletions

diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
index 12188526..43f10b8c 100644
--- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
+++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h
@@ -1002,8 +1002,6 @@ double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product(const
   }
 
   else{
-    GUDHI_CHECK(this->approx != second.approx || this->f != second.f, std::invalid_argument("Error: different values for representations"));
-
     int num_pts1 = this->d.size(); int num_pts2 = second.d.size(); double kernel_val = 0;
     for(int i = 0; i < num_pts1; i++)
       for(int j = 0; j < num_pts2; j++)
diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
index a3c0dc2f..6f67f7bc 100644
--- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
+++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
@@ -97,23 +97,7 @@ class Sliced_Wasserstein {
 
   // Compute the angle formed by two points of a PD
   double compute_angle(const Persistence_diagram & diag, int i, int j) const {
-    std::pair<double,double> vect; double x1,y1, x2,y2;
-    x1 = diag[i].first; y1 = diag[i].second;
-    x2 = diag[j].first; y2 = diag[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(vect.first*vect.first + vect.second*vect.second);
-    return asin(vect.second/norm);
+    if(diag[i].second == diag[j].second) return pi/2; else return atan((diag[j].first-diag[i].first)/(diag[i].second-diag[j].second));
   }
 
   // Compute the integral of |cos()| between alpha and beta, valid only if alpha is in [-pi,pi] and beta-alpha is in [0,pi]
@@ -145,10 +129,7 @@ class Sliced_Wasserstein {
   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 (diag1[p].first > diag2[q].first)
-      angle1 = theta1 - asin( (diag1[p].second-diag2[q].second)/norm  );
-    else
-      angle1 = theta1 - asin( (diag2[q].second-diag1[p].second)/norm  );
+    if (diag1[p].first == diag2[q].first) angle1 = theta1 - pi/2; else angle1 = theta1 - atan((diag1[p].second-diag2[q].second)/(diag1[p].first-diag2[q].first));
     double angle2 = angle1 + theta2 - theta1;
     double integral = compute_int_cos(angle1,angle2);
     return norm*integral;
@@ -164,24 +145,23 @@ class Sliced_Wasserstein {
     if(this->approx == -1){
 
       // Add projections onto diagonal.
-      int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double max_ordinate = std::numeric_limits<double>::lowest();
+      int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double min_ordinate = std::numeric_limits<double>::max(); double min_abscissa = std::numeric_limits<double>::max();
       for (int i = 0; i < n2; i++){
-        max_ordinate = std::max(max_ordinate, diagram2[i].second);
+        min_ordinate = std::min(min_ordinate, diagram2[i].second); min_abscissa = std::min(min_abscissa, diagram2[i].first);
         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);
+        min_ordinate = std::min(min_ordinate, diagram1[i].second); min_abscissa = std::min(min_abscissa, diagram1[i].first);
         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);
+      double thresh_y = pow(10,log10(min_ordinate)-5); double thresh_x = pow(10,log10(min_abscissa)-5);
       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);
+        diagram1[i].first += thresh_x*(1.0-2.0*rand()/RAND_MAX); diagram1[i].second += thresh_y*(1.0-2.0*rand()/RAND_MAX);
+        diagram2[i].first += thresh_x*(1.0-2.0*rand()/RAND_MAX); diagram2[i].second += thresh_y*(1.0-2.0*rand()/RAND_MAX);
       }
 
       // Compute all angles in both PDs.
@@ -201,8 +181,8 @@ class Sliced_Wasserstein {
       // 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 < 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); } );
+      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<int> order1(num_pts_dgm); std::vector<int> order2(num_pts_dgm);
@@ -274,6 +254,7 @@ class Sliced_Wasserstein {
  public:
 
   /** \brief Sliced Wasserstein kernel constructor.
+   * \implements Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_scalar_product
    * \ingroup Sliced_Wasserstein
    *
    * @param[in] _diagram  persistence diagram.
@@ -282,7 +263,6 @@ class Sliced_Wasserstein {
    *                      points on all directions are stored in memory to reduce computation time.
    *
    */
-  // This class implements the following concepts: Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_scalar_product
   Sliced_Wasserstein(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 10):diagram(_diagram), approx(_approx), sigma(_sigma) {build_rep();}
 
   /** \brief Evaluation of the kernel on a pair of diagrams.