summaryrefslogtreecommitdiff
path: root/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
diff options
context:
space:
mode:
authormcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb>2018-04-23 15:22:13 +0000
committermcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb>2018-04-23 15:22:13 +0000
commit541284f6f1bf7d4a76daac8a52850c7162a765cb (patch)
tree2ebecae35daf30dbcfb9c683f9587b93e117f443 /src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h
parent5e24206f945f66575c7c179d74e9661cf60ca3df (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.h243
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_