diff options
Diffstat (limited to 'src/Witness_complex/example')
-rw-r--r-- | src/Witness_complex/example/Landmark_choice_random_knn.h | 142 | ||||
-rw-r--r-- | src/Witness_complex/example/Landmark_choice_sparsification.h | 230 | ||||
-rw-r--r-- | src/Witness_complex/example/a0_persistence.cpp | 638 | ||||
-rw-r--r-- | src/Witness_complex/example/bench_rwit.cpp | 709 | ||||
-rw-r--r-- | src/Witness_complex/example/color-picker.cpp | 37 | ||||
-rw-r--r-- | src/Witness_complex/example/off_writer.h | 11 | ||||
-rw-r--r-- | src/Witness_complex/example/relaxed_delaunay.cpp | 180 | ||||
-rw-r--r-- | src/Witness_complex/example/relaxed_witness_complex_sphere.cpp | 461 | ||||
-rw-r--r-- | src/Witness_complex/example/relaxed_witness_persistence.cpp | 267 | ||||
-rw-r--r-- | src/Witness_complex/example/strange_sphere.cpp | 219 | ||||
-rw-r--r-- | src/Witness_complex/example/witness_complex_from_file.cpp | 36 |
11 files changed, 20 insertions, 2910 deletions
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 <http://www.gnu.org/licenses/>. - */ - -#ifndef LANDMARK_CHOICE_BY_RANDOM_KNN_H_ -#define LANDMARK_CHOICE_BY_RANDOM_KNN_H_ - -#include <utility> // for pair<> -#include <vector> -#include <cstddef> // for ptrdiff_t type - -//#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -//#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_incremental_neighbor_search.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> -//#include <CGAL/Kernel_d/Vector_d.h> -#include <CGAL/Random.h> - - -namespace Gudhi { - -namespace witness_complex { - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> 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<Traits_base> 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 <typename KNearestNeighbours, - typename Point_random_access_range, - typename Distance_matrix> - 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<Point_d> landmarks; - std::vector<int> 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<std::ptrdiff_t,Point_d*,Euclidean_distance> adapter(&(landmarks[0])); - Euclidean_distance ed; - Tree landmark_tree(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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 <http://www.gnu.org/licenses/>. - */ - -#ifndef LANDMARK_CHOICE_BY_SPARSIFICATION_H_ -#define LANDMARK_CHOICE_BY_SPARSIFICATION_H_ - -#include <utility> // for pair<> -#include <vector> -#include <cstddef> // for ptrdiff_t type -#include <algorithm> - -//#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -//#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_incremental_neighbor_search.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> -//#include <CGAL/Kernel_d/Vector_d.h> -#include <CGAL/Random.h> -#include <CGAL/Fuzzy_sphere.h> - -namespace Gudhi { - -namespace witness_complex { - - typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> 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<Traits_base> 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<STraits> Fuzzy_sphere; - - /** \brief Landmark selection function by as a sub-epsilon-net of the - * given set of points. - */ - template <typename Point_random_access_range> - 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<std::ptrdiff_t,Point_d*,Euclidean_distance> points_adapter(&(points[0])); - std::vector<bool> dropped_points(nbP, false); - - Tree witness_tree(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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<int> close_neighbors; - witness_tree.search(std::insert_iterator<std::vector<int>>(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 <typename KNearestNeighbours, - typename Point_random_access_range, - typename Distance_matrix> - 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<std::ptrdiff_t,Point_d*,Euclidean_distance> adapter(&(landmarks[0])); - Euclidean_distance ed; - Tree landmark_tree(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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 <typename Kernel, typename Point_container> - std::vector<typename Point_container::value_type> - 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<Kernel, Point_container> 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<typename Point_container::value_type> output; - - Points_ds points_ds(input_pts); - - std::vector<bool> 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 <http://www.gnu.org/licenses/>. - */ - -#include <gudhi/Simplex_tree.h> -#include <gudhi/A0_complex.h> -#include <gudhi/Relaxed_witness_complex.h> -#include <gudhi/Dim_lists.h> -#include <gudhi/reader_utils.h> -#include <gudhi/Persistent_cohomology.h> -#include "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <queue> -#include <iterator> -#include <string> - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -#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_d> 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<std::vector< int > > knn; - std::vector<std::vector< FT > > distances; - start = clock(); - //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); - - std::vector<Point_d> 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<double>(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<double>(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<double>(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<Simplex_tree<>> simplices(simplex_tree, limD); - - // simplices.collapse(); - // simplices.output_simplices(); - - // Simplex_tree<> collapsed_tree; - // for (auto sh: simplices) { - // std::vector<int> vertices; - // for (int v: collapsed_tree.simplex_vertex_range(sh)) - // vertices.push_back(v); - // collapsed_tree.insert_simplex(vertices); - // } - std::vector<int> 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<double> Point_t; - typedef Simplex_tree<Simplex_tree_options_fast_persistence> 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<Point_t>); - // 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<double>(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<ST, Field_Zp > 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<std::vector< int > > knn; - std::vector<std::vector< FT > > distances; - std::vector<Point_d> 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<int> & 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_endpoint> 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<int> 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<Simplex_tree<>> 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<std::vector< int > > const & knn, - std::vector<std::vector< FT > > const & distances, - unsigned nbL, - double alpha2, - std::vector<int>& 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<double>(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<double>(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<int> 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<Point_d> landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector<std::vector< int > > knn; - std::vector<std::vector< FT > > 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<int> 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<Point_d> landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector<std::vector< int > > knn; - std::vector<std::vector< FT > > 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 <http://www.gnu.org/licenses/>. - */ - -#include <gudhi/Simplex_tree.h> -#include <gudhi/A0_complex.h> -#include <gudhi/Relaxed_witness_complex.h> -#include <gudhi/Dim_lists.h> -#include <gudhi/reader_utils.h> -#include <gudhi/Persistent_cohomology.h> -#include <gudhi/Good_links.h> -#include "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <queue> -#include <iterator> -#include <string> - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -#include <boost/program_options.hpp> - -#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_d> Point_Vector; -//typedef Simplex_tree<Simplex_tree_options_fast_persistence> 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<int> & desired_homology - , double & min_persistence) { - namespace po = boost::program_options; - po::options_description hidden("Hidden options"); - hidden.add_options() - ("option", po::value<int>(& experiment_number), - "Experiment id.") - ("input-file", po::value<std::string>(&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<std::string>(&experiment_name)->default_value("witness"), - "The prefix of all the output files. Default is 'witness'") - ("landmarks,L", po::value<int>(&nbL)->default_value(0), - "Number of landmarks.") - ( "landmark-file,l", po::value<std::string>(&landmark_file), - "Name of a fike containing landmarks") - ("alpha2_s,A", po::value<double>(&alpha2_s)->default_value(0), - "Relaxation parameter for the strong complex.") - ("alpha2_w,a", po::value<double>(&alpha2_w)->default_value(0), - "Relaxation parameter for the weak complex.") - ("mu_epsilon,e", po::value<double>(&mu_epsilon)->default_value(0), - "Sparsification parameter.") - ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), - "Maximal dimension of the Witness complex we want to compute.") - - ("homology,H", po::value<std::vector<int>>(&desired_homology)->multitoken(), - "The desired Betti numbers.") - ("min-persistence,m", po::value<Filtration_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<std::vector<FT> >); - // 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<std::vector< int > > knn; - std::vector<std::vector< FT > > distances; - start = clock(); - //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); - - std::vector<Point_d> 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<double>(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<double>(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<double>(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<STree> simplices(simplex_tree, limD); - - // std::vector<Simplex_handle> simplices; - std::cout << "Starting collapses...\n"; - simplices.collapse(); - simplices.output_simplices(); - - STree collapsed_tree; - for (auto sh: simplices) { - std::vector<int> vertices; - for (int v: collapsed_tree.simplex_vertex_range(sh)) - vertices.push_back(v); - collapsed_tree.insert_simplex(vertices); - } - std::vector<int> 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<STree> 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<double>(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<ST, Field_Zp > 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<int> & 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_endpoint> 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<int> 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<STree> 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<std::vector< int > > const & knn, - std::vector<std::vector< FT > > const & distances, - Point_Vector & points, - unsigned nbL, - unsigned limD, - double alpha2_s, - double alpha2_w, - std::vector<int>& 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<double>(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<double>(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<int> desired_homology = {1}; - std::vector<Point_d> 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<std::vector< int > > knn; - std::vector<std::vector< FT > > 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<int> 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<Point_d> landmarks; - - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector<std::vector< int > > knn; - std::vector<std::vector< FT > > 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<int> 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<Point_d> landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - - std::vector<std::vector< int > > knn; - std::vector<std::vector< FT > > 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 <iostream> -#include <assert.h> - -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 <http://www.gnu.org/licenses/>. - */ - -#include <gudhi/Simplex_tree.h> -#include <gudhi/Relaxed_witness_complex.h> -#include <gudhi/Dim_lists.h> -#include <gudhi/reader_utils.h> -#include <gudhi/Persistent_cohomology.h> -#include "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <queue> -#include <iterator> -#include <string> - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -#include <CGAL/Epick_d.h> -#include <CGAL/Delaunay_triangulation.h> -//#include <CGAL/Sphere_d.h> - -#include "generators.h" -#include "output.h" - -using namespace Gudhi; -using namespace Gudhi::witness_complex; -using namespace Gudhi::persistent_cohomology; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -typedef K::Sphere_d Sphere_d; -typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; - -typedef std::vector<Point_d> 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<std::vector<int>> knn; - std::vector<std::vector<double>> 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<Simplex_tree<>> rwc(distances, - knn, - simplex_tree, - point_vector.size(), - alpha2, - limD); - std::vector<int> 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<int> 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<int> 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 <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <queue> -#include <iterator> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Relaxed_witness_complex.h" -#include "gudhi/reader_utils.h" -#include "gudhi/Collapse/Collapse.h" -//#include <boost/filesystem.hpp> - -//#include <CGAL/Delaunay_triangulation.h> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_incremental_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> - -#include <CGAL/Kernel_d/Vector_d.h> -#include <CGAL/point_generators_d.h> -#include <CGAL/constructions_d.h> -#include <CGAL/Fuzzy_sphere.h> -#include <CGAL/Random.h> - - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> 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<Traits_base> Euclidean_distance; - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> 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<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_incremental_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; - -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef std::vector<Point_d> Point_Vector; - -//typedef K::Equal_d Equal_d; -typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; -typedef CGAL::Random_points_in_ball_d<Point_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<Point_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 <int> > & 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 <std::vector<int>::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<Point_d> 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<FT> 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<Point_d> 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<int>& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - //std::vector<Point_d> landmarks; - int chosen_landmark; - //std::pair<Point_etiquette_map::iterator,bool> 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<int>& landmarks_ind, FT alpha) -{ - //********************Preface: origin point - unsigned D = W[0].size(); - std::vector<FT> 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 <int> > WL(nbP); - std::vector< std::vector< typename std::vector<int>::iterator > > ope_limits(nbP); - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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<int>::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark - Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(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<int>::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<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(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<double>({0.1,0.1})), Point_d(std::vector<double>({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<int> 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 <http://www.gnu.org/licenses/>. - */ - -#include <gudhi/Simplex_tree.h> -#include <gudhi/Relaxed_witness_complex.h> -#include <gudhi/Dim_lists.h> -#include <gudhi/reader_utils.h> -#include <gudhi/Persistent_cohomology.h> -#include "Landmark_choice_random_knn.h" -#include "Landmark_choice_sparsification.h" - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <queue> -#include <iterator> -#include <string> - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -#include "generators.h" -#include "output.h" - -using namespace Gudhi; -using namespace Gudhi::witness_complex; -using namespace Gudhi::persistent_cohomology; - -typedef std::vector<Point_d> 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<std::vector< int > > knn; - std::vector<std::vector< FT > > distances; - start = clock(); - //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); - std::vector<Point_d> 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<double>(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<double>(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<double>(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 <http://www.gnu.org/licenses/>. - */ -#define BOOST_PARAMETER_MAX_ARITY 12 - - -#include <sys/types.h> -#include <sys/stat.h> - -#include <gudhi/Simplex_tree.h> -#include <gudhi/Strange_witness_complex.h> -#include <gudhi/Landmark_choice_by_random_point.h> -#include <gudhi/reader_utils.h> -#include "output.h" - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <string> -#include <vector> - -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_incremental_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> - -#include <CGAL/Kernel_d/Vector_d.h> -#include <CGAL/point_generators_d.h> -#include <CGAL/constructions_d.h> -#include <CGAL/Fuzzy_sphere.h> -#include <CGAL/Random.h> - - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - - -#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<CGAL::Dynamic_dimension_tag> 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<Traits_base> Euclidean_distance; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -typedef CGAL::Orthogonal_incremental_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; -typedef std::vector<Point_d> 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<int>& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - //std::vector<Point_d> landmarks; - int chosen_landmark; - //std::pair<Point_etiquette_map::iterator,bool> 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<FT> 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<int> landmarks_ind; - landmark_choice(W, nbP, nbL, landmarks, landmarks_ind); - - - STraits traits(&(landmarks[0])); - Euclidean_distance ed; - std::vector< std::vector <int> > WL(nbP); - std::vector< std::vector< typename std::vector<int>::iterator > > ope_limits(nbP); - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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<int>::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark - Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(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<int, double> > 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<std::vector< int > > knn; - - landmarks_to_witness_complex(point_vector, nbL); - - end = clock(); - double time = static_cast<double>(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 <gudhi/Simplex_tree.h> #include <gudhi/Witness_complex.h> -#include <gudhi/Construct_closest_landmark_table.h> #include <gudhi/pick_n_random_points.h> #include <gudhi/reader_utils.h> +#include <CGAL/Epick_d.h> + #include <iostream> #include <fstream> #include <ctime> #include <string> #include <vector> +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef typename K::Point_d Point_d; +typedef typename Gudhi::witness_complex::Witness_complex<K> Witness_complex; typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::vector< std::vector <double> > Point_Vector; +typedef std::vector< Point_d > Point_vector; /** * \brief Customized version of read_points @@ -44,7 +48,7 @@ typedef std::vector< std::vector <double> > 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<std::vector< int > > 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<double>(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<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; } |