From be131d6f74a9264e15a0b1c1e72fa8967c4518bd Mon Sep 17 00:00:00 2001 From: mcarrier Date: Thu, 11 Jan 2018 16:31:30 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3129 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 567185cf3691252c2685c450d65787d824442bb2 --- src/Kernels/example/kernel_basic_example.cpp | 5 ++ src/Kernels/include/gudhi/kernel.h | 95 +++++++++++++--------------- 2 files changed, 49 insertions(+), 51 deletions(-) (limited to 'src') diff --git a/src/Kernels/example/kernel_basic_example.cpp b/src/Kernels/example/kernel_basic_example.cpp index 46e42c9d..85ce36d4 100644 --- a/src/Kernels/example/kernel_basic_example.cpp +++ b/src/Kernels/example/kernel_basic_example.cpp @@ -21,6 +21,11 @@ */ #include +#include +#include +#include +#include + void usage(int nbArgs, char *const progName) { std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; diff --git a/src/Kernels/include/gudhi/kernel.h b/src/Kernels/include/gudhi/kernel.h index 30864476..900db092 100644 --- a/src/Kernels/include/gudhi/kernel.h +++ b/src/Kernels/include/gudhi/kernel.h @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carrière * - * Copyright (C) 2017 INRIA (France) + * Copyright (C) 2018 INRIA (France) * * 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 @@ -23,45 +23,32 @@ #ifndef KERNEL_H_ #define KERNEL_H_ -#include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include #include #include -#include -#include #include -#include -#include -#include -#include -//#include -//#include #include -#include -#include +#include -using PD = std::vector >; -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 kernel { +using PD = std::vector >; +double pi = boost::math::constants::pi(); + +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;} double pss_weight(std::pair P){ if(P.second > P.first) return 1; else return -1; } +double arctan_weight(std::pair P){ + return atan(P.second - P.first); +} + @@ -78,7 +65,8 @@ double pss_weight(std::pair P){ * @param[in] weight weight function for the points in the diagrams. * */ -double lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ +template) > +double lpwgk(const PD & PD1, const PD & PD2, double sigma, Weight weight = arctan_weight){ int num_pts1 = PD1.size(); int num_pts2 = PD2.size(); double k = 0; for(int i = 0; i < num_pts1; i++) for(int j = 0; j < num_pts2; j++) @@ -94,10 +82,10 @@ double lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair(PD1[i].second,PD1[i].first)); PD pd2 = PD2; numpts = PD2.size(); for(int i = 0; i < numpts; i++) pd2.push_back(std::pair(PD2[i].second,PD2[i].first)); - return lpwgk(pd1, pd2, 2*sqrt(sigma), &pss_weight) / (2*8*3.14159265359*sigma); + return lpwgk(pd1, pd2, 2*sqrt(sigma), &pss_weight) / (2*8*pi*sigma); } /** \brief Computes the Gaussian Persistence Weighted Gaussian Kernel between two persistence diagrams. @@ -110,7 +98,8 @@ double pssk(PD PD1, PD PD2, double sigma){ * @param[in] weight weight function for the points in the diagrams. * */ -double gpwgk(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ +template) > +double gpwgk(const PD & PD1, const PD & PD2, double sigma, double tau, Weight weight = arctan_weight){ double k1 = lpwgk(PD1,PD1,sigma,weight); double k2 = lpwgk(PD2,PD2,sigma,weight); double k3 = lpwgk(PD1,PD2,sigma,weight); @@ -126,7 +115,8 @@ double gpwgk(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pai * @param[in] weight weight function for the points in the diagrams. * */ -double dpwg(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}){ +template) > +double dpwg(const PD & PD1, const PD & PD2, double sigma, Weight weight = arctan_weight){ double k1 = lpwgk(PD1,PD1,sigma,weight); double k2 = lpwgk(PD2,PD2,sigma,weight); double k3 = lpwgk(PD1,PD2,sigma,weight); @@ -158,24 +148,24 @@ double compute_angle(const PD & PersDiag, const int & i, const int & j){ // 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 <= 3.14159265359) || (alpha >= -3.14159265359 && alpha <= 0)); - if (alpha >= 0 && alpha <= 3.14159265359){ + //assert((alpha >= 0 && alpha <= pi) || (alpha >= -pi && alpha <= 0)); + if (alpha >= 0 && alpha <= pi){ if (cos(alpha) >= 0){ - if(3.14159265359/2 <= beta){res = 2-sin(alpha)-sin(beta);} + if(pi/2 <= beta){res = 2-sin(alpha)-sin(beta);} else{res = sin(beta)-sin(alpha);} } else{ - if(1.5*3.14159265359 <= beta){res = 2+sin(alpha)+sin(beta);} + if(1.5*pi <= beta){res = 2+sin(alpha)+sin(beta);} else{res = sin(alpha)-sin(beta);} } } - if (alpha >= -3.14159265359 && alpha <= 0){ + if (alpha >= -pi && alpha <= 0){ if (cos(alpha) <= 0){ - if(-3.14159265359/2 <= beta){res = 2+sin(alpha)+sin(beta);} + if(-pi/2 <= beta){res = 2+sin(alpha)+sin(beta);} else{res = sin(alpha)-sin(beta);} } else{ - if(3.14159265359/2 <= beta){res = 2-sin(alpha)-sin(beta);} + if(pi/2 <= beta){res = 2-sin(alpha)-sin(beta);} else{res = sin(beta)-sin(alpha);} } } @@ -200,9 +190,9 @@ double compute_sw(const std::vector > > & V1, 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 = -3.14159265359/2; + 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 != 3.14159265359/2){ + while(theta1 != pi/2){ 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, PD1, PD2); @@ -212,7 +202,7 @@ double compute_sw(const std::vector > > & V1, theta2 = std::min(U[ku].second, V[kv].second); } } - return sw/3.14159265359; + return sw/pi; } /** \brief Computes the Sliced Wasserstein distance between two persistence diagrams. @@ -234,7 +224,7 @@ double compute_sw(const std::vector > > & V1, 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) ) ); } - int N = PD1.size(); assert(N==PD2.size()); + int N = PD1.size(); // Slightly perturb the points so that the PDs are in generic positions. int mag = 0; while(max_ordinate > 10){mag++; max_ordinate/=10;} @@ -290,8 +280,8 @@ double compute_sw(const std::vector > > & V1, } for (int i = 0; i < N; i++){ - anglePerm1[order1[i]].push_back(std::pair(i,3.14159265359/2)); - anglePerm2[order2[i]].push_back(std::pair(i,3.14159265359/2)); + anglePerm1[order1[i]].push_back(std::pair(i,pi/2)); + anglePerm2[order2[i]].push_back(std::pair(i,pi/2)); } // Compute the SW distance with the list of inversions. @@ -344,7 +334,8 @@ double approx_dpwg_Fourier(const std::vector >& B1, con return std::sqrt((1.0/M)*(d1+d2)-2*d3); } -std::vector > Fourier_feat(PD D, std::vector > Z, double (*weight)(std::pair) = [](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; @@ -379,7 +370,8 @@ std::vector > random_Fourier(double sigma, int M = 1000 * @param[in] M number of Fourier features. * */ -double approx_lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}, int M = 1000){ +template) > +double approx_lpwgk(const PD & PD1, const PD & PD2, double sigma, Weight weight = arctan_weight, int M = 1000){ std::vector > Z = random_Fourier(sigma, M); std::vector > B1 = Fourier_feat(PD1,Z,weight); std::vector > B2 = Fourier_feat(PD2,Z,weight); @@ -395,10 +387,10 @@ double approx_lpwgk(PD PD1, PD PD2, double sigma, double (*weight)(std::pair(PD1[i].second,PD1[i].first)); PD pd2 = PD2; numpts = PD2.size(); for(int i = 0; i < numpts; i++) pd2.push_back(std::pair(PD2[i].second,PD2[i].first)); - return approx_lpwgk(pd1, pd2, 2*sqrt(sigma), &pss_weight, M) / (2*8*3.14159265359*sigma); + return approx_lpwgk(pd1, pd2, 2*sqrt(sigma), &pss_weight, M) / (2*8*pi*sigma); } @@ -413,7 +405,8 @@ double approx_pssk(PD PD1, PD PD2, double sigma, int M = 1000){ * @param[in] M number of Fourier features. * */ -double approx_gpwgk(PD PD1, PD PD2, double sigma, double tau, double (*weight)(std::pair) = [](std::pair P){return atan(P.second - P.first);}, int M = 1000){ +template) > +double approx_gpwgk(const PD & PD1, const PD & PD2, double sigma, double tau, Weight weight = arctan_weight, int M = 1000){ std::vector > Z = random_Fourier(sigma, M); std::vector > B1 = Fourier_feat(PD1,Z,weight); std::vector > B2 = Fourier_feat(PD2,Z,weight); @@ -431,7 +424,7 @@ double approx_gpwgk(PD PD1, PD PD2, double sigma, double tau, double (*weight)(s */ double approx_sw(PD PD1, PD PD2, int N = 100){ - double step = 3.14159265359/N; double sw = 0; + double step = pi/N; double sw = 0; // Add projections onto diagonal. int n1, n2; n1 = PD1.size(); n2 = PD2.size(); @@ -446,14 +439,14 @@ double approx_sw(PD PD1, PD PD2, int N = 100){ for (int i = 0; i < N; i++){ std::vector > L1, L2; for (int j = 0; j < n; j++){ - L1.push_back( std::pair(j, PD1[j].first*cos(-3.14159265359/2+i*step) + PD1[j].second*sin(-3.14159265359/2+i*step)) ); - L2.push_back( std::pair(j, PD2[j].first*cos(-3.14159265359/2+i*step) + PD2[j].second*sin(-3.14159265359/2+i*step)) ); + L1.push_back( std::pair(j, PD1[j].first*cos(-pi/2+i*step) + PD1[j].second*sin(-pi/2+i*step)) ); + L2.push_back( std::pair(j, PD2[j].first*cos(-pi/2+i*step) + PD2[j].second*sin(-pi/2+i*step)) ); } std::sort(L1.begin(),L1.end(), myComp); std::sort(L2.begin(),L2.end(), myComp); double f = 0; for (int j = 0; j < n; j++) f += std::abs(L1[j].second - L2[j].second); sw += f*step; } - return sw/3.14159265359; + return sw/pi; } /** \brief Computes an approximation of the Sliced Wasserstein Kernel between two persistence diagrams. -- cgit v1.2.3