diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h index e3ed2f6a20..163ed8ebbc 100644 --- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -22,6 +22,8 @@ #include // for std::abs, std::sqrt #include // for std::invalid_argument #include // for std::random_device +#include // for nextafter +#include namespace Gudhi { namespace Persistence_representations { @@ -64,6 +66,7 @@ class Sliced_Wasserstein { Persistence_diagram diagram; int approx; double sigma; + double genericity_factor; std::vector > projections, projections_diagonal; // ********************************** @@ -72,26 +75,22 @@ class Sliced_Wasserstein { void build_rep() { if (approx > 0) { - double step = pi / this->approx; int n = diagram.size(); - + double step = pi / this->approx; for (int i = 0; i < this->approx; i++) { std::vector l, l_diag; for (int j = 0; j < n; j++) { double px = diagram[j].first; double py = diagram[j].second; double proj_diag = (px + py) / 2; - l.push_back(px * cos(-pi / 2 + i * step) + py * sin(-pi / 2 + i * step)); l_diag.push_back(proj_diag * cos(-pi / 2 + i * step) + proj_diag * sin(-pi / 2 + i * step)); } - std::sort(l.begin(), l.end()); std::sort(l_diag.begin(), l_diag.end()); projections.push_back(std::move(l)); projections_diagonal.push_back(std::move(l_diag)); } - diagram.clear(); } } @@ -104,55 +103,75 @@ class Sliced_Wasserstein { 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] - double compute_int_cos(double alpha, double beta) const { - double res = 0; - if (alpha >= 0 && alpha <= pi) { - if (cos(alpha) >= 0) { - if (pi / 2 <= beta) { - res = 2 - sin(alpha) - sin(beta); - } else { - res = sin(beta) - sin(alpha); + // Compute the angle formed by two points of different PDs, as well as the sign of the difference of their ordinates + std::pair compute_inter_angle(const Persistence_diagram& diag1, const Persistence_diagram& diag2, int i, int j) const { + if (diag1[i].second == diag2[j].second) + return std::pair(pi / 2, 1); + else + if (diag1[i].second > diag2[j].second) + return std::pair(atan((diag2[j].first - diag1[i].first) / (diag1[i].second - diag2[j].second)), 1); + else + return std::pair(atan((diag2[j].first - diag1[i].first) / (diag1[i].second - diag2[j].second)), -1); + } + + // compute the integral of between theta1 and theta2 + double close_form_integral(double theta1, double theta2, double diff_x, double diff_y) const { + double int_cos = sin(theta2) - sin(theta1); + double int_sin = cos(theta1) - cos(theta2); + return diff_x * int_cos + diff_y * int_sin; + } + + double compute_integral(double theta1, double theta2, std::vector orderp1, std::vector orderp2, std::vector > > angles12, const Persistence_diagram& diag1, + const Persistence_diagram& diag2) const { + double integral = 0; + for (int i = 0; i < orderp1.size(); i++){ + + int dgm1_idx = orderp1[i]; + int dgm2_idx = orderp2[i]; + double diff_x = diag1[dgm1_idx].first - diag2[dgm2_idx].first; + double diff_y = diag1[dgm1_idx].second - diag2[dgm2_idx].second; + double angle12 = angles12[dgm1_idx][dgm2_idx].first; + int sign12 = angles12[dgm1_idx][dgm2_idx].second; + double small_integral; + + //std::cout << theta1 << " " << theta2 << " " << sign12 << " " << angle12 << std::endl; + + if (sign12 == 1){ + if (angle12 < theta1){ // || = within [theta1, theta2] + small_integral = this->close_form_integral(theta1, theta2, diff_x, diff_y); } - } else { - if (1.5 * pi <= beta) { - res = 2 + sin(alpha) + sin(beta); - } else { - res = sin(alpha) - sin(beta); + else{ + if (theta2 < angle12){ // || = - within [theta1, theta2] + small_integral = -this->close_form_integral(theta1, theta2, diff_x, diff_y); + } + else{ // || = - within [theta1, angle12] and || = within [angle12, theta2] + //std::cout << "breakdown: " << -this->close_form_integral(theta1, angle12, diff_x, diff_y) << " + " << this->close_form_integral(angle12, theta2, diff_x, diff_y) << std::endl; + small_integral = -this->close_form_integral(theta1, angle12, diff_x, diff_y) + this->close_form_integral(angle12, theta2, diff_x, diff_y); + } } } - } - if (alpha >= -pi && alpha <= 0) { - if (cos(alpha) <= 0) { - if (-pi / 2 <= beta) { - res = 2 + sin(alpha) + sin(beta); - } else { - res = sin(alpha) - sin(beta); + else{ + if (angle12 < theta1){ // || = - within [theta1, theta2] + small_integral = -this->close_form_integral(theta1, theta2, diff_x, diff_y); } - } else { - if (pi / 2 <= beta) { - res = 2 - sin(alpha) - sin(beta); - } else { - res = sin(beta) - sin(alpha); + else{ + if (theta2 < angle12){ // || = within [theta1, theta2] + small_integral = this->close_form_integral(theta1, theta2, diff_x, diff_y); + } + else{ // || = within [theta1, angle12] and || = - within [angle12, theta2] + //std::cout << "breakdown: " << this->close_form_integral(theta1, angle12, diff_x, diff_y) << " + " << - this->close_form_integral(angle12, theta2, diff_x, diff_y) << std::endl; + small_integral = this->close_form_integral(theta1, angle12, diff_x, diff_y) - this->close_form_integral(angle12, theta2, diff_x, diff_y); + } } } + + //std::cout << "small int = " << small_integral << std::endl; + integral += small_integral; + } - return res; - } - 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 - 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; + return integral; + } // Evaluation of the Sliced Wasserstein Distance between a pair of diagrams. @@ -162,68 +181,80 @@ class Sliced_Wasserstein { Persistence_diagram diagram1 = this->diagram; Persistence_diagram diagram2 = second.diagram; + double sw = 0; if (this->approx == -1) { - // Add projections onto diagonal. + int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); - double min_ordinate = std::numeric_limits::max(); - double min_abscissa = std::numeric_limits::max(); - double max_ordinate = std::numeric_limits::lowest(); - double max_abscissa = std::numeric_limits::lowest(); + + // Put diagrams in generic positions. + for (int i = 0; i < n1; i++) {diagram1[i].first += (2*i)*genericity_factor*DBL_EPSILON; diagram1[i].second += (2*i+1)*genericity_factor*DBL_EPSILON;} + for (int i = n1; i < n1+n2; i++) {diagram2[i-n1].first += (2*i)*genericity_factor*DBL_EPSILON; diagram2[i-n1].second += (2*i+1)*genericity_factor*DBL_EPSILON;} +// std::cout << (2*(n1+n2-1)+1)*10*DBL_EPSILON << std::endl; + + // Add projections onto diagonal. for (int i = 0; i < n2; i++) { - min_ordinate = std::min(min_ordinate, diagram2[i].second); - min_abscissa = std::min(min_abscissa, diagram2[i].first); - max_ordinate = std::max(max_ordinate, diagram2[i].second); - max_abscissa = std::max(max_abscissa, diagram2[i].first); - diagram1.emplace_back((diagram2[i].first + diagram2[i].second) / 2, - (diagram2[i].first + diagram2[i].second) / 2); + double proj_i = (diagram2[i].first + diagram2[i].second) / 2; + diagram1.emplace_back(proj_i,proj_i); } for (int i = 0; i < n1; i++) { - min_ordinate = std::min(min_ordinate, diagram1[i].second); - min_abscissa = std::min(min_abscissa, diagram1[i].first); - max_ordinate = std::max(max_ordinate, diagram1[i].second); - max_abscissa = std::max(max_abscissa, diagram1[i].first); - diagram2.emplace_back((diagram1[i].first + diagram1[i].second) / 2, - (diagram1[i].first + diagram1[i].second) / 2); + double proj_i = (diagram1[i].first + diagram1[i].second) / 2; + diagram2.emplace_back(proj_i, proj_i); } int num_pts_dgm = diagram1.size(); - // Slightly perturb the points so that the PDs are in generic positions. - double epsilon = 0.0001; - double thresh_y = (max_ordinate - min_ordinate) * epsilon; - double thresh_x = (max_abscissa - min_abscissa) * epsilon; - std::random_device rd; - std::default_random_engine re(rd()); - std::uniform_real_distribution uni(-1, 1); - for (int i = 0; i < num_pts_dgm; i++) { - double u = uni(re); - diagram1[i].first += u * thresh_x; - diagram1[i].second += u * thresh_y; - diagram2[i].first += u * thresh_x; - diagram2[i].second += u * thresh_y; - } - - // Compute all angles in both PDs. - std::vector > > angles1, angles2; + // Compute all intra-PD angles with the norms of the corresponding points (for sorting these angles later). + std::vector, std::pair > > > 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(i, j)); - angles2.emplace_back(theta2, std::pair(i, j)); + double p1ix = diagram1[i].first; double p1iy = diagram1[i].second; double p1jx = diagram1[j].first; double p1jy = diagram1[j].second; + double p2ix = diagram2[i].first; double p2iy = diagram2[i].second; double p2jx = diagram2[j].first; double p2jy = diagram2[j].second; + double n1i = p1ix*p1ix + p1iy*p1iy; double n1j = p1jx*p1jx + p1jy*p1jy; double n2i = p2ix*p2ix + p2iy*p2iy; double n2j = p2jx*p2jx + p2jy*p2jy; + double m1 = std::min(n1i,n1j); double M1 = std::max(n1i,n1j); double m2 = std::min(n2i,n2j); double M2 = std::max(n2i,n2j); + std::pair iM1(i,M1); std::pair jm1(j,m1); std::pair iM2(i,M2); std::pair jm2(j,m2); + angles1.emplace_back(theta1, std::pair, std::pair>(iM1, jm1)); + angles2.emplace_back(theta2, std::pair, std::pair>(iM2, jm2)); + } + } + int num_angles = angles1.size(); + + // Compute all inter-PD angles. + std::vector > > angles12(num_pts_dgm, std::vector >(num_pts_dgm)); + for (int i = 0; i < num_pts_dgm; i++) { + for (int j = 0; j < num_pts_dgm; j++) { + std::pair theta12 = compute_inter_angle(diagram1, diagram2, i, j); + angles12[i][j] = theta12; } } // Sort angles. std::sort(angles1.begin(), angles1.end(), - [](const std::pair >& p1, - const std::pair >& p2) { return (p1.first < p2.first); }); + [](const std::pair, std::pair> >& p1, + const std::pair, std::pair> >& p2) { + if (p1.first == p2.first){ + if (p1.second.first.second == p2.second.first.second) + return (p1.second.second.second < p2.second.second.second); + else + return (p1.second.first.second < p2.second.first.second); + } else + return (p1.first < p2.first); + }); std::sort(angles2.begin(), angles2.end(), - [](const std::pair >& p1, - const std::pair >& p2) { return (p1.first < p2.first); }); + [](const std::pair, std::pair> >& p1, + const std::pair, std::pair> >& p2) { + if (p1.first == p2.first){ + if (p1.second.first.second == p2.second.first.second) + return (p1.second.second.second < p2.second.second.second); + else + return (p1.second.first.second < p2.second.first.second); + } else + return (p1.first < p2.first); + }); // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2). std::vector orderp1, orderp2; @@ -233,15 +264,15 @@ class Sliced_Wasserstein { } std::sort(orderp1.begin(), orderp1.end(), [&](int i, int j) { if (diagram1[i].second != diagram1[j].second) - return (diagram1[i].second < diagram1[j].second); + return (diagram1[i].second > diagram1[j].second); else - return (diagram1[i].first > diagram1[j].first); + 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); + return (diagram2[i].second > diagram2[j].second); else - return (diagram2[i].first > diagram2[j].first); + return (diagram2[i].first < diagram2[j].first); }); // Find the inverses of the orders. @@ -252,62 +283,97 @@ class Sliced_Wasserstein { order2[orderp2[i]] = i; } - // Record all inversions of points in the orders as theta varies along the positive half-disk. - std::vector > > anglePerm1(num_pts_dgm); - std::vector > > 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; - } + double theta1 = -pi/2; + double theta2; + int try_idx1 = 0; + int try_idx2 = 0; + int order_to_swap, i_to_swap, j_to_swap; + + while( (try_idx1 < num_angles) || (try_idx2 < num_angles)){ + + // Find theta2 + next indices to swap + if ((try_idx1 < num_angles) && (try_idx2 < num_angles)){ + if ((angles1[try_idx1].first < angles2[try_idx2].first) && (try_idx1 < num_angles)){ + theta2 = angles1[try_idx1].first; + order_to_swap = 1; + i_to_swap = angles1[try_idx1].second.first.first; + j_to_swap = angles1[try_idx1].second.second.first; + try_idx1 += 1; + } + else{ + theta2 = angles2[try_idx2].first; + order_to_swap = 2; + i_to_swap = angles2[try_idx2].second.first.first; + j_to_swap = angles2[try_idx2].second.second.first; + try_idx2 += 1; + } + } + else{ + if (try_idx1 == num_angles){ + theta2 = angles2[try_idx2].first; + order_to_swap = 2; + i_to_swap = angles2[try_idx2].second.first.first; + j_to_swap = angles2[try_idx2].second.second.first; + try_idx2 += 1; + } + else{ + theta2 = angles1[try_idx1].first; + order_to_swap = 1; + i_to_swap = angles1[try_idx1].second.first.first; + j_to_swap = angles1[try_idx1].second.second.first; + try_idx1 += 1; + } + } - 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; - } + // compute integral + double integral_value = compute_integral(theta1, theta2, orderp1, orderp2, angles12, diagram1, diagram2); + sw += integral_value; + +// std::cout << theta1 << " " << theta2 << std::endl; +// std::cout << "orderp1 = [ "; +// for (int i = 0; i < orderp1.size(); i++) std::cout << orderp1[i] << " "; +// std::cout << "]" << std::endl; +// std::cout << "orderp2 = [ "; +// for (int i = 0; i < orderp2.size(); i++) std::cout << orderp2[i] << " "; +// std::cout << "]" << std::endl; +// std::cout << "int = " << integral_value << std::endl; + + // swap order + if (order_to_swap == 1){ + orderp1[order1[i_to_swap]] = j_to_swap; + orderp1[order1[j_to_swap]] = i_to_swap; + int tmp_pos = order1[i_to_swap]; + order1[i_to_swap] = order1[j_to_swap]; + order1[j_to_swap] = tmp_pos; + + } + else{ + orderp2[order2[i_to_swap]] = j_to_swap; + orderp2[order2[j_to_swap]] = i_to_swap; + int tmp_pos = order2[i_to_swap]; + order2[i_to_swap] = order2[j_to_swap]; + order2[j_to_swap] = tmp_pos; + } - 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); - } + // update theta1 + theta1 = theta2; - // Compute the SW distance with the list of inversions. - for (int i = 0; i < num_pts_dgm; i++) { - std::vector > 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); - } } + + // Last integral between last angle and pi / 2 + theta2 = pi / 2; + double integral_value = compute_integral(theta1, theta2, orderp1, orderp2, angles12, diagram1, diagram2); + sw += integral_value; + +// std::cout << theta1 << " " << theta2 << std::endl; +// std::cout << "orderp1 = [ "; +// for (int i = 0; i < orderp1.size(); i++) std::cout << orderp1[i] << " "; +// std::cout << "]" << std::endl; +// std::cout << "orderp2 = [ "; +// for (int i = 0; i < orderp2.size(); i++) std::cout << orderp2[i] << " "; +// std::cout << "]" << std::endl; +// std::cout << "int = " << integral_value << std::endl; + } else { double step = pi / this->approx; std::vector v1, v2; @@ -327,6 +393,7 @@ class Sliced_Wasserstein { } return sw / pi; + } public: @@ -341,8 +408,8 @@ class Sliced_Wasserstein { * directions are stored in memory to reduce computation time. * */ - Sliced_Wasserstein(const Persistence_diagram& _diagram, double _sigma = 1.0, int _approx = 10) - : diagram(_diagram), approx(_approx), sigma(_sigma) { + Sliced_Wasserstein(const Persistence_diagram& _diagram, double _sigma = 1.0, int _approx = 10, double _genericity_factor=10) + : diagram(_diagram), approx(_approx), sigma(_sigma), genericity_factor(_genericity_factor) { build_rep(); } @@ -373,6 +440,10 @@ class Sliced_Wasserstein { 2 * this->compute_scalar_product(second)); } + double sw_distance(const Sliced_Wasserstein& second) const { + return this->compute_sliced_wasserstein_distance(second); + } + }; // class Sliced_Wasserstein } // namespace Persistence_representations } // namespace Gudhi diff --git a/src/Persistence_representations/test/kernels.cpp b/src/Persistence_representations/test/kernels.cpp index 63089538e6..fed7dc7ef3 100644 --- a/src/Persistence_representations/test/kernels.cpp +++ b/src/Persistence_representations/test/kernels.cpp @@ -43,7 +43,55 @@ BOOST_AUTO_TEST_CASE(check_PWG) { BOOST_AUTO_TEST_CASE(check_SW) { Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); - SW sw1(v1, 1.0, 100); SW swex1(v1, 1.0, -1); - SW sw2(v2, 1.0, 100); SW swex2(v2, 1.0, -1); - BOOST_CHECK(std::abs(sw1.compute_scalar_product(sw2) - swex1.compute_scalar_product(swex2)) <= 1e-1); + SW sw1(v1, 1.0, 1000); SW swex1(v1, 1.0, -1); + SW sw2(v2, 1.0, 1000); SW swex2(v2, 1.0, -1); + std::cout << sw1.sw_distance(sw2) << std::endl; + std::cout << swex1.sw_distance(swex2) << std::endl; + std::cout << sw1.compute_scalar_product(sw2) << std::endl; + std::cout << swex1.compute_scalar_product(swex2) << std::endl; + BOOST_CHECK(std::abs(sw1.compute_scalar_product(sw2) - swex1.compute_scalar_product(swex2)) <= 1e-3); +} + +BOOST_AUTO_TEST_CASE(check_multiplicity_SW) { + Persistence_diagram v1, v2; + v1.emplace_back(0,1); v2.emplace_back(0,2); + v1.emplace_back(2,5); v2.emplace_back(1,4); + v1.emplace_back(3,5); v2.emplace_back(1,3); + v1.emplace_back(0,1); v2.emplace_back(1,3); + SW sw1(v1, 1.0, 10000); SW swex1(v1, 1.0, -1); + SW sw2(v2, 1.0, 10000); SW swex2(v2, 1.0, -1); + std::cout << sw1.sw_distance(sw2) << std::endl; + std::cout << swex1.sw_distance(swex2) << std::endl; + std::cout << sw1.compute_scalar_product(sw2) << std::endl; + std::cout << swex1.compute_scalar_product(swex2) << std::endl; + BOOST_CHECK(std::abs(sw1.compute_scalar_product(sw2) - swex1.compute_scalar_product(swex2)) <= 1e-3); +} + + +BOOST_AUTO_TEST_CASE(check_generic_SW) { + Persistence_diagram v1, v2; + double lower_bound = 0; + double upper_bound = 100; + std::random_device rd; // Will be used to obtain a seed for the random number engine + std::mt19937 gen(rd()); + std::uniform_real_distribution unif(lower_bound,upper_bound); + for (int i = 0; i < 10; i++){ + double ax_random_double = unif(gen); + double ay_random_double = unif(gen); + //std::cout << ax_random_double << " " << ax_random_double+ay_random_double << std::endl; + v1.emplace_back(ax_random_double,ax_random_double+ay_random_double); + } + for (int i = 0; i < 10; i++){ + double ax_random_double = unif(gen); + double ay_random_double = unif(gen); + //std::cout << ax_random_double << " " << ax_random_double+ay_random_double << std::endl; + v2.emplace_back(ax_random_double,ax_random_double+ay_random_double); + } + SW sw1(v1, 100.0, 10000); SW swex1(v1, 100.0, -1); + SW sw2(v2, 100.0, 10000); SW swex2(v2, 100.0, -1); + std::cout << sw1.sw_distance(sw2) << std::endl; + std::cout << swex1.sw_distance(swex2) << std::endl; + std::cout << sw1.compute_scalar_product(sw2) << std::endl; + std::cout << swex1.compute_scalar_product(swex2) << std::endl; + BOOST_CHECK(std::abs(sw1.compute_scalar_product(sw2) - swex1.compute_scalar_product(swex2)) <= 1e-3); }