From 93654d5708654a6071c1775580f625da625a08a8 Mon Sep 17 00:00:00 2001 From: skachano Date: Wed, 5 Oct 2016 15:02:50 +0000 Subject: Junk out git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1646 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: ce2dfd7033a6a8366f9861ced4ed64ac41dfb82e --- .../example/Landmark_choice_random_knn.h | 142 ----- .../example/Landmark_choice_sparsification.h | 230 ------- src/Witness_complex/example/a0_persistence.cpp | 638 ------------------ src/Witness_complex/example/bench_rwit.cpp | 709 --------------------- src/Witness_complex/example/color-picker.cpp | 37 -- src/Witness_complex/example/off_writer.h | 11 - src/Witness_complex/example/relaxed_delaunay.cpp | 180 ------ .../example/relaxed_witness_complex_sphere.cpp | 461 -------------- .../example/relaxed_witness_persistence.cpp | 267 -------- src/Witness_complex/example/strange_sphere.cpp | 219 ------- .../example/witness_complex_from_file.cpp | 36 +- src/Witness_complex/include/gudhi/A0_complex.h | 337 ---------- .../gudhi/Construct_closest_landmark_table.h | 90 --- .../include/gudhi/Dim_list_iterator.h | 155 ----- src/Witness_complex/include/gudhi/Dim_lists.h | 195 ------ src/Witness_complex/include/gudhi/Good_links.h | 449 ------------- .../gudhi/Landmark_choice_by_furthest_point.h | 105 --- .../gudhi/Landmark_choice_by_random_point.h | 96 --- .../include/gudhi/Relaxed_witness_complex.h | 379 ----------- .../include/gudhi/Strange_witness_complex.h | 232 ------- .../include/gudhi/Strong_relaxed_witness_complex.h | 337 ---------- .../include/gudhi/Weak_relaxed_witness_complex.h | 379 ----------- .../include/gudhi/Witness_complex.h | 7 +- 23 files changed, 21 insertions(+), 5670 deletions(-) delete mode 100644 src/Witness_complex/example/Landmark_choice_random_knn.h delete mode 100644 src/Witness_complex/example/Landmark_choice_sparsification.h delete mode 100644 src/Witness_complex/example/a0_persistence.cpp delete mode 100644 src/Witness_complex/example/bench_rwit.cpp delete mode 100644 src/Witness_complex/example/color-picker.cpp delete mode 100644 src/Witness_complex/example/off_writer.h delete mode 100644 src/Witness_complex/example/relaxed_delaunay.cpp delete mode 100644 src/Witness_complex/example/relaxed_witness_complex_sphere.cpp delete mode 100644 src/Witness_complex/example/relaxed_witness_persistence.cpp delete mode 100644 src/Witness_complex/example/strange_sphere.cpp delete mode 100644 src/Witness_complex/include/gudhi/A0_complex.h delete mode 100644 src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h delete mode 100644 src/Witness_complex/include/gudhi/Dim_list_iterator.h delete mode 100644 src/Witness_complex/include/gudhi/Dim_lists.h delete mode 100644 src/Witness_complex/include/gudhi/Good_links.h delete mode 100644 src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h delete mode 100644 src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h delete mode 100644 src/Witness_complex/include/gudhi/Relaxed_witness_complex.h delete mode 100644 src/Witness_complex/include/gudhi/Strange_witness_complex.h delete mode 100644 src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h delete mode 100644 src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h (limited to 'src/Witness_complex') diff --git a/src/Witness_complex/example/Landmark_choice_random_knn.h b/src/Witness_complex/example/Landmark_choice_random_knn.h deleted file mode 100644 index fb91e116..00000000 --- a/src/Witness_complex/example/Landmark_choice_random_knn.h +++ /dev/null @@ -1,142 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef LANDMARK_CHOICE_BY_RANDOM_KNN_H_ -#define LANDMARK_CHOICE_BY_RANDOM_KNN_H_ - -#include // for pair<> -#include -#include // for ptrdiff_t type - -//#include -#include -#include -//#include -#include -#include -#include -#include -#include -//#include -#include - - -namespace Gudhi { - -namespace witness_complex { - -typedef CGAL::Epick_d K; -typedef K::FT FT; -typedef K::Point_d Point_d; -typedef CGAL::Search_traits< FT, - Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d > Traits_base; -typedef CGAL::Euclidean_distance Euclidean_distance; -typedef CGAL::Search_traits_adapter< std::ptrdiff_t, - Point_d*, - Traits_base> STraits; -typedef CGAL::Distance_adapter< std::ptrdiff_t, - Point_d*, - Euclidean_distance > Distance_adapter; -typedef CGAL::Orthogonal_incremental_neighbor_search< STraits, - Distance_adapter > Neighbor_search; - typedef CGAL::Orthogonal_k_neighbor_search< STraits > Neighbor_search2; -typedef Neighbor_search::Tree Tree; - - - /** \brief Landmark choice strategy by taking random vertices for landmarks. - * \details It chooses nbL distinct landmarks from a random access range `points` - * and outputs a matrix {witness}*{closest landmarks} in knn. - */ - template - void landmark_choice_by_random_knn(Point_random_access_range const & points, - int nbL, - FT alpha, - unsigned limD, - KNearestNeighbours & knn, - Distance_matrix & distances) { - int nbP = points.end() - points.begin(); - assert(nbP >= nbL); - std::vector landmarks; - std::vector landmarks_ind; - Point_d p; - int chosen_landmark; - CGAL::Random rand; - // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety - int current_number_of_landmarks = 0; // counter for landmarks - for (; current_number_of_landmarks != nbL; current_number_of_landmarks++) { - do chosen_landmark = rand.get_int(0,nbP); - while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end()); - p = points[chosen_landmark]; - landmarks.push_back(p); - landmarks_ind.push_back(chosen_landmark); - } - // std::cout << "Choice finished!" << std::endl; - - //int dim = points.begin()->size(); - knn = KNearestNeighbours(nbP); - distances = Distance_matrix(nbP); - STraits traits(&(landmarks[0])); - CGAL::Distance_adapter adapter(&(landmarks[0])); - Euclidean_distance ed; - Tree landmark_tree(boost::counting_iterator(0), - boost::counting_iterator(nbL), - typename Tree::Splitter(), - traits); - for (int points_i = 0; points_i < nbP; points_i++) { - Point_d const & w = points[points_i]; - Neighbor_search search(landmark_tree, - w, - FT(0), - true, - adapter); - Neighbor_search::iterator search_it = search.begin(); - // Neighbor_search2 search(landmark_tree, - // w, limD+1, - // FT(0), - // true, - // adapter); - // Neighbor_search2::iterator search_it = search.begin(); - - while (knn[points_i].size() < limD) { - distances[points_i].push_back(sqrt(search_it->second)); - knn[points_i].push_back((search_it++)->first); - } - FT dtow = distances[points_i][limD-1]; - - if (alpha != 0) - while (search_it != search.end() && search_it->second < dtow + alpha) { - distances[points_i].push_back(sqrt(search_it->second)); - knn[points_i].push_back((search_it++)->first); - } - //std::cout << "k = " << knn[points_i].size() << std::endl; - } - } - -} // namespace witness_complex - -} // namespace Gudhi - -#endif // LANDMARK_CHOICE_BY_RANDOM_POINT_H_ diff --git a/src/Witness_complex/example/Landmark_choice_sparsification.h b/src/Witness_complex/example/Landmark_choice_sparsification.h deleted file mode 100644 index d5d5fba1..00000000 --- a/src/Witness_complex/example/Landmark_choice_sparsification.h +++ /dev/null @@ -1,230 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef LANDMARK_CHOICE_BY_SPARSIFICATION_H_ -#define LANDMARK_CHOICE_BY_SPARSIFICATION_H_ - -#include // for pair<> -#include -#include // for ptrdiff_t type -#include - -//#include -#include -#include -//#include -#include -#include -#include -#include -#include -//#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - typedef CGAL::Epick_d K; - typedef K::FT FT; - typedef K::Point_d Point_d; - typedef CGAL::Search_traits< FT, - Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d > Traits_base; - typedef CGAL::Euclidean_distance Euclidean_distance; - typedef CGAL::Search_traits_adapter< std::ptrdiff_t, - Point_d*, - Traits_base> STraits; - typedef CGAL::Distance_adapter< std::ptrdiff_t, - Point_d*, - Euclidean_distance > Distance_adapter; - typedef CGAL::Orthogonal_incremental_neighbor_search< STraits, - Distance_adapter > Neighbor_search; - typedef CGAL::Orthogonal_k_neighbor_search< STraits > Neighbor_search2; - typedef Neighbor_search::Tree Tree; - typedef CGAL::Fuzzy_sphere Fuzzy_sphere; - - /** \brief Landmark selection function by as a sub-epsilon-net of the - * given set of points. - */ - template - void landmark_choice_by_sparsification(Point_random_access_range & points, - unsigned nbL, - FT mu_epsilon, - Point_random_access_range & landmarks) - { - unsigned nbP = points.end() - points.begin(); - assert(nbP >= nbL); - CGAL::Random rand; - // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety - STraits points_traits(&(points[0])); - CGAL::Distance_adapter points_adapter(&(points[0])); - std::vector dropped_points(nbP, false); - - Tree witness_tree(boost::counting_iterator(0), - boost::counting_iterator(nbP), - typename Tree::Splitter(), - points_traits); - - for (unsigned points_i = 0; points_i < nbP; points_i++) { - if (dropped_points[points_i]) - continue; - Point_d & w = points[points_i]; - Fuzzy_sphere fs(w, mu_epsilon, 0, points_traits); - std::vector close_neighbors; - witness_tree.search(std::insert_iterator>(close_neighbors,close_neighbors.begin()),fs); - for (int i: close_neighbors) - dropped_points[i] = true; - } - - for (unsigned points_i = 0; points_i < nbP; points_i++) { - if (dropped_points[points_i]) - landmarks.push_back(points[points_i]); - } - - if (nbL < landmarks.size()) { - std::random_shuffle(landmarks.begin(), landmarks.end()); - landmarks.resize(nbL); - } - } - - - - - /** \brief Landmark choice strategy by taking random vertices for landmarks. - * \details It chooses nbL distinct landmarks from a random access range `points` - * and outputs a matrix {witness}*{closest landmarks} in knn. - */ - template - void build_distance_matrix(Point_random_access_range const & points, - Point_random_access_range & landmarks, - FT alpha, - unsigned limD, - KNearestNeighbours & knn, - Distance_matrix & distances) - { - int nbP = points.end() - points.begin(); - knn = KNearestNeighbours(nbP); - distances = Distance_matrix(nbP); - STraits traits(&(landmarks[0])); - CGAL::Distance_adapter adapter(&(landmarks[0])); - Euclidean_distance ed; - Tree landmark_tree(boost::counting_iterator(0), - boost::counting_iterator(landmarks.size()), - typename Tree::Splitter(), - traits); - for (int points_i = 0; points_i < nbP; points_i++) { - Point_d const & w = points[points_i]; - Neighbor_search search(landmark_tree, - w, - FT(0), - true, - adapter); - Neighbor_search::iterator search_it = search.begin(); - // Neighbor_search2 search(landmark_tree, - // w, limD+1, - // FT(0), - // true, - // adapter); - // Neighbor_search2::iterator search_it = search.begin(); - - while (knn[points_i].size() <= limD) { - distances[points_i].push_back(search_it->second); //!sq_dist - knn[points_i].push_back((search_it++)->first); - } - FT dtow = distances[points_i][limD]; - - while (search_it != search.end() && search_it->second < dtow + alpha) { - distances[points_i].push_back(search_it->second); - knn[points_i].push_back((search_it++)->first); - } - //std::cout << "k = " << knn[points_i].size() << std::endl; - } - } - - /* - template - std::vector - sparsify_point_set(const Kernel &k, - Point_container const& input_pts, - typename Kernel::FT min_squared_dist) - { - typedef typename CGAL::Tangential_complex_::Point_cloud_data_structure Points_ds; - typedef typename Points_ds::INS_iterator INS_iterator; - typedef typename Points_ds::INS_range INS_range; - - typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object(); - - // Create the output container - std::vector output; - - Points_ds points_ds(input_pts); - - std::vector dropped_points(input_pts.size(), false); - - // Parse the input points, and add them if they are not too close to - // the other points - std::size_t pt_idx = 0; - for (typename Point_container::const_iterator it_pt = input_pts.begin() ; - it_pt != input_pts.end(); - ++it_pt, ++pt_idx) - { - if (dropped_points[pt_idx]) - continue; - - output.push_back(*it_pt); - - INS_range ins_range = points_ds.query_incremental_ANN(*it_pt); - - // If another point Q is closer that min_squared_dist, mark Q to be dropped - for (INS_iterator nn_it = ins_range.begin() ; - nn_it != ins_range.end() ; - ++nn_it) - { - std::size_t neighbor_point_idx = nn_it->first; - // If the neighbor is too close, we drop the neighbor - if (nn_it->second < min_squared_dist) - { - // N.B.: If neighbor_point_idx < pt_idx, - // dropped_points[neighbor_point_idx] is already true but adding a - // test doesn't make things faster, so why bother? - dropped_points[neighbor_point_idx] = true; - } - else - break; - } - } - - return output; -} - */ - - -} // namespace witness_complex - -} // namespace Gudhi - -#endif // LANDMARK_CHOICE_BY_RANDOM_POINT_H_ diff --git a/src/Witness_complex/example/a0_persistence.cpp b/src/Witness_complex/example/a0_persistence.cpp deleted file mode 100644 index 422185ca..00000000 --- a/src/Witness_complex/example/a0_persistence.cpp +++ /dev/null @@ -1,638 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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 - * 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 -#include "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include "generators.h" -#include "output.h" -#include "output_tikz.h" - -using namespace Gudhi; -using namespace Gudhi::witness_complex; -using namespace Gudhi::persistent_cohomology; - -typedef std::vector Point_Vector; -typedef A0_complex< Simplex_tree<> > A0Complex; -typedef Simplex_tree<>::Simplex_handle Simplex_handle; - -typedef A0_complex< Simplex_tree<> > SRWit; -typedef Relaxed_witness_complex< Simplex_tree<> > WRWit; - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { - std::ifstream in_file(file_name.c_str(), std::ios::in); - if (!in_file.is_open()) { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while (getline(in_file, line)) { - std::vector< double > point; - std::istringstream iss(line); - while (iss >> x) { - point.push_back(x); - } - if (point.size() != 1) - points.push_back(point); - } - in_file.close(); -} - -void output_experiment_information(char * const file_name) -{ - std::cout << "Enter a valid experiment number. Usage: " - << file_name << " exp_no options\n"; - std::cout << "Experiment description:\n" - << "0 nbP nbL dim alpha limD mu_epsilon: " - << "Build persistence diagram on relaxed witness complex " - << "built from a point cloud on (dim-1)-dimensional sphere " - << "consisting of nbP witnesses and nbL landmarks. " - << "The maximal relaxation is alpha and the limit on simplicial complex " - << "dimension is limD.\n"; - std::cout << "1 file_name nbL alpha limD: " - << "Build persistence diagram on relaxed witness complex " - << "build from a point cloud stored in a file and nbL landmarks. " - << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n"; -} - -void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT mu_epsilon = 0.1) -{ - clock_t start, end; - Simplex_tree<> simplex_tree; - - // Choose landmarks - std::vector > knn; - std::vector > distances; - start = clock(); - //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); - - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - limD, - knn, - distances); - end = clock(); - double time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Choice of " << nbL << " landmarks took " - << time << " s. \n"; - // Compute witness complex - start = clock(); - A0Complex rw(distances, - knn, - simplex_tree, - nbL, - alpha2, - limD); - end = clock(); - time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Witness complex for " << nbL << " landmarks took " - << time << " s. \n"; - std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; - int streedim = 0; - for (auto s: simplex_tree.complex_simplex_range()) - if (simplex_tree.dimension(s) > streedim) - streedim = simplex_tree.dimension(s); - std::cout << "Dimension of the simplicial complex is " << streedim << std::endl; - - //std::cout << simplex_tree << "\n"; - - // Compute the persistence diagram of the complex - simplex_tree.set_dimension(limD); - persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, true); - int p = 3; - pcoh.init_coefficients( p ); //initilizes the coefficient field for homology - start = clock(); - pcoh.compute_persistent_cohomology( alpha2/10 ); - end = clock(); - time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Persistence diagram took " - << time << " s. \n"; - pcoh.output_diagram(); - - int chi = 0; - for (auto sh: simplex_tree.complex_simplex_range()) - chi += 1-2*(simplex_tree.dimension(sh)%2); - std::cout << "Euler characteristic is " << chi << std::endl; - - // Gudhi::witness_complex::Dim_lists> simplices(simplex_tree, limD); - - // simplices.collapse(); - // simplices.output_simplices(); - - // Simplex_tree<> collapsed_tree; - // for (auto sh: simplices) { - // std::vector vertices; - // for (int v: collapsed_tree.simplex_vertex_range(sh)) - // vertices.push_back(v); - // collapsed_tree.insert_simplex(vertices); - // } - std::vector landmarks_ind(nbL); - for (unsigned i = 0; i != distances.size(); ++i) { - if (distances[i][0] == 0) - landmarks_ind[knn[i][0]] = i; - } - //write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplices, false, true); - write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, "witness_before_collapse.mesh"); - - // collapsed_tree.set_dimension(limD); - // persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh2(collapsed_tree, true); - // pcoh2.init_coefficients( p ); //initilizes the coefficient field for homology - // pcoh2.compute_persistent_cohomology( alpha2/10 ); - // pcoh2.output_diagram(); - - // chi = 0; - // for (auto sh: simplices) - // chi += 1-2*(simplex_tree.dimension(sh)%2); - // std::cout << "Euler characteristic is " << chi << std::endl; - // write_witness_mesh(point_vector, landmarks_ind, collapsed_tree, collapsed_tree.complex_simplex_range(), false, true, "witness_after_collapse.mesh"); -} - -void rips_experiment(Point_Vector & points, double threshold, int dim_max) -{ - typedef std::vector Point_t; - typedef Simplex_tree ST; - clock_t start, end; - ST st; - - // Compute the proximity graph of the points - start = clock(); - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance); - // Construct the Rips complex in a Simplex Tree - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); - end = clock(); - - double time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Rips complex took " - << time << " s. \n"; - std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; - //std::cout << " and has dimension " << st.dimension() << " \n"; - - // Sort the simplices in the order of the filtration - st.initialize_filtration(); - - // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology pcoh(st); - // initializes the coefficient field for homology - int p = 3; - double min_persistence = -1; //threshold/5; - pcoh.init_coefficients(p); - pcoh.compute_persistent_cohomology(min_persistence); - pcoh.output_diagram(); -} - - -int experiment0 (int argc, char * const argv[]) -{ - if (argc != 8) { - std::cerr << "Usage: " << argv[0] - << " 0 nbP nbL dim alpha limD mu_epsilon\n"; - return 0; - } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - int nbP = atoi(argv[2]); - int nbL = atoi(argv[3]); - int dim = atoi(argv[4]); - double alpha = atof(argv[5]); - int limD = atoi(argv[6]); - double mu_epsilon = atof(argv[7]); - - // Read the point file - Point_Vector point_vector; - generate_points_sphere(point_vector, nbP, dim); - std::cout << "Successfully generated " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - rw_experiment(point_vector, nbL, alpha, limD); - return 0; -} - -// int experiment1 (int argc, char * const argv[]) -// { -// if (argc != 3) { -// std::cerr << "Usage: " << argv[0] -// << " 1 file_name\n"; -// return 0; -// } -// /* -// boost::filesystem::path p; -// for (; argc > 2; --argc, ++argv) -// p /= argv[1]; -// */ - -// std::string file_name = argv[2]; - -// // Read the point file -// Point_Vector point_vector; -// read_points_cust(file_name, point_vector); -// std::cout << "The file contains " << point_vector.size() << " points.\n"; -// std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - -// bool ok = false; -// int nbL, limD; -// double alpha; -// while (!ok) { -// std::cout << "Relaxed witness complex: parameters nbL, alpha, limD.\n"; -// std::cout << "Enter nbL: "; -// std::cin >> nbL; -// std::cout << "Enter alpha: "; -// std::cin >> alpha; -// std::cout << "Enter limD: "; -// std::cin >> limD; -// std::cout << "Start relaxed witness complex...\n"; -// rw_experiment(point_vector, nbL, alpha, limD); -// std::cout << "Is the result correct? [y/n]: "; -// char answer; -// std::cin >> answer; -// switch (answer) { -// case 'n': -// ok = false; break; -// default : -// ok = true; break; -// } -// } -// // ok = false; -// // while (!ok) { -// // std::cout << "Rips complex: parameters threshold, limD.\n"; -// // std::cout << "Enter threshold: "; -// // std::cin >> alpha; -// // std::cout << "Enter limD: "; -// // std::cin >> limD; -// // std::cout << "Start Rips complex...\n"; -// // rips_experiment(point_vector, alpha, limD); -// // std::cout << "Is the result correct? [y/n]: "; -// // char answer; -// // std::cin >> answer; -// // switch (answer) { -// // case 'n': -// // ok = false; break; -// // default : -// // ok = true; break; -// // } -// // } -// return 0; -// } - -int experiment1 (int argc, char * const argv[]) -{ - if (argc != 8) { - std::cerr << "Usage: " << argv[0] - << " 1 file_name nbL alpha mu_epsilon limD experiment_name\n"; - return 0; - } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - std::string file_name = argv[2]; - int nbL = atoi(argv[3]), limD = atoi(argv[6]); - double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]); - std::string experiment_name = argv[7]; - - // Read the point file - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - std::cout << "The file contains " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - Simplex_tree<> simplex_tree; - std::vector > knn; - std::vector > distances; - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - limD, - knn, - distances); - - rw_experiment(point_vector, nbL, alpha2, limD, mu_epsilon); - - // ok = false; - // while (!ok) { - // std::cout << "Rips complex: parameters threshold, limD.\n"; - // std::cout << "Enter threshold: "; - // std::cin >> alpha; - // std::cout << "Enter limD: "; - // std::cin >> limD; - // std::cout << "Start Rips complex...\n"; - // rips_experiment(point_vector, alpha, limD); - // std::cout << "Is the result correct? [y/n]: "; - // char answer; - // std::cin >> answer; - // switch (answer) { - // case 'n': - // ok = false; break; - // default : - // ok = true; break; - // } - // } - return 0; -} - - -/******************************************************************************************** - * Length of the good interval experiment - *******************************************************************************************/ - -struct Pers_endpoint { - double alpha; - bool start; - int dim; - Pers_endpoint(double alpha_, bool start_, int dim_) - : alpha(alpha_), start(start_), dim(dim_) - {} -}; - -/* -struct less_than_key { - inline bool operator() (const MyStruct& struct1, const MyStruct& struct2) { - return (struct1.key < struct2.key); - } -}; -*/ - -double good_interval_length(const std::vector & desired_homology, Simplex_tree<> & simplex_tree, double alpha2) -{ - int nbL = simplex_tree.num_vertices(); - int p = 3; - persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, true); - pcoh.init_coefficients( p ); //initilizes the coefficient field for homology - pcoh.compute_persistent_cohomology( -1 ); - std::ofstream out_stream("pers_diag.tmp"); - pcoh.output_diagram(out_stream); - out_stream.close(); - std::ifstream in_stream("pers_diag.tmp", std::ios::in); - std::string line; - std::vector pers_endpoints; - while (getline(in_stream, line)) { - int p, dim; - double alpha_start, alpha_end; - std::istringstream iss(line); - iss >> p >> dim >> alpha_start >> alpha_end; - if (alpha_start != alpha_end) { - if (alpha_end < alpha_start) - alpha_end = alpha2; - pers_endpoints.push_back(Pers_endpoint(alpha_start, true, dim)); - pers_endpoints.push_back(Pers_endpoint(alpha_end, false, dim)); - } - } - std::cout << "Pers_endpoints.size = " << pers_endpoints.size() << std::endl; - in_stream.close(); - std::sort(pers_endpoints.begin(), - pers_endpoints.end(), - [](const Pers_endpoint & p1, const Pers_endpoint & p2){ - return p1.alpha < p2.alpha;} - ); - write_barcodes("pers_diag.tmp", alpha2); - /* - for (auto p: pers_endpoints) { - std::cout << p.alpha << " " << p.dim << " " << p.start << "\n"; - } - */ - std::vector current_homology(nbL-1,0); - current_homology[0] = 1; // for the compulsary "0 0 inf" entry - double good_start = 0, good_end = 0; - double sum_intervals = 0; - int num_pieces = 0; - bool interval_in_process = (desired_homology == current_homology); - for (auto p: pers_endpoints) { - /* - std::cout << "Treating " << p.alpha << " " << p.dim << " " << p.start - << " ["; - for (int v: current_homology) - std::cout << v << " "; - std::cout << "]\n"; - */ - if (p.start) - current_homology[p.dim]++; - else - current_homology[p.dim]--; - if (interval_in_process) { - good_end = p.alpha; - sum_intervals += good_end - good_start; - std::cout << "good_start = " << good_start - << ", good_end = " << good_end << "\n"; - - Gudhi::witness_complex::Dim_lists> simplices(simplex_tree, nbL-1, (good_end - good_start)/2); - simplices.collapse(); - simplices.output_simplices(); - interval_in_process = false; - //break; - } - else if (desired_homology == current_homology) { - interval_in_process = true; - good_start = p.alpha; - num_pieces++; - } - } - std::cout << "Number of good homology intervals: " << num_pieces << "\n"; - return sum_intervals; -} - -void run_comparison(std::vector > const & knn, - std::vector > const & distances, - unsigned nbL, - double alpha2, - std::vector& desired_homology) -{ - clock_t start, end; - Simplex_tree<> simplex_tree; - - start = clock(); - SRWit srwit(distances, - knn, - simplex_tree, - nbL, - alpha2, - nbL-1); - end = clock(); - std::cout << "SRWit.size = " << simplex_tree.num_simplices() << std::endl; - simplex_tree.set_dimension(nbL-1); - - std::cout << "Good homology interval length for SRWit is " - << good_interval_length(desired_homology, simplex_tree, alpha2) << "\n"; - std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; - - /* - Simplex_tree<> simplex_tree2; - start = clock(); - WRWit wrwit(distances, - knn, - simplex_tree2, - nbL, - alpha2, - nbL-1); - end = clock(); - std::cout << "WRWit.size = " << simplex_tree2.num_simplices() << std::endl; - simplex_tree.set_dimension(nbL-1); - - std::cout << "Good homology interval length for WRWit is " - << good_interval_length(desired_homology, simplex_tree2, alpha2) << "\n"; - std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; - */ -} - -int experiment2(int argc, char * const argv[]) -{ - for (unsigned d = 2; d < 2; d++) { - // Sphere S^d - Point_Vector point_vector; - unsigned N = 1; - double alpha2 = 2.4 - 0.4*d; - switch (d) { - case 1: alpha2 = 2.2; break; - case 2: alpha2 = 1.8; break; - case 3: alpha2 = 1.5; break; - case 4: alpha2 = 1.4; break; - default: alpha2 = 1.4; break; - } - std::cout << "alpha2 = " << alpha2 << "\n"; - unsigned nbL = 20; - std::vector desired_homology(nbL-1,0); - desired_homology[0] = 1; desired_homology[d] = 1; - - - for (unsigned i = 1; i <= N; ++i) { - unsigned nbW = 1000*i;//, nbL = 20; - double mu_epsilon = 1/sqrt(nbL); - std::cout << "Running test S"<< d <<", |W|=" << nbW << ", |L|=" << nbL << std::endl; - generate_points_sphere(point_vector, i*1000, d+1); - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector > knn; - std::vector > distances; - - std::cout << "|L| after sparsification: " << landmarks.size() << "\n"; - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - nbL-1, - knn, - distances); - run_comparison(knn, distances, nbL, alpha2, desired_homology); - } - } - { - // SO(3) - Point_Vector point_vector; - double alpha2 = 0.6; - std::cout << "alpha2 = " << alpha2 << "\n"; - unsigned nbL = 150; - std::vector desired_homology(nbL-1,0); - desired_homology[0] = 1; desired_homology[1] = 1; desired_homology[2] = 1; //Kl - // desired_homology[0] = 1; desired_homology[3] = 1; //SO3 - - double mu_epsilon = 1/sqrt(nbL); - if (argc < 3) std::cerr << "No file name indicated!\n"; - read_points_cust(argv[2], point_vector); - int nbW = point_vector.size(); - std::cout << "Running test SO(3), |W|=" << nbW << ", |L|=" << nbL << std::endl; - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector > knn; - std::vector > distances; - - std::cout << "|L| after sparsification: " << landmarks.size() << "\n"; - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - nbL-1, - knn, - distances); - run_comparison(knn, distances, nbL, alpha2, desired_homology); - } - return 0; -} - -int experiment3(int argc, char * const argv[]) -{ - // COLLAPSES EXPERIMENT - - return 0; -} - -int main (int argc, char * const argv[]) -{ - if (argc == 1) { - output_experiment_information(argv[0]); - return 1; - } - switch (atoi(argv[1])) { - case 0 : - return experiment0(argc, argv); - break; - case 1 : - return experiment1(argc, argv); - break; - case 2 : - return experiment2(argc, argv); - break; - default : - output_experiment_information(argv[0]); - return 1; - } -} diff --git a/src/Witness_complex/example/bench_rwit.cpp b/src/Witness_complex/example/bench_rwit.cpp deleted file mode 100644 index 2d3a009c..00000000 --- a/src/Witness_complex/example/bench_rwit.cpp +++ /dev/null @@ -1,709 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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 - * 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 -#include -#include "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include - -#include "generators.h" -#include "output.h" -#include "output_tikz.h" - -using namespace Gudhi; -using namespace Gudhi::witness_complex; -using namespace Gudhi::persistent_cohomology; - -typedef std::vector Point_Vector; -//typedef Simplex_tree STree; -typedef Simplex_tree<> STree; -typedef A0_complex< STree> A0Complex; -typedef STree::Simplex_handle Simplex_handle; - -typedef A0_complex< STree > SRWit; -typedef Relaxed_witness_complex< STree > WRWit; - -/** Program options *************************************************************** -*********************************************************************************** -*********************************************************************************** -*********************************************************************************** -**********************************************************************************/ - - -void program_options(int argc, char * const argv[] - , int & experiment_number - , std::string & filepoints - , std::string & landmark_file - , std::string & experiment_name - , int & nbL - , double & alpha2_s - , double & alpha2_w - , double & mu_epsilon - , int & dim_max - , std::vector & desired_homology - , double & min_persistence) { - namespace po = boost::program_options; - po::options_description hidden("Hidden options"); - hidden.add_options() - ("option", po::value(& experiment_number), - "Experiment id.") - ("input-file", po::value(&filepoints), - "Name of file containing a point set. Format is one point per line: X1 ... Xd "); - - po::options_description visible("Allowed options", 100); - visible.add_options() - ("help,h", "produce help message") - ("output-file,o", po::value(&experiment_name)->default_value("witness"), - "The prefix of all the output files. Default is 'witness'") - ("landmarks,L", po::value(&nbL)->default_value(0), - "Number of landmarks.") - ( "landmark-file,l", po::value(&landmark_file), - "Name of a fike containing landmarks") - ("alpha2_s,A", po::value(&alpha2_s)->default_value(0), - "Relaxation parameter for the strong complex.") - ("alpha2_w,a", po::value(&alpha2_w)->default_value(0), - "Relaxation parameter for the weak complex.") - ("mu_epsilon,e", po::value(&mu_epsilon)->default_value(0), - "Sparsification parameter.") - ("cpx-dimension,d", po::value(&dim_max)->default_value(1), - "Maximal dimension of the Witness complex we want to compute.") - - ("homology,H", po::value>(&desired_homology)->multitoken(), - "The desired Betti numbers.") - ("min-persistence,m", po::value(&min_persistence), - "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); - - po::positional_options_description pos; - pos.add("option", 1); - pos.add("input-file", 2); - - po::options_description all; - all.add(visible).add(hidden); - - po::variables_map vm; - po::store(po::command_line_parser(argc, argv). - options(all).positional(pos).run(), vm); - po::notify(vm); - - if (vm.count("help") || !vm.count("input-file")) { - std::cout << std::endl; - std::cout << "Compute the persistent homology with coefficient field Z/3Z \n"; - std::cout << "of a Strong relaxed witness complex defined on a set of input points.\n \n"; - std::cout << "The output diagram contains one bar per line, written with the convention: \n"; - std::cout << " p dim b d \n"; - std::cout << "where dim is the dimension of the homological feature,\n"; - std::cout << "b and d are respectively the birth and death of the feature and \n"; - std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; - - std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; - std::cout << visible << std::endl; - std::abort(); - } -} - - - - - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { - std::ifstream in_file(file_name.c_str(), std::ios::in); - if (!in_file.is_open()) { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while (getline(in_file, line)) { - std::vector< double > point; - std::istringstream iss(line); - while (iss >> x) { - point.push_back(x); - } - if (point.size() != 1) - points.push_back(point); - } - in_file.close(); -} - -void rips(Point_Vector & points, double alpha2, int dim_max, STree& st) -{ - Graph_t prox_graph = compute_proximity_graph(points, sqrt(alpha2) - , euclidean_distance >); - // Construct the Rips complex in a Simplex Tree - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); -} - -void output_experiment_information(char * const file_name) -{ - std::cout << "Enter a valid experiment number. Usage: " - << file_name << " exp_no options\n"; - std::cout << "Experiment description:\n" - << "0 nbP nbL dim alpha limD mu_epsilon: " - << "Build persistence diagram on relaxed witness complex " - << "built from a point cloud on (dim-1)-dimensional sphere " - << "consisting of nbP witnesses and nbL landmarks. " - << "The maximal relaxation is alpha and the limit on simplicial complex " - << "dimension is limD.\n"; - std::cout << "1 file_name nbL alpha limD: " - << "Build persistence diagram on relaxed witness complex " - << "build from a point cloud stored in a file and nbL landmarks. " - << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n"; -} - -void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT mu_epsilon = 0.1, - std::string mesh_filename = "witness") -{ - clock_t start, end; - STree simplex_tree; - - // Choose landmarks - std::vector > knn; - std::vector > distances; - start = clock(); - //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); - - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - limD, - knn, - distances); - end = clock(); - double time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Choice of " << nbL << " landmarks took " - << time << " s. \n"; - // Compute witness complex - start = clock(); - A0Complex rw(distances, - knn, - simplex_tree, - nbL, - alpha2, - limD); - end = clock(); - time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Witness complex for " << nbL << " landmarks took " - << time << " s. \n"; - std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; - //std::cout << simplex_tree << "\n"; - - // Compute the persistence diagram of the complex - simplex_tree.set_dimension(limD); - persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh(simplex_tree, true); - int p = 3; - pcoh.init_coefficients( p ); //initilizes the coefficient field for homology - start = clock(); - pcoh.compute_persistent_cohomology( alpha2/10 ); - end = clock(); - time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Persistence diagram took " - << time << " s. \n"; - pcoh.output_diagram(); - - int chi = 0; - for (auto sh: simplex_tree.complex_simplex_range()) - chi += 1-2*(simplex_tree.dimension(sh)%2); - std::cout << "Euler characteristic is " << chi << std::endl; - - Gudhi::witness_complex::Dim_lists simplices(simplex_tree, limD); - - // std::vector simplices; - std::cout << "Starting collapses...\n"; - simplices.collapse(); - simplices.output_simplices(); - - STree collapsed_tree; - for (auto sh: simplices) { - std::vector vertices; - for (int v: collapsed_tree.simplex_vertex_range(sh)) - vertices.push_back(v); - collapsed_tree.insert_simplex(vertices); - } - std::vector landmarks_ind(nbL); - for (unsigned i = 0; i != distances.size(); ++i) { - if (distances[i][0] == 0) - landmarks_ind[knn[i][0]] = i; - } - //write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplices, false, true); - write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, mesh_filename+"_before_collapse.mesh"); - - collapsed_tree.set_dimension(limD); - persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh2(collapsed_tree, true); - pcoh2.init_coefficients( p ); //initilizes the coefficient field for homology - pcoh2.compute_persistent_cohomology( alpha2/10 ); - pcoh2.output_diagram(); - - chi = 0; - for (auto sh: simplices) - chi += 1-2*(simplex_tree.dimension(sh)%2); - std::cout << "Euler characteristic is " << chi << std::endl; - write_witness_mesh(point_vector, landmarks_ind, collapsed_tree, collapsed_tree.complex_simplex_range(), false, true, mesh_filename+"_after_collapse.mesh"); - Gudhi::Good_links gl(collapsed_tree); - if (gl.complex_is_pseudomanifold()) - std::cout << "Collapsed complex is a pseudomanifold.\n"; - else - std::cout << "Collapsed complex is NOT a pseudomanifold.\n"; - bool good = true; - for (auto v: collapsed_tree.complex_vertex_range()) - if (!gl.has_good_link(v)) { - std::cout << "Bad link around " << v << std::endl; - good = false; - } - if (good) - std::cout << "All links are good.\n"; - else - std::cout << "There are bad links.\n"; -} - -void rips_experiment(Point_Vector & points, double threshold, int dim_max) -{ - typedef STree ST; - clock_t start, end; - ST st; - - // Compute the proximity graph of the points - start = clock(); - rips(points, threshold, dim_max, st); - end = clock(); - - double time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Rips complex took " - << time << " s. \n"; - std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; - //std::cout << " and has dimension " << st.dimension() << " \n"; - - // Sort the simplices in the order of the filtration - st.initialize_filtration(); - - // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology pcoh(st); - // initializes the coefficient field for homology - int p = 3; - double min_persistence = -1; //threshold/5; - pcoh.init_coefficients(p); - pcoh.compute_persistent_cohomology(min_persistence); - pcoh.output_diagram(); -} - - -int experiment0 (int argc, char * const argv[]) -{ - if (argc != 8) { - std::cerr << "Usage: " << argv[0] - << " 0 nbP nbL dim alpha limD mu_epsilon\n"; - return 0; - } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - int nbP = atoi(argv[2]); - int nbL = atoi(argv[3]); - int dim = atoi(argv[4]); - double alpha = atof(argv[5]); - int limD = atoi(argv[6]); - double mu_epsilon = atof(argv[7]); - - // Read the point file - Point_Vector point_vector; - generate_points_sphere(point_vector, nbP, dim); - std::cout << "Successfully generated " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - rw_experiment(point_vector, nbL, alpha, limD); - return 0; -} - - -/******************************************************************************************** - * Length of the good interval experiment - *******************************************************************************************/ - -struct Pers_endpoint { - double alpha; - bool start; - int dim; - Pers_endpoint(double alpha_, bool start_, int dim_) - : alpha(alpha_), start(start_), dim(dim_) - {} -}; - -/* -struct less_than_key { - inline bool operator() (const MyStruct& struct1, const MyStruct& struct2) { - return (struct1.key < struct2.key); - } -}; -*/ - -double good_interval_length(const std::vector & desired_homology, STree & simplex_tree, double alpha2) -{ - int nbL = simplex_tree.num_vertices(); - int p = 3; - persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh(simplex_tree, true); - pcoh.init_coefficients( p ); //initilizes the coefficient field for homology - pcoh.compute_persistent_cohomology( -1 ); - std::ofstream out_stream("pers_diag.tmp"); - pcoh.output_diagram(out_stream); - out_stream.close(); - std::ifstream in_stream("pers_diag.tmp", std::ios::in); - std::string line; - std::vector pers_endpoints; - while (getline(in_stream, line)) { - unsigned p, dim; - double alpha_start, alpha_end; - std::istringstream iss(line); - iss >> p >> dim >> alpha_start >> alpha_end; - if (iss.fail()) - alpha_end = alpha2; - //std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n"; - //if (dim < desired_homology.size()+1) - if (alpha_start != alpha_end) { - // if (alpha_end < alpha_start) - // alpha_end = alpha2; - pers_endpoints.push_back(Pers_endpoint(alpha_start, true, dim)); - pers_endpoints.push_back(Pers_endpoint(alpha_end, false, dim)); - std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n"; - } - } - std::cout << "desired_homology.size() = " << desired_homology.size() << "\n"; - for (auto nd: desired_homology) - std::cout << nd << std::endl; - std::cout << "Pers_endpoints.size = " << pers_endpoints.size() << std::endl; - in_stream.close(); - std::sort(pers_endpoints.begin(), - pers_endpoints.end(), - [](const Pers_endpoint & p1, const Pers_endpoint & p2){ - return p1.alpha < p2.alpha;} - ); - write_barcodes("pers_diag.tmp", alpha2); - /* - for (auto p: pers_endpoints) { - std::cout << p.alpha << " " << p.dim << " " << p.start << "\n"; - } - */ - std::vector current_homology(desired_homology.size(),0); - //current_homology[0] = 1; // for the compulsary "0 0 inf" entry - double good_start = 0, good_end = 0; - double sum_intervals = 0; - int num_pieces = 0; - bool interval_in_process = (desired_homology == current_homology); - for (auto p: pers_endpoints) { - /* - std::cout << p.alpha << " " << p.dim << " "; - if (p.start) - std::cout << "s\n"; - else - std::cout << "e\n"; - */ - /* - std::cout << "Treating " << p.alpha << " " << p.dim << " " << p.start - << " ["; - for (int v: current_homology) - std::cout << v << " "; - std::cout << "]\n"; - */ - if (p.start) - current_homology[p.dim]++; - else - current_homology[p.dim]--; - if (interval_in_process) { - good_end = p.alpha; - sum_intervals += good_end - good_start; - std::cout << "good_start = " << good_start - << ", good_end = " << good_end << "\n"; - - Gudhi::witness_complex::Dim_lists simplices(simplex_tree, nbL-1, (good_end - good_start)/2); - //simplices.collapse(); - //simplices.output_simplices(); - interval_in_process = false; - //break; - } - else if (desired_homology == current_homology) { - interval_in_process = true; - good_start = p.alpha; - num_pieces++; - } - } - std::cout << "Number of good homology intervals: " << num_pieces << "\n"; - return sum_intervals; -} - - - -void run_comparison(std::vector > const & knn, - std::vector > const & distances, - Point_Vector & points, - unsigned nbL, - unsigned limD, - double alpha2_s, - double alpha2_w, - std::vector& desired_homology) -{ - clock_t start, end; - STree simplex_tree; - - //std::cout << "alpha2 = " << alpha2_s << "\n"; - start = clock(); - SRWit srwit(distances, - knn, - simplex_tree, - nbL, - alpha2_s, - limD); - end = clock(); - std::cout << "SRWit.size = " << simplex_tree.num_simplices() << std::endl; - simplex_tree.set_dimension(desired_homology.size()); - - std::cout << "Good homology interval length for SRWit is " - << good_interval_length(desired_homology, simplex_tree, alpha2_s) << "\n"; - std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; - int chi = 0; - for (auto sh: simplex_tree.complex_simplex_range()) - chi += 1-2*(simplex_tree.dimension(sh)%2); - std::cout << "Euler characteristic is " << chi << std::endl; - - - STree simplex_tree2; - std::cout << "alpha2 = " << alpha2_w << "\n"; - start = clock(); - WRWit wrwit(distances, - knn, - simplex_tree2, - nbL, - alpha2_w, - limD); - end = clock(); - std::cout << "WRWit.size = " << simplex_tree2.num_simplices() << std::endl; - simplex_tree2.set_dimension(nbL-1); - - std::cout << "Good homology interval length for WRWit is " - << good_interval_length(desired_homology, simplex_tree2, alpha2_w) << "\n"; - std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; - chi = 0; - for (auto sh: simplex_tree2.complex_simplex_range()) - chi += 1-2*(simplex_tree2.dimension(sh)%2); - std::cout << "Euler characteristic is " << chi << std::endl; - - //write_witness_mesh(points, landmarks_ind, simplex_tree2, simplex_tree2.complex_simplex_range(), false, true, "wrwit.mesh"); - - -} - -int experiment1 (int argc, char * const argv[]) -{ - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - // std::string file_name = argv[2]; - // int nbL = atoi(argv[3]), limD = atoi(argv[6]); - // double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]); - // std::string experiment_name = argv[7]; - - int option = 1; - std::string file_name, landmark_file; - int nbL = 0, limD; - double alpha2_s, alpha2_w, mu_epsilon, min_pers; - std::string experiment_name; - std::vector desired_homology = {1}; - std::vector landmarks; - - program_options(argc, argv, option, file_name, landmark_file, experiment_name, nbL, alpha2_s, alpha2_w, mu_epsilon, limD, desired_homology, min_pers); - - // Read the point file - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - //std::cout << "The file contains " << point_vector.size() << " points.\n"; - //std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - //std::cout << "Limit dimension for the complexes is " << limD << ".\n"; - - if (landmark_file == "") - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha2_s, limD, knn, distances); - else - read_points_cust(landmark_file, landmarks); - nbL = landmarks.size(); - STree simplex_tree; - std::vector > knn; - std::vector > distances; - - //Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2_s, - limD, - knn, - distances); - - run_comparison(knn, distances, point_vector, nbL, limD, alpha2_s, alpha2_w, desired_homology); - return 0; -} - - -int experiment2(int argc, char * const argv[]) -{ - for (unsigned d = 3; d < 4; d++) { - // Sphere S^d - Point_Vector point_vector; - unsigned N = 1; - double alpha2 = 2.4 - 0.4*d; - switch (d) { - case 1: alpha2 = 2.2; break; - case 2: alpha2 = 1.7; break; - case 3: alpha2 = 1.5; break; - case 4: alpha2 = 1.4; break; - default: alpha2 = 1.4; break; - } - unsigned nbL = 20; - std::vector desired_homology(nbL-1,0); - desired_homology[0] = 1; desired_homology[d] = 1; - - - for (unsigned i = 1; i <= N; ++i) { - unsigned nbW = 1000*i;//, nbL = 20; - double mu_epsilon = 1/sqrt(nbL); - std::cout << "Running test S"<< d <<", |W|=" << nbW << ", |L|=" << nbL << std::endl; - generate_points_sphere(point_vector, i*1000, d+1); - std::vector landmarks; - - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector > knn; - std::vector > distances; - - - std::cout << "|L| after sparsification: " << landmarks.size() << "\n"; - - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - nbL-1, - knn, - distances); - run_comparison(knn, distances, point_vector, nbL, nbL-1, alpha2, alpha2, desired_homology); - } - } - /* - { - // SO(3) - Point_Vector point_vector; - double alpha2 = 0.6; - std::cout << "alpha2 = " << alpha2 << "\n"; - unsigned nbL = 150; - std::vector desired_homology(nbL-1,0); - desired_homology[0] = 1; desired_homology[1] = 1; desired_homology[2] = 1; //Kl - // desired_homology[0] = 1; desired_homology[3] = 1; //SO3 - - double mu_epsilon = 1/sqrt(nbL); - if (argc < 3) std::cerr << "No file name indicated!\n"; - read_points_cust(argv[2], point_vector); - int nbW = point_vector.size(); - std::cout << "Running test SO(3), |W|=" << nbW << ", |L|=" << nbL << std::endl; - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector > knn; - std::vector > distances; - - std::cout << "|L| after sparsification: " << landmarks.size() << "\n"; - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - nbL-1, - knn, - distances); - run_comparison(knn, distances, point_vector, nbL, alpha2, desired_homology); - } - */ - return 0; -} - -int experiment3(int argc, char * const argv[]) -{ - // Both witnesses and landmarks are given as input - - - return 0; -} - -int main (int argc, char * const argv[]) -{ - if (argc == 1) { - output_experiment_information(argv[0]); - return 1; - } - switch (atoi(argv[1])) { - case 0 : - return experiment0(argc, argv); - break; - case 1 : - return experiment1(argc, argv); - break; - case 2 : - return experiment2(argc, argv); - break; - case 3 : - return experiment3(argc, argv); - break; - default : - output_experiment_information(argv[0]); - return 1; - } -} diff --git a/src/Witness_complex/example/color-picker.cpp b/src/Witness_complex/example/color-picker.cpp deleted file mode 100644 index fcea033e..00000000 --- a/src/Witness_complex/example/color-picker.cpp +++ /dev/null @@ -1,37 +0,0 @@ -#include -#include - -struct Color { - double r, g, b; - Color() : r(0), g(0), b(0) { } - Color(double r_, double g_, double b_) - : r(r_), g(g_), b(b_) { } -}; - -// Convert a double in [0,1] to a color -template< typename ColorType > -void double_to_color(double t, ColorType & color) { - assert(t >= 0 && t <= 1); - double s = 1.0/6; - if (t >= 0 && t <= s) - color = ColorType(1, 6*t, 0); - else if (t > s && t <= 2*s) - color = ColorType(1-6*(t-s), 1, 0); - else if (t > 2*s && t <= 3*s) - color = ColorType(0, 1, 6*(t-2*s)); - else if (t > 3*s && t <= 4*s) - color = ColorType(0, 1-6*(t-3*s), 1); - else if (t > 4*s && t <= 5*s) - color = ColorType(6*(t-4*s), 0, 1); - else - color = ColorType(1, 0, 1-6*(t-5*s)); -} - -int main (int argc, char * const argv[]) { - double number_to_convert = 0; - std::cout << "Enter a real number in [0,1]: "; - std::cin >> number_to_convert; - Color color; - double_to_color(number_to_convert, color); - std::cout << "The corresponding color in rgb is: (" << color.r << ", " << color.g << ", " << color.b << ")\n"; -} diff --git a/src/Witness_complex/example/off_writer.h b/src/Witness_complex/example/off_writer.h deleted file mode 100644 index 8cafe9e2..00000000 --- a/src/Witness_complex/example/off_writer.h +++ /dev/null @@ -1,11 +0,0 @@ -#ifndef OFF_WRITER_H_ -#define OFF_WRITER_H_ - - -class Off_Writer { - -private: - -} - -#endif diff --git a/src/Witness_complex/example/relaxed_delaunay.cpp b/src/Witness_complex/example/relaxed_delaunay.cpp deleted file mode 100644 index b0190446..00000000 --- a/src/Witness_complex/example/relaxed_delaunay.cpp +++ /dev/null @@ -1,180 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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 - * 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 "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -//#include - -#include "generators.h" -#include "output.h" - -using namespace Gudhi; -using namespace Gudhi::witness_complex; -using namespace Gudhi::persistent_cohomology; - -typedef CGAL::Epick_d K; -typedef K::Point_d Point_d; -typedef K::Sphere_d Sphere_d; -typedef CGAL::Delaunay_triangulation Delaunay_triangulation; - -typedef std::vector Point_Vector; -typedef Relaxed_witness_complex< Simplex_tree<> > RelaxedWitnessComplex; -typedef Simplex_tree<>::Simplex_handle Simplex_handle; - - - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { - std::ifstream in_file(file_name.c_str(), std::ios::in); - if (!in_file.is_open()) { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while (getline(in_file, line)) { - std::vector< double > point; - std::istringstream iss(line); - while (iss >> x) { - point.push_back(x); - } - if (point.size() != 1) - points.push_back(point); - } - in_file.close(); -} - -int main (int argc, char * const argv[]) -{ - if (argc != 4) { - std::cerr << "Usage: " << argv[0] - << " 1 file_name alpha limD\n"; - return 0; - } - std::string file_name = argv[1]; - double alpha2 = atof(argv[2]); - int limD = atoi(argv[3]); - - // Read points - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - generate_points_random_box(point_vector, 200, 2); - write_points(file_name, point_vector); - - std::cout << "The file contains " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - // 1. Compute Delaunay centers - Delaunay_triangulation delaunay(point_vector[0].size()); - delaunay.insert(point_vector.begin(), point_vector.end()); - Point_Vector del_centers; - for (auto f_it = delaunay.full_cells_begin(); f_it != delaunay.full_cells_end(); ++f_it) { - if (delaunay.is_infinite(f_it)) - continue; - Point_Vector vertices; - for (auto v_it = f_it->vertices_begin(); v_it != f_it->vertices_end(); ++v_it) - vertices.push_back((*v_it)->point()); - Sphere_d sphere(vertices.begin(), vertices.end()); - del_centers.push_back(sphere.center()); - } - std::cout << "Delaunay center count: " << del_centers.size() << ".\n"; - - // 2. Build Relaxed Witness Complex - std::vector> knn; - std::vector> distances; - Gudhi::witness_complex::build_distance_matrix(del_centers, // aka witnesses - point_vector, // aka landmarks - alpha2, - limD, - knn, - distances); - - write_wl("wl_distances.txt", distances); - Simplex_tree<> simplex_tree; - Gudhi::witness_complex::Relaxed_witness_complex> rwc(distances, - knn, - simplex_tree, - point_vector.size(), - alpha2, - limD); - std::vector dim_simplices(limD+1); - for (auto sh: simplex_tree.complex_simplex_range()) { - dim_simplices[simplex_tree.dimension(sh)]++; - } - for (unsigned i =0; i != dim_simplices.size(); ++i) - std::cout << "dim[" << i << "]: " << dim_simplices[i] << " simplices.\n"; - - std::vector landmarks_ind; - for (unsigned i = 0; i < point_vector.size(); ++i) - landmarks_ind.push_back(i); - write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), true, true, "relaxed_delaunay.mesh"); - - // 3. Check if the thing is Relaxed Delaunay - for (auto sh: simplex_tree.complex_simplex_range()) { - Point_Vector vertices; - for (auto v: simplex_tree.simplex_vertex_range(sh)) - vertices.push_back(point_vector[v]); - Sphere_d sphere(vertices.begin(), vertices.end()); - Point_d center = sphere.center(); - double r2 = sphere.squared_radius(); - typename K::Squared_distance_d dist2; - std::vector v_inds; - for (auto v: simplex_tree.simplex_vertex_range(sh)) - v_inds.push_back(v); - auto range_begin = std::begin(v_inds); - auto range_end = std::end(v_inds); - if (simplex_tree.dimension(sh) == (int)point_vector[0].size()) - for (auto v: simplex_tree.complex_vertex_range()) - if (std::find(range_begin, range_end, v) == range_end) { - if (dist2(point_vector[v], center) < r2 - alpha2) - std::cout << "WARNING! The vertex " << point_vector[v] << " is inside the (r2-alpha2)-ball (" << center << ", " << r2 << ") distance is " << dist2(point_vector[v], center) << "\n"; - } - } -} diff --git a/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp deleted file mode 100644 index 067321ce..00000000 --- a/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp +++ /dev/null @@ -1,461 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 -#include -#include - -#include -#include -//#include - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Relaxed_witness_complex.h" -#include "gudhi/reader_utils.h" -#include "gudhi/Collapse/Collapse.h" -//#include - -//#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - - -#include -#include -#include -#include - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef CGAL::Epick_d K; -typedef K::FT FT; -typedef K::Point_d Point_d; -typedef CGAL::Search_traits< - FT, Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d> Traits_base; -typedef CGAL::Euclidean_distance Euclidean_distance; - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter Euclidean_adapter; -//typedef CGAL::Kd_tree Kd_tree; -typedef CGAL::Orthogonal_incremental_neighbor_search> Neighbor_search; -typedef Neighbor_search::Tree Tree; -typedef Neighbor_search::Distance Distance; -typedef Neighbor_search::iterator KNS_iterator; -typedef Neighbor_search::iterator KNS_range; -typedef boost::container::flat_map Point_etiquette_map; -typedef CGAL::Kd_tree Tree2; - -typedef CGAL::Fuzzy_sphere Fuzzy_sphere; - -typedef std::vector Point_Vector; - -//typedef K::Equal_d Equal_d; -typedef CGAL::Random_points_in_cube_d Random_cube_iterator; -typedef CGAL::Random_points_in_ball_d Random_point_iterator; - -bool toric=false; - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust ( std::string file_name , Point_Vector & points) -{ - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector< double > point; - std::istringstream iss( line ); - while(iss >> x) { point.push_back(x); } - Point_d p(point.begin(), point.end()); - if (point.size() != 1) - points.push_back(p); - } - in_file.close(); -} - - -void generate_points_sphere(Point_Vector& W, int nbP, int dim) -{ - CGAL::Random_points_on_sphere_d rp(dim,1); - for (int i = 0; i < nbP; i++) - W.push_back(*rp++); -} - - -void write_wl( std::string file_name, std::vector< std::vector > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} - -void write_rl( std::string file_name, std::vector< std::vector ::iterator> > & rl) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : rl) - { - for (auto l: w) - ofs << *l << " "; - ofs << "\n"; - } - ofs.close(); -} - -std::vector convert_to_torus(std::vector< Point_d>& points) -{ - std::vector< Point_d > points_torus; - for (auto p: points) - { - FT theta = M_PI*p[0]; - FT phi = M_PI*p[1]; - std::vector p_torus; - p_torus.push_back((1+0.2*cos(theta))*cos(phi)); - p_torus.push_back((1+0.2*cos(theta))*sin(phi)); - p_torus.push_back(0.2*sin(theta)); - points_torus.push_back(Point_d(p_torus)); - } - return points_torus; -} - - -void write_points_torus( std::string file_name, std::vector< Point_d > & points) -{ - std::ofstream ofs (file_name, std::ofstream::out); - std::vector points_torus = convert_to_torus(points); - for (auto w : points_torus) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); -} - - -void write_points( std::string file_name, std::vector< Point_d > & points) -{ - if (toric) write_points_torus(file_name, points); - else - { - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : points) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); - } -} - - -void write_edges_torus(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - Point_Vector l_torus = convert_to_torus(landmarks); - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = l_torus[u].cartesian_begin(); it != l_torus[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = l_torus[v].cartesian_begin(); it != l_torus[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); -} - -void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - if (toric) write_edges_torus(file_name, witness_complex, landmarks); - else - { - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); - } -} - - -/** Function that chooses landmarks from W and place it in the kd-tree L. - * Note: nbL hould be removed if the code moves to Witness_complex - */ -void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - //std::vector landmarks; - int chosen_landmark; - //std::pair res = std::make_pair(L_i.begin(),false); - Point_d* p; - CGAL::Random rand; - for (int i = 0; i < nbL; i++) - { - // while (!res.second) - // { - do chosen_landmark = rand.get_int(0,nbP); - while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end()); - //rand++; - //std::cout << "Chose " << chosen_landmark << std::endl; - p = &W[chosen_landmark]; - //L_i.emplace(chosen_landmark,i); - // } - landmarks.push_back(*p); - landmarks_ind.push_back(chosen_landmark); - //std::cout << "Added landmark " << chosen_landmark << std::endl; - } - } - - -void landmarks_to_witness_complex(Point_Vector &W, Point_Vector& landmarks, std::vector& landmarks_ind, FT alpha) -{ - //********************Preface: origin point - unsigned D = W[0].size(); - std::vector orig_vector; - for (unsigned i = 0; i < D; i++) - orig_vector.push_back(0); - Point_d origin(orig_vector); - //Distance dist; - //dist.transformed_distance(0,1); - //******************** Constructing a WL matrix - int nbP = W.size(); - int nbL = landmarks.size(); - STraits traits(&(landmarks[0])); - Euclidean_distance ed; - std::vector< std::vector > WL(nbP); - std::vector< std::vector< typename std::vector::iterator > > ope_limits(nbP); - Tree L(boost::counting_iterator(0), - boost::counting_iterator(nbL), - typename Tree::Splitter(), - traits); - - std::cout << "Enter (D+1) nearest landmarks\n"; - //std::cout << "Size of the tree is " << L.size() << std::endl; - for (int i = 0; i < nbP; i++) - { - //std::cout << "Entered witness number " << i << std::endl; - Point_d& w = W[i]; - std::queue< typename std::vector::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark - Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter(&(landmarks[0]))); - Neighbor_search::iterator search_it = search.begin(); - - //Incremental search and filling WL - while (WL[i].size() < D) - WL[i].push_back((search_it++)->first); - FT dtow = ed.transformed_distance(w, landmarks[WL[i][D-1]]); - while (search_it->second < dtow + alpha) - WL[i].push_back((search_it++)->first); - - //Filling the (1+epsilon)-limits table - for (std::vector::iterator wl_it = WL[i].begin(); wl_it != WL[i].end(); ++wl_it) - { - ope_queue.push(wl_it); - FT d_to_curr_l = ed.transformed_distance(w, landmarks[*wl_it]); - //std::cout << "d_to_curr_l=" << d_to_curr_l << std::endl; - //std::cout << "d_to_front+alpha=" << d_to_curr_l << std::endl; - while (d_to_curr_l > alpha + ed.transformed_distance(w, landmarks[*(ope_queue.front())])) - { - ope_limits[i].push_back(wl_it); - ope_queue.pop(); - } - } - while (ope_queue.size() > 0) - { - ope_limits[i].push_back(WL[i].end()); - ope_queue.pop(); - } - //std::cout << "Safely constructed a point\n"; - ////Search D+1 nearest neighbours from the tree of landmarks L - /* - if (w[0]>0.95) - std::cout << i << std::endl; - */ - //K_neighbor_search search(L, w, D, FT(0), true, - // CGAL::Distance_adapter(&(landmarks[0])) ); - //std::cout << "Safely found nearest landmarks\n"; - /* - for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) - { - //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; - //Point_etiquette_map::iterator itm = L_i.find(it->first); - //assert(itm != L_i.end()); - //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; - WL[i].push_back(it->first); - //std::cout << "ITFIRST " << it->first << std::endl; - //std::cout << i << " " << it->first << ": " << it->second << std::endl; - } - */ - } - //std::cout << "\n"; - - //std::string out_file = "wl_result"; - write_wl("wl_result",WL); - write_rl("rl_result",ope_limits); - - //******************** Constructng a witness complex - std::cout << "Entered witness complex construction\n"; - Witness_complex<> witnessComplex; - witnessComplex.setNbL(nbL); - witnessComplex.relaxed_witness_complex(WL, ope_limits); - char buffer[100]; - int i = sprintf(buffer,"stree_result.txt"); - - if (i >= 0) - { - std::string out_file = (std::string)buffer; - std::ofstream ofs (out_file, std::ofstream::out); - witnessComplex.st_to_file(ofs); - ofs.close(); - } - write_edges("landmarks/edges", witnessComplex, landmarks); - std::cout << Distance().transformed_distance(Point_d(std::vector({0.1,0.1})), Point_d(std::vector({1.9,1.9}))) << std::endl; -} - - -int main (int argc, char * const argv[]) -{ - - if (argc != 5) - { - std::cerr << "Usage: " << argv[0] - << " nbP nbL dim alpha\n"; - return 0; - } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - int nbP = atoi(argv[1]); - int nbL = atoi(argv[2]); - int dim = atoi(argv[3]); - double alpha = atof(argv[4]); - //clock_t start, end; - //Construct the Simplex Tree - Witness_complex<> witnessComplex; - - std::cout << "Let the carnage begin!\n"; - Point_Vector point_vector; - //read_points_cust(file_name, point_vector); - generate_points_sphere(point_vector, nbP, dim); - /* - for (auto &p: point_vector) - { - assert(std::count(point_vector.begin(),point_vector.end(),p) == 1); - } - */ - //std::cout << "Successfully read the points\n"; - //witnessComplex.setNbL(nbL); - Point_Vector L; - std::vector chosen_landmarks; - landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); - //start = clock(); - - write_points("landmarks/initial_pointset",point_vector); - write_points("landmarks/initial_landmarks",L); - - landmarks_to_witness_complex(point_vector, L, chosen_landmarks, alpha); - //end = clock(); - - /* - std::cout << "Landmark choice took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - start = clock(); - witnessComplex.witness_complex(WL); - // - end = clock(); - std::cout << "Howdy world! The process took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - */ - - /* - out_file = "output/"+file_name+"_"+argv[2]+".stree"; - std::ofstream ofs (out_file, std::ofstream::out); - witnessComplex.st_to_file(ofs); - ofs.close(); - - out_file = "output/"+file_name+"_"+argv[2]+".badlinks"; - std::ofstream ofs2(out_file, std::ofstream::out); - witnessComplex.write_bad_links(ofs2); - ofs2.close(); - */ -} diff --git a/src/Witness_complex/example/relaxed_witness_persistence.cpp b/src/Witness_complex/example/relaxed_witness_persistence.cpp deleted file mode 100644 index b4837594..00000000 --- a/src/Witness_complex/example/relaxed_witness_persistence.cpp +++ /dev/null @@ -1,267 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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 - * 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 "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include "generators.h" -#include "output.h" - -using namespace Gudhi; -using namespace Gudhi::witness_complex; -using namespace Gudhi::persistent_cohomology; - -typedef std::vector Point_Vector; -typedef Relaxed_witness_complex< Simplex_tree<> > RelaxedWitnessComplex; -typedef Simplex_tree<>::Simplex_handle Simplex_handle; - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { - std::ifstream in_file(file_name.c_str(), std::ios::in); - if (!in_file.is_open()) { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while (getline(in_file, line)) { - std::vector< double > point; - std::istringstream iss(line); - while (iss >> x) { - point.push_back(x); - } - if (point.size() != 1) - points.push_back(point); - } - in_file.close(); -} - -void output_experiment_information(char * const file_name) -{ - std::cout << "Enter a valid experiment number. Usage: " - << file_name << " exp_no options\n"; - std::cout << "Experiment description:\n" - << "0 nbP nbL dim alpha limD mu_epsilon: " - << "Build persistence diagram on relaxed witness complex " - << "built from a point cloud on (dim-1)-dimensional sphere " - << "consisting of nbP witnesses and nbL landmarks. " - << "The maximal relaxation is alpha and the limit on simplicial complex " - << "dimension is limD.\n"; - std::cout << "1 file_name nbL alpha limD: " - << "Build persistence diagram on relaxed witness complex " - << "build from a point cloud stored in a file and nbL landmarks. " - << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n"; -} - -void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha, int limD, FT mu_epsilon = 0.1) -{ - clock_t start, end; - Simplex_tree<> simplex_tree; - - // Choose landmarks - std::vector > knn; - std::vector > distances; - start = clock(); - //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha, - limD, - knn, - distances); - end = clock(); - double time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Choice of " << nbL << " landmarks took " - << time << " s. \n"; - // Compute witness complex - start = clock(); - RelaxedWitnessComplex rw(distances, - knn, - simplex_tree, - nbL, - alpha, - limD); - end = clock(); - time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Relaxed witness complex for " << nbL << " landmarks took " - << time << " s. \n"; - std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; - // std::cout << simplex_tree << "\n"; - - // Compute the persistence diagram of the complex - simplex_tree.set_dimension(limD); - persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, true); - int p = 3; - pcoh.init_coefficients( p ); //initilizes the coefficient field for homology - start = clock(); - pcoh.compute_persistent_cohomology( alpha/1000 ); - end = clock(); - time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Persistence diagram took " - << time << " s. \n"; - pcoh.output_diagram(); - - int chi = 0; - for (auto sh: simplex_tree.complex_simplex_range()) - chi += 1-2*(simplex_tree.dimension(sh)%2); - std::cout << "Euler characteristic is " << chi << std::endl; - -} - - -int experiment0 (int argc, char * const argv[]) -{ - if (argc != 8) { - std::cerr << "Usage: " << argv[0] - << " 0 nbP nbL dim alpha limD mu_epsilon\n"; - return 0; - } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - int nbP = atoi(argv[2]); - int nbL = atoi(argv[3]); - int dim = atoi(argv[4]); - double alpha = atof(argv[5]); - int limD = atoi(argv[6]); - double mu_epsilon = atof(argv[7]); - - // Read the point file - Point_Vector point_vector; - generate_points_sphere(point_vector, nbP, dim); - std::cout << "Successfully generated " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - rw_experiment(point_vector, nbL, alpha, limD); -} - -int experiment1 (int argc, char * const argv[]) -{ - if (argc != 3) { - std::cerr << "Usage: " << argv[0] - << " 1 file_name\n"; - return 0; - } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - std::string file_name = argv[2]; - - // Read the point file - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - std::cout << "The file contains " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - bool ok = false; - int nbL, limD; - double alpha; - while (!ok) { - std::cout << "Relaxed witness complex: parameters nbL, alpha, limD.\n"; - std::cout << "Enter nbL: "; - std::cin >> nbL; - std::cout << "Enter alpha: "; - std::cin >> alpha; - std::cout << "Enter limD: "; - std::cin >> limD; - std::cout << "Start relaxed witness complex...\n"; - rw_experiment(point_vector, nbL, alpha, limD); - std::cout << "Is the result correct? [y/n]: "; - char answer; - std::cin >> answer; - switch (answer) { - case 'n': - ok = false; break; - default : - ok = true; break; - } - } - // ok = false; - // while (!ok) { - // std::cout << "Rips complex: parameters threshold, limD.\n"; - // std::cout << "Enter threshold: "; - // std::cin >> alpha; - // std::cout << "Enter limD: "; - // std::cin >> limD; - // std::cout << "Start Rips complex...\n"; - // rips_experiment(point_vector, alpha, limD); - // std::cout << "Is the result correct? [y/n]: "; - // char answer; - // std::cin >> answer; - // switch (answer) { - // case 'n': - // ok = false; break; - // default : - // ok = true; break; - // } - // } -} - -int main (int argc, char * const argv[]) -{ - if (argc == 1) { - output_experiment_information(argv[0]); - return 1; - } - switch (atoi(argv[1])) { - case 0 : - return experiment0(argc, argv); - case 1 : - return experiment1(argc, argv); - default : - output_experiment_information(argv[0]); - return 1; - } -} diff --git a/src/Witness_complex/example/strange_sphere.cpp b/src/Witness_complex/example/strange_sphere.cpp deleted file mode 100644 index 9188bd8e..00000000 --- a/src/Witness_complex/example/strange_sphere.cpp +++ /dev/null @@ -1,219 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ -#define BOOST_PARAMETER_MAX_ARITY 12 - - -#include -#include - -#include -#include -#include -#include -#include "output.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 "generators.h" - -using namespace Gudhi; -using namespace Gudhi::witness_complex; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef Strange_witness_complex< Simplex_tree<> > WitnessComplex; - -typedef CGAL::Epick_d K; -typedef K::FT FT; -typedef K::Point_d Point_d; -typedef CGAL::Search_traits< - FT, Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d> Traits_base; -typedef CGAL::Euclidean_distance Euclidean_distance; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -typedef CGAL::Orthogonal_incremental_neighbor_search> Neighbor_search; -typedef Neighbor_search::Tree Tree; -typedef Neighbor_search::Distance Distance; -typedef Neighbor_search::iterator KNS_iterator; -typedef Neighbor_search::iterator KNS_range; -typedef boost::container::flat_map Point_etiquette_map; -typedef CGAL::Kd_tree Tree2; -typedef CGAL::Fuzzy_sphere Fuzzy_sphere; -typedef std::vector Point_Vector; - -/** Function that chooses landmarks from W and place it in the kd-tree L. - * Note: nbL hould be removed if the code moves to Witness_complex - */ -void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - //std::vector landmarks; - int chosen_landmark; - //std::pair res = std::make_pair(L_i.begin(),false); - Point_d* p; - CGAL::Random rand; - for (int i = 0; i < nbL; i++) - { - // while (!res.second) - // { - do chosen_landmark = rand.get_int(0,nbP); - while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end()); - //rand++; - //std::cout << "Chose " << chosen_landmark << std::endl; - p = &W[chosen_landmark]; - //L_i.emplace(chosen_landmark,i); - // } - landmarks.push_back(*p); - landmarks_ind.push_back(chosen_landmark); - //std::cout << "Added landmark " << chosen_landmark << std::endl; - } - } - -void landmarks_to_witness_complex(Point_Vector &W, int nbL) -{ - //********************Preface: origin point - unsigned D = W[0].size(); - std::vector orig_vector; - for (unsigned i = 0; i < D; i++) - orig_vector.push_back(0); - Point_d origin(orig_vector); - //Distance dist; - //dist.transformed_distance(0,1); - //******************** Constructing a WL matrix - int nbP = W.size(); - Point_Vector landmarks; - std::vector landmarks_ind; - landmark_choice(W, nbP, nbL, landmarks, landmarks_ind); - - - STraits traits(&(landmarks[0])); - Euclidean_distance ed; - std::vector< std::vector > WL(nbP); - std::vector< std::vector< typename std::vector::iterator > > ope_limits(nbP); - Tree L(boost::counting_iterator(0), - boost::counting_iterator(nbL), - typename Tree::Splitter(), - traits); - - std::cout << "Enter (D+1) nearest landmarks\n"; - //std::cout << "Size of the tree is " << L.size() << std::endl; - for (int i = 0; i < nbP; i++) { - //std::cout << "Entered witness number " << i << std::endl; - Point_d& w = W[i]; - std::queue< typename std::vector::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark - Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter(&(landmarks[0]))); - Neighbor_search::iterator search_it = search.begin(); - - //Incremental search and filling WL - while (WL[i].size() < D) - WL[i].push_back((search_it++)->first); - } - write_wl("WL.txt", WL); //! - //******************** Constructng a witness complex - std::cout << "Entered witness complex construction\n"; - - Simplex_tree<> simplex_tree; - WitnessComplex(WL, simplex_tree, nbL, W[0].size()-1); - //std::cout << simplex_tree << "\n"; //! - write_witness_mesh(W, landmarks_ind, simplex_tree, false, true); -} - - - -/** Write a gnuplot readable file. - * Data range is a random access range of pairs (arg, value) - */ -template < typename Data_range > -void write_data(Data_range & data, std::string filename) { - std::ofstream ofs(filename, std::ofstream::out); - for (auto entry : data) - ofs << entry.first << ", " << entry.second << "\n"; - ofs.close(); -} - - - -int main(int argc, char * const argv[]) { - if (argc != 4) { - std::cerr << "Usage: " << argv[0] - << " nbP nbL dim\n"; - return 0; - } - - int nbP = atoi(argv[1]); - int nbL = atoi(argv[2]); - int dim = atoi(argv[3]); - clock_t start, end; - - // Construct the Simplex Tree - Simplex_tree<> simplex_tree; - - std::vector< std::pair > l_time; - - // Read the point file - Point_Vector point_vector; - generate_points_sphere(point_vector, nbP, dim); - std::cout << "Successfully generated " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - // Choose landmarks - start = clock(); - std::vector > knn; - - landmarks_to_witness_complex(point_vector, nbL); - - end = clock(); - double time = static_cast(end - start) / CLOCKS_PER_SEC; - std::cout << "Witness complex for " << nbL << " landmarks took " - << time << " s. \n"; - l_time.push_back(std::make_pair(nbP, time)); - //write_data(l_time, "w_time.dat"); -} diff --git a/src/Witness_complex/example/witness_complex_from_file.cpp b/src/Witness_complex/example/witness_complex_from_file.cpp index 5e9f0e81..e5859b2a 100644 --- a/src/Witness_complex/example/witness_complex_from_file.cpp +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -25,18 +25,22 @@ #include #include -#include #include #include +#include + #include #include #include #include #include +typedef CGAL::Epick_d K; +typedef typename K::Point_d Point_d; +typedef typename Gudhi::witness_complex::Witness_complex Witness_complex; typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::vector< std::vector > Point_Vector; +typedef std::vector< Point_d > Point_vector; /** * \brief Customized version of read_points @@ -44,7 +48,7 @@ typedef std::vector< std::vector > Point_Vector; * */ inline void -read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { +read_points_cust(std::string file_name, Point_vector & points) { std::ifstream in_file(file_name.c_str(), std::ios::in); if (!in_file.is_open()) { std::cerr << "Unable to open file " << file_name << std::endl; @@ -59,44 +63,44 @@ read_points_cust(std::string file_name, std::vector< std::vector< double > > & p point.push_back(x); } if (point.size() != 1) - points.push_back(point); + points.push_back(Point_d(point)); } in_file.close(); } int main(int argc, char * const argv[]) { - if (argc != 3) { + if (argc != 4) { std::cerr << "Usage: " << argv[0] - << " path_to_point_file nbL \n"; + << " path_to_point_file nbL alpha^2\n"; return 0; } std::string file_name = argv[1]; int nbL = atoi(argv[2]); + double alpha2 = atof(argv[3]); clock_t start, end; // Construct the Simplex Tree Gudhi::Simplex_tree<> simplex_tree; // Read the point file - Point_Vector point_vector, landmarks; + Point_vector point_vector, landmarks; read_points_cust(file_name, point_vector); std::cout << "Successfully read " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n"; // Choose landmarks - start = clock(); - std::vector > knn; - Gudhi::subsampling::pick_n_random_points(point_vector, 100, std::back_inserter(landmarks)); - Gudhi::witness_complex::construct_closest_landmark_table(point_vector, landmarks, knn); - end = clock(); - std::cout << "Landmark choice for " << nbL << " landmarks took " - << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; + Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); // Compute witness complex start = clock(); - Gudhi::witness_complex::witness_complex(knn, nbL, point_vector[0].size(), simplex_tree); + Witness_complex witness_complex(landmarks.begin(), + landmarks.end(), + point_vector.begin(), + point_vector.end()); + witness_complex.create_complex(simplex_tree, alpha2); end = clock(); std::cout << "Witness complex took " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; } diff --git a/src/Witness_complex/include/gudhi/A0_complex.h b/src/Witness_complex/include/gudhi/A0_complex.h deleted file mode 100644 index 2018e1c8..00000000 --- a/src/Witness_complex/include/gudhi/A0_complex.h +++ /dev/null @@ -1,337 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef A0_COMPLEX_H_ -#define A0_COMPLEX_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" -#include -#include -#include -#include -#include -#include -#include -#include - -// Needed for nearest neighbours -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** \addtogroup simplex_tree - * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. - */ -template< class Simplicial_complex > -class A0_complex { - -private: - struct Active_witness { - int witness_id; - int landmark_id; - - Active_witness(int witness_id_, int landmark_id_) - : witness_id(witness_id_), - landmark_id(landmark_id_) { } - }; - -private: - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - typedef typename Simplicial_complex::Filtration_value FT; - - typedef std::vector< double > Point_t; - typedef std::vector< Point_t > Point_Vector; - - // typedef typename Simplicial_complex::Filtration_value Filtration_value; - typedef std::vector< Vertex_handle > typeVectorVertex; - typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; - typedef std::pair< Simplex_handle, bool > typePairSimplexBool; - - typedef int Witness_id; - typedef int Landmark_id; - typedef std::list< Vertex_handle > ActiveWitnessList; - - private: - int nbL; // Number of landmarks - Simplicial_complex& sc; // Simplicial complex - - public: - /** @name Constructor - */ - - //@{ - - /** - * \brief Iterative construction of the relaxed witness complex. - * \details The witness complex is written in sc_ basing on a matrix knn - * of k nearest neighbours of the form {witnesses}x{landmarks} and - * and a matrix distances of distances to these landmarks from witnesses. - * The parameter alpha defines relaxation and - * limD defines the - * - * The line lengths in one given matrix can differ, - * however both matrices have the same corresponding line lengths. - * - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges. - * - * Constructor takes into account at most (dim+1) - * first landmarks from each landmark range to construct simplices. - * - * Landmarks are supposed to be in [0,nbL_-1] - */ - - template< typename KNearestNeighbours > - A0_complex(std::vector< std::vector > const & distances, - KNearestNeighbours const & knn, - Simplicial_complex & sc_, - int nbL_, - double alpha2, - unsigned limD) : nbL(nbL_), sc(sc_) { - int nbW = knn.size(); - typeVectorVertex vv; - //int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - ActiveWitnessList active_w;// = new ActiveWitnessList(); - for (int i = 0; i != nbL; ++i) { - // initial fill of 0-dimensional simplices - // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore - //counter++; - vv = {i}; - sc.insert_simplex(vv, Filtration_value(0.0)); - /* TODO Error if not inserted : normally no need here though*/ - } - for (int i=0; i != nbW; ++i) { - // int i_end = limD+1; - // if (knn[i].size() < limD+1) - // i_end = knn[i].size(); - // double dist_wL = *(distances[i].begin()); - // while (distances[i][i_end] > dist_wL + alpha2) - // i_end--; - // add_all_witnessed_faces(distances[i].begin(), - // knn[i].begin(), - // knn[i].begin() + i_end + 1); - unsigned j_end = 0; - while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) { - std::vector simplex; - for (unsigned j = 0; j <= j_end; ++j) - simplex.push_back(knn[i][j]); - assert(distances[i][j_end] - distances[i][0] >= 0); - sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]); - j_end++; - } - } - sc.set_dimension(limD); - } - - //@} - -private: - /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness - * Iterator is needed to know until how far we can take landmarks to form simplexes - * simplex is the prefix of the simplexes to insert - * The output value indicates if the witness rests active or not - */ - void add_all_witnessed_faces(std::vector::const_iterator curr_d, - std::vector::const_iterator curr_l, - std::vector::const_iterator end) - { - std::vector simplex; - std::vector::const_iterator l_end = curr_l; - for (; l_end != end; ++l_end) { - std::vector::const_iterator l_it = curr_l; - std::vector::const_iterator d_it = curr_d; - simplex = {}; - for (; l_it != l_end; ++l_it, ++d_it) - simplex.push_back(*l_it); - sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d); - } - } - - /** \brief Check if the facets of the k-dimensional simplex witnessed - * by witness witness_id are already in the complex. - * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id - */ - bool all_faces_in(std::vector& simplex, double* filtration_value) - { - std::vector< int > facet; - for (std::vector::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) - { - facet.clear(); - for (std::vector::iterator it = simplex.begin(); it != simplex.end(); ++it) - if (it != not_it) - facet.push_back(*it); - Simplex_handle facet_sh = sc.find(facet); - if (facet_sh == sc.null_simplex()) - return false; - else if (sc.filtration(facet_sh) > *filtration_value) - *filtration_value = sc.filtration(facet_sh); - } - return true; - } - - bool is_face(Simplex_handle face, Simplex_handle coface) - { - // vertex range is sorted in decreasing order - auto fvr = sc.simplex_vertex_range(face); - auto cfvr = sc.simplex_vertex_range(coface); - auto fv_it = fvr.begin(); - auto cfv_it = cfvr.begin(); - while (fv_it != fvr.end() && cfv_it != cfvr.end()) { - if (*fv_it < *cfv_it) - ++cfv_it; - else if (*fv_it == *cfv_it) { - ++cfv_it; - ++fv_it; - } - else - return false; - - } - return (fv_it == fvr.end()); - } - - // void erase_simplex(Simplex_handle sh) - // { - // auto siblings = sc.self_siblings(sh); - // auto oncles = siblings->oncles(); - // int prev_vertex = siblings->parent(); - // siblings->members().erase(sh->first); - // if (siblings->members().empty()) { - // typename typedef Simplicial_complex::Siblings Siblings; - // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); - // assert(!sc.has_children(oncles->members().find(prev_vertex))); - // //delete siblings; - // } - - // } - - void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) - { - erase_simplex(coface_sh); - erase_simplex(face_sh); - } - -public: - // void collapse(std::vector& simplices) - // { - // // Get a vector of simplex handles ordered by filtration value - // std::cout << sc << std::endl; - // //std::vector simplices; - // for (Simplex_handle sh: sc.filtration_simplex_range()) - // simplices.push_back(sh); - // // std::sort(simplices.begin(), - // // simplices.end(), - // // [&](Simplex_handle sh1, Simplex_handle sh2) - // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); - // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); - // // Double iteration - // auto face_it = simplices.rbegin(); - // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { - // int coface_count = 0; - // auto reduced_coface = simplices.rbegin(); - // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) - // if (face_it != coface_it && is_face(*face_it, *coface_it)) { - // coface_count++; - // if (coface_count == 1) - // reduced_coface = coface_it; - // else - // break; - // } - // if (coface_count == 1) { - // std::cout << "Erase ( "; - // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) - // std::cout << v << " "; - - // simplices.erase(--(reduced_coface.base())); - // //elementary_collapse(*face_it, *reduced_coface); - // std::cout << ") and then ( "; - // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) - // std::cout << v << " "; - // std::cout << ")\n"; - // simplices.erase(--((face_it++).base())); - // //face_it = simplices.rbegin(); - // //std::cout << "Size of vector: " << simplices.size() << "\n"; - // } - // else - // face_it++; - // } - // sc.initialize_filtration(); - // //std::cout << sc << std::endl; - // } - - template - void collapse(Dim_lists& dim_lists) - { - dim_lists.collapse(); - } - -private: - - /** Collapse recursively boundary faces of the given simplex - * if its filtration is bigger than alpha_lim. - */ - void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) - { - for (Simplex_handle face_it : sc.boundary_simplex_range()) { - - } - - } - -}; //class Relaxed_witness_complex - -} // namespace witness_complex - -} // namespace Gudhi - -#endif diff --git a/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h b/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h deleted file mode 100644 index ef711c34..00000000 --- a/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h +++ /dev/null @@ -1,90 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_ -#define CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_ - -#include - -#include // for priority_queue<> -#include // for pair<> -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** - * \ingroup witness_complex - * \brief Construct the closest landmark tables for all witnesses. - * \details Output a table 'knn', each line of which represents a witness and - * consists of landmarks sorted by - * euclidean distance from the corresponding witness. - * - * The type WitnessContainer is a random access range and - * the type LandmarkContainer is a range. - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges and - * Vertex_handle is the label type of a vertex in a simplicial complex. - * Closest_landmark_range needs to have push_back operation. - */ - - template - void construct_closest_landmark_table(WitnessContainer const &points, - LandmarkContainer const &landmarks, - KNearestNeighbours &knn) { - int nbP = boost::size(points); - assert(nbP >= boost::size(landmarks)); - - int dim = boost::size(*std::begin(points)); - typedef std::pair dist_i; - typedef bool (*comp)(dist_i, dist_i); - knn = KNearestNeighbours(nbP); - for (int points_i = 0; points_i < nbP; points_i++) { - std::priority_queue, comp> l_heap([](dist_i j1, dist_i j2) { - return j1.first > j2.first; - }); - typename LandmarkContainer::const_iterator landmarks_it; - int landmarks_i = 0; - for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end(); - ++landmarks_it, landmarks_i++) { - dist_i dist = std::make_pair(euclidean_distance(points[points_i], *landmarks_it), landmarks_i); - l_heap.push(dist); - } - for (int i = 0; i < dim + 1; i++) { - dist_i dist = l_heap.top(); - knn[points_i].push_back(dist.second); - l_heap.pop(); - } - } - } - -} // namespace witness_complex - -} // namespace Gudhi - -#endif // CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_ diff --git a/src/Witness_complex/include/gudhi/Dim_list_iterator.h b/src/Witness_complex/include/gudhi/Dim_list_iterator.h deleted file mode 100644 index 3f23e8c9..00000000 --- a/src/Witness_complex/include/gudhi/Dim_list_iterator.h +++ /dev/null @@ -1,155 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef DIM_LISTS_ITERATOR_H_ -#define DIM_LISTS_ITERATOR_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** \addtogroup simplex_tree - * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. - */ -template< class Simplex_handle > -class Dim_lists_iterator { - -private: - - typedef typename std::list::iterator List_iterator; - typedef typename std::vector>::reverse_iterator Vector_iterator; - typedef typename Gudhi::witness_complex::Dim_lists_iterator Iterator; - - - List_iterator sh_; - Vector_iterator curr_list_; - typename std::vector>& table_; - -public: - - Dim_lists_iterator(List_iterator sh, - Vector_iterator curr_list, - typename std::vector>& table) - : sh_(sh), curr_list_(curr_list), table_(table) - { - } - - Simplex_handle operator*() - { - return *sh_; - } - - Iterator operator++() - { - increment(); - return (*this); - } - - Iterator operator++(int) - { - Iterator prev_it(sh_, curr_list_, table_); - increment(); - return prev_it; - } - - Iterator dim_begin() - { - return Iterator(curr_list_->begin(), curr_list_, table_); - } - - Iterator dim_end() - { - return Iterator(curr_list_->end(), curr_list_, table_); - } - - Iterator dimp1_begin() - { - return Iterator((curr_list_-1)->begin(), curr_list_-1, table_); - } - - Iterator dimp1_end() - { - return Iterator((curr_list_-1)->end(), curr_list_-1, table_); - } - - bool operator==(const Iterator& it2) const - { - return (sh_ == it2.sh_); - } - - bool operator!=(const Iterator& it2) const - { - return (sh_ != it2.sh_); - } - - void remove_incr() - { - - } - -private: - - void increment() - { - if (++sh_ == curr_list_->end()) - if (++curr_list_ != table_.rend()) - sh_ = curr_list_->begin(); - // The iterator of the end of the table is the end of the last list - } - - -}; //class Dim_lists_iterator - -} // namespace witness_complex - -} // namespace Gudhi - -#endif diff --git a/src/Witness_complex/include/gudhi/Dim_lists.h b/src/Witness_complex/include/gudhi/Dim_lists.h deleted file mode 100644 index 1d1b03c3..00000000 --- a/src/Witness_complex/include/gudhi/Dim_lists.h +++ /dev/null @@ -1,195 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef DIM_LISTS_H_ -#define DIM_LISTS_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** \addtogroup simplex_tree - * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. - */ -template< class Simplicial_complex > -class Dim_lists { - -private: - - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Gudhi::witness_complex::Dim_lists_iterator Iterator; - - std::vector> table_; - Simplicial_complex& sc_; - -public: - - Dim_lists(Simplicial_complex & sc, int dim, double alpha_max = 100) - : sc_(sc) - { - table_ = std::vector>(dim+1); - for (auto sh: sc.filtration_simplex_range()) { - if (sc_.filtration(sh) < alpha_max) - table_[sc.dimension(sh)].push_front(sh); - } - auto t_it = table_.rbegin(); - while (t_it->empty()) { - t_it++; - table_.pop_back(); - } - } - - Iterator begin() - { - return Iterator(table_.rbegin()->begin(), table_.rbegin(), table_); - } - - Iterator end() - { - return Iterator(table_[0].end(), table_.rend(), table_); - } - - unsigned size() - { - unsigned curr_size = 0; - for (auto l: table_) - curr_size += l.size(); - return curr_size; - } - - bool is_face(Simplex_handle face, Simplex_handle coface) - { - // vertex range is sorted in decreasing order - auto fvr = sc_.simplex_vertex_range(face); - auto cfvr = sc_.simplex_vertex_range(coface); - auto fv_it = fvr.begin(); - auto cfv_it = cfvr.begin(); - while (fv_it != fvr.end() && cfv_it != cfvr.end()) { - if (*fv_it < *cfv_it) - ++cfv_it; - else if (*fv_it == *cfv_it) { - ++cfv_it; - ++fv_it; - } - else - return false; - - } - return (fv_it == fvr.end()); - } - - void output_simplices() { - std::cout << "Size of vector: " << size() << std::endl; - for (auto line_it = table_.rbegin(); line_it != table_.rend(); ++line_it) - for (auto sh: *line_it) { - std::cout << sc_.dimension(sh) << " "; - for (auto v : sc_.simplex_vertex_range(sh)) - std::cout << v << " "; - std::cout << sc_.filtration(sh) << "\n"; - } - } - - void collapse() - { - auto coface_list_it = table_.rbegin(); - auto face_list_it = table_.rbegin()+1; - for ( ; - face_list_it != table_.rend(); - ++face_list_it) { - auto face_it = face_list_it->begin(); - while (face_it != face_list_it->end() && sc_.filtration(*face_it) != 0) { - int coface_count = 0; - auto reduced_coface = coface_list_it->begin(); - for (auto coface_it = coface_list_it->begin(); coface_it != coface_list_it->end() && sc_.filtration(*coface_it) != 0; ++coface_it) - if (is_face(*face_it, *coface_it)) { - coface_count++; - if (coface_count == 1) - reduced_coface = coface_it; - else - break; - } - if (coface_count == 1) { - /* - std::cout << "Erase ( "; - for (auto v: sc_.simplex_vertex_range(*reduced_coface)) - std::cout << v << " "; - */ - coface_list_it->erase(reduced_coface); - /* - std::cout << ") and then ( "; - for (auto v: sc_.simplex_vertex_range(*face_it)) - std::cout << v << " "; - std::cout << ")\n"; - */ - face_list_it->erase(face_it); - face_it = face_list_it->begin(); - } - else - face_it++; - } - if ((coface_list_it++)->empty()) - table_.pop_back(); - } - } - - bool is_pseudomanifold() - { - - return true; - } - -}; //class Dim_lists - -} // namespace witness_complex - -} // namespace Gudhi - -#endif diff --git a/src/Witness_complex/include/gudhi/Good_links.h b/src/Witness_complex/include/gudhi/Good_links.h deleted file mode 100644 index c5e993e5..00000000 --- a/src/Witness_complex/include/gudhi/Good_links.h +++ /dev/null @@ -1,449 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef GOOD_LINKS_H_ -#define GOOD_LINKS_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - -template < typename Simplicial_complex > -class Good_links { - - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - typedef std::vector Vertex_vector; - - // graph is an adjacency list - typedef typename boost::adjacency_list Adj_graph; - // map that gives to a certain simplex its node in graph and its dimension - //typedef std::pair Reference; - typedef boost::graph_traits::vertex_descriptor Vertex_t; - typedef boost::graph_traits::edge_descriptor Edge_t; - typedef boost::graph_traits::adjacency_iterator Adj_it; - typedef std::pair Out_edge_it; - - typedef boost::container::flat_map Graph_map; - typedef boost::container::flat_map Inv_graph_map; - -public: - Good_links(Simplicial_complex& sc): sc_(sc) - { - int dim = 0; - for (auto sh: sc_.complex_simplex_range()) - if (sc_.dimension(sh) > dim) - dim = sc_.dimension(sh); - dimension = dim; - count_good = Vertex_vector(dim); - count_bad = Vertex_vector(dim); - } - -private: - - Simplicial_complex& sc_; - unsigned dimension; - std::vector count_good; - std::vector count_bad; - - void add_vertices_to_link_graph(Vertex_vector& star_vertices, - typename Vertex_vector::iterator curr_v, - Adj_graph& adj_graph, - Graph_map& d_map, - Graph_map& f_map, - Vertex_vector& curr_simplex, - int curr_d, - int link_dimension) - { - Simplex_handle sh; - Vertex_t vert; - typename Vertex_vector::iterator it; - //std::pair resPair; - //typename Graph_map::iterator resPair; - //Add vertices - //std::cout << "Entered add vertices\n"; - for (it = curr_v; it != star_vertices.end(); ++it) { - curr_simplex.push_back(*it); //push next vertex in question - curr_simplex.push_back(star_vertices[0]); //push the center of the star - /* - std::cout << "Searching for "; - print_vector(curr_simplex); - std::cout << " curr_dim " << curr_d << " d " << dimension << ""; - */ - Vertex_vector curr_simplex_copy(curr_simplex); - sh = sc_.find(curr_simplex_copy); //a simplex of the star - curr_simplex.pop_back(); //pop the center of the star - curr_simplex_copy = Vertex_vector(curr_simplex); - if (sh != sc_.null_simplex()) { - //std::cout << " added\n"; - if (curr_d == link_dimension) { - sh = sc_.find(curr_simplex_copy); //a simplex of the link - assert(sh != sc_.null_simplex()); //ASSERT! - vert = boost::add_vertex(adj_graph); - d_map.emplace(sh,vert); - } - else { - if (curr_d == link_dimension-1) { - sh = sc_.find(curr_simplex_copy); //a simplex of the link - assert(sh != sc_.null_simplex()); - vert = boost::add_vertex(adj_graph); - f_map.emplace(sh,vert); - } - - //delete (&curr_simplex_copy); //Just so it doesn't stack - add_vertices_to_link_graph(star_vertices, - it+1, - adj_graph, - d_map, - f_map, - curr_simplex, - curr_d+1, link_dimension); - } - } - /* - else - std::cout << "\n"; - */ - curr_simplex.pop_back(); //pop the vertex in question - } - } - - void add_edges_to_link_graph(Adj_graph& adj_graph, - Graph_map& d_map, - Graph_map& f_map) - { - Simplex_handle sh; - // Add edges - //std::cout << "Entered add edges:\n"; - typename Graph_map::iterator map_it; - for (auto d_map_pair : d_map) { - //std::cout << "*"; - sh = d_map_pair.first; - Vertex_t d_vert = d_map_pair.second; - for (auto facet_sh : sc_.boundary_simplex_range(sh)) - //for (auto f_map_it : f_map) - { - //std::cout << "'"; - map_it = f_map.find(facet_sh); - //We must have all the facets in the graph at this point - assert(map_it != f_map.end()); - Vertex_t f_vert = map_it->second; - //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n"; - boost::add_edge(d_vert,f_vert,adj_graph); - } - } - } - - - - /* \brief Verifies if the simplices formed by vertices given by link_vertices - * form a pseudomanifold. - * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional - * faces and edges represent adjacency between them. - */ - bool link_is_pseudomanifold(Vertex_vector& star_vertices, - int dimension) - { - Adj_graph adj_graph; - Graph_map d_map, f_map; - // d_map = map for d-dimensional simplices - // f_map = map for its facets - Vertex_vector init_vector = {}; - add_vertices_to_link_graph(star_vertices, - star_vertices.begin()+1, - adj_graph, - d_map, - f_map, - init_vector, - 0, dimension); - //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; - //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; - add_edges_to_link_graph(adj_graph, - d_map, - f_map); - for (auto f_map_it : f_map) { - //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; - if (boost::out_degree(f_map_it.second, adj_graph) != 2){ - count_bad[dimension]++; - return false; - } - } - // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices - // What is left is to check the connexity - //std::vector components(boost::num_vertices(adj_graph)); - return true; //Forget the connexity - //return (boost::connected_components(adj_graph, &components[0]) == 1); - } - - int star_dim(Vertex_vector& star_vertices, - typename Vertex_vector::iterator curr_v, - int curr_d, - Vertex_vector& curr_simplex, - typename std::vector< int >::iterator curr_dc) - { - //std::cout << "Entered star_dim for " << *(curr_v-1) << "\n"; - Simplex_handle sh; - int final_d = curr_d; - typename Vertex_vector::iterator it; - typename Vertex_vector::iterator dc_it; - //std::cout << "Current vertex is " << - for (it = curr_v, dc_it = curr_dc; it != star_vertices.end(); ++it, ++dc_it) - { - curr_simplex.push_back(*it); - Vertex_vector curr_simplex_copy(curr_simplex); - /* - std::cout << "Searching for "; - print_vector(curr_simplex); - std::cout << " curr_dim " << curr_d << " final_dim " << final_d; - */ - sh = sc_.find(curr_simplex_copy); //Need a copy because find sorts the vector and I want star center to be the first - if (sh != sc_.null_simplex()) - { - //std::cout << " -> " << *it << "\n"; - int d = star_dim(star_vertices, - it+1, - curr_d+1, - curr_simplex, - dc_it); - if (d >= final_d) - { - final_d = d; - //std::cout << d << " "; - //print_vector(curr_simplex); - //std::cout << std::endl; - } - if (d >= *dc_it) - *dc_it = d; - } - /* - else - std::cout << "\n"; - */ - curr_simplex.pop_back(); - } - return final_d; - } - -public: - - /** \brief Returns true if the link is good - */ - bool has_good_link(Vertex_handle v) - { - typedef Vertex_vector typeVectorVertex; - typeVectorVertex star_vertices; - // Fill star_vertices - star_vertices.push_back(v); - for (auto u: sc_.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u != v && sc_.find(edge) != sc_.null_simplex()) - star_vertices.push_back(u); - } - // Find the dimension - typeVectorVertex init_simplex = {star_vertices[0]}; - bool is_pure = true; - std::vector dim_coface(star_vertices.size(), 1); - int d = star_dim(star_vertices, - star_vertices.begin()+1, - 0, - init_simplex, - dim_coface.begin()+1) - 1; //link_dim = star_dim - 1 - - assert(init_simplex.size() == 1); - // if (!is_pure) - // std::cout << "Found an impure star around " << v << "\n"; - for (int dc: dim_coface) - is_pure = (dc == dim_coface[0]); - /* - if (d == count_good.size()) - { - std::cout << "Found a star of dimension " << (d+1) << " around " << v << "\nThe star is "; - print_vector(star_vertices); std::cout << std::endl; - } - */ - //if (d == -1) count_bad[0]++; - bool b= (is_pure && link_is_pseudomanifold(star_vertices,d)); - if (d != -1) {if (b) count_good[d]++; else count_bad[d]++;} - if (!is_pure) count_bad[0]++; - return (d != -1 && b && is_pure); - } - -private: - void add_max_simplices_to_graph(Vertex_vector& star_vertices, - typename Vertex_vector::iterator curr_v, - Adj_graph& adj_graph, - Graph_map& d_map, - Graph_map& f_map, - Inv_graph_map& inv_d_map, - Vertex_vector& curr_simplex, - int curr_d, - int link_dimension) - { - Simplex_handle sh; - Vertex_t vert; - typename Vertex_vector::iterator it; - //std::pair resPair; - //typename Graph_map::iterator resPair; - //Add vertices - //std::cout << "Entered add vertices\n"; - for (it = curr_v; it != star_vertices.end(); ++it) { - curr_simplex.push_back(*it); //push next vertex in question - //curr_simplex.push_back(star_vertices[0]); //push the center of the star - /* - std::cout << "Searching for "; - print_vector(curr_simplex); - std::cout << " curr_dim " << curr_d << " d " << dimension << ""; - */ - Vertex_vector curr_simplex_copy(curr_simplex); - sh = sc_.find(curr_simplex_copy); //a simplex of the star - //curr_simplex.pop_back(); //pop the center of the star - curr_simplex_copy = Vertex_vector(curr_simplex); - if (sh != sc_.null_simplex()) { - //std::cout << " added\n"; - if (curr_d == link_dimension) { - sh = sc_.find(curr_simplex_copy); //a simplex of the link - assert(sh != sc_.null_simplex()); //ASSERT! - vert = boost::add_vertex(adj_graph); - d_map.emplace(sh,vert); - inv_d_map.emplace(vert,sh); - } - else { - if (curr_d == link_dimension-1) { - sh = sc_.find(curr_simplex_copy); //a simplex of the link - assert(sh != sc_.null_simplex()); - vert = boost::add_vertex(adj_graph); - f_map.emplace(sh,vert); - } - //delete (&curr_simplex_copy); //Just so it doesn't stack - add_max_simplices_to_graph(star_vertices, - it+1, - adj_graph, - d_map, - f_map, - inv_d_map, - curr_simplex, - curr_d+1, - link_dimension); - } - } - /* - else - std::cout << "\n"; - */ - curr_simplex.pop_back(); //pop the vertex in question - } - } - -public: - bool complex_is_pseudomanifold() - { - Adj_graph adj_graph; - Graph_map d_map, f_map; - // d_map = map for d-dimensional simplices - // f_map = map for its facets - Inv_graph_map inv_d_map; - Vertex_vector init_vector = {}; - std::vector star_vertices; - for (int v: sc_.complex_vertex_range()) - star_vertices.push_back(v); - add_max_simplices_to_graph(star_vertices, - star_vertices.begin(), - adj_graph, - d_map, - f_map, - inv_d_map, - init_vector, - 0, dimension); - //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; - //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; - add_edges_to_link_graph(adj_graph, - d_map, - f_map); - for (auto f_map_it : f_map) { - //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; - if (boost::out_degree(f_map_it.second, adj_graph) != 2) { - // if (boost::out_degree(f_map_it.second, adj_graph) >= 3) { - // // std::cout << "This simplex has 3+ cofaces: "; - // // for(auto v : sc_.simplex_vertex_range(f_map_it.first)) - // // std::cout << v << " "; - // // std::cout << std::endl; - // Adj_it ai, ai_end; - // for (std::tie(ai, ai_end) = boost::adjacent_vertices(f_map_it.second, adj_graph); ai != ai_end; ++ai) { - // auto it = inv_d_map.find(*ai); - // assert (it != inv_d_map.end()); - // typename Simplicial_complex::Simplex_handle sh = it->second; - // for(auto v : sc_.simplex_vertex_range(sh)) - // std::cout << v << " "; - // std::cout << std::endl; - // } - // } - count_bad[dimension]++; - return false; - } - } - // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices - // What is left is to check the connexity - //std::vector components(boost::num_vertices(adj_graph)); - return true; //Forget the connexity - //return (boost::connected_components(adj_graph, &components[0]) == 1); - } - - int number_good_links(int dim) - { - return count_good[dim]; - } - - int number_bad_links(int dim) - { - return count_bad[dim]; - } - -}; - -} // namespace Gudhi - -#endif diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h deleted file mode 100644 index df93155b..00000000 --- a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h +++ /dev/null @@ -1,105 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ -#define LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ - -#include - -#include // for numeric_limits<> -#include -#include // for sort -#include - -namespace Gudhi { - -namespace witness_complex { - - typedef std::vector typeVectorVertex; - - /** - * \ingroup witness_complex - * \brief Landmark choice strategy by iteratively adding the furthest witness from the - * current landmark set as the new landmark. - * \details It chooses nbL landmarks from a random access range `points` and - * writes {witness}*{closest landmarks} matrix in `knn`. - * - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges - * - */ - - template - void landmark_choice_by_furthest_point(Point_random_access_range const &points, - int nbL, - KNearestNeighbours &knn) { - int nb_points = boost::size(points); - assert(nb_points >= nbL); - // distance matrix witness x landmarks - std::vector> wit_land_dist(nb_points, std::vector()); - // landmark list - typeVectorVertex chosen_landmarks; - - knn = KNearestNeighbours(nb_points, std::vector()); - int current_number_of_landmarks = 0; // counter for landmarks - double curr_max_dist = 0; // used for defining the furhest point from L - const double infty = std::numeric_limits::infinity(); // infinity (see next entry) - std::vector< double > dist_to_L(nb_points, infty); // vector of current distances to L from points - - // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety - // or better yet std::uniform_int_distribution - int rand_int = rand() % nb_points; - int curr_max_w = rand_int; // For testing purposes a pseudo-random number is used here - - for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { - // curr_max_w at this point is the next landmark - chosen_landmarks.push_back(curr_max_w); - unsigned i = 0; - for (auto& p : points) { - double curr_dist = euclidean_distance(p, *(std::begin(points) + chosen_landmarks[current_number_of_landmarks])); - wit_land_dist[i].push_back(curr_dist); - knn[i].push_back(current_number_of_landmarks); - if (curr_dist < dist_to_L[i]) - dist_to_L[i] = curr_dist; - ++i; - } - curr_max_dist = 0; - for (i = 0; i < dist_to_L.size(); i++) - if (dist_to_L[i] > curr_max_dist) { - curr_max_dist = dist_to_L[i]; - curr_max_w = i; - } - } - for (int i = 0; i < nb_points; ++i) - std::sort(std::begin(knn[i]), - std::end(knn[i]), - [&wit_land_dist, i](int a, int b) { - return wit_land_dist[i][a] < wit_land_dist[i][b]; }); - } - -} // namespace witness_complex - -} // namespace Gudhi - -#endif // LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h deleted file mode 100644 index ebf6aad1..00000000 --- a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h +++ /dev/null @@ -1,96 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef LANDMARK_CHOICE_BY_RANDOM_POINT_H_ -#define LANDMARK_CHOICE_BY_RANDOM_POINT_H_ - -#include - -#include // for priority_queue<> -#include // for pair<> -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** - * \ingroup witness_complex - * \brief Landmark choice strategy by taking random vertices for landmarks. - * \details It chooses nbL distinct landmarks from a random access range `points` - * and outputs a matrix {witness}*{closest landmarks} in knn. - * - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges and - * Vertex_handle is the label type of a vertex in a simplicial complex. - * Closest_landmark_range needs to have push_back operation. - */ - - template - void landmark_choice_by_random_point(Point_random_access_range const &points, - int nbL, - KNearestNeighbours &knn) { - int nbP = boost::size(points); - assert(nbP >= nbL); - std::set landmarks; - int current_number_of_landmarks = 0; // counter for landmarks - - // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety - int chosen_landmark = rand() % nbP; - for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { - while (landmarks.find(chosen_landmark) != landmarks.end()) - chosen_landmark = rand() % nbP; - landmarks.insert(chosen_landmark); - } - - int dim = boost::size(*std::begin(points)); - typedef std::pair dist_i; - typedef bool (*comp)(dist_i, dist_i); - knn = KNearestNeighbours(nbP); - for (int points_i = 0; points_i < nbP; points_i++) { - std::priority_queue, comp> l_heap([](dist_i j1, dist_i j2) { - return j1.first > j2.first; - }); - std::set::iterator landmarks_it; - int landmarks_i = 0; - for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end(); - ++landmarks_it, landmarks_i++) { - dist_i dist = std::make_pair(euclidean_distance(points[points_i], points[*landmarks_it]), landmarks_i); - l_heap.push(dist); - } - for (int i = 0; i < dim + 1; i++) { - dist_i dist = l_heap.top(); - knn[points_i].push_back(dist.second); - l_heap.pop(); - } - } - } - -} // namespace witness_complex - -} // namespace Gudhi - -#endif // LANDMARK_CHOICE_BY_RANDOM_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h deleted file mode 100644 index 2971149a..00000000 --- a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h +++ /dev/null @@ -1,379 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef RELAXED_WITNESS_COMPLEX_H_ -#define RELAXED_WITNESS_COMPLEX_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" -#include -#include -#include -#include -#include -#include -#include -#include - -// Needed for nearest neighbours -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** \addtogroup simplex_tree - * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. - */ -template< class Simplicial_complex > -class Relaxed_witness_complex { - -private: - struct Active_witness { - int witness_id; - int landmark_id; - - Active_witness(int witness_id_, int landmark_id_) - : witness_id(witness_id_), - landmark_id(landmark_id_) { } - }; - -private: - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - typedef typename Simplicial_complex::Filtration_value FT; - - typedef std::vector< double > Point_t; - typedef std::vector< Point_t > Point_Vector; - - // typedef typename Simplicial_complex::Filtration_value Filtration_value; - typedef std::vector< Vertex_handle > typeVectorVertex; - typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; - typedef std::pair< Simplex_handle, bool > typePairSimplexBool; - - typedef int Witness_id; - typedef int Landmark_id; - typedef std::list< Vertex_handle > ActiveWitnessList; - - private: - int nbL; // Number of landmarks - Simplicial_complex& sc; // Simplicial complex - - public: - /** @name Constructor - */ - - //@{ - - /** - * \brief Iterative construction of the relaxed witness complex. - * \details The witness complex is written in sc_ basing on a matrix knn - * of k nearest neighbours of the form {witnesses}x{landmarks} and - * and a matrix distances of distances to these landmarks from witnesses. - * The parameter alpha defines relaxation and - * limD defines the - * - * The line lengths in one given matrix can differ, - * however both matrices have the same corresponding line lengths. - * - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges. - * - * Constructor takes into account at most (dim+1) - * first landmarks from each landmark range to construct simplices. - * - * Landmarks are supposed to be in [0,nbL_-1] - */ - - template< typename KNearestNeighbours > - Relaxed_witness_complex(std::vector< std::vector > const & distances, - KNearestNeighbours const & knn, - Simplicial_complex & sc_, - int nbL_, - double alpha, - int limD) : nbL(nbL_), sc(sc_) { - int nbW = knn.size(); - typeVectorVertex vv; - //int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - ActiveWitnessList active_w;// = new ActiveWitnessList(); - for (int i = 0; i != nbL; ++i) { - // initial fill of 0-dimensional simplices - // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore - //counter++; - vv = {i}; - sc.insert_simplex(vv, Filtration_value(0.0)); - /* TODO Error if not inserted : normally no need here though*/ - } - int k = 1; /* current dimension in iterative construction */ - for (int i=0; i != nbW; ++i) - active_w.push_back(i); - while (!active_w.empty() && k <= limD && k < nbL ) - { - typename ActiveWitnessList::iterator aw_it = active_w.begin(); - while (aw_it != active_w.end()) - { - std::vector simplex; - bool ok = add_all_faces_of_dimension(k, - alpha, - std::numeric_limits::infinity(), - distances[*aw_it].begin(), - knn[*aw_it].begin(), - simplex, - knn[*aw_it].end()); - if (!ok) - active_w.erase(aw_it++); //First increase the iterator and then erase the previous element - else - aw_it++; - } - k++; - } - sc.set_dimension(limD); - } - - //@} - -private: - /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness - * Iterator is needed to know until how far we can take landmarks to form simplexes - * simplex is the prefix of the simplexes to insert - * The output value indicates if the witness rests active or not - */ - bool add_all_faces_of_dimension(int dim, - double alpha, - double norelax_dist, - std::vector::const_iterator curr_d, - std::vector::const_iterator curr_l, - std::vector& simplex, - std::vector::const_iterator end) - { - if (curr_l == end) - return false; - bool will_be_active = false; - std::vector::const_iterator l_it = curr_l; - std::vector::const_iterator d_it = curr_d; - if (dim > 0) - for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { - simplex.push_back(*l_it); - if (sc.find(simplex) != sc.null_simplex()) - will_be_active = add_all_faces_of_dimension(dim-1, - alpha, - norelax_dist, - d_it+1, - l_it+1, - simplex, - end) || will_be_active; - simplex.pop_back(); - // If norelax_dist is infinity, change to first omitted distance - if (*d_it <= norelax_dist) - norelax_dist = *d_it; - will_be_active = add_all_faces_of_dimension(dim-1, - alpha, - norelax_dist, - d_it+1, - l_it+1, - simplex, - end) || will_be_active; - } - else if (dim == 0) - for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { - simplex.push_back(*l_it); - double filtration_value = 0; - // if norelax_dist is infinite, relaxation is 0. - if (*d_it > norelax_dist) - filtration_value = *d_it - norelax_dist; - if (all_faces_in(simplex, &filtration_value)) { - will_be_active = true; - sc.insert_simplex(simplex, filtration_value); - } - simplex.pop_back(); - // If norelax_dist is infinity, change to first omitted distance - if (*d_it < norelax_dist) - norelax_dist = *d_it; - } - return will_be_active; - } - - /** \brief Check if the facets of the k-dimensional simplex witnessed - * by witness witness_id are already in the complex. - * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id - */ - bool all_faces_in(std::vector& simplex, double* filtration_value) - { - std::vector< int > facet; - for (std::vector::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) - { - facet.clear(); - for (std::vector::iterator it = simplex.begin(); it != simplex.end(); ++it) - if (it != not_it) - facet.push_back(*it); - Simplex_handle facet_sh = sc.find(facet); - if (facet_sh == sc.null_simplex()) - return false; - else if (sc.filtration(facet_sh) > *filtration_value) - *filtration_value = sc.filtration(facet_sh); - } - return true; - } - - bool is_face(Simplex_handle face, Simplex_handle coface) - { - // vertex range is sorted in decreasing order - auto fvr = sc.simplex_vertex_range(face); - auto cfvr = sc.simplex_vertex_range(coface); - auto fv_it = fvr.begin(); - auto cfv_it = cfvr.begin(); - while (fv_it != fvr.end() && cfv_it != cfvr.end()) { - if (*fv_it < *cfv_it) - ++cfv_it; - else if (*fv_it == *cfv_it) { - ++cfv_it; - ++fv_it; - } - else - return false; - - } - return (fv_it == fvr.end()); - } - - // void erase_simplex(Simplex_handle sh) - // { - // auto siblings = sc.self_siblings(sh); - // auto oncles = siblings->oncles(); - // int prev_vertex = siblings->parent(); - // siblings->members().erase(sh->first); - // if (siblings->members().empty()) { - // typename typedef Simplicial_complex::Siblings Siblings; - // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); - // assert(!sc.has_children(oncles->members().find(prev_vertex))); - // //delete siblings; - // } - - // } - - void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) - { - erase_simplex(coface_sh); - erase_simplex(face_sh); - } - -public: - // void collapse(std::vector& simplices) - // { - // // Get a vector of simplex handles ordered by filtration value - // std::cout << sc << std::endl; - // //std::vector simplices; - // for (Simplex_handle sh: sc.filtration_simplex_range()) - // simplices.push_back(sh); - // // std::sort(simplices.begin(), - // // simplices.end(), - // // [&](Simplex_handle sh1, Simplex_handle sh2) - // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); - // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); - // // Double iteration - // auto face_it = simplices.rbegin(); - // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { - // int coface_count = 0; - // auto reduced_coface = simplices.rbegin(); - // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) - // if (face_it != coface_it && is_face(*face_it, *coface_it)) { - // coface_count++; - // if (coface_count == 1) - // reduced_coface = coface_it; - // else - // break; - // } - // if (coface_count == 1) { - // std::cout << "Erase ( "; - // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) - // std::cout << v << " "; - - // simplices.erase(--(reduced_coface.base())); - // //elementary_collapse(*face_it, *reduced_coface); - // std::cout << ") and then ( "; - // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) - // std::cout << v << " "; - // std::cout << ")\n"; - // simplices.erase(--((face_it++).base())); - // //face_it = simplices.rbegin(); - // //std::cout << "Size of vector: " << simplices.size() << "\n"; - // } - // else - // face_it++; - // } - // sc.initialize_filtration(); - // //std::cout << sc << std::endl; - // } - - template - void collapse(Dim_lists& dim_lists) - { - dim_lists.collapse(); - } - -private: - - /** Collapse recursively boundary faces of the given simplex - * if its filtration is bigger than alpha_lim. - */ - void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) - { - for (Simplex_handle face_it : sc.boundary_simplex_range()) { - - } - - } - -}; //class Relaxed_witness_complex - -} // namespace witness_complex - -} // namespace Gudhi - -#endif diff --git a/src/Witness_complex/include/gudhi/Strange_witness_complex.h b/src/Witness_complex/include/gudhi/Strange_witness_complex.h deleted file mode 100644 index c1c2912b..00000000 --- a/src/Witness_complex/include/gudhi/Strange_witness_complex.h +++ /dev/null @@ -1,232 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef WITNESS_COMPLEX_H_ -#define WITNESS_COMPLEX_H_ - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - -/** - \class Witness_complex - \brief Constructs the witness complex for the given set of witnesses and landmarks. - \ingroup witness_complex - */ -template< class Simplicial_complex> -class Strange_witness_complex { - private: - struct Active_witness { - int witness_id; - int landmark_id; - - Active_witness(int witness_id_, int landmark_id_) - : witness_id(witness_id_), - landmark_id(landmark_id_) { } - }; - - private: - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - - typedef std::vector< double > Point_t; - typedef std::vector< Point_t > Point_Vector; - - // typedef typename Simplicial_complex::Filtration_value Filtration_value; - typedef std::vector< Vertex_handle > typeVectorVertex; - typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; - typedef std::pair< Simplex_handle, bool > typePairSimplexBool; - - typedef int Witness_id; - typedef int Landmark_id; - typedef std::list< Vertex_handle > ActiveWitnessList; - - private: - int nbL; // Number of landmarks - Simplicial_complex& sc; // Simplicial complex - - public: - ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /** @name Constructor - */ - - //@{ - - // Witness_range> - - /** - * \brief Iterative construction of the witness complex. - * \details The witness complex is written in sc_ basing on a matrix knn of - * nearest neighbours of the form {witnesses}x{landmarks}. - * - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges. - * - * Constructor takes into account at most (dim+1) - * first landmarks from each landmark range to construct simplices. - * - * Landmarks are supposed to be in [0,nbL_-1] - */ - template< typename KNearestNeighbors > - Strange_witness_complex(KNearestNeighbors const & knn, - Simplicial_complex & sc_, - int nbL_, - int dim) : nbL(nbL_), sc(sc_) { - // Construction of the active witness list - int nbW = knn.size(); - typeVectorVertex vv; - int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - for (int i = 0; i != nbL; ++i) { - // initial fill of 0-dimensional simplices - // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore - counter++; - vv = {i}; - sc.insert_simplex(vv); - // TODO(SK) Error if not inserted : normally no need here though - } - for (int i = 0; i != nbW; ++i) { - std::vector simplex; - simplex.reserve(dim+1); - for (int j = 0; j <= dim; j++) - simplex.push_back(knn[i][j]); - sc.insert_simplex_and_subfaces(simplex); - } - } - - //@} - - private: - /** \brief Check if the facets of the k-dimensional simplex witnessed - * by witness witness_id are already in the complex. - * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id - */ - template - bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) { - // std::cout << "All face in with the landmark " << inserted_vertex << std::endl; - std::vector< Vertex_handle > facet; - // Vertex_handle curr_vh = curr_sh->first; - // CHECK ALL THE FACETS - for (int i = 0; i != k + 1; ++i) { - facet = {}; - for (int j = 0; j != k + 1; ++j) { - if (j != i) { - facet.push_back(knn[witness_id][j]); - } - } // endfor - if (sc.find(facet) == sc.null_simplex()) - return false; - // std::cout << "++++ finished loop safely\n"; - } // endfor - return true; - } - - template - static void print_vector(const std::vector& v) { - std::cout << "["; - if (!v.empty()) { - std::cout << *(v.begin()); - for (auto it = v.begin() + 1; it != v.end(); ++it) { - std::cout << ","; - std::cout << *it; - } - } - std::cout << "]"; - } - - public: - /** - * \brief Verification if every simplex in the complex is witnessed by witnesses in knn. - * \param print_output =true will print the witnesses for each simplex - * \remark Added for debugging purposes. - */ - template< class KNearestNeighbors > - bool is_witness_complex(KNearestNeighbors const & knn, bool print_output) { - // bool final_result = true; - for (Simplex_handle sh : sc.complex_simplex_range()) { - bool is_witnessed = false; - typeVectorVertex simplex; - int nbV = 0; // number of verticed in the simplex - for (int v : sc.simplex_vertex_range(sh)) - simplex.push_back(v); - nbV = simplex.size(); - for (typeVectorVertex w : knn) { - bool has_vertices = true; - for (int v : simplex) - if (std::find(w.begin(), w.begin() + nbV, v) == w.begin() + nbV) { - has_vertices = false; - // break; - } - if (has_vertices) { - is_witnessed = true; - if (print_output) { - std::cout << "The simplex "; - print_vector(simplex); - std::cout << " is witnessed by the witness "; - print_vector(w); - std::cout << std::endl; - } - break; - } - } - if (!is_witnessed) { - if (print_output) { - std::cout << "The following simplex is not witnessed "; - print_vector(simplex); - std::cout << std::endl; - } - assert(is_witnessed); - return false; - } - } - return true; - } -}; - -} // namespace witness_complex - -} // namespace Gudhi - -#endif // WITNESS_COMPLEX_H_ diff --git a/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h deleted file mode 100644 index ee863a42..00000000 --- a/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h +++ /dev/null @@ -1,337 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef STRONG_RELAXED_WITNESS_COMPLEX_H_ -#define STRONG_RELAXED_WITNESS_COMPLEX_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" -#include -#include -#include -#include -#include -#include -#include -#include - -// Needed for nearest neighbours -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** \addtogroup simplex_tree - * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. - */ -template< class Simplicial_complex > -class Strong_relaxed_witness_complex { - -private: - struct Active_witness { - int witness_id; - int landmark_id; - - Active_witness(int witness_id_, int landmark_id_) - : witness_id(witness_id_), - landmark_id(landmark_id_) { } - }; - -private: - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - typedef typename Simplicial_complex::Filtration_value FT; - - typedef std::vector< double > Point_t; - typedef std::vector< Point_t > Point_Vector; - - // typedef typename Simplicial_complex::Filtration_value Filtration_value; - typedef std::vector< Vertex_handle > typeVectorVertex; - typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; - typedef std::pair< Simplex_handle, bool > typePairSimplexBool; - - typedef int Witness_id; - typedef int Landmark_id; - typedef std::list< Vertex_handle > ActiveWitnessList; - - private: - int nbL; // Number of landmarks - Simplicial_complex& sc; // Simplicial complex - - public: - /** @name Constructor - */ - - //@{ - - /** - * \brief Iterative construction of the relaxed witness complex. - * \details The witness complex is written in sc_ basing on a matrix knn - * of k nearest neighbours of the form {witnesses}x{landmarks} and - * and a matrix distances of distances to these landmarks from witnesses. - * The parameter alpha defines relaxation and - * limD defines the - * - * The line lengths in one given matrix can differ, - * however both matrices have the same corresponding line lengths. - * - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges. - * - * Constructor takes into account at most (dim+1) - * first landmarks from each landmark range to construct simplices. - * - * Landmarks are supposed to be in [0,nbL_-1] - */ - - template< typename KNearestNeighbours > - Strong_relaxed_witness_complex(std::vector< std::vector > const & distances, - KNearestNeighbours const & knn, - Simplicial_complex & sc_, - int nbL_, - double alpha2, - unsigned limD) : nbL(nbL_), sc(sc_) { - int nbW = knn.size(); - typeVectorVertex vv; - //int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - ActiveWitnessList active_w;// = new ActiveWitnessList(); - for (int i = 0; i != nbL; ++i) { - // initial fill of 0-dimensional simplices - // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore - //counter++; - vv = {i}; - sc.insert_simplex(vv, Filtration_value(0.0)); - /* TODO Error if not inserted : normally no need here though*/ - } - for (int i=0; i != nbW; ++i) { - // int i_end = limD+1; - // if (knn[i].size() < limD+1) - // i_end = knn[i].size(); - // double dist_wL = *(distances[i].begin()); - // while (distances[i][i_end] > dist_wL + alpha2) - // i_end--; - // add_all_witnessed_faces(distances[i].begin(), - // knn[i].begin(), - // knn[i].begin() + i_end + 1); - unsigned j_end = 0; - while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) { - std::vector simplex; - for (unsigned j = 0; j <= j_end; ++j) - simplex.push_back(knn[i][j]); - assert(distances[i][j_end] - distances[i][0] >= 0); - sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]); - j_end++; - } - } - sc.set_dimension(limD); - } - - //@} - -private: - /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness - * Iterator is needed to know until how far we can take landmarks to form simplexes - * simplex is the prefix of the simplexes to insert - * The output value indicates if the witness rests active or not - */ - void add_all_witnessed_faces(std::vector::const_iterator curr_d, - std::vector::const_iterator curr_l, - std::vector::const_iterator end) - { - std::vector simplex; - std::vector::const_iterator l_end = curr_l; - for (; l_end != end; ++l_end) { - std::vector::const_iterator l_it = curr_l; - std::vector::const_iterator d_it = curr_d; - simplex = {}; - for (; l_it != l_end; ++l_it, ++d_it) - simplex.push_back(*l_it); - sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d); - } - } - - /** \brief Check if the facets of the k-dimensional simplex witnessed - * by witness witness_id are already in the complex. - * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id - */ - bool all_faces_in(std::vector& simplex, double* filtration_value) - { - std::vector< int > facet; - for (std::vector::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) - { - facet.clear(); - for (std::vector::iterator it = simplex.begin(); it != simplex.end(); ++it) - if (it != not_it) - facet.push_back(*it); - Simplex_handle facet_sh = sc.find(facet); - if (facet_sh == sc.null_simplex()) - return false; - else if (sc.filtration(facet_sh) > *filtration_value) - *filtration_value = sc.filtration(facet_sh); - } - return true; - } - - bool is_face(Simplex_handle face, Simplex_handle coface) - { - // vertex range is sorted in decreasing order - auto fvr = sc.simplex_vertex_range(face); - auto cfvr = sc.simplex_vertex_range(coface); - auto fv_it = fvr.begin(); - auto cfv_it = cfvr.begin(); - while (fv_it != fvr.end() && cfv_it != cfvr.end()) { - if (*fv_it < *cfv_it) - ++cfv_it; - else if (*fv_it == *cfv_it) { - ++cfv_it; - ++fv_it; - } - else - return false; - - } - return (fv_it == fvr.end()); - } - - // void erase_simplex(Simplex_handle sh) - // { - // auto siblings = sc.self_siblings(sh); - // auto oncles = siblings->oncles(); - // int prev_vertex = siblings->parent(); - // siblings->members().erase(sh->first); - // if (siblings->members().empty()) { - // typename typedef Simplicial_complex::Siblings Siblings; - // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); - // assert(!sc.has_children(oncles->members().find(prev_vertex))); - // //delete siblings; - // } - - // } - - void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) - { - erase_simplex(coface_sh); - erase_simplex(face_sh); - } - -public: - // void collapse(std::vector& simplices) - // { - // // Get a vector of simplex handles ordered by filtration value - // std::cout << sc << std::endl; - // //std::vector simplices; - // for (Simplex_handle sh: sc.filtration_simplex_range()) - // simplices.push_back(sh); - // // std::sort(simplices.begin(), - // // simplices.end(), - // // [&](Simplex_handle sh1, Simplex_handle sh2) - // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); - // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); - // // Double iteration - // auto face_it = simplices.rbegin(); - // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { - // int coface_count = 0; - // auto reduced_coface = simplices.rbegin(); - // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) - // if (face_it != coface_it && is_face(*face_it, *coface_it)) { - // coface_count++; - // if (coface_count == 1) - // reduced_coface = coface_it; - // else - // break; - // } - // if (coface_count == 1) { - // std::cout << "Erase ( "; - // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) - // std::cout << v << " "; - - // simplices.erase(--(reduced_coface.base())); - // //elementary_collapse(*face_it, *reduced_coface); - // std::cout << ") and then ( "; - // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) - // std::cout << v << " "; - // std::cout << ")\n"; - // simplices.erase(--((face_it++).base())); - // //face_it = simplices.rbegin(); - // //std::cout << "Size of vector: " << simplices.size() << "\n"; - // } - // else - // face_it++; - // } - // sc.initialize_filtration(); - // //std::cout << sc << std::endl; - // } - - template - void collapse(Dim_lists& dim_lists) - { - dim_lists.collapse(); - } - -private: - - /** Collapse recursively boundary faces of the given simplex - * if its filtration is bigger than alpha_lim. - */ - void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) - { - for (Simplex_handle face_it : sc.boundary_simplex_range()) { - - } - - } - -}; //class Strong_relaxed_witness_complex - -} // namespace witness_complex - -} // namespace Gudhi - -#endif diff --git a/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h deleted file mode 100644 index a1aedd8e..00000000 --- a/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h +++ /dev/null @@ -1,379 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 - * 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 . - */ - -#ifndef WEAK_RELAXED_WITNESS_COMPLEX_H_ -#define WEAK_RELAXED_WITNESS_COMPLEX_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" -#include -#include -#include -#include -#include -#include -#include -#include - -// Needed for nearest neighbours -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - -namespace witness_complex { - - /** \addtogroup simplex_tree - * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. - */ -template< class Simplicial_complex > -class Weak_relaxed_witness_complex { - -private: - struct Active_witness { - int witness_id; - int landmark_id; - - Active_witness(int witness_id_, int landmark_id_) - : witness_id(witness_id_), - landmark_id(landmark_id_) { } - }; - -private: - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - typedef typename Simplicial_complex::Filtration_value FT; - - typedef std::vector< double > Point_t; - typedef std::vector< Point_t > Point_Vector; - - // typedef typename Simplicial_complex::Filtration_value Filtration_value; - typedef std::vector< Vertex_handle > typeVectorVertex; - typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; - typedef std::pair< Simplex_handle, bool > typePairSimplexBool; - - typedef int Witness_id; - typedef int Landmark_id; - typedef std::list< Vertex_handle > ActiveWitnessList; - - private: - int nbL; // Number of landmarks - Simplicial_complex& sc; // Simplicial complex - - public: - /** @name Constructor - */ - - //@{ - - /** - * \brief Iterative construction of the relaxed witness complex. - * \details The witness complex is written in sc_ basing on a matrix knn - * of k nearest neighbours of the form {witnesses}x{landmarks} and - * and a matrix distances of distances to these landmarks from witnesses. - * The parameter alpha defines relaxation and - * limD defines the - * - * The line lengths in one given matrix can differ, - * however both matrices have the same corresponding line lengths. - * - * The type KNearestNeighbors can be seen as - * Witness_range>, where - * Witness_range and Closest_landmark_range are random access ranges. - * - * Constructor takes into account at most (dim+1) - * first landmarks from each landmark range to construct simplices. - * - * Landmarks are supposed to be in [0,nbL_-1] - */ - - template< typename KNearestNeighbours > - Weak_relaxed_witness_complex(std::vector< std::vector > const & distances, - KNearestNeighbours const & knn, - Simplicial_complex & sc_, - int nbL_, - double alpha, - int limD) : nbL(nbL_), sc(sc_) { - int nbW = knn.size(); - typeVectorVertex vv; - //int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - ActiveWitnessList active_w;// = new ActiveWitnessList(); - for (int i = 0; i != nbL; ++i) { - // initial fill of 0-dimensional simplices - // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore - //counter++; - vv = {i}; - sc.insert_simplex(vv, Filtration_value(0.0)); - /* TODO Error if not inserted : normally no need here though*/ - } - int k = 1; /* current dimension in iterative construction */ - for (int i=0; i != nbW; ++i) - active_w.push_back(i); - while (!active_w.empty() && k <= limD && k < nbL ) - { - typename ActiveWitnessList::iterator aw_it = active_w.begin(); - while (aw_it != active_w.end()) - { - std::vector simplex; - bool ok = add_all_faces_of_dimension(k, - alpha, - std::numeric_limits::infinity(), - distances[*aw_it].begin(), - knn[*aw_it].begin(), - simplex, - knn[*aw_it].end()); - if (!ok) - active_w.erase(aw_it++); //First increase the iterator and then erase the previous element - else - aw_it++; - } - k++; - } - sc.set_dimension(limD); - } - - //@} - -private: - /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness - * Iterator is needed to know until how far we can take landmarks to form simplexes - * simplex is the prefix of the simplexes to insert - * The output value indicates if the witness rests active or not - */ - bool add_all_faces_of_dimension(int dim, - double alpha, - double norelax_dist, - std::vector::const_iterator curr_d, - std::vector::const_iterator curr_l, - std::vector& simplex, - std::vector::const_iterator end) - { - if (curr_l == end) - return false; - bool will_be_active = false; - std::vector::const_iterator l_it = curr_l; - std::vector::const_iterator d_it = curr_d; - if (dim > 0) - for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { - simplex.push_back(*l_it); - if (sc.find(simplex) != sc.null_simplex()) - will_be_active = add_all_faces_of_dimension(dim-1, - alpha, - norelax_dist, - d_it+1, - l_it+1, - simplex, - end) || will_be_active; - simplex.pop_back(); - // If norelax_dist is infinity, change to first omitted distance - if (*d_it <= norelax_dist) - norelax_dist = *d_it; - will_be_active = add_all_faces_of_dimension(dim-1, - alpha, - norelax_dist, - d_it+1, - l_it+1, - simplex, - end) || will_be_active; - } - else if (dim == 0) - for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { - simplex.push_back(*l_it); - double filtration_value = 0; - // if norelax_dist is infinite, relaxation is 0. - if (*d_it > norelax_dist) - filtration_value = *d_it - norelax_dist; - if (all_faces_in(simplex, &filtration_value)) { - will_be_active = true; - sc.insert_simplex(simplex, filtration_value); - } - simplex.pop_back(); - // If norelax_dist is infinity, change to first omitted distance - if (*d_it < norelax_dist) - norelax_dist = *d_it; - } - return will_be_active; - } - - /** \brief Check if the facets of the k-dimensional simplex witnessed - * by witness witness_id are already in the complex. - * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id - */ - bool all_faces_in(std::vector& simplex, double* filtration_value) - { - std::vector< int > facet; - for (std::vector::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) - { - facet.clear(); - for (std::vector::iterator it = simplex.begin(); it != simplex.end(); ++it) - if (it != not_it) - facet.push_back(*it); - Simplex_handle facet_sh = sc.find(facet); - if (facet_sh == sc.null_simplex()) - return false; - else if (sc.filtration(facet_sh) > *filtration_value) - *filtration_value = sc.filtration(facet_sh); - } - return true; - } - - bool is_face(Simplex_handle face, Simplex_handle coface) - { - // vertex range is sorted in decreasing order - auto fvr = sc.simplex_vertex_range(face); - auto cfvr = sc.simplex_vertex_range(coface); - auto fv_it = fvr.begin(); - auto cfv_it = cfvr.begin(); - while (fv_it != fvr.end() && cfv_it != cfvr.end()) { - if (*fv_it < *cfv_it) - ++cfv_it; - else if (*fv_it == *cfv_it) { - ++cfv_it; - ++fv_it; - } - else - return false; - - } - return (fv_it == fvr.end()); - } - - // void erase_simplex(Simplex_handle sh) - // { - // auto siblings = sc.self_siblings(sh); - // auto oncles = siblings->oncles(); - // int prev_vertex = siblings->parent(); - // siblings->members().erase(sh->first); - // if (siblings->members().empty()) { - // typename typedef Simplicial_complex::Siblings Siblings; - // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); - // assert(!sc.has_children(oncles->members().find(prev_vertex))); - // //delete siblings; - // } - - // } - - void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) - { - erase_simplex(coface_sh); - erase_simplex(face_sh); - } - -public: - // void collapse(std::vector& simplices) - // { - // // Get a vector of simplex handles ordered by filtration value - // std::cout << sc << std::endl; - // //std::vector simplices; - // for (Simplex_handle sh: sc.filtration_simplex_range()) - // simplices.push_back(sh); - // // std::sort(simplices.begin(), - // // simplices.end(), - // // [&](Simplex_handle sh1, Simplex_handle sh2) - // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); - // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); - // // Double iteration - // auto face_it = simplices.rbegin(); - // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { - // int coface_count = 0; - // auto reduced_coface = simplices.rbegin(); - // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) - // if (face_it != coface_it && is_face(*face_it, *coface_it)) { - // coface_count++; - // if (coface_count == 1) - // reduced_coface = coface_it; - // else - // break; - // } - // if (coface_count == 1) { - // std::cout << "Erase ( "; - // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) - // std::cout << v << " "; - - // simplices.erase(--(reduced_coface.base())); - // //elementary_collapse(*face_it, *reduced_coface); - // std::cout << ") and then ( "; - // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) - // std::cout << v << " "; - // std::cout << ")\n"; - // simplices.erase(--((face_it++).base())); - // //face_it = simplices.rbegin(); - // //std::cout << "Size of vector: " << simplices.size() << "\n"; - // } - // else - // face_it++; - // } - // sc.initialize_filtration(); - // //std::cout << sc << std::endl; - // } - - template - void collapse(Dim_lists& dim_lists) - { - dim_lists.collapse(); - } - -private: - - /** Collapse recursively boundary faces of the given simplex - * if its filtration is bigger than alpha_lim. - */ - void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) - { - for (Simplex_handle face_it : sc.boundary_simplex_range()) { - - } - - } - -}; //class Weak_relaxed_witness_complex - -} // namespace witness_complex - -} // namespace Gudhi - -#endif diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index e732cb18..a16f9270 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -173,8 +173,6 @@ private: else aw_it++; } - std::cout << "Active witnesses after dim=" << k << " is finished: " << active_witnesses.size() << "\n"; - std::cout << complex << "\n"; k++; } return true; @@ -233,11 +231,8 @@ private: simplex.push_back(l_it->first); double filtration_value = 0; // if norelax_dist is infinite, relaxation is 0. - //std::cout << "landmark_id=" << l_it->first << " distance=" << l_it->second << "\n"; - // std::size_t landmark_id = l_it->first; - // double distance = l_it->second; if (l_it->second > norelax_dist2) - filtration_value = l_it->second - norelax_dist2; + filtration_value = l_it->second - norelax_dist2; if (all_faces_in(simplex, &filtration_value, sc)) { will_be_active = true; sc.insert_simplex(simplex, filtration_value); -- cgit v1.2.3