From 230eac46395aeb406ca79280a8d62fc35b6f41a3 Mon Sep 17 00:00:00 2001 From: mcarrier Date: Tue, 30 Oct 2018 21:07:52 +0000 Subject: removed template Kernel2D git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3965 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 0cfdece6eb19cb8601371bc96e1dbb66ae41a086 --- .../example/persistence_heat_maps.cpp | 16 +-- .../include/gudhi/Persistence_heat_maps.h | 119 +++++++++++---------- .../include/gudhi/Sliced_Wasserstein.h | 7 +- src/Persistence_representations/test/kernels.cpp | 10 +- 4 files changed, 84 insertions(+), 68 deletions(-) (limited to 'src/Persistence_representations') diff --git a/src/Persistence_representations/example/persistence_heat_maps.cpp b/src/Persistence_representations/example/persistence_heat_maps.cpp index c177c009..0beb1052 100644 --- a/src/Persistence_representations/example/persistence_heat_maps.cpp +++ b/src/Persistence_representations/example/persistence_heat_maps.cpp @@ -27,9 +27,13 @@ #include #include -using Gaussian_kernel = Gudhi::Persistence_representations::Gaussian_kernel_2D; + +std::function, std::pair)> Gaussian_function(double sigma){ + return [=](std::pair p, std::pair q){return std::exp( -( (p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second) )/(sigma) );}; +} + using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; -using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; +using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; int main(int argc, char** argv) { // create two simple vectors with birth--death pairs: @@ -79,10 +83,10 @@ int main(int argc, char** argv) { std::cout << "Scalar product is : " << hm1.compute_scalar_product(hm2) << std::endl; // Mathieu's code ************************************************************************************************************ - Persistence_heat_maps hm1k(persistence1); - Persistence_heat_maps hm2k(persistence2); - Persistence_heat_maps hm1i(persistence1, 20, 20, 0, 11, 0, 11); - Persistence_heat_maps hm2i(persistence2, 20, 20, 0, 11, 0, 11); + Persistence_heat_maps hm1k(persistence1, Gaussian_function(1.0)); + Persistence_heat_maps hm2k(persistence2, Gaussian_function(1.0)); + Persistence_heat_maps hm1i(persistence1, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11); + Persistence_heat_maps hm2i(persistence2, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11); std::cout << "Scalar product computed with exact 2D kernel on grid is : " << hm1i.compute_scalar_product(hm2i) << std::endl; std::cout << "Scalar product computed with exact 2D kernel is : " << hm1k.compute_scalar_product(hm2k) << std::endl; // *************************************************************************************************************************** diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h index b2354432..db51cc14 100644 --- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h +++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h @@ -175,11 +175,11 @@ class weight_by_setting_maximal_interval_to_have_length_one { }; // Mathieu's code *********************************************************************************************************************** -class Gaussian_kernel_2D { - public: - double operator()(const std::pair& p, const std::pair& q) const { - return (1.0/(std::sqrt(2*pi))) * std::exp(-((p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second)) / 2); } -}; +//class Gaussian_kernel_2D { +// public: +// double operator()(const std::pair& p, const std::pair& q) const { +// return (1.0/(std::sqrt(2*pi))) * std::exp(-((p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second)) / 2); } +//}; // ************************************************************************************************************************************** /** @@ -191,7 +191,7 @@ class Gaussian_kernel_2D { // This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, // Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product -template +template class Persistence_heat_maps { public: /** @@ -220,7 +220,7 @@ class Persistence_heat_maps { *std::numeric_limits::max(), in which case the program compute the values based on the data. **/ Persistence_heat_maps(const std::vector >& interval, - std::vector > filter, //= create_Gaussian_filter(5, 1), Mathieu: I had to comment this default parameter otherwise there is an ambiguity with constructor defined line 263 + std::vector > filter = create_Gaussian_filter(5, 1), bool erase_below_diagonal = false, size_t number_of_pixels = 1000, double min_ = std::numeric_limits::max(), double max_ = std::numeric_limits::max()); @@ -256,13 +256,16 @@ class Persistence_heat_maps { /** * Construction that takes as inputs (1) the diagram, and (2) the grid parameters (min, max and number of samples for x and y axes) **/ - Persistence_heat_maps(const Persistence_diagram & interval, size_t number_of_x_pixels, size_t number_of_y_pixels, + Persistence_heat_maps(const std::vector >& interval, + const std::function, std::pair)>& kernel, + size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x = 0, double max_x = 1, double min_y = 0, double max_y = 1); /** * Construction that takes only the diagram as input (weight and 2D kernel are template parameters) **/ - Persistence_heat_maps(const Persistence_diagram & interval); + Persistence_heat_maps(const std::vector >& interval, + const std::function, std::pair)>& kernel); // ************************************************************************************************************************************** @@ -545,8 +548,8 @@ class Persistence_heat_maps { // Mathieu's code *********************************************************************************************************************** bool discrete = true; // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false) - Kernel2D k; - Persistence_diagram d; + std::function, std::pair)> kernel; + std::vector > interval; std::vector weights; // ************************************************************************************************************************************** @@ -561,34 +564,36 @@ class Persistence_heat_maps { // Mathieu's code *********************************************************************************************************************** -template -Persistence_heat_maps::Persistence_heat_maps(const Persistence_diagram & diagram, - size_t number_of_x_pixels, size_t number_of_y_pixels, - double min_x, double max_x, - double min_y, double max_y) { +template +Persistence_heat_maps::Persistence_heat_maps(const std::vector >& interval, + const std::function, std::pair)>& kernel, + size_t number_of_x_pixels, size_t number_of_y_pixels, + double min_x, double max_x, + double min_y, double max_y) { this->discrete = true; this->min_ = min_x; this->max_ = max_x; this->heat_map.resize(number_of_y_pixels); double step_x = (max_x - min_x)/(number_of_x_pixels - 1); double step_y = (max_y - min_y)/(number_of_y_pixels - 1); - int num_pts = diagram.size(); + int num_pts = interval.size(); for(size_t i = 0; i < number_of_y_pixels; i++){ double y = min_y + i*step_y; this->heat_map[i].reserve(number_of_x_pixels); for(size_t j = 0; j < number_of_x_pixels; j++){ double x = min_x + j*step_x; std::pair grid_point(x,y); double pixel_value = 0; - for(int k = 0; k < num_pts; k++) pixel_value += this->f(diagram[k]) * this->k(diagram[k], grid_point); + for(int k = 0; k < num_pts; k++) pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point); this->heat_map[i].push_back(pixel_value); } } this->set_up_parameters_for_basic_classes(); } -template -Persistence_heat_maps::Persistence_heat_maps(const Persistence_diagram& diagram) { - this->discrete = false; this->d = diagram; - int num_pts = this->d.size(); - for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->d[i])); +template +Persistence_heat_maps::Persistence_heat_maps(const std::vector >& interval, + const std::function, std::pair)>& kernel) { + this->discrete = false; this->interval = interval; this->kernel = kernel; + int num_pts = this->interval.size(); + for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->interval[i])); this->set_up_parameters_for_basic_classes(); } // ************************************************************************************************************************************** @@ -596,8 +601,8 @@ Persistence_heat_maps::Persistence_heat_maps(cons // if min_ == max_, then the program is requested to set up the values itself based on persistence intervals -template -void Persistence_heat_maps::construct(const std::vector >& intervals_, +template +void Persistence_heat_maps::construct(const std::vector >& intervals_, std::vector > filter, bool erase_below_diagonal, size_t number_of_pixels, double min_, double max_) { @@ -701,16 +706,16 @@ void Persistence_heat_maps::construct(const std:: } } // construct -template -Persistence_heat_maps::Persistence_heat_maps( +template +Persistence_heat_maps::Persistence_heat_maps( const std::vector >& interval, std::vector > filter, bool erase_below_diagonal, size_t number_of_pixels, double min_, double max_) { this->construct(interval, filter, erase_below_diagonal, number_of_pixels, min_, max_); this->set_up_parameters_for_basic_classes(); } -template -Persistence_heat_maps::Persistence_heat_maps(const char* filename, +template +Persistence_heat_maps::Persistence_heat_maps(const char* filename, std::vector > filter, bool erase_below_diagonal, size_t number_of_pixels, double min_, double max_, unsigned dimension) { @@ -724,8 +729,8 @@ Persistence_heat_maps::Persistence_heat_maps(cons this->set_up_parameters_for_basic_classes(); } -template -std::vector > Persistence_heat_maps::check_and_initialize_maps( +template +std::vector > Persistence_heat_maps::check_and_initialize_maps( const std::vector& maps) { // checking if all the heat maps are of the same size: for (size_t i = 0; i != maps.size(); ++i) { @@ -746,8 +751,8 @@ std::vector > Persistence_heat_maps -void Persistence_heat_maps::compute_median(const std::vector& maps) { +template +void Persistence_heat_maps::compute_median(const std::vector& maps) { std::vector > heat_maps = this->check_and_initialize_maps(maps); std::vector to_compute_median(maps.size()); @@ -766,8 +771,8 @@ void Persistence_heat_maps::compute_median(const this->max_ = maps[0]->max_; } -template -void Persistence_heat_maps::compute_mean(const std::vector& maps) { +template +void Persistence_heat_maps::compute_mean(const std::vector& maps) { std::vector > heat_maps = this->check_and_initialize_maps(maps); for (size_t i = 0; i != heat_maps.size(); ++i) { for (size_t j = 0; j != heat_maps[i].size(); ++j) { @@ -783,8 +788,8 @@ void Persistence_heat_maps::compute_mean(const st this->max_ = maps[0]->max_; } -template -void Persistence_heat_maps::compute_percentage_of_active( +template +void Persistence_heat_maps::compute_percentage_of_active( const std::vector& maps, size_t cutoff) { std::vector > heat_maps = this->check_and_initialize_maps(maps); @@ -806,8 +811,8 @@ void Persistence_heat_maps::compute_percentage_of this->max_ = maps[0]->max_; } -template -void Persistence_heat_maps::plot(const char* filename) const { +template +void Persistence_heat_maps::plot(const char* filename) const { std::ofstream out; std::stringstream gnuplot_script; gnuplot_script << filename << "_GnuplotScript"; @@ -825,8 +830,8 @@ void Persistence_heat_maps::plot(const char* file << gnuplot_script.str().c_str() << "\'\"" << std::endl; } -template -void Persistence_heat_maps::print_to_file(const char* filename) const { +template +void Persistence_heat_maps::print_to_file(const char* filename) const { std::ofstream out; out.open(filename); @@ -841,8 +846,8 @@ void Persistence_heat_maps::print_to_file(const c out.close(); } -template -void Persistence_heat_maps::load_from_file(const char* filename) { +template +void Persistence_heat_maps::load_from_file(const char* filename) { bool dbg = false; std::ifstream in; @@ -890,8 +895,8 @@ void Persistence_heat_maps::load_from_file(const } // Concretizations of virtual methods: -template -std::vector Persistence_heat_maps::vectorize(int number_of_function) const { +template +std::vector Persistence_heat_maps::vectorize(int number_of_function) const { std::vector result; if(!discrete){std::cout << "No vectorize method in case of infinite dimensional vectorization" << std::endl; return result;} @@ -913,8 +918,8 @@ std::vector Persistence_heat_maps::vector return result; } -template -double Persistence_heat_maps::distance(const Persistence_heat_maps& second, double power) const { +template +double Persistence_heat_maps::distance(const Persistence_heat_maps& second, double power) const { if(this->discrete){ // first we need to check if (*this) and second are defined on the same domain and have the same dimensions: if (!this->check_if_the_same(second)) { @@ -950,8 +955,8 @@ double Persistence_heat_maps::distance(const Pers } } -template -double Persistence_heat_maps::project_to_R(int number_of_function) const { +template +double Persistence_heat_maps::project_to_R(int number_of_function) const { double result = 0; for (size_t i = 0; i != this->heat_map.size(); ++i) { for (size_t j = 0; j != this->heat_map[i].size(); ++j) { @@ -961,14 +966,14 @@ double Persistence_heat_maps::project_to_R(int nu return result; } -template -void Persistence_heat_maps::compute_average( +template +void Persistence_heat_maps::compute_average( const std::vector& to_average) { this->compute_mean(to_average); } -template -double Persistence_heat_maps::compute_scalar_product(const Persistence_heat_maps& second) const { +template +double Persistence_heat_maps::compute_scalar_product(const Persistence_heat_maps& second) const { if(discrete){ // first we need to check if (*this) and second are defined on the same domain and have the same dimensions: @@ -991,12 +996,12 @@ double Persistence_heat_maps::compute_scalar_prod // Mathieu's code *********************************************************************************************************************** else{ - int num_pts1 = this->d.size(); int num_pts2 = second.d.size(); double kernel_val = 0; + int num_pts1 = this->interval.size(); int num_pts2 = second.interval.size(); double kernel_val = 0; for(int i = 0; i < num_pts1; i++){ - std::pair pi = this->d[i]; + std::pair pi = this->interval[i]; for(int j = 0; j < num_pts2; j++){ - std::pair pj = second.d[j]; - kernel_val += this->weights[i] * second.weights[j] * this->k(pi, pj); + std::pair pj = second.interval[j]; + kernel_val += this->weights[i] * second.weights[j] * this->kernel(pi, pj); } } return kernel_val; diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h index 79fe8690..a0191dd7 100644 --- a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -163,10 +163,11 @@ class Sliced_Wasserstein { // Slightly perturb the points so that the PDs are in generic positions. double thresh_y = max_ordinate * 0.00001; double thresh_x = max_abscissa * 0.00001; - srand(time(NULL)); + std::random_device rd; std::default_random_engine re(rd()); std::uniform_real_distribution uni(-1,1); + double epsilon = uni(re); for (int i = 0; i < num_pts_dgm; i++){ - diagram1[i].first += thresh_x*(1.0-2.0*rand()/RAND_MAX); diagram1[i].second += thresh_y*(1.0-2.0*rand()/RAND_MAX); - diagram2[i].first += thresh_x*(1.0-2.0*rand()/RAND_MAX); diagram2[i].second += thresh_y*(1.0-2.0*rand()/RAND_MAX); + diagram1[i].first += thresh_x*epsilon; diagram1[i].second += thresh_y*epsilon; + diagram2[i].first += thresh_x*epsilon; diagram2[i].second += thresh_y*epsilon; } // Compute all angles in both PDs. diff --git a/src/Persistence_representations/test/kernels.cpp b/src/Persistence_representations/test/kernels.cpp index c95e8086..b8d02d4c 100644 --- a/src/Persistence_representations/test/kernels.cpp +++ b/src/Persistence_representations/test/kernels.cpp @@ -40,10 +40,16 @@ using SW = Gudhi::Persistence_representations::Sliced_Wasserstein; using PWG = Gudhi::Persistence_representations::Persistence_heat_maps; using Persistence_diagram = std::vector >; +std::function, std::pair)> Gaussian_function(double sigma){ + return [=](std::pair p, std::pair q){ + return (1/std::sqrt(2*Gudhi::Persistence_representations::pi)*sigma) * std::exp( -( (p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second) )/(2*sigma) ); + }; +} + BOOST_AUTO_TEST_CASE(check_PWG) { Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); - PWG pwg1(v1, Gudhi::Persistence_representations::Gaussian_kernel(1.0)); - PWG pwg2(v2, Gudhi::Persistence_representations::Gaussian_kernel(1.0)); + PWG pwg1(v1, Gaussian_function(1.0)); + PWG pwg2(v2, Gaussian_function(1.0)); BOOST_CHECK(std::abs(pwg1.compute_scalar_product(pwg2) - std::exp(-0.5)/(std::sqrt(2*Gudhi::Persistence_representations::pi))) <= 1e-3); } -- cgit v1.2.3