diff options
author | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2018-04-23 15:22:13 +0000 |
---|---|---|
committer | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2018-04-23 15:22:13 +0000 |
commit | 541284f6f1bf7d4a76daac8a52850c7162a765cb (patch) | |
tree | 2ebecae35daf30dbcfb9c683f9587b93e117f443 /src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h | |
parent | 5e24206f945f66575c7c179d74e9661cf60ca3df (diff) |
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3387 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 3fe2ae4af0c7cadf507fc5148c05dcf664c5e151
Diffstat (limited to 'src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h')
-rw-r--r-- | src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h | 243 |
1 files changed, 116 insertions, 127 deletions
diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h index 6a9a607e..235918fe 100644 --- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -40,19 +40,17 @@ #include <utility> #include <functional> -using PD = std::vector<std::pair<double,double> >; - namespace Gudhi { namespace Persistence_representations { /** * \class Sliced_Wasserstein gudhi/Sliced_Wasserstein.h - * \brief A class implementing the Sliced Wasserstein Kernel. + * \brief A class implementing the Sliced Wasserstein kernel. * * \ingroup Persistence_representations * * \details - * The Sliced Wasserstein Kernel is defined as a Gaussian-like Kernel between persistence diagrams, where the distance used for + * 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: * @@ -65,15 +63,14 @@ namespace Persistence_representations { * * \f$ k(D_1,D_2) = {\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right).\f$ * - * For more details, please consult <i>Sliced Wasserstein Kernel for Persistence Diagrams</i>\cite pmlr-v70-carriere17a . - * It implements the following concepts: Topological_data_with_distances, Topological_data_with_scalar_product. + * For more details, please see \cite pmlr-v70-carriere17a . * **/ class Sliced_Wasserstein { protected: - PD diagram; + Persistence_diagram diagram; int approx; double sigma; std::vector<std::vector<double> > projections, projections_diagonal; @@ -107,7 +104,7 @@ class Sliced_Wasserstein { } - /** \brief Sliced Wasserstein Kernel constructor. + /** \brief Sliced Wasserstein kernel constructor. * \ingroup Sliced_Wasserstein * * @param[in] _diagram persistence diagram. @@ -115,21 +112,14 @@ class Sliced_Wasserstein { * @param[in] _approx number of directions used to approximate the integral in the Sliced Wasserstein distance, set to -1 for exact computation. * */ - Sliced_Wasserstein(PD _diagram, double _sigma = 1.0, int _approx = 100){diagram = _diagram; approx = _approx; sigma = _sigma; build_rep();} - - PD get_diagram() const {return this->diagram;} - int get_approx() const {return this->approx;} - double get_sigma() const {return this->sigma;} - - - + Sliced_Wasserstein(const Persistence_diagram & _diagram, double _sigma = 1.0, int _approx = 100){diagram = _diagram; approx = _approx; sigma = _sigma; build_rep();} // ********************************** // Utils. // ********************************** // Compute the angle formed by two points of a PD - double compute_angle(PD diag, int i, int j){ + 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; @@ -150,7 +140,7 @@ class Sliced_Wasserstein { } // 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 compute_int_cos(double alpha, double beta) const { double res = 0; if (alpha >= 0 && alpha <= pi){ if (cos(alpha) >= 0){ @@ -175,13 +165,13 @@ class Sliced_Wasserstein { return res; } - 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( (PD1[p].first-PD2[q].first)*(PD1[p].first-PD2[q].first) + (PD1[p].second-PD2[q].second)*(PD1[p].second-PD2[q].second) ); + 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 (PD1[p].first > PD2[q].first) - angle1 = theta1 - asin( (PD1[p].second-PD2[q].second)/norm ); + if (diag1[p].first > diag2[q].first) + angle1 = theta1 - asin( (diag1[p].second-diag2[q].second)/norm ); else - angle1 = theta1 - asin( (PD2[q].second-PD1[p].second)/norm ); + angle1 = theta1 - asin( (diag2[q].second-diag1[p].second)/norm ); double angle2 = angle1 + theta2 - theta1; double integral = compute_int_cos(angle1,angle2); return norm*integral; @@ -197,134 +187,133 @@ class Sliced_Wasserstein { /** \brief Evaluation of the Sliced Wasserstein Distance between a pair of diagrams. * \ingroup Sliced_Wasserstein * + * @pre approx attribute needs to be the same for both instances. * @param[in] second other instance of class Sliced_Wasserstein. - * For warning in red: - * @warning approx parameter needs to be the same for both instances. + * * */ - double compute_sliced_wasserstein_distance(Sliced_Wasserstein second) { + double compute_sliced_wasserstein_distance(const Sliced_Wasserstein & second) const { - GUDHI_CHECK(this->approx != second.approx, std::invalid_argument("Error: different approx values for representations")); + GUDHI_CHECK(this->approx != second.approx, std::invalid_argument("Error: different approx values for representations")); - PD diagram1 = this->diagram; PD diagram2 = second.diagram; double sw = 0; + Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram; double sw = 0; - if(this->approx == -1){ + 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(); - for (int i = 0; i < n2; i++){ - max_ordinate = std::max(max_ordinate, diagram2[i].second); - 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); - 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); - 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); - } + // Add projections onto diagonal. + int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); double max_ordinate = std::numeric_limits<double>::lowest(); + for (int i = 0; i < n2; i++){ + max_ordinate = std::max(max_ordinate, diagram2[i].second); + 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); + 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); + 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); + } - // Compute all angles in both PDs. - std::vector<std::pair<double, std::pair<int,int> > > angles1, angles2; - for (int i = 0; i < num_pts_dgm; i++){ - for (int j = i+1; j < num_pts_dgm; j++){ - double theta1 = compute_angle(diagram1,i,j); double theta2 = compute_angle(diagram2,i,j); - angles1.emplace_back(theta1, std::pair<int,int>(i,j)); - angles2.emplace_back(theta2, std::pair<int,int>(i,j)); - } + // Compute all angles in both PDs. + std::vector<std::pair<double, std::pair<int,int> > > angles1, angles2; + for (int i = 0; i < num_pts_dgm; i++){ + for (int j = i+1; j < num_pts_dgm; j++){ + double theta1 = compute_angle(diagram1,i,j); double theta2 = compute_angle(diagram2,i,j); + angles1.emplace_back(theta1, std::pair<int,int>(i,j)); + angles2.emplace_back(theta2, std::pair<int,int>(i,j)); } + } - // Sort angles. - std::sort(angles1.begin(), angles1.end(), [=](std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);}); - std::sort(angles2.begin(), angles2.end(), [=](std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);}); - - // 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); } ); - - // Find the inverses of the orders. - std::vector<int> order1(num_pts_dgm); std::vector<int> order2(num_pts_dgm); - for(int i = 0; i < num_pts_dgm; i++) for (int j = 0; j < num_pts_dgm; j++) if(orderp1[j] == i){ order1[i] = j; break; } - for(int i = 0; i < num_pts_dgm; i++) for (int j = 0; j < num_pts_dgm; 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(num_pts_dgm); - std::vector<std::vector<std::pair<int,double> > > anglePerm2(num_pts_dgm); - - 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]].emplace_back(p,theta); - anglePerm1[order1[q]].emplace_back(q,theta); - int a = order1[p]; int b = order1[q]; order1[p] = b; order1[q] = a; - } + // Sort angles. + std::sort(angles1.begin(), angles1.end(), [=](std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);}); + std::sort(angles2.begin(), angles2.end(), [=](std::pair<double, std::pair<int,int> >& p1, const std::pair<double, std::pair<int,int> >& p2){return (p1.first < p2.first);}); + + // 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); } ); + + // Find the inverses of the orders. + std::vector<int> order1(num_pts_dgm); std::vector<int> order2(num_pts_dgm); + for(int i = 0; i < num_pts_dgm; i++) for (int j = 0; j < num_pts_dgm; j++) if(orderp1[j] == i){ order1[i] = j; break; } + for(int i = 0; i < num_pts_dgm; i++) for (int j = 0; j < num_pts_dgm; 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(num_pts_dgm); + std::vector<std::vector<std::pair<int,double> > > anglePerm2(num_pts_dgm); + + 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]].emplace_back(p,theta); + anglePerm1[order1[q]].emplace_back(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]].emplace_back(p,theta); - anglePerm2[order2[q]].emplace_back(q,theta); - int a = order2[p]; int b = order2[q]; order2[p] = b; order2[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]].emplace_back(p,theta); + anglePerm2[order2[q]].emplace_back(q,theta); + int a = order2[p]; int b = order2[q]; order2[p] = b; order2[q] = a; + } - for (int i = 0; i < num_pts_dgm; i++){ - anglePerm1[order1[i]].emplace_back(i,pi/2); - anglePerm2[order2[i]].emplace_back(i,pi/2); - } + for (int i = 0; i < num_pts_dgm; i++){ + anglePerm1[order1[i]].emplace_back(i,pi/2); + anglePerm2[order2[i]].emplace_back(i,pi/2); + } - // Compute the SW distance with the list of inversions. - for (int i = 0; i < num_pts_dgm; i++){ - std::vector<std::pair<int,double> > u,v; u = anglePerm1[i]; v = anglePerm2[i]; - double theta1, theta2; theta1 = -pi/2; - unsigned int ku, kv; ku = 0; kv = 0; theta2 = std::min(u[ku].second,v[kv].second); - while(theta1 != pi/2){ - if(diagram1[u[ku].first].first != diagram2[v[kv].first].first || diagram1[u[ku].first].second != diagram2[v[kv].first].second) - if(theta1 != theta2) - sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2); - 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); - } + // Compute the SW distance with the list of inversions. + for (int i = 0; i < num_pts_dgm; i++){ + std::vector<std::pair<int,double> > u,v; u = anglePerm1[i]; v = anglePerm2[i]; + double theta1, theta2; theta1 = -pi/2; + unsigned int ku, kv; ku = 0; kv = 0; theta2 = std::min(u[ku].second,v[kv].second); + while(theta1 != pi/2){ + if(diagram1[u[ku].first].first != diagram2[v[kv].first].first || diagram1[u[ku].first].second != diagram2[v[kv].first].second) + if(theta1 != theta2) + sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2); + 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); } } + } - else{ - - double step = pi/this->approx; - - for (int i = 0; i < this->approx; i++){ + else{ - std::vector<double> v1; std::vector<double> l1 = this->projections[i]; std::vector<double> l1bis = second.projections_diagonal[i]; std::merge(l1.begin(), l1.end(), l1bis.begin(), l1bis.end(), std::back_inserter(v1)); - std::vector<double> v2; std::vector<double> l2 = second.projections[i]; std::vector<double> l2bis = this->projections_diagonal[i]; std::merge(l2.begin(), l2.end(), l2bis.begin(), l2bis.end(), std::back_inserter(v2)); - int n = v1.size(); double f = 0; - for (int j = 0; j < n; j++) f += std::abs(v1[j] - v2[j]); - sw += f*step; + double step = pi/this->approx; + for (int i = 0; i < this->approx; i++){ - } + std::vector<double> v1; std::vector<double> l1 = this->projections[i]; std::vector<double> l1bis = second.projections_diagonal[i]; std::merge(l1.begin(), l1.end(), l1bis.begin(), l1bis.end(), std::back_inserter(v1)); + std::vector<double> v2; std::vector<double> l2 = second.projections[i]; std::vector<double> l2bis = this->projections_diagonal[i]; std::merge(l2.begin(), l2.end(), l2bis.begin(), l2bis.end(), std::back_inserter(v2)); + int n = v1.size(); double f = 0; + for (int j = 0; j < n; j++) f += std::abs(v1[j] - v2[j]); + sw += f*step; } + } - return sw/pi; + return sw/pi; } /** \brief Evaluation of the kernel on a pair of diagrams. * \ingroup Sliced_Wasserstein * - * @param[in] second other instance of class Sliced_Wasserstein. Warning: sigma and approx parameters need to be the same for both instances!!! + * @pre approx and sigma attributes need to be the same for both instances. + * @param[in] second other instance of class Sliced_Wasserstein. * */ - double compute_scalar_product(Sliced_Wasserstein second){ + double compute_scalar_product(const Sliced_Wasserstein & second) const { GUDHI_CHECK(this->sigma != second.sigma, std::invalid_argument("Error: different sigma values for representations")); return std::exp(-compute_sliced_wasserstein_distance(second)/(2*this->sigma*this->sigma)); } @@ -332,10 +321,11 @@ class Sliced_Wasserstein { /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. * \ingroup Sliced_Wasserstein * - * @param[in] second other instance of class Sliced_Wasserstein. Warning: sigma and approx parameters need to be the same for both instances!!! + * @pre approx and sigma attributes need to be the same for both instances. + * @param[in] second other instance of class Sliced_Wasserstein. * */ - double distance(Sliced_Wasserstein second) { + double distance(const Sliced_Wasserstein & second) const { GUDHI_CHECK(this->sigma != second.sigma, std::invalid_argument("Error: different sigma values for representations")); return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), 0.5); } @@ -343,9 +333,8 @@ class Sliced_Wasserstein { -}; - -} // namespace Sliced_Wasserstein +}; // class Sliced_Wasserstein +} // namespace Persistence_representations } // namespace Gudhi #endif // SLICED_WASSERSTEIN_H_ |