From ff0dc023588e3b33bc4bc7f26ce1f68c647ae441 Mon Sep 17 00:00:00 2001 From: mcarrier Date: Fri, 16 Feb 2018 15:43:29 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3251 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 80f084fc990df6e5c6b60ac83514220aba2ceb5c --- .../example/CMakeLists.txt | 10 +++ .../example/persistence_weighted_gaussian.cpp | 34 +++++----- .../example/sliced_wasserstein.cpp | 16 +++-- .../include/gudhi/Persistence_weighted_gaussian.h | 78 ++++++++++++---------- .../include/gudhi/Sliced_Wasserstein.h | 26 +++++--- 5 files changed, 98 insertions(+), 66 deletions(-) (limited to 'src/Persistence_representations') diff --git a/src/Persistence_representations/example/CMakeLists.txt b/src/Persistence_representations/example/CMakeLists.txt index 54d719ac..79d39c4d 100644 --- a/src/Persistence_representations/example/CMakeLists.txt +++ b/src/Persistence_representations/example/CMakeLists.txt @@ -27,3 +27,13 @@ add_test(NAME Persistence_representations_example_heat_maps COMMAND $) install(TARGETS Persistence_representations_example_heat_maps DESTINATION bin) +add_executable ( Sliced_Wasserstein sliced_wasserstein.cpp ) +add_test(NAME Sliced_Wasserstein + COMMAND $) +install(TARGETS Sliced_Wasserstein DESTINATION bin) + +add_executable ( Persistence_weighted_gaussian persistence_weighted_gaussian.cpp ) +add_test(NAME Persistence_weighted_gaussian + COMMAND $) +install(TARGETS Persistence_weighted_gaussian DESTINATION bin) + diff --git a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp index e95b9445..a0e820ea 100644 --- a/src/Persistence_representations/example/persistence_weighted_gaussian.cpp +++ b/src/Persistence_representations/example/persistence_weighted_gaussian.cpp @@ -44,21 +44,24 @@ int main(int argc, char** argv) { persistence2.push_back(std::make_pair(3, 5)); persistence2.push_back(std::make_pair(6, 10)); - PWG PWG1(persistence1); - PWG PWG2(persistence2); double sigma = 1; double tau = 1; - int m = 1000; + int m = 10000; + PWG PWG1(persistence1, sigma, m, PWG::arctan_weight); + PWG PWG2(persistence2, sigma, m, PWG::arctan_weight); + + PWG PWGex1(persistence1, sigma, -1, PWG::arctan_weight); + PWG PWGex2(persistence2, sigma, -1, PWG::arctan_weight); // Linear PWG - std::cout << PWG1.compute_scalar_product (PWG2, sigma, PWG::arctan_weight, m) << std::endl; - std::cout << PWG1.compute_scalar_product (PWG2, sigma, PWG::arctan_weight, -1) << std::endl; + std::cout << PWG1.compute_scalar_product (PWG2) << std::endl; + std::cout << PWGex1.compute_scalar_product (PWGex2) << std::endl; - std::cout << PWG1.distance (PWG2, sigma, PWG::arctan_weight, m) << std::endl; - std::cout << PWG1.distance (PWG2, sigma, PWG::arctan_weight, -1) << std::endl; + std::cout << PWG1.distance (PWG2) << std::endl; + std::cout << PWGex1.distance (PWGex2) << std::endl; @@ -68,8 +71,8 @@ int main(int argc, char** argv) { // Gaussian PWG - std::cout << std::exp( -PWG1.distance (PWG2, sigma, PWG::arctan_weight, m, 2) ) / (2*tau*tau) << std::endl; - std::cout << std::exp( -PWG1.distance (PWG2, sigma, PWG::arctan_weight, -1, 2) ) / (2*tau*tau) << std::endl; + std::cout << std::exp( -PWG1.distance (PWG2, 2) ) / (2*tau*tau) << std::endl; + std::cout << std::exp( -PWGex1.distance (PWGex2, 2) ) / (2*tau*tau) << std::endl; @@ -82,14 +85,15 @@ int main(int argc, char** argv) { PD pd1 = persistence1; int numpts = persistence1.size(); for(int i = 0; i < numpts; i++) pd1.emplace_back(persistence1[i].second,persistence1[i].first); PD pd2 = persistence2; numpts = persistence2.size(); for(int i = 0; i < numpts; i++) pd2.emplace_back(persistence2[i].second,persistence2[i].first); - PWG pwg1(pd1); - PWG pwg2(pd2); + PWG pwg1(pd1, 2*std::sqrt(sigma), m, PWG::pss_weight); + PWG pwg2(pd2, 2*std::sqrt(sigma), m, PWG::pss_weight); + + PWG pwgex1(pd1, 2*std::sqrt(sigma), -1, PWG::pss_weight); + PWG pwgex2(pd2, 2*std::sqrt(sigma), -1, PWG::pss_weight); - std::cout << pwg1.compute_scalar_product (pwg2, 2*std::sqrt(sigma), PWG::pss_weight, m) / (16*pi*sigma) << std::endl; - std::cout << pwg1.compute_scalar_product (pwg2, 2*std::sqrt(sigma), PWG::pss_weight, -1) / (16*pi*sigma) << std::endl; + std::cout << pwg1.compute_scalar_product (pwg2) / (16*pi*sigma) << std::endl; + std::cout << pwgex1.compute_scalar_product (pwgex2) / (16*pi*sigma) << std::endl; - std::cout << pwg1.distance (pwg2, 2*std::sqrt(sigma), PWG::pss_weight, m) / (16*pi*sigma) << std::endl; - std::cout << pwg1.distance (pwg2, 2*std::sqrt(sigma), PWG::pss_weight, -1) / (16*pi*sigma) << std::endl; return 0; diff --git a/src/Persistence_representations/example/sliced_wasserstein.cpp b/src/Persistence_representations/example/sliced_wasserstein.cpp index 673d8474..f153fbe8 100644 --- a/src/Persistence_representations/example/sliced_wasserstein.cpp +++ b/src/Persistence_representations/example/sliced_wasserstein.cpp @@ -43,13 +43,17 @@ int main(int argc, char** argv) { persistence2.push_back(std::make_pair(3, 5)); persistence2.push_back(std::make_pair(6, 10)); - SW SW1(persistence1); - SW SW2(persistence2); - std::cout << SW1.compute_sliced_wasserstein_distance(SW2,100) << std::endl; - std::cout << SW1.compute_sliced_wasserstein_distance(SW2,-1) << std::endl; - std::cout << SW1.compute_scalar_product(SW2,1,100) << std::endl; - std::cout << SW1.distance(SW2,1,100,1) << std::endl; + SW sw1(persistence1, 1, 100); + SW sw2(persistence2, 1, 100); + + SW swex1(persistence1, 1, -1); + SW swex2(persistence2, 1, -1); + + std::cout << sw1.compute_sliced_wasserstein_distance(sw2) << std::endl; + std::cout << swex1.compute_sliced_wasserstein_distance(swex2) << std::endl; + std::cout << sw1.compute_scalar_product(sw2) << std::endl; + std::cout << swex1.distance(swex2) << std::endl; return 0; } diff --git a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h index 2884885c..2b25b9a8 100644 --- a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h +++ b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h @@ -45,6 +45,7 @@ double pi = boost::math::constants::pi(); using PD = std::vector >; +using Weight = std::function) >; namespace Gudhi { namespace Persistence_representations { @@ -53,11 +54,18 @@ class Persistence_weighted_gaussian{ protected: PD diagram; + Weight weight; + double sigma; + int approx; public: - Persistence_weighted_gaussian(PD _diagram){diagram = _diagram;} + Persistence_weighted_gaussian(PD _diagram){diagram = _diagram; sigma = 1.0; approx = 1000; weight = arctan_weight;} + Persistence_weighted_gaussian(PD _diagram, double _sigma, int _approx, Weight _weight){diagram = _diagram; sigma = _sigma; approx = _approx; weight = _weight;} PD get_diagram(){return this->diagram;} + double get_sigma(){return this->sigma;} + int get_approx(){return this->approx;} + Weight get_weight(){return this->weight;} // ********************************** @@ -65,38 +73,37 @@ class Persistence_weighted_gaussian{ // ********************************** - static double pss_weight(std::pair P){ - if(P.second > P.first) return 1; + static double pss_weight(std::pair p){ + if(p.second > p.first) return 1; else return -1; } - static double arctan_weight(std::pair P){ - return atan(P.second - P.first); + static double arctan_weight(std::pair p){ + return atan(p.second - p.first); } - template) > > - std::vector > Fourier_feat(PD D, std::vector > Z, Weight weight = arctan_weight){ - int m = D.size(); std::vector > B; int M = Z.size(); - for(int i = 0; i < M; i++){ - double d1 = 0; double d2 = 0; double zx = Z[i].first; double zy = Z[i].second; - for(int j = 0; j < m; j++){ - double x = D[j].first; double y = D[j].second; - d1 += weight(D[j])*cos(x*zx + y*zy); - d2 += weight(D[j])*sin(x*zx + y*zy); + std::vector > Fourier_feat(PD diag, std::vector > z, Weight weight = arctan_weight){ + int md = diag.size(); std::vector > b; int mz = z.size(); + for(int i = 0; i < mz; i++){ + double d1 = 0; double d2 = 0; double zx = z[i].first; double zy = z[i].second; + for(int j = 0; j < md; j++){ + double x = diag[j].first; double y = diag[j].second; + d1 += weight(diag[j])*cos(x*zx + y*zy); + d2 += weight(diag[j])*sin(x*zx + y*zy); } - B.emplace_back(d1,d2); + b.emplace_back(d1,d2); } - return B; + return b; } - std::vector > random_Fourier(double sigma, int M = 1000){ - std::normal_distribution distrib(0,1); std::vector > Z; std::random_device rd; - for(int i = 0; i < M; i++){ + std::vector > random_Fourier(double sigma, int m = 1000){ + std::normal_distribution distrib(0,1); std::vector > z; std::random_device rd; + for(int i = 0; i < m; i++){ std::mt19937 e1(rd()); std::mt19937 e2(rd()); double zx = distrib(e1); double zy = distrib(e2); - Z.emplace_back(zx/sigma,zy/sigma); + z.emplace_back(zx/sigma,zy/sigma); } - return Z; + return z; } @@ -106,32 +113,33 @@ class Persistence_weighted_gaussian{ // ********************************** - template) > > - double compute_scalar_product(Persistence_weighted_gaussian second, double sigma, Weight weight = arctan_weight, int m = 1000){ + double compute_scalar_product(Persistence_weighted_gaussian second){ PD diagram1 = this->diagram; PD diagram2 = second.diagram; - if(m == -1){ + if(this->approx == -1){ int num_pts1 = diagram1.size(); int num_pts2 = diagram2.size(); double k = 0; for(int i = 0; i < num_pts1; i++) for(int j = 0; j < num_pts2; j++) - k += weight(diagram1[i])*weight(diagram2[j])*exp(-((diagram1[i].first - diagram2[j].first) * (diagram1[i].first - diagram2[j].first) + - (diagram1[i].second - diagram2[j].second) * (diagram1[i].second - diagram2[j].second)) - /(2*sigma*sigma)); + k += this->weight(diagram1[i])*this->weight(diagram2[j])*exp(-((diagram1[i].first - diagram2[j].first) * (diagram1[i].first - diagram2[j].first) + + (diagram1[i].second - diagram2[j].second) * (diagram1[i].second - diagram2[j].second)) + /(2*this->sigma*this->sigma)); return k; } else{ - std::vector > z = random_Fourier(sigma, m); - std::vector > b1 = Fourier_feat(diagram1,z,weight); - std::vector > b2 = Fourier_feat(diagram2,z,weight); - double d = 0; for(int i = 0; i < m; i++) d += b1[i].first*b2[i].first + b1[i].second*b2[i].second; - return d/m; + std::vector > z = random_Fourier(this->sigma, this->approx); + std::vector > b1 = Fourier_feat(diagram1,z,this->weight); + std::vector > b2 = Fourier_feat(diagram2,z,this->weight); + double d = 0; for(int i = 0; i < this->approx; i++) d += b1[i].first*b2[i].first + b1[i].second*b2[i].second; + return d/this->approx; } } - template) > > - double distance(Persistence_weighted_gaussian second, double sigma, Weight weight = arctan_weight, int m = 1000, double power = 1) { - return std::pow(this->compute_scalar_product(*this, sigma, weight, m) + second.compute_scalar_product(second, sigma, weight, m)-2*this->compute_scalar_product(second, sigma, weight, m), power/2.0); + double distance(Persistence_weighted_gaussian second, double power = 1) { + if(this->sigma != second.get_sigma() || this->approx != second.get_approx()){ + std::cout << "Error: different representations!" << std::endl; return 0; + } + else return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), power/2.0); } diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h index 4fa6151f..ad1a6c42 100644 --- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -53,11 +53,16 @@ class Sliced_Wasserstein { protected: PD diagram; + int approx; + double sigma; public: - Sliced_Wasserstein(PD _diagram){diagram = _diagram;} + Sliced_Wasserstein(PD _diagram){diagram = _diagram; approx = 100; sigma = 0.001;} + Sliced_Wasserstein(PD _diagram, double _sigma, int _approx){diagram = _diagram; approx = _approx; sigma = _sigma;} PD get_diagram(){return this->diagram;} + int get_approx(){return this->approx;} + double get_sigma(){return this->sigma;} // ********************************** @@ -130,11 +135,11 @@ class Sliced_Wasserstein { // Scalar product + distance. // ********************************** - double compute_sliced_wasserstein_distance(Sliced_Wasserstein second, int approx) { + double compute_sliced_wasserstein_distance(Sliced_Wasserstein second) { PD diagram1 = this->diagram; PD diagram2 = second.diagram; double sw = 0; - if(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::lowest(); @@ -226,7 +231,7 @@ class Sliced_Wasserstein { else{ - double step = pi/approx; + double step = pi/this->approx; // Add projections onto diagonal. int n1, n2; n1 = diagram1.size(); n2 = diagram2.size(); @@ -238,7 +243,7 @@ class Sliced_Wasserstein { // Sort and compare all projections. #ifdef GUDHI_USE_TBB - tbb::parallel_for(0, approx, [&](int i){ + tbb::parallel_for(0, this->approx, [&](int i){ std::vector > l1, l2; for (int j = 0; j < n; j++){ l1.emplace_back( j, diagram1[j].first*cos(-pi/2+i*step) + diagram1[j].second*sin(-pi/2+i*step) ); @@ -250,7 +255,7 @@ class Sliced_Wasserstein { sw += f*step; }); #else - for (int i = 0; i < approx; i++){ + for (int i = 0; i < this->approx; i++){ std::vector > l1, l2; for (int j = 0; j < n; j++){ l1.emplace_back( j, diagram1[j].first*cos(-pi/2+i*step) + diagram1[j].second*sin(-pi/2+i*step) ); @@ -268,12 +273,13 @@ class Sliced_Wasserstein { } - double compute_scalar_product(Sliced_Wasserstein second, double sigma, int approx = 100) { - return std::exp(-compute_sliced_wasserstein_distance(second, approx)/(2*sigma*sigma)); + double compute_scalar_product(Sliced_Wasserstein second){ + return std::exp(-compute_sliced_wasserstein_distance(second)/(2*this->sigma*this->sigma)); } - double distance(Sliced_Wasserstein second, double sigma, int approx = 100, double power = 1) { - return std::pow(this->compute_scalar_product(*this, sigma, approx) + second.compute_scalar_product(second, sigma, approx)-2*this->compute_scalar_product(second, sigma, approx), power/2.0); + double distance(Sliced_Wasserstein second, double power = 1) { + if(this->sigma != second.sigma || this->approx != second.approx){std::cout << "Error: different representations!" << std::endl; return 0;} + else return std::pow(this->compute_scalar_product(*this) + second.compute_scalar_product(second)-2*this->compute_scalar_product(second), power/2.0); } -- cgit v1.2.3