From 21e6cb713dd3db78e68b4140ab2d69508dad01af Mon Sep 17 00:00:00 2001 From: mcarrier Date: Fri, 22 Dec 2017 15:53:56 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3103 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f7ac7ffce4fbfb04b099797a02999d4a93e3f22b --- src/Kernels/example/kernel_basic_example.cpp | 49 +++++++ src/Kernels/include/gudhi/PSS.h | 4 +- src/Kernels/include/gudhi/SW.h | 198 +++++++++++++++------------ src/Kernels/include/gudhi/figtree-0.9.3.zip | Bin 0 -> 1229617 bytes src/cmake/modules/GUDHI_modules.cmake | 4 +- 5 files changed, 166 insertions(+), 89 deletions(-) create mode 100644 src/Kernels/example/kernel_basic_example.cpp create mode 100644 src/Kernels/include/gudhi/figtree-0.9.3.zip diff --git a/src/Kernels/example/kernel_basic_example.cpp b/src/Kernels/example/kernel_basic_example.cpp new file mode 100644 index 00000000..0a8d83b3 --- /dev/null +++ b/src/Kernels/example/kernel_basic_example.cpp @@ -0,0 +1,49 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * Copyright (C) 2015 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +//#include +//#include + +#include +#include +#include // for pair +#include // for numeric_limits + +int main() { + std::vector< std::pair > v1, v2; + + v1.emplace_back(2.7, 3.7); + v1.emplace_back(9.6, 14.); + v1.emplace_back(34.2, 34.974); + + v2.emplace_back(2.8, 4.45); + v2.emplace_back(9.5, 14.1); + + + double b1 = Gudhi::sliced_wasserstein::compute_approximate_SW (v1, v2); + double b2 = Gudhi::sliced_wasserstein::compute_exact_SW (v1, v2); + + std::cout << "Approximate Sliced Wasserstein distance = " << b1 << std::endl; + std::cout << "Exact Sliced Wasserstein distance = " << b2 << std::endl; + +} diff --git a/src/Kernels/include/gudhi/PSS.h b/src/Kernels/include/gudhi/PSS.h index 70743c47..5111a09f 100644 --- a/src/Kernels/include/gudhi/PSS.h +++ b/src/Kernels/include/gudhi/PSS.h @@ -56,8 +56,8 @@ #include #include -#include "../../figtree-0.9.3/include/figtree.h" -#include "../../figtree-0.9.3/external/ann_1.1.1/include/ANN/ANN.h" +#include "figtree.h" +#include "ANN.h" using PD = std::vector >; diff --git a/src/Kernels/include/gudhi/SW.h b/src/Kernels/include/gudhi/SW.h index 6871d990..0b041252 100644 --- a/src/Kernels/include/gudhi/SW.h +++ b/src/Kernels/include/gudhi/SW.h @@ -55,47 +55,36 @@ using PD = std::vector >; -std::vector > PDi, PDj; - -bool compOri(const int& p, const int& q){ - if(PDi[p].second != PDi[q].second) - return (PDi[p].second < PDi[q].second); - else - return (PDi[p].first > PDi[q].first); -} - -bool compOrj(const int& p, const int& q){ - if(PDj[p].second != PDj[q].second) - return (PDj[p].second < PDj[q].second); - else - return (PDj[p].first > PDj[q].first); -} - -bool sortAngle(const std::pair >& p1, const std::pair >& p2){ - return p1.first < p2.first; -} - +bool sortAngle(const std::pair >& p1, const std::pair >& p2){return (p1.first < p2.first);} bool myComp(const std::pair & P1, const std::pair & P2){return P1.second < P2.second;} namespace Gudhi { namespace sliced_wasserstein { +// ******************************************************************** +// Approximate computation. +// ******************************************************************** + +/** \brief Computes an approximation of the Sliced Wasserstein distance between two persistence diagrams + * + * @param[in] N number of points sampled on the circle. + * + */ + double compute_approximate_SW(PD PD1, PD PD2, int N = 100){ double step = NUMPI/N; double sw = 0; // Add projections onto diagonal. - // ****************************** int n1, n2; n1 = PD1.size(); n2 = PD2.size(); for (int i = 0; i < n2; i++) - PD1.push_back(std::pair( (PD2[i].first+PD2[i].second)/2, (PD2[i].first+PD2[i].second)/2) ); + PD1.push_back(std::pair( (PD2[i].first + PD2[i].second)/2, (PD2[i].first + PD2[i].second)/2) ); for (int i = 0; i < n1; i++) - PD2.push_back(std::pair( (PD1[i].first+PD1[i].second)/2, (PD1[i].first+PD1[i].second)/2) ); + PD2.push_back(std::pair( (PD1[i].first + PD1[i].second)/2, (PD1[i].first + PD1[i].second)/2) ); int n = PD1.size(); // Sort and compare all projections. - // ********************************* //#pragma omp parallel for for (int i = 0; i < N; i++){ std::vector > L1, L2; @@ -110,8 +99,55 @@ double compute_approximate_SW(PD PD1, PD PD2, int N = 100){ return sw/NUMPI; } -double compute_int_cos(const double& alpha, const double& beta){ // Valid only if alpha is in [-pi,pi] and beta-alpha is in [0,pi] - double res; + + + + + + + + + + + + + + + + + + +// ******************************************************************** +// Exact computation. +// ******************************************************************** + + + +// Compute the angle formed by two points of a PD +double compute_angle(const PD & PersDiag, const int & i, const int & j){ + std::pair vect; double x1,y1, x2,y2; + x1 = PersDiag[i].first; y1 = PersDiag[i].second; + x2 = PersDiag[j].first; y2 = PersDiag[j].second; + if (y1 - y2 > 0){ + vect.first = y1 - y2; + vect.second = x2 - x1;} + else{ + if(y1 - y2 < 0){ + vect.first = y2 - y1; + vect.second = x1 - x2; + } + else{ + vect.first = 0; + vect.second = abs(x1 - x2);} + } + double norm = std::sqrt(pow(vect.first,2) + pow(vect.second,2)); + return asin(vect.second/norm); +} + +// 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 res = 0; assert((alpha >= 0 && alpha <= NUMPI) || (alpha >= -NUMPI && alpha <= 0)); if (alpha >= 0 && alpha <= NUMPI){ if (cos(alpha) >= 0){ @@ -136,75 +172,76 @@ double compute_int_cos(const double& alpha, const double& beta){ // Valid only i return res; } -double compute_int(const double& theta1, const double& theta2, const int& p, const int& q){ - double norm = std::sqrt(pow(PDi[p].first-PDj[q].first,2) + pow(PDi[p].second-PDj[q].second,2)); +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(pow(PD1[p].first-PD2[q].first,2) + pow(PD1[p].second-PD2[q].second,2)); double angle1; - if (PDi[p].first > PDj[q].first) - angle1 = theta1 - asin( (PDi[p].second-PDj[q].second)/norm ); + if (PD1[p].first > PD2[q].first) + angle1 = theta1 - asin( (PD1[p].second-PD2[q].second)/norm ); else - angle1 = theta1 - asin( (PDj[q].second-PDi[p].second)/norm ); - double angle2 = angle1+theta2-theta1; + angle1 = theta1 - asin( (PD2[q].second-PD1[p].second)/norm ); + double angle2 = angle1 + theta2 - theta1; double integral = compute_int_cos(angle1,angle2); return norm*integral; } -double compute_sw(const std::vector > >& V1, \ - const std::vector > >& V2){ + + +double compute_sw(const std::vector > > & V1, const std::vector > > & V2, const PD & PD1, const PD & PD2){ int N = V1.size(); double sw = 0; for (int i = 0; i < N; i++){ std::vector > U,V; U = V1[i]; V = V2[i]; double theta1, theta2; theta1 = -NUMPI/2; - int ku, kv; ku = 0; kv = 0; theta2 = std::min(U[ku].second,V[kv].second); + unsigned int ku, kv; ku = 0; kv = 0; theta2 = std::min(U[ku].second,V[kv].second); while(theta1 != NUMPI/2){ - if(PDi[U[ku].first].first != PDj[V[kv].first].first || PDi[U[ku].first].second != PDj[V[kv].first].second) + if(PD1[U[ku].first].first != PD2[V[kv].first].first || PD1[U[ku].first].second != PD2[V[kv].first].second) if(theta1 != theta2) - sw += compute_int(theta1,theta2,U[ku].first,V[kv].first); + sw += Gudhi::sliced_wasserstein::compute_int(theta1, theta2, U[ku].first, V[kv].first, PD1, PD2); theta1 = theta2; - if ( (theta2 == U[ku].second) && ku < U.size()-1 ){ku++;} - if ( (theta2 == V[kv].second) && kv < V.size()-1 ){kv++;} + 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); } } return sw/NUMPI; } -double compute_angle(const PD& PersDiag, const int& i, const int& j){ - std::pair vect; double x1,y1, x2,y2; - x1 = PersDiag[i].first; y1 = PersDiag[i].second; - x2 = PersDiag[j].first; y2 = PersDiag[j].second; - if (y1 - y2 > 0){ - vect.first = y1 - y2; - vect.second = x2 - x1;} - else{ - if(y1 - y2 < 0){ - vect.first = y2 - y1; - vect.second = x1 - x2; - } - else{ - vect.first = 0; - vect.second = abs(x1 - x2);} - } - double norm = std::sqrt(pow(vect.first,2) + pow(vect.second,2)); - return asin(vect.second/norm); -} + + + + + + + + + + + + + + + + + + +/** \brief Computes the Sliced Wasserstein distance between two persistence diagrams + * + */ double compute_exact_SW(PD PD1, PD PD2){ // Add projections onto diagonal. - // ****************************** - int n1, n2; n1 = PD1.size(); n2 = PD2.size(); double max_ordinate = std::numeric_limits::min(); + int n1, n2; n1 = PD1.size(); n2 = PD2.size(); double max_ordinate = std::numeric_limits::lowest(); for (int i = 0; i < n2; i++){ max_ordinate = std::max(max_ordinate, PD2[i].second); - PD1.push_back(std::pair( ((PD2[i].first+PD2[i].second)/2), ((PD2[i].first+PD2[i].second)/2)) ); + PD1.push_back( std::pair( ((PD2[i].first+PD2[i].second)/2), ((PD2[i].first+PD2[i].second)/2) ) ); } for (int i = 0; i < n1; i++){ max_ordinate = std::max(max_ordinate, PD1[i].second); - PD2.push_back(std::pair( ((PD1[i].first+PD1[i].second)/2), ((PD1[i].first+PD1[i].second)/2)) ); + PD2.push_back( std::pair( ((PD1[i].first+PD1[i].second)/2), ((PD1[i].first+PD1[i].second)/2) ) ); } int N = PD1.size(); assert(N==PD2.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)); @@ -214,43 +251,30 @@ double compute_exact_SW(PD PD1, PD PD2){ } // Compute all angles in both PDs. - // ******************************* std::vector > > angles1, angles2; for (int i = 0; i < N; i++){ for (int j = i+1; j < N; j++){ - double theta1 = compute_angle(PD1,i,j); double theta2 = compute_angle(PD2,i,j); + double theta1 = Gudhi::sliced_wasserstein::compute_angle(PD1,i,j); double theta2 = Gudhi::sliced_wasserstein::compute_angle(PD2,i,j); angles1.push_back(std::pair >(theta1, std::pair(i,j))); angles2.push_back(std::pair >(theta2, std::pair(i,j))); } } // Sort angles. - // ************ std::sort(angles1.begin(), angles1.end(), sortAngle); std::sort(angles2.begin(), angles2.end(), sortAngle); - // Initialize orders of the points of both PD (given by ordinates when theta = -pi/2). - // *********************************************************************************** - PDi = PD1; PDj = PD2; + // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2). std::vector orderp1, orderp2; - for (int i = 0; i < N; i++){orderp1.push_back(i); orderp2.push_back(i);} - std::sort(orderp1.begin(),orderp1.end(),compOri); std::sort(orderp2.begin(),orderp2.end(),compOrj); + for (int i = 0; i < N; i++){ orderp1.push_back(i); orderp2.push_back(i); } + std::sort( orderp1.begin(), orderp1.end(), [=](int i, int j){ if(PD1[i].second != PD1[j].second) return (PD1[i].second < PD1[j].second); else return (PD1[i].first > PD1[j].first); } ); + std::sort( orderp2.begin(), orderp2.end(), [=](int i, int j){ if(PD2[i].second != PD2[j].second) return (PD2[i].second < PD2[j].second); else return (PD2[i].first > PD2[j].first); } ); // Find the inverses of the orders. - // ******************************** std::vector order1(N); std::vector order2(N); - for(int i = 0; i < N; i++){ - for (int j = 0; j < N; j++) - if(orderp1[j] == i) - order1[i] = j; - } - for(int i = 0; i < N; i++){ - for (int j = 0; j < N; j++) - if(orderp2[j] == i) - order2[i] = j; - } + for(int i = 0; i < N; i++) for (int j = 0; j < N; j++) if(orderp1[j] == i){ order1[i] = j; break; } + for(int i = 0; i < N; i++) for (int j = 0; j < N; 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 > > anglePerm1(N); std::vector > > anglePerm2(N); @@ -276,11 +300,15 @@ double compute_exact_SW(PD PD1, PD PD2){ } // Compute the SW distance with the list of inversions. - // **************************************************** - return compute_sw(anglePerm1,anglePerm2); + return Gudhi::sliced_wasserstein::compute_sw(anglePerm1, anglePerm2, PD1, PD2); } + + + + + }} #endif diff --git a/src/Kernels/include/gudhi/figtree-0.9.3.zip b/src/Kernels/include/gudhi/figtree-0.9.3.zip new file mode 100644 index 00000000..a9468274 Binary files /dev/null and b/src/Kernels/include/gudhi/figtree-0.9.3.zip differ diff --git a/src/cmake/modules/GUDHI_modules.cmake b/src/cmake/modules/GUDHI_modules.cmake index f95d0c34..205ee8a1 100644 --- a/src/cmake/modules/GUDHI_modules.cmake +++ b/src/cmake/modules/GUDHI_modules.cmake @@ -16,8 +16,8 @@ function(add_gudhi_module file_path) endfunction(add_gudhi_module) -option(WITH_GUDHI_BENCHMARK "Activate/desactivate benchmark compilation" OFF) -option(WITH_GUDHI_EXAMPLE "Activate/desactivate examples compilation and installation" OFF) +option(WITH_GUDHI_BENCHMARK "Activate/desactivate benchmark compilation" ON) +option(WITH_GUDHI_EXAMPLE "Activate/desactivate examples compilation and installation" ON) option(WITH_GUDHI_PYTHON "Activate/desactivate python module compilation and installation" ON) option(WITH_GUDHI_TEST "Activate/desactivate examples compilation and installation" ON) option(WITH_GUDHI_UTILITIES "Activate/desactivate utilities compilation and installation" ON) -- cgit v1.2.3