diff options
Diffstat (limited to 'src')
38 files changed, 6272 insertions, 104 deletions
diff --git a/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h index 4c9c04d9..5963caa3 100644 --- a/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h +++ b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h @@ -63,7 +63,7 @@ namespace cubical_complex { * For further details and theory of cubical complexes, please consult \cite kaczynski2004computational as well as the * following paper \cite peikert2012topological . * - * \section datastructure Data structure. + * \section cubicalcomplexdatastructure Data structure. * * The implementation of Cubical complex provides a representation of complexes that occupy a rectangular region in * \f$\mathbb{R}^n\f$. This extra assumption allows for a memory efficient way of storing cubical complexes in a form diff --git a/src/Contraction/include/gudhi/Edge_contraction.h b/src/Contraction/include/gudhi/Edge_contraction.h index 5af13c3e..61f2d945 100644 --- a/src/Contraction/include/gudhi/Edge_contraction.h +++ b/src/Contraction/include/gudhi/Edge_contraction.h @@ -41,7 +41,7 @@ namespace contraction { \author David Salinas -\section Introduction +\section edgecontractionintroduction Introduction The purpose of this package is to offer a user-friendly interface for edge contraction simplification of huge simplicial complexes. It uses the \ref skbl data-structure whose size remains small during simplification diff --git a/src/Doxyfile b/src/Doxyfile index 949838b8..51950b3d 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -782,7 +782,9 @@ RECURSIVE = YES EXCLUDE = data/ \ example/ \ - GudhUI/ + GudhUI/ \ + cmake/ \ + debian/ # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or # directories that are symbolic links (a Unix file system feature) are excluded @@ -1803,18 +1805,6 @@ GENERATE_XML = NO XML_OUTPUT = xml -# The XML_SCHEMA tag can be used to specify a XML schema, which can be used by a -# validating XML parser to check the syntax of the XML files. -# This tag requires that the tag GENERATE_XML is set to YES. - -XML_SCHEMA = - -# The XML_DTD tag can be used to specify a XML DTD, which can be used by a -# validating XML parser to check the syntax of the XML files. -# This tag requires that the tag GENERATE_XML is set to YES. - -XML_DTD = - # If the XML_PROGRAMLISTING tag is set to YES doxygen will dump the program # listings (including syntax highlighting and cross-referencing information) to # the XML output. Note that enabling this will significantly increase the size diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index 0cba6361..433cfd3e 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -139,7 +139,7 @@ namespace persistent_cohomology { by increasing filtration values (breaking ties so as a simplex appears after its subsimplices of same filtration value) provides an indexing scheme. -\section Examples +\section pcohexamples Examples We provide several example files: run these examples with -h for details on their use, and read the README file. diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index d97d1b63..96d3e73a 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -45,6 +45,7 @@ if(GMP_FOUND) target_link_libraries(rips_multifield_persistence ${TBB_LIBRARIES}) target_link_libraries(performance_rips_persistence ${TBB_LIBRARIES}) endif(TBB_FOUND) + message("CGAL_LIBRARY = ${CGAL_LIBRARY}") add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.2 -d 3 -p 2 -q 71 -m 100) endif(GMPXX_FOUND) diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h index bd907131..32fe411c 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h @@ -42,7 +42,7 @@ namespace skeleton_blocker { \author David Salinas -\section Introduction +\section skblintroduction Introduction The Skeleton-Blocker data-structure proposes a light encoding for simplicial complexes by storing only an *implicit* representation of its simplices \cite socg_blockers_2011,\cite blockers2012. @@ -53,7 +53,7 @@ This data-structure handles all simplicial complexes operations such as are operations that do not require simplex enumeration such as edge iteration, link computation or simplex contraction. -\section Definitions +\section skbldefinitions Definitions We recall briefly classical definitions of simplicial complexes \cite Munkres-elementsalgtop1984. @@ -108,7 +108,7 @@ and point access in addition. -\subsection Visitor +\subsection skblvisitor Visitor The class Skeleton_blocker_complex has a visitor that is called when usual operations such as adding an edge or remove a vertex are called. You may want to use this visitor to compute statistics or to update another data-structure (for instance this visitor is heavily used in the \ref contr package). @@ -116,7 +116,7 @@ You may want to use this visitor to compute statistics or to update another data -\section Example +\section skblexample Example \subsection Iterating Iterating through vertices, edges, blockers and simplices diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h index b0ee35f5..2c9602fa 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h @@ -147,7 +147,7 @@ struct Trie { } void remove_leaf() { - assert(is_leaf); + assert(is_leaf()); if (!is_root()) parent_->childs.erase(this); } diff --git a/src/Skeleton_blocker/test/test_skeleton_blocker_complex.cpp b/src/Skeleton_blocker/test/test_skeleton_blocker_complex.cpp index 7b303935..4f9888ba 100644 --- a/src/Skeleton_blocker/test/test_skeleton_blocker_complex.cpp +++ b/src/Skeleton_blocker/test/test_skeleton_blocker_complex.cpp @@ -358,7 +358,8 @@ BOOST_AUTO_TEST_CASE(test_skeleton_iterator_blockers) { num_blockers = 0; for (auto blockers : complex.blocker_range()) { -#ifndef _WIN64 +// If not windows - _WIN32 is for windows 32 and 64 bits +#ifndef _WIN32 for (auto block_ptr = myBlockers.begin(); block_ptr < myBlockers.end(); block_ptr++) if (*block_ptr == *blockers) myBlockers.erase(block_ptr); @@ -366,7 +367,8 @@ BOOST_AUTO_TEST_CASE(test_skeleton_iterator_blockers) { num_blockers++; } BOOST_CHECK(num_blockers == 4); -#ifndef _WIN64 +// If not windows - _WIN32 is for windows 32 and 64 bits +#ifndef _WIN32 BOOST_CHECK(myBlockers.empty()); #endif } diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 4d67e0d0..0054775d 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -10,7 +10,8 @@ if(CGAL_FOUND) if (EIGEN3_FOUND) add_executable ( witness_complex_sphere witness_complex_sphere.cpp ) target_link_libraries(witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - add_test( witness_complex_sphere_10 ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere 10) + add_test( witness_complex_sphere_10 ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere 10) + add_executable ( relaxed_witness_persistence relaxed_witness_persistence.cpp ) endif(EIGEN3_FOUND) endif (NOT CGAL_VERSION VERSION_LESS 4.6.0) endif() diff --git a/src/Witness_complex/example/Landmark_choice_random_knn.h b/src/Witness_complex/example/Landmark_choice_random_knn.h new file mode 100644 index 00000000..fb91e116 --- /dev/null +++ b/src/Witness_complex/example/Landmark_choice_random_knn.h @@ -0,0 +1,142 @@ +/* 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 new file mode 100644 index 00000000..d5d5fba1 --- /dev/null +++ b/src/Witness_complex/example/Landmark_choice_sparsification.h @@ -0,0 +1,230 @@ +/* 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 new file mode 100644 index 00000000..422185ca --- /dev/null +++ b/src/Witness_complex/example/a0_persistence.cpp @@ -0,0 +1,638 @@ +/* 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 new file mode 100644 index 00000000..2d3a009c --- /dev/null +++ b/src/Witness_complex/example/bench_rwit.cpp @@ -0,0 +1,709 @@ +/* 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 new file mode 100644 index 00000000..fcea033e --- /dev/null +++ b/src/Witness_complex/example/color-picker.cpp @@ -0,0 +1,37 @@ +#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/generators.h b/src/Witness_complex/example/generators.h index ac445261..731a52b0 100644 --- a/src/Witness_complex/example/generators.h +++ b/src/Witness_complex/example/generators.h @@ -25,10 +25,12 @@ #include <CGAL/Epick_d.h> #include <CGAL/point_generators_d.h> +#include <CGAL/Random.h> #include <fstream> #include <string> #include <vector> +#include <cmath> typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; typedef K::FT FT; @@ -144,4 +146,21 @@ void generate_points_sphere(Point_Vector& W, int nbP, int dim) { W.push_back(*rp++); } +/** \brief Generate nbP points on a (flat) d-torus embedded in R^{2d} + * + */ +void generate_points_torus(Point_Vector& W, int nbP, int dim) { + CGAL::Random rand; + const double pi = std::acos(-1); + for (int i = 0; i < nbP; i++) { + std::vector<FT> point; + for (int j = 0; j < dim; j++) { + double alpha = rand.uniform_real((double)0, 2*pi); + point.push_back(sin(alpha)); + point.push_back(cos(alpha)); + } + W.push_back(Point_d(point)); + } +} + #endif // EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ diff --git a/src/Witness_complex/example/off_writer.h b/src/Witness_complex/example/off_writer.h new file mode 100644 index 00000000..8cafe9e2 --- /dev/null +++ b/src/Witness_complex/example/off_writer.h @@ -0,0 +1,11 @@ +#ifndef OFF_WRITER_H_ +#define OFF_WRITER_H_ + + +class Off_Writer { + +private: + +} + +#endif diff --git a/src/Witness_complex/example/output.h b/src/Witness_complex/example/output.h new file mode 100644 index 00000000..d3f534af --- /dev/null +++ b/src/Witness_complex/example/output.h @@ -0,0 +1,308 @@ +#ifndef OUTPUT_H +#define OUTPUT_H + +#include <fstream> +#include <vector> +#include <string> + +#include <gudhi/Simplex_tree.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Delaunay_triangulation.h> + +//typename Gudhi::Witness_complex<> Witness_complex; + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::Point_d Point_d; +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; + +/** \brief Write the table of the nearest landmarks to each witness + * to a file. + */ +template <class Value> +void write_wl( std::string file_name, std::vector< std::vector <Value> > & WL) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : WL) + { + for (auto l: w) + ofs << l << " "; + ofs << "\n"; + } + ofs.close(); +} + +/** \brief Write the coordinates of points in points to a file. + * + */ +void write_points( std::string file_name, std::vector< Point_d > & points) +{ + 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(); +} + +/** Write edges of a witness complex in a file. + * The format of an edge is coordinates of u \n coordinates of v \n\n\n + * This format is compatible with gnuplot + */ +template< typename STree > +void write_edges(std::string file_name, STree& witness_complex, Point_Vector& landmarks) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto u: witness_complex.complex_vertex_range()) + for (auto v: witness_complex.complex_vertex_range()) + { + std::vector<int> 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(); +} + +/** \brief Write triangles (tetrahedra in 3d) of a simplicial complex in a file, compatible with medit. + * `landmarks_ind` represents the set of landmark indices in W + * `st` is the Simplex_tree to be visualized, + * `shr` is the Simplex_handle_range of simplices in `st` to be visualized + * `is2d` should be true if the simplicial complex is 2d, false if 3d + * `l_is_v` = landmark is vertex + */ +template <typename SimplexHandleRange, + typename STree > +void write_witness_mesh(Point_Vector& W, std::vector<int>& landmarks_ind, STree& st, SimplexHandleRange const & shr, bool is2d, bool l_is_v, std::string file_name = "witness.mesh") +{ + std::ofstream ofs (file_name, std::ofstream::out); + if (is2d) + ofs << "MeshVersionFormatted 1\nDimension 2\n"; + else + ofs << "MeshVersionFormatted 1\nDimension 3\n"; + + if (!l_is_v) + ofs << "Vertices\n" << W.size() << "\n"; + else + ofs << "Vertices\n" << landmarks_ind.size() << "\n"; + + if (l_is_v) + for (auto p_it : landmarks_ind) { + for (auto coord = W[p_it].cartesian_begin(); coord != W[p_it].cartesian_end() && coord != W[p_it].cartesian_begin()+3 ; ++coord) + ofs << *coord << " "; + ofs << "508\n"; + } + else + for (auto p_it : W) { + for (auto coord = p_it.cartesian_begin(); coord != p_it.cartesian_end() && coord != p_it.cartesian_begin()+3 ; ++coord) + ofs << *coord << " "; + ofs << "508\n"; + } + + // int num_triangles = W.size(), num_tetrahedra = 0; + int num_edges = 0, num_triangles = 0, num_tetrahedra = 0; + if (!l_is_v) { + for (auto sh_it : shr) + if (st.dimension(sh_it) == 1) + num_edges++; + else if (st.dimension(sh_it) == 2) + num_triangles++; + else if (st.dimension(sh_it) == 3) + num_tetrahedra++; + ofs << "Edges " << num_edges << "\n"; + for (auto sh_it : shr) { + if (st.dimension(sh_it) == 1) { + for (auto v_it : st.simplex_vertex_range(sh_it)) + ofs << landmarks_ind[v_it]+1 << " "; + ofs << "200\n"; + } + } + ofs << "Triangles " << num_triangles << "\n"; + for (unsigned i = 0; i < W.size(); ++i) + ofs << i << " " << i << " " << i << " " << "508\n"; + for (auto sh_it : shr) + { + if (st.dimension(sh_it) == 2) { + for (auto v_it : st.simplex_vertex_range(sh_it)) + ofs << landmarks_ind[v_it]+1 << " "; + ofs << "508\n"; + } + } + ofs << "Tetrahedra " << num_tetrahedra << "\n"; + for (auto sh_it : shr) + { + if (st.dimension(sh_it) == 3) { + for (auto v_it : st.simplex_vertex_range(sh_it)) + ofs << landmarks_ind[v_it]+1 << " "; + ofs << "250\n"; + } + } + } + else { + for (auto sh_it : shr) + if (st.dimension(sh_it) == 1) + num_edges++; + else if (st.dimension(sh_it) == 2) + num_triangles++; + else if (st.dimension(sh_it) == 3) + num_tetrahedra++; + ofs << "Edges " << num_edges << "\n"; + for (auto sh_it : shr) { + if (st.dimension(sh_it) == 1) { + for (auto v_it : st.simplex_vertex_range(sh_it)) + ofs << v_it+1 << " "; + ofs << "200\n"; + } + } + ofs << "Triangles " << num_triangles << "\n"; + for (auto sh_it : shr) + { + if (st.dimension(sh_it) == 2) { + for (auto v_it : st.simplex_vertex_range(sh_it)) + ofs << v_it+1 << " "; + ofs << "508\n"; + } + } + ofs << "Tetrahedra " << num_tetrahedra << "\n"; + for (auto sh_it : shr) + { + if (st.dimension(sh_it) == 3) { + for (auto v_it : st.simplex_vertex_range(sh_it)) + ofs << v_it+1 << " "; + ofs << "250\n"; + } + } + } + + ofs << "End\n"; + /* + else + { + ofs << "Tetrahedra " << t.number_of_finite_full_cells()+1 << "\n"; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) + ofs << index_of_vertex[*vh_it] << " "; + ofs << "508\n"; + } + ofs << nbV << " " << nbV << " " << nbV << " " << nbV << " " << 208 << "\n"; + ofs << "End\n"; + } + */ + ofs.close(); +} + +void write_witness_mesh(Point_Vector& W, std::vector<int>& landmarks_ind, Gudhi::Simplex_tree<>& st, bool is2d, bool l_is_v, std::string file_name = "witness.mesh") +{ + write_witness_mesh(W, landmarks_ind, st, st.complex_simplex_range(), is2d, l_is_v, file_name); +} + +/** \brief Write triangles (tetrahedra in 3d) of a Delaunay + * triangulation in a file, compatible with medit. + */ +void write_delaunay_mesh(Delaunay_triangulation& t, const Point_d& p, bool is2d) +{ + std::ofstream ofs ("delaunay.mesh", std::ofstream::out); + int nbV = t.number_of_vertices()+1; + if (is2d) + ofs << "MeshVersionFormatted 1\nDimension 2\n"; + else + ofs << "MeshVersionFormatted 1\nDimension 3\n"; + ofs << "Vertices\n" << nbV << "\n"; + int ind = 1; //index of a vertex + std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + { + if (t.is_infinite(v_it)) + continue; + // Add maximum 3 coordinates + for (auto coord = v_it->point().cartesian_begin(); coord != v_it->point().cartesian_end() && coord != v_it->point().cartesian_begin()+3; ++coord) + ofs << *coord << " "; + ofs << "508\n"; + index_of_vertex[v_it] = ind++; + } + for (auto coord = p.cartesian_begin(); coord != p.cartesian_end(); ++coord) + ofs << *coord << " "; + ofs << "208\n"; + if (is2d) + { + ofs << "Triangles " << t.number_of_finite_full_cells()+1 << "\n"; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) + ofs << index_of_vertex[*vh_it] << " "; + ofs << "508\n"; + } + ofs << nbV << " " << nbV << " " << nbV << " " << 208 << "\n"; + ofs << "End\n"; + } + else if (p.size() == 3) + { + ofs << "Tetrahedra " << t.number_of_finite_full_cells()+1 << "\n"; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) + ofs << index_of_vertex[*vh_it] << " "; + ofs << "508\n"; + } + ofs << nbV << " " << nbV << " " << nbV << " " << nbV << " " << 208 << "\n"; + ofs << "End\n"; + } + else if (p.size() == 4) + { + ofs << "Tetrahedra " << 5*(t.number_of_finite_full_cells())+1 << "\n"; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) + { + for (auto vh_it2 = fc_it->vertices_begin(); vh_it2 != fc_it->vertices_end(); ++vh_it2) + if (vh_it != vh_it2) + ofs << index_of_vertex[*vh_it2] << " "; + ofs << "508\n"; + } + } + ofs << nbV << " " << nbV << " " << nbV << " " << nbV << " " << 208 << "\n"; + ofs << "End\n"; + } + ofs.close(); +} + +/////////////////////////////////////////////////////////////////////// +// PRINT VECTOR +/////////////////////////////////////////////////////////////////////// + +template <typename T> +void print_vector(std::vector<T> v) +{ + std::cout << "["; + if (!v.empty()) + { + std::cout << *(v.begin()); + for (auto it = v.begin()+1; it != v.end(); ++it) + { + std::cout << ","; + std::cout << *it; + } + } + std::cout << "]"; +} + + +#endif diff --git a/src/Witness_complex/example/output_tikz.h b/src/Witness_complex/example/output_tikz.h new file mode 100644 index 00000000..281c65b7 --- /dev/null +++ b/src/Witness_complex/example/output_tikz.h @@ -0,0 +1,178 @@ +#ifndef OUTPUT_TIKZ_H +#define OUTPUT_TIKZ_H + +#include <vector> +#include <string> +#include <algorithm> +#include <fstream> +#include <cmath> + +typedef double FT; + + +//////////////////AUX///////////////////// + +void write_preamble(std::ofstream& ofs) +{ + ofs << "\\documentclass{standalone}\n" + << "\\usepackage{tikz}\n\n" + << "\\begin{document}\n" + << "\\begin{tikzpicture}\n"; +} + +void write_end(std::ofstream& ofs) +{ + ofs << "\\end{tikzpicture}\n" + << "\\end{document}"; +} + + +/////////////////MAIN////////////////////// + +void write_tikz_plot(std::vector<FT> data, std::string filename) +{ + int n = data.size(); + FT vmax = *(std::max_element(data.begin(), data.end())); + //std::cout << std::log10(vmax) << " " << std::floor(std::log10(vmax)); + + FT order10 = pow(10,std::floor(std::log10(vmax))); + int digit = std::floor( vmax / order10) + 1; + if (digit == 4 || digit == 6) digit = 5; + if (digit > 6) digit = 10; + FT plot_max = digit*order10; + std::cout << plot_max << " " << vmax; + FT hstep = 10.0/(n-1); + FT wstep = 10.0 / plot_max; + + std::cout << "(eps_max-eps_min)/(N-48) = " << (vmax-*data.begin())/(data.size()-48) << "\n"; + std::ofstream ofs(filename, std::ofstream::out); + + ofs << + "\\documentclass{standalone}\n" << + "\\usepackage[utf8]{inputenc}\n" << + "\\usepackage{amsmath}\n" << + "\\usepackage{tikz}\n\n" << + "\\begin{document}\n" << + "\\begin{tikzpicture}\n"; + + ofs << "\\draw[->] (0,0) -- (0,11);" << std::endl << + "\\draw[->] (0,0) -- (11,0);" << std::endl << + "\\foreach \\i in {1,...,10}" << std::endl << + "\\draw (0,\\i) -- (-0.05,\\i);" << std::endl << + "\\foreach \\i in {1,...,10}" << std::endl << + "\\draw (\\i,0) -- (\\i,-0.05);" << std::endl << std::endl << + + "\\foreach \\i in {1,...,10}" << std::endl << + "\\draw[dashed] (-0.05,\\i) -- (11,\\i);" << std::endl << std::endl << + + "\\node at (-0.5,11) {$*$}; " << std::endl << + "\\node at (11,-0.5) {$*$}; " << std::endl << + "\\node at (-0.5,-0.5) {0}; " << std::endl << + "\\node at (-0.5,10) {" << plot_max << "}; " << std::endl << + "%\\node at (10,-0.5) {2}; " << std::endl; + + ofs << "\\draw[red] (0," << wstep*data[0] << ")"; + for (int i = 1; i < n; ++i) + ofs << " -- (" << hstep*i << "," << wstep*data[i] << ")"; + ofs << ";\n"; + + ofs << + "\\end{tikzpicture}\n" << + "\\end{document}"; + + ofs.close(); +} + + +// A little script to make a tikz histogram of epsilon distribution +// Returns the average epsilon +void write_histogram(std::vector<double> histo, std::string file_name = "histogram.tikz", std::string xaxis = "$\\epsilon/\\epsilon_{max}$", std::string yaxis = "$\\epsilon$", FT max_x = 1) +{ + int n = histo.size(); + + std::ofstream ofs (file_name, std::ofstream::out); + FT barwidth = 20.0/n; + FT max_value = *(std::max_element(histo.begin(), histo.end())); + std::cout << max_value << std::endl; + FT ten_power = pow(10, ceil(log10(max_value))); + FT max_histo = ten_power; + if (max_value/ten_power > 1) { + if (max_value/ten_power < 2) + max_histo = 0.2*ten_power; + else if (max_value/ten_power < 5) + max_histo = 0.5*ten_power; + } + std::cout << ceil(log10(max_value)) << std::endl << max_histo << std::endl; + FT unitht = max_histo/10.0; + write_preamble(ofs); + + ofs << "\\draw[->] (0,0) -- (0,11);\n" << + "\\draw[->] (0,0) -- (21,0);\n" << + "\\foreach \\i in {1,...,10}\n" << + "\\draw (0,\\i) -- (-0.1,\\i);\n" << + "\\foreach \\i in {1,...,20}\n" << + "\\draw (\\i,0) -- (\\i,-0.1);\n" << + + "\\node at (-1,11) {" << yaxis << "};\n" << + "\\node at (22,-1) {" << xaxis << "};\n" << + "\\node at (-0.5,-0.5) {0};\n" << + "\\node at (-0.5,10) {" << max_histo << "};\n" << + "\\node at (20,-0.5) {" << max_x << "};\n"; + + for (int i = 0; i < n; ++i) + ofs << "\\draw (" << barwidth*i << "," << histo[i]/unitht << ") -- (" + << barwidth*(i+1) << "," << histo[i]/unitht << ") -- (" + << barwidth*(i+1) << ",0) -- (" << barwidth*i << ",0) -- cycle;\n"; + + write_end(ofs); + ofs.close(); +} + +struct Pers_interval { + double alpha_start, alpha_end; + int dim; + Pers_interval(double alpha_start_, double alpha_end_, int dim_) + : alpha_start(alpha_start_), alpha_end(alpha_end_), dim(dim_) + {} +}; + +void write_barcodes(std::string in_file, double alpha2, std::string out_file = "barcodes.tikz.tex") +{ + std::ifstream ifs(in_file, std::ios::in); + std::string line; + std::vector<Pers_interval> pers_intervals; + while (getline(ifs, 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_intervals.push_back(Pers_interval(alpha_start, alpha_end, dim)); + } + } + ifs.close(); + std::ofstream ofs (out_file, std::ofstream::out); + write_preamble(ofs); + double barwidth = 0.01; + int i = 0; + for (auto interval: pers_intervals) { + std::string color = "black"; + switch (interval.dim) { + case 0: color = "orange"; break; + case 1: color = "red"; break; + case 2: color = "blue"; break; + case 3: color = "green"; break; + case 4: color = "yellow"; break; + default: color = "orange"; break; + } + ofs << "\\fill[" << color << "] (" << interval.alpha_start << "," << barwidth*i << ") rectangle (" + << interval.alpha_end << "," << barwidth*(i+1) <<");\n"; + i++; + } + write_end(ofs); + ofs.close(); +} + +#endif diff --git a/src/Witness_complex/example/relaxed_delaunay.cpp b/src/Witness_complex/example/relaxed_delaunay.cpp new file mode 100644 index 00000000..b0190446 --- /dev/null +++ b/src/Witness_complex/example/relaxed_delaunay.cpp @@ -0,0 +1,180 @@ +/* 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 new file mode 100644 index 00000000..067321ce --- /dev/null +++ b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp @@ -0,0 +1,461 @@ +/* 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 new file mode 100644 index 00000000..b4837594 --- /dev/null +++ b/src/Witness_complex/example/relaxed_witness_persistence.cpp @@ -0,0 +1,267 @@ +/* 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 new file mode 100644 index 00000000..9188bd8e --- /dev/null +++ b/src/Witness_complex/example/strange_sphere.cpp @@ -0,0 +1,219 @@ +/* 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_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index e6f88274..3d4a54aa 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -31,6 +31,8 @@ #include <gudhi/pick_n_random_points.h> #include <gudhi/reader_utils.h> +#include <CGAL/Epick_d.h> + #include <iostream> #include <fstream> #include <ctime> @@ -80,7 +82,12 @@ int main(int argc, char * const argv[]) { Gudhi::witness_complex::construct_closest_landmark_table(point_vector, landmarks, knn); // Compute witness complex - Gudhi::witness_complex::witness_complex(knn, number_of_landmarks, point_vector[0].size(), simplex_tree); + // Gudhi::witness_complex::witness_complex(knn, number_of_landmarks, point_vector[0].size(), simplex_tree); + Gudhi::witness_complex::Witness_complex<Gudhi::Simplex_tree<>, CGAL::Epick_d<CGAL::Dynamic_dimension_tag>>(landmarks.begin(), + landmarks.end(), + point_vector.begin(), + point_vector.end(), + simplex_tree); end = clock(); double time = static_cast<double>(end - start) / CLOCKS_PER_SEC; std::cout << "Witness complex for " << number_of_landmarks << " landmarks took " diff --git a/src/Witness_complex/include/gudhi/A0_complex.h b/src/Witness_complex/include/gudhi/A0_complex.h new file mode 100644 index 00000000..2018e1c8 --- /dev/null +++ b/src/Witness_complex/include/gudhi/A0_complex.h @@ -0,0 +1,337 @@ +/* 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 A0_COMPLEX_H_ +#define A0_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#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_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + +namespace witness_complex { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class A0_complex { + +private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } + }; + +private: + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + // typedef typename Simplicial_complex::Filtration_value Filtration_value; + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL; // Number of landmarks + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Constructor takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + + template< typename KNearestNeighbours > + A0_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha2, + unsigned limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + typeVectorVertex vv; + //int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + ActiveWitnessList active_w;// = new ActiveWitnessList(); + for (int i = 0; i != nbL; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + //counter++; + vv = {i}; + sc.insert_simplex(vv, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ + } + for (int i=0; i != nbW; ++i) { + // int i_end = limD+1; + // if (knn[i].size() < limD+1) + // i_end = knn[i].size(); + // double dist_wL = *(distances[i].begin()); + // while (distances[i][i_end] > dist_wL + alpha2) + // i_end--; + // add_all_witnessed_faces(distances[i].begin(), + // knn[i].begin(), + // knn[i].begin() + i_end + 1); + unsigned j_end = 0; + while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) { + std::vector<int> simplex; + for (unsigned j = 0; j <= j_end; ++j) + simplex.push_back(knn[i][j]); + assert(distances[i][j_end] - distances[i][0] >= 0); + sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]); + j_end++; + } + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness + * Iterator is needed to know until how far we can take landmarks to form simplexes + * simplex is the prefix of the simplexes to insert + * The output value indicates if the witness rests active or not + */ + void add_all_witnessed_faces(std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>::const_iterator end) + { + std::vector<int> simplex; + std::vector<int>::const_iterator l_end = curr_l; + for (; l_end != end; ++l_end) { + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + simplex = {}; + for (; l_it != l_end; ++l_it, ++d_it) + simplex.push_back(*l_it); + sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d); + } + } + + /** \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Dim_list_iterator.h b/src/Witness_complex/include/gudhi/Dim_list_iterator.h new file mode 100644 index 00000000..3f23e8c9 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Dim_list_iterator.h @@ -0,0 +1,155 @@ +/* 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 DIM_LISTS_ITERATOR_H_ +#define DIM_LISTS_ITERATOR_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + +namespace witness_complex { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplex_handle > +class Dim_lists_iterator { + +private: + + typedef typename std::list<Simplex_handle>::iterator List_iterator; + typedef typename std::vector<std::list<Simplex_handle>>::reverse_iterator Vector_iterator; + typedef typename Gudhi::witness_complex::Dim_lists_iterator<Simplex_handle> Iterator; + + + List_iterator sh_; + Vector_iterator curr_list_; + typename std::vector<std::list<Simplex_handle>>& table_; + +public: + + Dim_lists_iterator(List_iterator sh, + Vector_iterator curr_list, + typename std::vector<std::list<Simplex_handle>>& table) + : sh_(sh), curr_list_(curr_list), table_(table) + { + } + + Simplex_handle operator*() + { + return *sh_; + } + + Iterator operator++() + { + increment(); + return (*this); + } + + Iterator operator++(int) + { + Iterator prev_it(sh_, curr_list_, table_); + increment(); + return prev_it; + } + + Iterator dim_begin() + { + return Iterator(curr_list_->begin(), curr_list_, table_); + } + + Iterator dim_end() + { + return Iterator(curr_list_->end(), curr_list_, table_); + } + + Iterator dimp1_begin() + { + return Iterator((curr_list_-1)->begin(), curr_list_-1, table_); + } + + Iterator dimp1_end() + { + return Iterator((curr_list_-1)->end(), curr_list_-1, table_); + } + + bool operator==(const Iterator& it2) const + { + return (sh_ == it2.sh_); + } + + bool operator!=(const Iterator& it2) const + { + return (sh_ != it2.sh_); + } + + void remove_incr() + { + + } + +private: + + void increment() + { + if (++sh_ == curr_list_->end()) + if (++curr_list_ != table_.rend()) + sh_ = curr_list_->begin(); + // The iterator of the end of the table is the end of the last list + } + + +}; //class Dim_lists_iterator + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Dim_lists.h b/src/Witness_complex/include/gudhi/Dim_lists.h new file mode 100644 index 00000000..1d1b03c3 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Dim_lists.h @@ -0,0 +1,195 @@ +/* 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 DIM_LISTS_H_ +#define DIM_LISTS_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include <gudhi/Dim_list_iterator.h> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + +namespace witness_complex { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Dim_lists { + +private: + + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Gudhi::witness_complex::Dim_lists_iterator<Simplex_handle> Iterator; + + std::vector<std::list<Simplex_handle>> table_; + Simplicial_complex& sc_; + +public: + + Dim_lists(Simplicial_complex & sc, int dim, double alpha_max = 100) + : sc_(sc) + { + table_ = std::vector<std::list<Simplex_handle>>(dim+1); + for (auto sh: sc.filtration_simplex_range()) { + if (sc_.filtration(sh) < alpha_max) + table_[sc.dimension(sh)].push_front(sh); + } + auto t_it = table_.rbegin(); + while (t_it->empty()) { + t_it++; + table_.pop_back(); + } + } + + Iterator begin() + { + return Iterator(table_.rbegin()->begin(), table_.rbegin(), table_); + } + + Iterator end() + { + return Iterator(table_[0].end(), table_.rend(), table_); + } + + unsigned size() + { + unsigned curr_size = 0; + for (auto l: table_) + curr_size += l.size(); + return curr_size; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc_.simplex_vertex_range(face); + auto cfvr = sc_.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + void output_simplices() { + std::cout << "Size of vector: " << size() << std::endl; + for (auto line_it = table_.rbegin(); line_it != table_.rend(); ++line_it) + for (auto sh: *line_it) { + std::cout << sc_.dimension(sh) << " "; + for (auto v : sc_.simplex_vertex_range(sh)) + std::cout << v << " "; + std::cout << sc_.filtration(sh) << "\n"; + } + } + + void collapse() + { + auto coface_list_it = table_.rbegin(); + auto face_list_it = table_.rbegin()+1; + for ( ; + face_list_it != table_.rend(); + ++face_list_it) { + auto face_it = face_list_it->begin(); + while (face_it != face_list_it->end() && sc_.filtration(*face_it) != 0) { + int coface_count = 0; + auto reduced_coface = coface_list_it->begin(); + for (auto coface_it = coface_list_it->begin(); coface_it != coface_list_it->end() && sc_.filtration(*coface_it) != 0; ++coface_it) + if (is_face(*face_it, *coface_it)) { + coface_count++; + if (coface_count == 1) + reduced_coface = coface_it; + else + break; + } + if (coface_count == 1) { + /* + std::cout << "Erase ( "; + for (auto v: sc_.simplex_vertex_range(*reduced_coface)) + std::cout << v << " "; + */ + coface_list_it->erase(reduced_coface); + /* + std::cout << ") and then ( "; + for (auto v: sc_.simplex_vertex_range(*face_it)) + std::cout << v << " "; + std::cout << ")\n"; + */ + face_list_it->erase(face_it); + face_it = face_list_it->begin(); + } + else + face_it++; + } + if ((coface_list_it++)->empty()) + table_.pop_back(); + } + } + + bool is_pseudomanifold() + { + + return true; + } + +}; //class Dim_lists + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Good_links.h b/src/Witness_complex/include/gudhi/Good_links.h new file mode 100644 index 00000000..c5e993e5 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Good_links.h @@ -0,0 +1,449 @@ +/* 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 GOOD_LINKS_H_ +#define GOOD_LINKS_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include <gudhi/Dim_list_iterator.h> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + +template < typename Simplicial_complex > +class Good_links { + + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef std::vector<Vertex_handle> Vertex_vector; + + // graph is an adjacency list + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph; + // map that gives to a certain simplex its node in graph and its dimension + //typedef std::pair<boost::vecS,Color> Reference; + typedef boost::graph_traits<Adj_graph>::vertex_descriptor Vertex_t; + typedef boost::graph_traits<Adj_graph>::edge_descriptor Edge_t; + typedef boost::graph_traits<Adj_graph>::adjacency_iterator Adj_it; + typedef std::pair<Adj_it, Adj_it> Out_edge_it; + + typedef boost::container::flat_map<Simplex_handle, Vertex_t> Graph_map; + typedef boost::container::flat_map<Vertex_t, Simplex_handle> Inv_graph_map; + +public: + Good_links(Simplicial_complex& sc): sc_(sc) + { + int dim = 0; + for (auto sh: sc_.complex_simplex_range()) + if (sc_.dimension(sh) > dim) + dim = sc_.dimension(sh); + dimension = dim; + count_good = Vertex_vector(dim); + count_bad = Vertex_vector(dim); + } + +private: + + Simplicial_complex& sc_; + unsigned dimension; + std::vector<int> count_good; + std::vector<int> count_bad; + + void add_vertices_to_link_graph(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + Vertex_vector& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename Vertex_vector::iterator it; + //std::pair<typename Graph_map::iterator,bool> resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) { + curr_simplex.push_back(*it); //push next vertex in question + curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + Vertex_vector curr_simplex_copy(curr_simplex); + sh = sc_.find(curr_simplex_copy); //a simplex of the star + curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = Vertex_vector(curr_simplex); + if (sh != sc_.null_simplex()) { + //std::cout << " added\n"; + if (curr_d == link_dimension) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + } + else { + if (curr_d == link_dimension-1) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_vertices_to_link_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + curr_simplex, + curr_d+1, link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } + + void add_edges_to_link_graph(Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map) + { + Simplex_handle sh; + // Add edges + //std::cout << "Entered add edges:\n"; + typename Graph_map::iterator map_it; + for (auto d_map_pair : d_map) { + //std::cout << "*"; + sh = d_map_pair.first; + Vertex_t d_vert = d_map_pair.second; + for (auto facet_sh : sc_.boundary_simplex_range(sh)) + //for (auto f_map_it : f_map) + { + //std::cout << "'"; + map_it = f_map.find(facet_sh); + //We must have all the facets in the graph at this point + assert(map_it != f_map.end()); + Vertex_t f_vert = map_it->second; + //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n"; + boost::add_edge(d_vert,f_vert,adj_graph); + } + } + } + + + + /* \brief Verifies if the simplices formed by vertices given by link_vertices + * form a pseudomanifold. + * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional + * faces and edges represent adjacency between them. + */ + bool link_is_pseudomanifold(Vertex_vector& star_vertices, + int dimension) + { + Adj_graph adj_graph; + Graph_map d_map, f_map; + // d_map = map for d-dimensional simplices + // f_map = map for its facets + Vertex_vector init_vector = {}; + add_vertices_to_link_graph(star_vertices, + star_vertices.begin()+1, + adj_graph, + d_map, + f_map, + init_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, + d_map, + f_map); + for (auto f_map_it : f_map) { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2){ + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector<int> components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + int star_dim(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + int curr_d, + Vertex_vector& curr_simplex, + typename std::vector< int >::iterator curr_dc) + { + //std::cout << "Entered star_dim for " << *(curr_v-1) << "\n"; + Simplex_handle sh; + int final_d = curr_d; + typename Vertex_vector::iterator it; + typename Vertex_vector::iterator dc_it; + //std::cout << "Current vertex is " << + for (it = curr_v, dc_it = curr_dc; it != star_vertices.end(); ++it, ++dc_it) + { + curr_simplex.push_back(*it); + Vertex_vector curr_simplex_copy(curr_simplex); + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " final_dim " << final_d; + */ + sh = sc_.find(curr_simplex_copy); //Need a copy because find sorts the vector and I want star center to be the first + if (sh != sc_.null_simplex()) + { + //std::cout << " -> " << *it << "\n"; + int d = star_dim(star_vertices, + it+1, + curr_d+1, + curr_simplex, + dc_it); + if (d >= final_d) + { + final_d = d; + //std::cout << d << " "; + //print_vector(curr_simplex); + //std::cout << std::endl; + } + if (d >= *dc_it) + *dc_it = d; + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); + } + return final_d; + } + +public: + + /** \brief Returns true if the link is good + */ + bool has_good_link(Vertex_handle v) + { + typedef Vertex_vector typeVectorVertex; + typeVectorVertex star_vertices; + // Fill star_vertices + star_vertices.push_back(v); + for (auto u: sc_.complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u != v && sc_.find(edge) != sc_.null_simplex()) + star_vertices.push_back(u); + } + // Find the dimension + typeVectorVertex init_simplex = {star_vertices[0]}; + bool is_pure = true; + std::vector<int> dim_coface(star_vertices.size(), 1); + int d = star_dim(star_vertices, + star_vertices.begin()+1, + 0, + init_simplex, + dim_coface.begin()+1) - 1; //link_dim = star_dim - 1 + + assert(init_simplex.size() == 1); + // if (!is_pure) + // std::cout << "Found an impure star around " << v << "\n"; + for (int dc: dim_coface) + is_pure = (dc == dim_coface[0]); + /* + if (d == count_good.size()) + { + std::cout << "Found a star of dimension " << (d+1) << " around " << v << "\nThe star is "; + print_vector(star_vertices); std::cout << std::endl; + } + */ + //if (d == -1) count_bad[0]++; + bool b= (is_pure && link_is_pseudomanifold(star_vertices,d)); + if (d != -1) {if (b) count_good[d]++; else count_bad[d]++;} + if (!is_pure) count_bad[0]++; + return (d != -1 && b && is_pure); + } + +private: + void add_max_simplices_to_graph(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + Inv_graph_map& inv_d_map, + Vertex_vector& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename Vertex_vector::iterator it; + //std::pair<typename Graph_map::iterator,bool> resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) { + curr_simplex.push_back(*it); //push next vertex in question + //curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + Vertex_vector curr_simplex_copy(curr_simplex); + sh = sc_.find(curr_simplex_copy); //a simplex of the star + //curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = Vertex_vector(curr_simplex); + if (sh != sc_.null_simplex()) { + //std::cout << " added\n"; + if (curr_d == link_dimension) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + inv_d_map.emplace(vert,sh); + } + else { + if (curr_d == link_dimension-1) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_max_simplices_to_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + inv_d_map, + curr_simplex, + curr_d+1, + link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } + +public: + bool complex_is_pseudomanifold() + { + Adj_graph adj_graph; + Graph_map d_map, f_map; + // d_map = map for d-dimensional simplices + // f_map = map for its facets + Inv_graph_map inv_d_map; + Vertex_vector init_vector = {}; + std::vector<int> star_vertices; + for (int v: sc_.complex_vertex_range()) + star_vertices.push_back(v); + add_max_simplices_to_graph(star_vertices, + star_vertices.begin(), + adj_graph, + d_map, + f_map, + inv_d_map, + init_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, + d_map, + f_map); + for (auto f_map_it : f_map) { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2) { + // if (boost::out_degree(f_map_it.second, adj_graph) >= 3) { + // // std::cout << "This simplex has 3+ cofaces: "; + // // for(auto v : sc_.simplex_vertex_range(f_map_it.first)) + // // std::cout << v << " "; + // // std::cout << std::endl; + // Adj_it ai, ai_end; + // for (std::tie(ai, ai_end) = boost::adjacent_vertices(f_map_it.second, adj_graph); ai != ai_end; ++ai) { + // auto it = inv_d_map.find(*ai); + // assert (it != inv_d_map.end()); + // typename Simplicial_complex::Simplex_handle sh = it->second; + // for(auto v : sc_.simplex_vertex_range(sh)) + // std::cout << v << " "; + // std::cout << std::endl; + // } + // } + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector<int> components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + int number_good_links(int dim) + { + return count_good[dim]; + } + + int number_bad_links(int dim) + { + return count_bad[dim]; + } + +}; + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h new file mode 100644 index 00000000..df93155b --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h @@ -0,0 +1,105 @@ +/* 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_FURTHEST_POINT_H_ +#define LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ + +#include <boost/range/size.hpp> + +#include <limits> // for numeric_limits<> +#include <iterator> +#include <algorithm> // for sort +#include <vector> + +namespace Gudhi { + +namespace witness_complex { + + typedef std::vector<int> typeVectorVertex; + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by iteratively adding the furthest witness from the + * current landmark set as the new landmark. + * \details It chooses nbL landmarks from a random access range `points` and + * writes {witness}*{closest landmarks} matrix in `knn`. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges + * + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_furthest_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nb_points = boost::size(points); + assert(nb_points >= nbL); + // distance matrix witness x landmarks + std::vector<std::vector<double>> wit_land_dist(nb_points, std::vector<double>()); + // landmark list + typeVectorVertex chosen_landmarks; + + knn = KNearestNeighbours(nb_points, std::vector<int>()); + int current_number_of_landmarks = 0; // counter for landmarks + double curr_max_dist = 0; // used for defining the furhest point from L + const double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nb_points, infty); // vector of current distances to L from points + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + // or better yet std::uniform_int_distribution + int rand_int = rand() % nb_points; + int curr_max_w = rand_int; // For testing purposes a pseudo-random number is used here + + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + // curr_max_w at this point is the next landmark + chosen_landmarks.push_back(curr_max_w); + unsigned i = 0; + for (auto& p : points) { + double curr_dist = euclidean_distance(p, *(std::begin(points) + chosen_landmarks[current_number_of_landmarks])); + wit_land_dist[i].push_back(curr_dist); + knn[i].push_back(current_number_of_landmarks); + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + ++i; + } + curr_max_dist = 0; + for (i = 0; i < dist_to_L.size(); i++) + if (dist_to_L[i] > curr_max_dist) { + curr_max_dist = dist_to_L[i]; + curr_max_w = i; + } + } + for (int i = 0; i < nb_points; ++i) + std::sort(std::begin(knn[i]), + std::end(knn[i]), + [&wit_land_dist, i](int a, int b) { + return wit_land_dist[i][a] < wit_land_dist[i][b]; }); + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h new file mode 100644 index 00000000..ebf6aad1 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h @@ -0,0 +1,96 @@ +/* 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_POINT_H_ +#define LANDMARK_CHOICE_BY_RANDOM_POINT_H_ + +#include <boost/range/size.hpp> + +#include <queue> // for priority_queue<> +#include <utility> // for pair<> +#include <iterator> +#include <vector> +#include <set> + +namespace Gudhi { + +namespace witness_complex { + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by taking random vertices for landmarks. + * \details It chooses nbL distinct landmarks from a random access range `points` + * and outputs a matrix {witness}*{closest landmarks} in knn. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges and + * Vertex_handle is the label type of a vertex in a simplicial complex. + * Closest_landmark_range needs to have push_back operation. + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_random_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nbP = boost::size(points); + assert(nbP >= nbL); + std::set<int> landmarks; + int current_number_of_landmarks = 0; // counter for landmarks + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + int chosen_landmark = rand() % nbP; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + while (landmarks.find(chosen_landmark) != landmarks.end()) + chosen_landmark = rand() % nbP; + landmarks.insert(chosen_landmark); + } + + int dim = boost::size(*std::begin(points)); + typedef std::pair<double, int> dist_i; + typedef bool (*comp)(dist_i, dist_i); + knn = KNearestNeighbours(nbP); + for (int points_i = 0; points_i < nbP; points_i++) { + std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([](dist_i j1, dist_i j2) { + return j1.first > j2.first; + }); + std::set<int>::iterator landmarks_it; + int landmarks_i = 0; + for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end(); + ++landmarks_it, landmarks_i++) { + dist_i dist = std::make_pair(euclidean_distance(points[points_i], points[*landmarks_it]), landmarks_i); + l_heap.push(dist); + } + for (int i = 0; i < dim + 1; i++) { + dist_i dist = l_heap.top(); + knn[points_i].push_back(dist.second); + l_heap.pop(); + } + } + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_RANDOM_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h new file mode 100644 index 00000000..2971149a --- /dev/null +++ b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h @@ -0,0 +1,379 @@ +/* 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 RELAXED_WITNESS_COMPLEX_H_ +#define RELAXED_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#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_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + +namespace witness_complex { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Relaxed_witness_complex { + +private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } + }; + +private: + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + // typedef typename Simplicial_complex::Filtration_value Filtration_value; + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL; // Number of landmarks + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Constructor takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + + template< typename KNearestNeighbours > + Relaxed_witness_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha, + int limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + typeVectorVertex vv; + //int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + ActiveWitnessList active_w;// = new ActiveWitnessList(); + for (int i = 0; i != nbL; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + //counter++; + vv = {i}; + sc.insert_simplex(vv, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ + } + int k = 1; /* current dimension in iterative construction */ + for (int i=0; i != nbW; ++i) + active_w.push_back(i); + while (!active_w.empty() && k <= limD && k < nbL ) + { + typename ActiveWitnessList::iterator aw_it = active_w.begin(); + while (aw_it != active_w.end()) + { + std::vector<int> simplex; + bool ok = add_all_faces_of_dimension(k, + alpha, + std::numeric_limits<double>::infinity(), + distances[*aw_it].begin(), + knn[*aw_it].begin(), + simplex, + knn[*aw_it].end()); + if (!ok) + active_w.erase(aw_it++); //First increase the iterator and then erase the previous element + else + aw_it++; + } + k++; + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness + * Iterator is needed to know until how far we can take landmarks to form simplexes + * simplex is the prefix of the simplexes to insert + * The output value indicates if the witness rests active or not + */ + bool add_all_faces_of_dimension(int dim, + double alpha, + double norelax_dist, + std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>& simplex, + std::vector<int>::const_iterator end) + { + if (curr_l == end) + return false; + bool will_be_active = false; + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + if (dim > 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + if (sc.find(simplex) != sc.null_simplex()) + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it <= norelax_dist) + norelax_dist = *d_it; + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + } + else if (dim == 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + double filtration_value = 0; + // if norelax_dist is infinite, relaxation is 0. + if (*d_it > norelax_dist) + filtration_value = *d_it - norelax_dist; + if (all_faces_in(simplex, &filtration_value)) { + will_be_active = true; + sc.insert_simplex(simplex, filtration_value); + } + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it < norelax_dist) + norelax_dist = *d_it; + } + return will_be_active; + } + + /** \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Strange_witness_complex.h b/src/Witness_complex/include/gudhi/Strange_witness_complex.h new file mode 100644 index 00000000..c1c2912b --- /dev/null +++ b/src/Witness_complex/include/gudhi/Strange_witness_complex.h @@ -0,0 +1,232 @@ +/* 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 WITNESS_COMPLEX_H_ +#define WITNESS_COMPLEX_H_ + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> + +#include <gudhi/distance_functions.h> + +#include <algorithm> +#include <utility> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <ctime> +#include <iostream> + +namespace Gudhi { + +namespace witness_complex { + +/** + \class Witness_complex + \brief Constructs the witness complex for the given set of witnesses and landmarks. + \ingroup witness_complex + */ +template< class Simplicial_complex> +class Strange_witness_complex { + private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } + }; + + private: + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + // typedef typename Simplicial_complex::Filtration_value Filtration_value; + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL; // Number of landmarks + Simplicial_complex& sc; // Simplicial complex + + public: + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @name Constructor + */ + + //@{ + + // Witness_range<Closest_landmark_range<Vertex_handle>> + + /** + * \brief Iterative construction of the witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn of + * nearest neighbours of the form {witnesses}x{landmarks}. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Constructor takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + template< typename KNearestNeighbors > + Strange_witness_complex(KNearestNeighbors const & knn, + Simplicial_complex & sc_, + int nbL_, + int dim) : nbL(nbL_), sc(sc_) { + // Construction of the active witness list + int nbW = knn.size(); + typeVectorVertex vv; + int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + for (int i = 0; i != nbL; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + counter++; + vv = {i}; + sc.insert_simplex(vv); + // TODO(SK) Error if not inserted : normally no need here though + } + for (int i = 0; i != nbW; ++i) { + std::vector<int> simplex; + simplex.reserve(dim+1); + for (int j = 0; j <= dim; j++) + simplex.push_back(knn[i][j]); + sc.insert_simplex_and_subfaces(simplex); + } + } + + //@} + + private: + /** \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + template <typename KNearestNeighbors> + bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) { + // std::cout << "All face in with the landmark " << inserted_vertex << std::endl; + std::vector< Vertex_handle > facet; + // Vertex_handle curr_vh = curr_sh->first; + // CHECK ALL THE FACETS + for (int i = 0; i != k + 1; ++i) { + facet = {}; + for (int j = 0; j != k + 1; ++j) { + if (j != i) { + facet.push_back(knn[witness_id][j]); + } + } // endfor + if (sc.find(facet) == sc.null_simplex()) + return false; + // std::cout << "++++ finished loop safely\n"; + } // endfor + return true; + } + + template <typename T> + static void print_vector(const std::vector<T>& v) { + std::cout << "["; + if (!v.empty()) { + std::cout << *(v.begin()); + for (auto it = v.begin() + 1; it != v.end(); ++it) { + std::cout << ","; + std::cout << *it; + } + } + std::cout << "]"; + } + + public: + /** + * \brief Verification if every simplex in the complex is witnessed by witnesses in knn. + * \param print_output =true will print the witnesses for each simplex + * \remark Added for debugging purposes. + */ + template< class KNearestNeighbors > + bool is_witness_complex(KNearestNeighbors const & knn, bool print_output) { + // bool final_result = true; + for (Simplex_handle sh : sc.complex_simplex_range()) { + bool is_witnessed = false; + typeVectorVertex simplex; + int nbV = 0; // number of verticed in the simplex + for (int v : sc.simplex_vertex_range(sh)) + simplex.push_back(v); + nbV = simplex.size(); + for (typeVectorVertex w : knn) { + bool has_vertices = true; + for (int v : simplex) + if (std::find(w.begin(), w.begin() + nbV, v) == w.begin() + nbV) { + has_vertices = false; + // break; + } + if (has_vertices) { + is_witnessed = true; + if (print_output) { + std::cout << "The simplex "; + print_vector(simplex); + std::cout << " is witnessed by the witness "; + print_vector(w); + std::cout << std::endl; + } + break; + } + } + if (!is_witnessed) { + if (print_output) { + std::cout << "The following simplex is not witnessed "; + print_vector(simplex); + std::cout << std::endl; + } + assert(is_witnessed); + return false; + } + } + return true; + } +}; + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // WITNESS_COMPLEX_H_ diff --git a/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h new file mode 100644 index 00000000..ee863a42 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h @@ -0,0 +1,337 @@ +/* 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 STRONG_RELAXED_WITNESS_COMPLEX_H_ +#define STRONG_RELAXED_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#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_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + +namespace witness_complex { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Strong_relaxed_witness_complex { + +private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } + }; + +private: + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + // typedef typename Simplicial_complex::Filtration_value Filtration_value; + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL; // Number of landmarks + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Constructor takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + + template< typename KNearestNeighbours > + Strong_relaxed_witness_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha2, + unsigned limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + typeVectorVertex vv; + //int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + ActiveWitnessList active_w;// = new ActiveWitnessList(); + for (int i = 0; i != nbL; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + //counter++; + vv = {i}; + sc.insert_simplex(vv, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ + } + for (int i=0; i != nbW; ++i) { + // int i_end = limD+1; + // if (knn[i].size() < limD+1) + // i_end = knn[i].size(); + // double dist_wL = *(distances[i].begin()); + // while (distances[i][i_end] > dist_wL + alpha2) + // i_end--; + // add_all_witnessed_faces(distances[i].begin(), + // knn[i].begin(), + // knn[i].begin() + i_end + 1); + unsigned j_end = 0; + while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) { + std::vector<int> simplex; + for (unsigned j = 0; j <= j_end; ++j) + simplex.push_back(knn[i][j]); + assert(distances[i][j_end] - distances[i][0] >= 0); + sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]); + j_end++; + } + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness + * Iterator is needed to know until how far we can take landmarks to form simplexes + * simplex is the prefix of the simplexes to insert + * The output value indicates if the witness rests active or not + */ + void add_all_witnessed_faces(std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>::const_iterator end) + { + std::vector<int> simplex; + std::vector<int>::const_iterator l_end = curr_l; + for (; l_end != end; ++l_end) { + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + simplex = {}; + for (; l_it != l_end; ++l_it, ++d_it) + simplex.push_back(*l_it); + sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d); + } + } + + /** \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Strong_relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h new file mode 100644 index 00000000..a1aedd8e --- /dev/null +++ b/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h @@ -0,0 +1,379 @@ +/* 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 WEAK_RELAXED_WITNESS_COMPLEX_H_ +#define WEAK_RELAXED_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#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_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + +namespace witness_complex { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Weak_relaxed_witness_complex { + +private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } + }; + +private: + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + // typedef typename Simplicial_complex::Filtration_value Filtration_value; + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL; // Number of landmarks + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Constructor takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + + template< typename KNearestNeighbours > + Weak_relaxed_witness_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha, + int limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + typeVectorVertex vv; + //int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + ActiveWitnessList active_w;// = new ActiveWitnessList(); + for (int i = 0; i != nbL; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + //counter++; + vv = {i}; + sc.insert_simplex(vv, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ + } + int k = 1; /* current dimension in iterative construction */ + for (int i=0; i != nbW; ++i) + active_w.push_back(i); + while (!active_w.empty() && k <= limD && k < nbL ) + { + typename ActiveWitnessList::iterator aw_it = active_w.begin(); + while (aw_it != active_w.end()) + { + std::vector<int> simplex; + bool ok = add_all_faces_of_dimension(k, + alpha, + std::numeric_limits<double>::infinity(), + distances[*aw_it].begin(), + knn[*aw_it].begin(), + simplex, + knn[*aw_it].end()); + if (!ok) + active_w.erase(aw_it++); //First increase the iterator and then erase the previous element + else + aw_it++; + } + k++; + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness + * Iterator is needed to know until how far we can take landmarks to form simplexes + * simplex is the prefix of the simplexes to insert + * The output value indicates if the witness rests active or not + */ + bool add_all_faces_of_dimension(int dim, + double alpha, + double norelax_dist, + std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>& simplex, + std::vector<int>::const_iterator end) + { + if (curr_l == end) + return false; + bool will_be_active = false; + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + if (dim > 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + if (sc.find(simplex) != sc.null_simplex()) + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it <= norelax_dist) + norelax_dist = *d_it; + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + } + else if (dim == 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + double filtration_value = 0; + // if norelax_dist is infinite, relaxation is 0. + if (*d_it > norelax_dist) + filtration_value = *d_it - norelax_dist; + if (all_faces_in(simplex, &filtration_value)) { + will_be_active = true; + sc.insert_simplex(simplex, filtration_value); + } + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it < norelax_dist) + norelax_dist = *d_it; + } + return will_be_active; + } + + /** \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Weak_relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 489cdf11..4c489e7a 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -31,6 +31,7 @@ #include <boost/range/size.hpp> #include <gudhi/distance_functions.h> +#include <gudhi/Kd_tree_search.h> #include <algorithm> #include <utility> @@ -42,8 +43,10 @@ #include <ctime> #include <iostream> -namespace Gudhi { +namespace gss = Gudhi::spatial_searching; +namespace Gudhi { + namespace witness_complex { // /* @@ -52,8 +55,19 @@ namespace witness_complex { // \brief Constructs the witness complex for the given set of witnesses and landmarks. // \ingroup witness_complex // */ -template< class SimplicialComplex> +template< class SimplicialComplex, + class Kernel_ > class Witness_complex { + typedef Kernel_ K; + typedef typename K::Point_d Point_d; + typedef typename K::FT FT; + typedef std::vector<Point_d> Point_range; + typedef gss::Kd_tree_search<Kernel_, Point_range> Kd_tree; + typedef typename Kd_tree::INS_range Nearest_landmark_range; + typedef typename std::vector<Nearest_landmark_range> Nearest_landmark_table; + typedef typename Nearest_landmark_table::const_iterator Nearest_landmark_table_iterator; + //typedef std::vector<std::pair<unsigned,FT>> Nearest_landmarks; + private: struct Active_witness { int witness_id; @@ -81,8 +95,11 @@ class Witness_complex { private: int nbL_; // Number of landmarks - SimplicialComplex& sc_; // Simplicial complex - + SimplicialComplex& sc_; // Simplicial complex + Point_range witnesses_, landmarks_; + Kd_tree landmark_tree_; + std::vector<Nearest_landmark_range> nearest_landmarks_; + public: ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /* @name Constructor @@ -93,7 +110,7 @@ class Witness_complex { // Witness_range<Closest_landmark_range<Vertex_handle>> /* - * \brief Iterative construction of the witness complex. + * \brief Iterative construction of the (weak) witness complex. * \details The witness complex is written in sc_ basing on a matrix knn of * nearest neighbours of the form {witnesses}x{landmarks}. * @@ -106,45 +123,58 @@ class Witness_complex { * * Landmarks are supposed to be in [0,nbL_-1] */ - template< typename KNearestNeighbors > - Witness_complex(KNearestNeighbors const & knn, - int nbL, - int dim, - SimplicialComplex & sc) : nbL_(nbL), sc_(sc) { - // Construction of the active witness list - int nbW = boost::size(knn); + template< typename InputIteratorLandmarks, + typename InputIteratorWitnesses > + Witness_complex(InputIteratorLandmarks landmarks_first, + InputIteratorLandmarks landmarks_last, + InputIteratorWitnesses witnesses_first, + InputIteratorWitnesses witnesses_last, + SimplicialComplex& sc) + : witnesses_(witnesses_first, witnesses_last), landmarks_(landmarks_first, landmarks_last), landmark_tree_(landmarks_), sc_(sc) + { + for (Point_d w : witnesses_) + nearest_landmarks_.push_back(landmark_tree_.query_incremental_nearest_neighbors(w)); + } + + Point_d get_point( std::size_t vertex ) const + { + return landmarks_[vertex]; + } + + template < typename SimplicialComplexForWitness > + bool create_complex(SimplicialComplexForWitness& complex, + FT max_alpha_square) const + { + unsigned nbW = witnesses_.size(), + nbL = landmarks_.size(); typeVectorVertex vv; - int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - ActiveWitnessList active_w; // = new ActiveWitnessList(); - for (Vertex_handle i = 0; i != nbL_; ++i) { + ActiveWitnessList active_w;// = new ActiveWitnessList(); + for (int i = 0; i != nbL; ++i) { // initial fill of 0-dimensional simplices // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore - counter++; + //counter++; vv = {i}; - sc_.insert_simplex(vv); - // TODO(SK) Error if not inserted : normally no need here though + complex.insert_simplex(vv, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ } - int k = 1; /* current dimension in iterative construction */ - for (int i = 0; i != nbW; ++i) + unsigned k = 1; /* current dimension in iterative construction */ + for (int i=0; i != nbW; ++i) active_w.push_back(i); - while (!active_w.empty() && k < dim) { - typename ActiveWitnessList::iterator it = active_w.begin(); - while (it != active_w.end()) { - typeVectorVertex simplex_vector; - /* THE INSERTION: Checking if all the subfaces are in the simplex tree*/ - bool ok = all_faces_in(knn, *it, k); - if (ok) { - for (int i = 0; i != k + 1; ++i) - simplex_vector.push_back(knn[*it][i]); - sc_.insert_simplex(simplex_vector); - // TODO(SK) Error if not inserted : normally no need here though - ++it; - } else { - active_w.erase(it++); // First increase the iterator and then erase the previous element - } + while (!active_w.empty() && k < nbL ) { + typename ActiveWitnessList::iterator aw_it = active_w.begin(); + while (aw_it != active_w.end()) { + std::vector<int> simplex; + bool ok = add_all_faces_of_dimension(k, + max_alpha_square, + std::numeric_limits<double>::infinity(), + nearest_landmarks_[*aw_it].begin(), + simplex, + complex, + nearest_landmarks_[*aw_it].end()); + if (!ok) + active_w.erase(aw_it++); //First increase the iterator and then erase the previous element + else + aw_it++; } k++; } @@ -153,40 +183,114 @@ class Witness_complex { //@} private: - /* \brief Check if the facets of the k-dimensional simplex witnessed - * by witness witness_id are already in the complex. - * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + + /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness + * Iterator is needed to know until how far we can take landmarks to form simplexes + * simplex is the prefix of the simplexes to insert + * The output value indicates if the witness rests active or not */ - template <typename KNearestNeighbors> - bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) { - std::vector< Vertex_handle > facet; - // CHECK ALL THE FACETS - for (int i = 0; i != k + 1; ++i) { - facet = {}; - for (int j = 0; j != k + 1; ++j) { - if (j != i) { - facet.push_back(knn[witness_id][j]); + template < typename SimplicialComplexForWitness > + bool add_all_faces_of_dimension(int dim, + double alpha2, + double norelax_dist2, + Nearest_landmark_table_iterator curr_l, + std::vector<int>& simplex, + SimplicialComplexForWitness& sc, + Nearest_landmark_table_iterator end) + { + if (curr_l == end) + return false; + bool will_be_active = false; + Nearest_landmark_table_iterator l_it = curr_l; + if (dim > 0) + for (; l_it->second - alpha2 <= norelax_dist2 && l_it != end; ++l_it) { + simplex.push_back(l_it->first); + if (sc.find(simplex) != sc.null_simplex()) + will_be_active = add_all_faces_of_dimension(dim-1, + alpha2, + norelax_dist2, + l_it+1, + simplex, + sc, + end) || will_be_active; + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (l_it->second <= norelax_dist2) + norelax_dist2 = l_it->second; + will_be_active = add_all_faces_of_dimension(dim-1, + alpha2, + norelax_dist2, + l_it+1, + simplex, + sc, + end) || will_be_active; + } + else if (dim == 0) + for (; l_it->second - alpha2 <= norelax_dist2 && l_it != end; ++l_it) { + simplex.push_back(l_it->second); + double filtration_value = 0; + // if norelax_dist is infinite, relaxation is 0. + if (l_it->second > norelax_dist2) + filtration_value = l_it->second - norelax_dist2; + if (all_faces_in(simplex, &filtration_value, sc)) { + will_be_active = true; + sc.insert_simplex(simplex, filtration_value); } - } // endfor - if (sc_.find(facet) == sc_.null_simplex()) - return false; - } // endfor - return true; + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (l_it->second < norelax_dist2) + norelax_dist2 = l_it->second; + } + return will_be_active; } - - template <typename T> - static void print_vector(const std::vector<T>& v) { - std::cout << "["; - if (!v.empty()) { - std::cout << *(v.begin()); - for (auto it = v.begin() + 1; it != v.end(); ++it) { - std::cout << ","; - std::cout << *it; + + /** \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + template < typename SimplicialComplexForWitness > + bool all_faces_in(std::vector<int>& simplex, + double* filtration_value, + SimplicialComplexForWitness& sc) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); } - } - std::cout << "]"; + return true; } + // bool is_face(Simplex_handle face, Simplex_handle coface) + // { + // // vertex range is sorted in decreasing order + // auto fvr = sc.simplex_vertex_range(face); + // auto cfvr = sc.simplex_vertex_range(coface); + // auto fv_it = fvr.begin(); + // auto cfv_it = cfvr.begin(); + // while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + // if (*fv_it < *cfv_it) + // ++cfv_it; + // else if (*fv_it == *cfv_it) { + // ++cfv_it; + // ++fv_it; + // } + // else + // return false; + + // } + // return (fv_it == fvr.end()); + // } + + public: // /* // * \brief Verification if every simplex in the complex is witnessed by witnesses in knn. @@ -250,13 +354,13 @@ class Witness_complex { * * Landmarks are supposed to be in [0,nbL_-1] */ - template <class KNearestNeighbors, class SimplicialComplexForWitness> - void witness_complex(KNearestNeighbors const & knn, - int nbL, - int dim, - SimplicialComplexForWitness & sc) { - Witness_complex<SimplicialComplexForWitness>(knn, nbL, dim, sc); - } + // template <class KNearestNeighbors, class SimplicialComplexForWitness> + // void witness_complex(KNearestNeighbors const & knn, + // int nbL, + // int dim, + // SimplicialComplexForWitness & sc) { + // Witness_complex<SimplicialComplexForWitness>(knn, nbL, dim, sc); + // } } // namespace witness_complex diff --git a/src/cmake/modules/GUDHI_user_version_target.txt b/src/cmake/modules/GUDHI_user_version_target.txt index 0ab36cfc..805f0a83 100644 --- a/src/cmake/modules/GUDHI_user_version_target.txt +++ b/src/cmake/modules/GUDHI_user_version_target.txt @@ -48,7 +48,7 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI) - set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Contraction;Hasse_complex;Persistent_cohomology;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Witness_complex") + set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Contraction;Hasse_complex;Persistent_cohomology;Simplex_tree;Skeleton_blocker;Witness_complex") foreach(GUDHI_MODULE ${GUDHI_MODULES}) # doc files diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 0983051d..21cf6925 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -315,8 +315,8 @@ make \endverbatim * @example Bitmap_cubical_complex/Bitmap_cubical_complex.cpp * @example Bitmap_cubical_complex/Bitmap_cubical_complex_periodic_boundary_conditions.cpp * @example Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp - * @example common/CGAL_3D_points_off_reader.cpp - * @example common/CGAL_points_off_reader.cpp + * @example common/example_CGAL_3D_points_off_reader.cpp + * @example common/example_CGAL_points_off_reader.cpp * @example Contraction/Garland_heckbert.cpp * @example Contraction/Rips_contraction.cpp * @example Persistent_cohomology/alpha_complex_3d_persistence.cpp diff --git a/src/common/include/gudhi/Points_3D_off_io.h b/src/common/include/gudhi/Points_3D_off_io.h index 2647f11e..b0d24998 100644 --- a/src/common/include/gudhi/Points_3D_off_io.h +++ b/src/common/include/gudhi/Points_3D_off_io.h @@ -132,12 +132,12 @@ class Points_3D_off_visitor_reader { * * @code template<class InputIterator > Point_3::Point_3(double x, double y, double z) @endcode * - * @section Example + * @section point3doffioexample Example * * This example loads points from an OFF file and builds a vector of CGAL points in dimension 3. * Then, it is asked to display the points. * - * @include common/CGAL_3D_points_off_reader.cpp + * @include common/example_CGAL_3D_points_off_reader.cpp * * When launching: * diff --git a/src/common/include/gudhi/Points_off_io.h b/src/common/include/gudhi/Points_off_io.h index 18b23e84..29af8a8a 100644 --- a/src/common/include/gudhi/Points_off_io.h +++ b/src/common/include/gudhi/Points_off_io.h @@ -114,7 +114,7 @@ class Points_off_visitor_reader { * * where d is the point dimension. * - * \section Example + * \section pointoffioexample Example * * This example loads points from an OFF file and builds a vector of points (vector of double). * Then, it is asked to display the points. |