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

Former-commit-id: 80f084fc990df6e5c6b60ac83514220aba2ceb5c

5 files changed, 98 insertions, 66 deletions

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 $<TARGET_FILE:Persistence_representations_example_heat_maps>)
 install(TARGETS Persistence_representations_example_heat_maps DESTINATION bin)
 
+add_executable ( Sliced_Wasserstein sliced_wasserstein.cpp )
+add_test(NAME Sliced_Wasserstein
+    COMMAND $<TARGET_FILE:Sliced_Wasserstein>)
+install(TARGETS Sliced_Wasserstein DESTINATION bin)
+
+add_executable ( Persistence_weighted_gaussian persistence_weighted_gaussian.cpp )
+add_test(NAME Persistence_weighted_gaussian
+    COMMAND $<TARGET_FILE:Persistence_weighted_gaussian>)
+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<double>();
 using PD = std::vector<std::pair<double,double> >;
+using Weight = std::function<double (std::pair<double,double>) >;
 
 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<double,double> P){
-    if(P.second > P.first)  return 1;
+  static double pss_weight(std::pair<double,double> p){
+    if(p.second > p.first)  return 1;
     else return -1;
   }
 
-  static double arctan_weight(std::pair<double,double> P){
-    return atan(P.second - P.first);
+  static double arctan_weight(std::pair<double,double> p){
+    return atan(p.second - p.first);
   }
 
-  template<class Weight = std::function<double (std::pair<double,double>) > >
-  std::vector<std::pair<double,double> > Fourier_feat(PD D, std::vector<std::pair<double,double> > Z, Weight weight = arctan_weight){
-    int m = D.size(); std::vector<std::pair<double,double> > 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<std::pair<double,double> > Fourier_feat(PD diag, std::vector<std::pair<double,double> > z, Weight weight = arctan_weight){
+    int md = diag.size(); std::vector<std::pair<double,double> > 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<std::pair<double,double> > random_Fourier(double sigma, int M = 1000){
-    std::normal_distribution<double> distrib(0,1); std::vector<std::pair<double,double> > Z; std::random_device rd;
-    for(int i = 0; i < M; i++){
+  std::vector<std::pair<double,double> > random_Fourier(double sigma, int m = 1000){
+    std::normal_distribution<double> distrib(0,1); std::vector<std::pair<double,double> > 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<class Weight = std::function<double (std::pair<double,double>) > >
-  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<std::pair<double,double> > z =  random_Fourier(sigma, m);
-      std::vector<std::pair<double,double> > b1 = Fourier_feat(diagram1,z,weight);
-      std::vector<std::pair<double,double> > 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<std::pair<double,double> > z =  random_Fourier(this->sigma, this->approx);
+      std::vector<std::pair<double,double> > b1 = Fourier_feat(diagram1,z,this->weight);
+      std::vector<std::pair<double,double> > 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<class Weight = std::function<double (std::pair<double,double>) > >
-  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<double>::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<std::pair<int,double> > 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<std::pair<int,double> > 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);
   }