diff options
14 files changed, 620 insertions, 227 deletions
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index f82e8ce9..969daba6 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -35,6 +35,7 @@ #include <algorithm> // for sort #include <vector> #include <numeric> // for iota +#include <cstddef> namespace Gudhi { @@ -61,7 +62,7 @@ class Bitmap_cubical_complex : public T { //*********************************************// // Typedefs and typenames //*********************************************// - typedef size_t Simplex_key; + typedef std::size_t Simplex_key; typedef typename T::filtration_type Filtration_value; typedef Simplex_key Simplex_handle; @@ -82,7 +83,7 @@ class Bitmap_cubical_complex : public T { if (globalDbg) { std::cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n"; } - for (size_t i = 0; i != this->total_number_of_cells; ++i) { + for (std::size_t i = 0; i != this->total_number_of_cells; ++i) { this->key_associated_to_simplex[i] = i; } // we initialize this only once, in each constructor, when the bitmap is constructed. @@ -99,7 +100,7 @@ class Bitmap_cubical_complex : public T { Bitmap_cubical_complex(const std::vector<unsigned>& dimensions, const std::vector<Filtration_value>& top_dimensional_cells) : T(dimensions, top_dimensional_cells), key_associated_to_simplex(this->total_number_of_cells + 1) { - for (size_t i = 0; i != this->total_number_of_cells; ++i) { + for (std::size_t i = 0; i != this->total_number_of_cells; ++i) { this->key_associated_to_simplex[i] = i; } // we initialize this only once, in each constructor, when the bitmap is constructed. @@ -120,7 +121,7 @@ class Bitmap_cubical_complex : public T { std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed) : T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed), key_associated_to_simplex(this->total_number_of_cells + 1) { - for (size_t i = 0; i != this->total_number_of_cells; ++i) { + for (std::size_t i = 0; i != this->total_number_of_cells; ++i) { this->key_associated_to_simplex[i] = i; } // we initialize this only once, in each constructor, when the bitmap is constructed. @@ -141,7 +142,7 @@ class Bitmap_cubical_complex : public T { /** * Returns number of all cubes in the complex. **/ - size_t num_simplices() const { return this->total_number_of_cells; } + std::size_t num_simplices() const { return this->total_number_of_cells; } /** * Returns a Simplex_handle to a cube that do not exist in this complex. @@ -156,7 +157,7 @@ class Bitmap_cubical_complex : public T { /** * Returns dimension of the complex. **/ - inline size_t dimension() const { return this->sizes.size(); } + inline std::size_t dimension() const { return this->sizes.size(); } /** * Return dimension of a cell pointed by the Simplex_handle. @@ -308,7 +309,7 @@ class Bitmap_cubical_complex : public T { private: Bitmap_cubical_complex<T>* b; - size_t position; + std::size_t position; }; /** @@ -379,7 +380,7 @@ class Bitmap_cubical_complex : public T { * Function needed for compatibility with Gudhi. Not useful for other purposes. **/ std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) { - std::vector<size_t> bdry = this->get_boundary_of_a_cell(sh); + std::vector<std::size_t> bdry = this->get_boundary_of_a_cell(sh); if (globalDbg) { std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n"; std::cerr << "bdry.size() : " << bdry.size() << std::endl; @@ -401,9 +402,9 @@ class Bitmap_cubical_complex : public T { // Iterator over all simplices of the complex in the order of the indexing scheme. // 'value_type' must be 'Simplex_handle'. public: - Skeleton_simplex_iterator(Bitmap_cubical_complex* b, size_t d) : b(b), dimension(d) { + Skeleton_simplex_iterator(Bitmap_cubical_complex* b, std::size_t d) : b(b), dimension(d) { if (globalDbg) { - std::cerr << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , size_t d )\n"; + std::cerr << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , std::size_t d )\n"; } // find the position of the first simplex of a dimension d this->position = 0; @@ -469,7 +470,7 @@ class Bitmap_cubical_complex : public T { private: Bitmap_cubical_complex<T>* b; - size_t position; + std::size_t position; unsigned dimension; }; @@ -519,8 +520,8 @@ class Bitmap_cubical_complex : public T { friend class is_before_in_filtration<T>; protected: - std::vector<size_t> key_associated_to_simplex; - std::vector<size_t> simplex_associated_to_key; + std::vector<std::size_t> key_associated_to_simplex; + std::vector<std::size_t> simplex_associated_to_key; }; // Bitmap_cubical_complex template <typename T> @@ -528,7 +529,7 @@ void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key() { if (globalDbg) { std::cerr << "void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration() \n"; } - this->simplex_associated_to_key = std::vector<size_t>(this->data.size()); + this->simplex_associated_to_key = std::vector<std::size_t>(this->data.size()); std::iota(std::begin(simplex_associated_to_key), std::end(simplex_associated_to_key), 0); #ifdef GUDHI_USE_TBB tbb::parallel_sort(simplex_associated_to_key.begin(), simplex_associated_to_key.end(), @@ -538,7 +539,7 @@ void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key() { #endif // we still need to deal here with a key_associated_to_simplex: - for (size_t i = 0; i != simplex_associated_to_key.size(); ++i) { + for (std::size_t i = 0; i != simplex_associated_to_key.size(); ++i) { this->key_associated_to_simplex[simplex_associated_to_key[i]] = i; } } @@ -558,8 +559,8 @@ class is_before_in_filtration { return fil1 < fil2; } // in this case they are on the same filtration level, so the dimension decide. - size_t dim1 = CC_->get_dimension_of_a_cell(sh1); - size_t dim2 = CC_->get_dimension_of_a_cell(sh2); + std::size_t dim1 = CC_->get_dimension_of_a_cell(sh1); + std::size_t dim2 = CC_->get_dimension_of_a_cell(sh2); if (dim1 != dim2) { return dim1 < dim2; } diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h index 4b072f10..705b68a0 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h @@ -25,6 +25,7 @@ #include <iostream> #include <vector> +#include <cstddef> namespace Gudhi { @@ -63,14 +64,14 @@ class counter { * If the value of the function is false, that means, that the counter have reached its end-value. **/ bool increment() { - size_t i = 0; + std::size_t i = 0; while ((i != this->end.size()) && (this->current[i] == this->end[i])) { ++i; } if (i == this->end.size())return false; ++this->current[i]; - for (size_t j = 0; j != i; ++j) { + for (std::size_t j = 0; j != i; ++j) { this->current[j] = this->begin[j]; } return true; @@ -80,7 +81,7 @@ class counter { * Function to check if we are at the end of counter. **/ bool isFinal() { - for (size_t i = 0; i != this->current.size(); ++i) { + for (std::size_t i = 0; i != this->current.size(); ++i) { if (this->current[i] == this->end[i])return true; } return false; @@ -93,7 +94,7 @@ class counter { **/ std::vector< unsigned > find_opposite(const std::vector< bool >& directionsForPeriodicBCond) { std::vector< unsigned > result; - for (size_t i = 0; i != this->current.size(); ++i) { + for (std::size_t i = 0; i != this->current.size(); ++i) { if ((this->current[i] == this->end[i]) && (directionsForPeriodicBCond[i] == true)) { result.push_back(this->begin[i]); } else { @@ -108,7 +109,7 @@ class counter { **/ std::vector< bool > directions_of_finals() { std::vector< bool > result; - for (size_t i = 0; i != this->current.size(); ++i) { + for (std::size_t i = 0; i != this->current.size(); ++i) { if (this->current[i] == this->end[i]) { result.push_back(true); } else { @@ -123,7 +124,7 @@ class counter { **/ friend std::ostream& operator<<(std::ostream& out, const counter& c) { // std::cerr << "c.current.size() : " << c.current.size() << endl; - for (size_t i = 0; i != c.current.size(); ++i) { + for (std::size_t i = 0; i != c.current.size(); ++i) { out << c.current[i] << " "; } return out; diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 4adadce6..bf257be1 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -34,6 +34,7 @@ #include <limits> #include <utility> #include <stdexcept> +#include <cstddef> namespace Gudhi { @@ -103,7 +104,7 @@ class Bitmap_cubical_complex_base { * The boundary elements are guaranteed to be returned so that the * incidence coefficients of boundary elements are alternating. */ - virtual inline std::vector<size_t> get_boundary_of_a_cell(size_t cell) const; + virtual inline std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const; /** * The functions get_coboundary_of_a_cell, get_coboundary_of_a_cell, * get_dimension_of_a_cell and get_cell_data are the basic @@ -118,7 +119,7 @@ class Bitmap_cubical_complex_base { * not guaranteed to be returned with alternating incidence numbers. * **/ - virtual inline std::vector<size_t> get_coboundary_of_a_cell(size_t cell) const; + virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const; /** * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of @@ -139,15 +140,15 @@ class Bitmap_cubical_complex_base { * @exception std::logic_error In case when the cube \f$B\f$ is not n-1 * dimensional face of a cube \f$A\f$. **/ - virtual int compute_incidence_between_cells(size_t coface, size_t face) const { + virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) const { // first get the counters for coface and face: std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface); std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face); // coface_counter and face_counter should agree at all positions except from one: int number_of_position_in_which_counters_do_not_agree = -1; - size_t number_of_full_faces_that_comes_before = 0; - for (size_t i = 0; i != coface_counter.size(); ++i) { + std::size_t number_of_full_faces_that_comes_before = 0; + for (std::size_t i = 0; i != coface_counter.size(); ++i) { if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) { ++number_of_full_faces_that_comes_before; } @@ -180,7 +181,7 @@ class Bitmap_cubical_complex_base { * To compute incidence between cells use compute_incidence_between_cells * procedure **/ - inline unsigned get_dimension_of_a_cell(size_t cell) const; + inline unsigned get_dimension_of_a_cell(std::size_t cell) const; /** * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position. @@ -188,7 +189,7 @@ class Bitmap_cubical_complex_base { * code do not check if we have a filtration or not. i.e. it do not check if the value of a filtration of a cell is * not smaller than the value of a filtration of its boundary and not greater than the value of its coboundary. **/ - inline T& get_cell_data(size_t cell); + inline T& get_cell_data(std::size_t cell); /** * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells. @@ -222,10 +223,10 @@ class Bitmap_cubical_complex_base { * equally distributed in the range of data. * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute * persistence gets worst. When dealing with this type of data, one may want to put different values on cells to - * some number of bins. The function put_data_to_bins( size_t number_of_bins ) is designed for that purpose. + * some number of bins. The function put_data_to_bins( std::size_t number_of_bins ) is designed for that purpose. * The parameter of the function is the number of bins (distinct values) we want to have in the cubical complex. **/ - void put_data_to_bins(size_t number_of_bins); + void put_data_to_bins(std::size_t number_of_bins); /** * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values @@ -285,11 +286,11 @@ class Bitmap_cubical_complex_base { * boundary and coboundary and dimension * and in function get_cell_data to get a filtration of a cell. */ - size_t operator*() { return this->counter; } + std::size_t operator*() { return this->counter; } friend class Bitmap_cubical_complex_base; protected: - size_t counter; + std::size_t counter; }; /** @@ -329,26 +330,26 @@ class Bitmap_cubical_complex_base { /** * Boundary_range class provides ranges for boundary iterators. **/ - typedef typename std::vector<size_t>::const_iterator Boundary_iterator; - typedef typename std::vector<size_t> Boundary_range; + typedef typename std::vector<std::size_t>::const_iterator Boundary_iterator; + typedef typename std::vector<std::size_t> Boundary_range; /** * boundary_simplex_range creates an object of a Boundary_simplex_range class * that provides ranges for the Boundary_simplex_iterator. **/ - Boundary_range boundary_range(size_t sh) { return this->get_boundary_of_a_cell(sh); } + Boundary_range boundary_range(std::size_t sh) { return this->get_boundary_of_a_cell(sh); } /** * Coboundary_range class provides ranges for boundary iterators. **/ - typedef typename std::vector<size_t>::const_iterator Coboundary_iterator; - typedef typename std::vector<size_t> Coboundary_range; + typedef typename std::vector<std::size_t>::const_iterator Coboundary_iterator; + typedef typename std::vector<std::size_t> Coboundary_range; /** * boundary_simplex_range creates an object of a Boundary_simplex_range class * that provides ranges for the Boundary_simplex_iterator. **/ - Coboundary_range coboundary_range(size_t sh) { return this->get_coboundary_of_a_cell(sh); } + Coboundary_range coboundary_range(std::size_t sh) { return this->get_coboundary_of_a_cell(sh); } /** * @brief Iterator through top dimensional cells of the complex. The cells appear in order they are stored @@ -357,18 +358,18 @@ class Bitmap_cubical_complex_base { class Top_dimensional_cells_iterator : std::iterator<std::input_iterator_tag, T> { public: Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) { - this->counter = std::vector<size_t>(b.dimension()); + this->counter = std::vector<std::size_t>(b.dimension()); // std::fill( this->counter.begin() , this->counter.end() , 0 ); } Top_dimensional_cells_iterator operator++() { // first find first element of the counter that can be increased: - size_t dim = 0; + std::size_t dim = 0; while ((dim != this->b.dimension()) && (this->counter[dim] == this->b.sizes[dim] - 1)) ++dim; if (dim != this->b.dimension()) { ++this->counter[dim]; - for (size_t i = 0; i != dim; ++i) { + for (std::size_t i = 0; i != dim; ++i) { this->counter[i] = 0; } } else { @@ -392,7 +393,7 @@ class Bitmap_cubical_complex_base { bool operator==(const Top_dimensional_cells_iterator& rhs) const { if (&this->b != &rhs.b) return false; if (this->counter.size() != rhs.counter.size()) return false; - for (size_t i = 0; i != this->counter.size(); ++i) { + for (std::size_t i = 0; i != this->counter.size(); ++i) { if (this->counter[i] != rhs.counter[i]) return false; } return true; @@ -407,25 +408,25 @@ class Bitmap_cubical_complex_base { * boundary and coboundary and dimension * and in function get_cell_data to get a filtration of a cell. */ - size_t operator*() { return this->compute_index_in_bitmap(); } + std::size_t operator*() { return this->compute_index_in_bitmap(); } - size_t compute_index_in_bitmap() const { - size_t index = 0; - for (size_t i = 0; i != this->counter.size(); ++i) { + std::size_t compute_index_in_bitmap() const { + std::size_t index = 0; + for (std::size_t i = 0; i != this->counter.size(); ++i) { index += (2 * this->counter[i] + 1) * this->b.multipliers[i]; } return index; } void print_counter() const { - for (size_t i = 0; i != this->counter.size(); ++i) { + for (std::size_t i = 0; i != this->counter.size(); ++i) { std::cout << this->counter[i] << " "; } } friend class Bitmap_cubical_complex_base; protected: - std::vector<size_t> counter; + std::vector<std::size_t> counter; Bitmap_cubical_complex_base& b; }; @@ -442,7 +443,7 @@ class Bitmap_cubical_complex_base { **/ Top_dimensional_cells_iterator top_dimensional_cells_iterator_end() { Top_dimensional_cells_iterator a(*this); - for (size_t i = 0; i != this->dimension(); ++i) { + for (std::size_t i = 0; i != this->dimension(); ++i) { a.counter[i] = this->sizes[i] - 1; } a.counter[0]++; @@ -471,7 +472,7 @@ class Bitmap_cubical_complex_base { //****************************************************************************************************************// //****************************************************************************************************************// - inline size_t number_cells() const { return this->total_number_of_cells; } + inline std::size_t number_cells() const { return this->total_number_of_cells; } //****************************************************************************************************************// //****************************************************************************************************************// @@ -482,11 +483,11 @@ class Bitmap_cubical_complex_base { std::vector<unsigned> sizes; std::vector<unsigned> multipliers; std::vector<T> data; - size_t total_number_of_cells; + std::size_t total_number_of_cells; void set_up_containers(const std::vector<unsigned>& sizes) { unsigned multiplier = 1; - for (size_t i = 0; i != sizes.size(); ++i) { + for (std::size_t i = 0; i != sizes.size(); ++i) { this->sizes.push_back(sizes[i]); this->multipliers.push_back(multiplier); multiplier *= 2 * sizes[i] + 1; @@ -495,18 +496,18 @@ class Bitmap_cubical_complex_base { this->total_number_of_cells = multiplier; } - size_t compute_position_in_bitmap(const std::vector<unsigned>& counter) { - size_t position = 0; - for (size_t i = 0; i != this->multipliers.size(); ++i) { + std::size_t compute_position_in_bitmap(const std::vector<unsigned>& counter) { + std::size_t position = 0; + for (std::size_t i = 0; i != this->multipliers.size(); ++i) { position += this->multipliers[i] * counter[i]; } return position; } - std::vector<unsigned> compute_counter_for_given_cell(size_t cell) const { + std::vector<unsigned> compute_counter_for_given_cell(std::size_t cell) const { std::vector<unsigned> counter; counter.reserve(this->sizes.size()); - for (size_t dim = this->sizes.size(); dim != 0; --dim) { + for (std::size_t dim = this->sizes.size(); dim != 0; --dim) { counter.push_back(cell / this->multipliers[dim - 1]); cell = cell % this->multipliers[dim - 1]; } @@ -523,40 +524,38 @@ class Bitmap_cubical_complex_base { }; template <typename T> -void Bitmap_cubical_complex_base<T>::put_data_to_bins(size_t number_of_bins) { - bool bdg = false; +void Bitmap_cubical_complex_base<T>::put_data_to_bins(std::size_t number_of_bins) { + bool dbg = false; std::pair<T, T> min_max = this->min_max_filtration(); T dx = (min_max.second - min_max.first) / (T)number_of_bins; // now put the data into the appropriate bins: - for (size_t i = 0; i != this->data.size(); ++i) { - if (bdg) { + for (std::size_t i = 0; i != this->data.size(); ++i) { + if (dbg) { std::cerr << "Before binning : " << this->data[i] << std::endl; } this->data[i] = min_max.first + dx * (this->data[i] - min_max.first) / number_of_bins; - if (bdg) { + if (dbg) { std::cerr << "After binning : " << this->data[i] << std::endl; - getchar(); } } } template <typename T> void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) { - bool bdg = false; + bool dbg = false; std::pair<T, T> min_max = this->min_max_filtration(); - size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin; + std::size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin; // now put the data into the appropriate bins: - for (size_t i = 0; i != this->data.size(); ++i) { - if (bdg) { + for (std::size_t i = 0; i != this->data.size(); ++i) { + if (dbg) { std::cerr << "Before binning : " << this->data[i] << std::endl; } this->data[i] = min_max.first + diameter_of_bin * (this->data[i] - min_max.first) / number_of_bins; - if (bdg) { + if (dbg) { std::cerr << "After binning : " << this->data[i] << std::endl; - getchar(); } } } @@ -564,7 +563,7 @@ void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) { template <typename T> std::pair<T, T> Bitmap_cubical_complex_base<T>::min_max_filtration() { std::pair<T, T> min_max(std::numeric_limits<T>::max(), std::numeric_limits<T>::min()); - for (size_t i = 0; i != this->data.size(); ++i) { + for (std::size_t i = 0; i != this->data.size(); ++i) { if (this->data[i] < min_max.first) min_max.first = this->data[i]; if (this->data[i] > min_max.second) min_max.second = this->data[i]; } @@ -590,23 +589,24 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) { this->set_up_containers(sizes_in_following_directions); - size_t number_of_top_dimensional_elements = 1; - for (size_t i = 0; i != sizes_in_following_directions.size(); ++i) { + std::size_t number_of_top_dimensional_elements = 1; + for (std::size_t i = 0; i != sizes_in_following_directions.size(); ++i) { number_of_top_dimensional_elements *= sizes_in_following_directions[i]; } if (number_of_top_dimensional_elements != top_dimensional_cells.size()) { - std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<size_t> sizes_in_following_directions" - << ", std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " - << "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector." + std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<std::size_t> " + << "sizes_in_following_directions, std::vector<T> top_dimensional_cells ). Number of top dimensional " + << "elements that follow from sizes_in_following_directions vector is different than the size of " + << "top_dimensional_cells vector." << std::endl; throw( - "Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions," + "Error in constructor Bitmap_cubical_complex_base( std::vector<std::size_t> sizes_in_following_directions," "std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector."); } Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); - size_t index = 0; + std::size_t index = 0; for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) { this->get_cell_data(*it) = top_dimensional_cells[index]; ++index; @@ -630,15 +630,17 @@ void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus if (dbg) { std::cerr << "dimensionOfData : " << dimensionOfData << std::endl; - getchar(); } std::vector<unsigned> sizes; sizes.reserve(dimensionOfData); - for (size_t i = 0; i != dimensionOfData; ++i) { + // all dimensions multiplied + std::size_t dimensions = 1; + for (std::size_t i = 0; i != dimensionOfData; ++i) { unsigned size_in_this_dimension; inFiltration >> size_in_this_dimension; sizes.push_back(size_in_this_dimension); + dimensions *= size_in_this_dimension; if (dbg) { std::cerr << "size_in_this_dimension : " << size_in_this_dimension << std::endl; } @@ -648,9 +650,11 @@ void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); it = this->top_dimensional_cells_iterator_begin(); - while (!inFiltration.eof()) { - T filtrationLevel; - inFiltration >> filtrationLevel; + T filtrationLevel; + for (std::size_t i = 0; i < dimensions; ++i) { + if (!(inFiltration >> filtrationLevel) || (inFiltration.eof())) { + throw std::ios_base::failure("Bad Perseus file format."); + } if (dbg) { std::cerr << "Cell of an index : " << it.compute_index_in_bitmap() << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) @@ -659,6 +663,7 @@ void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus this->get_cell_data(*it) = filtrationLevel; ++it; } + inFiltration.close(); this->impose_lower_star_filtration(); } @@ -697,15 +702,15 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_ } template <typename T> -std::vector<size_t> Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell) const { - std::vector<size_t> boundary_elements; +std::vector<std::size_t> Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(std::size_t cell) const { + std::vector<std::size_t> boundary_elements; // Speed traded of for memory. Check if it is better in practice. boundary_elements.reserve(this->dimension() * 2); - size_t sum_of_dimensions = 0; - size_t cell1 = cell; - for (size_t i = this->multipliers.size(); i != 0; --i) { + std::size_t sum_of_dimensions = 0; + std::size_t cell1 = cell; + for (std::size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; if (position % 2 == 1) { if (sum_of_dimensions % 2) { @@ -724,11 +729,11 @@ std::vector<size_t> Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_ } template <typename T> -std::vector<size_t> Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell) const { +std::vector<std::size_t> Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(std::size_t cell) const { std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell); - std::vector<size_t> coboundary_elements; - size_t cell1 = cell; - for (size_t i = this->multipliers.size(); i != 0; --i) { + std::vector<std::size_t> coboundary_elements; + std::size_t cell1 = cell; + for (std::size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; if (position % 2 == 0) { if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) { @@ -744,11 +749,11 @@ std::vector<size_t> Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(siz } template <typename T> -unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell) const { +unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(std::size_t cell) const { bool dbg = false; if (dbg) std::cerr << "\n\n\n Computing position o a cell of an index : " << cell << std::endl; unsigned dimension = 0; - for (size_t i = this->multipliers.size(); i != 0; --i) { + for (std::size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell / this->multipliers[i - 1]; if (dbg) { @@ -756,7 +761,6 @@ unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell) co std::cerr << "cell : " << cell << std::endl; std::cerr << "position : " << position << std::endl; std::cerr << "multipliers[" << i - 1 << "] = " << this->multipliers[i - 1] << std::endl; - getchar(); } if (position % 2 == 1) { @@ -769,7 +773,7 @@ unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell) co } template <typename T> -inline T& Bitmap_cubical_complex_base<T>::get_cell_data(size_t cell) { +inline T& Bitmap_cubical_complex_base<T>::get_cell_data(std::size_t cell) { return this->data[cell]; } @@ -780,12 +784,12 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { // this vector will be used to check which elements have already been taken care of in imposing lower star filtration std::vector<bool> is_this_cell_considered(this->data.size(), false); - size_t size_to_reserve = 1; - for (size_t i = 0; i != this->multipliers.size(); ++i) { - size_to_reserve *= (size_t)((this->multipliers[i] - 1) / 2); + std::size_t size_to_reserve = 1; + for (std::size_t i = 0; i != this->multipliers.size(); ++i) { + size_to_reserve *= (std::size_t)((this->multipliers[i] - 1) / 2); } - std::vector<size_t> indices_to_consider; + std::vector<std::size_t> indices_to_consider; indices_to_consider.reserve(size_to_reserve); // we assume here that we already have a filtration on the top dimensional cells and // we have to extend it to lower ones. @@ -797,27 +801,24 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { while (indices_to_consider.size()) { if (dbg) { std::cerr << "indices_to_consider in this iteration \n"; - for (size_t i = 0; i != indices_to_consider.size(); ++i) { + for (std::size_t i = 0; i != indices_to_consider.size(); ++i) { std::cout << indices_to_consider[i] << " "; } - getchar(); } - std::vector<size_t> new_indices_to_consider; - for (size_t i = 0; i != indices_to_consider.size(); ++i) { - std::vector<size_t> bd = this->get_boundary_of_a_cell(indices_to_consider[i]); - for (size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) { + std::vector<std::size_t> new_indices_to_consider; + for (std::size_t i = 0; i != indices_to_consider.size(); ++i) { + std::vector<std::size_t> bd = this->get_boundary_of_a_cell(indices_to_consider[i]); + for (std::size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) { if (dbg) { std::cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[bd[boundaryIt]] << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[indices_to_consider[i]] << std::endl; - getchar(); } if (this->data[bd[boundaryIt]] > this->data[indices_to_consider[i]]) { this->data[bd[boundaryIt]] = this->data[indices_to_consider[i]]; if (dbg) { std::cerr << "Setting the value of a cell : " << bd[boundaryIt] << " to : " << this->data[indices_to_consider[i]] << std::endl; - getchar(); } } if (is_this_cell_considered[bd[boundaryIt]] == false) { @@ -831,8 +832,8 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { } template <typename T> -bool compareFirstElementsOfTuples(const std::pair<std::pair<T, size_t>, char>& first, - const std::pair<std::pair<T, size_t>, char>& second) { +bool compareFirstElementsOfTuples(const std::pair<std::pair<T, std::size_t>, char>& first, + const std::pair<std::pair<T, std::size_t>, char>& second) { if (first.first.first < second.first.first) { return true; } else { diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h index e2f86f3b..4a0d1c74 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h @@ -29,6 +29,7 @@ #include <limits> // for numeric_limits<> #include <vector> #include <stdexcept> +#include <cstddef> namespace Gudhi { @@ -94,7 +95,7 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c * The boundary elements are guaranteed to be returned so that the * incidence coefficients are alternating. */ - virtual std::vector<size_t> get_boundary_of_a_cell(size_t cell) const; + virtual std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const; /** * A version of a function that return coboundary of a given cell for an object of @@ -104,7 +105,7 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c * To compute incidence between cells use compute_incidence_between_cells * procedure */ - virtual std::vector<size_t> get_coboundary_of_a_cell(size_t cell) const; + virtual std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const; /** * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of @@ -125,15 +126,15 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c * @exception std::logic_error In case when the cube \f$B\f$ is not n-1 * dimensional face of a cube \f$A\f$. **/ - virtual int compute_incidence_between_cells(size_t coface, size_t face) { + virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) { // first get the counters for coface and face: std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface); std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face); // coface_counter and face_counter should agree at all positions except from one: int number_of_position_in_which_counters_do_not_agree = -1; - size_t number_of_full_faces_that_comes_before = 0; - for (size_t i = 0; i != coface_counter.size(); ++i) { + std::size_t number_of_full_faces_that_comes_before = 0; + for (std::size_t i = 0; i != coface_counter.size(); ++i) { if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) { ++number_of_full_faces_that_comes_before; } @@ -165,7 +166,7 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c void set_up_containers(const std::vector<unsigned>& sizes) { unsigned multiplier = 1; - for (size_t i = 0; i != sizes.size(); ++i) { + for (std::size_t i = 0; i != sizes.size(); ++i) { this->sizes.push_back(sizes[i]); this->multipliers.push_back(multiplier); @@ -198,7 +199,7 @@ void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_comp this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed; this->set_up_containers(dimensions); - size_t i = 0; + std::size_t i = 0; for (auto it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) { this->get_cell_data(*it) = topDimensionalCells[i]; ++i; @@ -229,7 +230,7 @@ Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_comp std::vector<unsigned> sizes; sizes.reserve(dimensionOfData); - for (size_t i = 0; i != dimensionOfData; ++i) { + for (std::size_t i = 0; i != dimensionOfData; ++i) { int size_in_this_dimension; inFiltration >> size_in_this_dimension; if (size_in_this_dimension < 0) { @@ -285,19 +286,19 @@ Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_comp // ***********************Methods************************ // template <typename T> -std::vector<size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell( - size_t cell) const { +std::vector<std::size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell( + std::size_t cell) const { bool dbg = false; if (dbg) { std::cerr << "Computations of boundary of a cell : " << cell << std::endl; } - std::vector<size_t> boundary_elements; + std::vector<std::size_t> boundary_elements; boundary_elements.reserve(this->dimension() * 2); - size_t cell1 = cell; - size_t sum_of_dimensions = 0; + std::size_t cell1 = cell; + std::size_t sum_of_dimensions = 0; - for (size_t i = this->multipliers.size(); i != 0; --i) { + for (std::size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction. if (position % 2 == 1) { @@ -351,12 +352,12 @@ std::vector<size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>: } template <typename T> -std::vector<size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell( - size_t cell) const { +std::vector<std::size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell( + std::size_t cell) const { std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell); - std::vector<size_t> coboundary_elements; - size_t cell1 = cell; - for (size_t i = this->multipliers.size(); i != 0; --i) { + std::vector<std::size_t> coboundary_elements; + std::size_t cell1 = cell; + for (std::size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; // if the cell has zero length in this direction, then it will have cbd in this direction. if (position % 2 == 0) { diff --git a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp index 554eeba6..c1de0ef8 100644 --- a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp @@ -45,14 +45,7 @@ using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; using Vertex_handle = Simplex_tree::Vertex_handle; using Filtration_value = Simplex_tree::Filtration_value; -using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS -, boost::property < vertex_filtration_t, Filtration_value > -, boost::property < edge_filtration_t, Filtration_value > ->; -using Edge_t = std::pair< Vertex_handle, Vertex_handle >; - -template< typename InputPointRange, typename Distance > -Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance); +using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; @@ -81,8 +74,9 @@ int main(int argc, char * argv[]) { Points_off_reader off_reader(off_file_points); // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold - , Gudhi::Euclidean_distance()); + Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), + threshold, + Gudhi::Euclidean_distance()); // Construct the Rips complex in a Simplex Tree Simplex_tree st; @@ -170,48 +164,3 @@ void program_options(int argc, char * argv[] std::abort(); } } - -/** Output the proximity graph of the points. - * - * If points contains n elements, the proximity graph is the graph - * with n vertices, and an edge [u,v] iff the distance function between - * points u and v is smaller than threshold. - * - * The type PointCloud furnishes .begin() and .end() methods, that return - * iterators with value_type Point. - */ -template< typename InputPointRange, typename Distance > -Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance) { - std::vector< Edge_t > edges; - std::vector< Filtration_value > edges_fil; - - Vertex_handle idx_u, idx_v; - Filtration_value fil; - idx_u = 0; - for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { - idx_v = idx_u + 1; - for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) { - fil = distance(*it_u, *it_v); - if (fil <= threshold) { - edges.emplace_back(idx_u, idx_v); - edges_fil.push_back(fil); - } - } - ++idx_u; - } - - Graph_t skel_graph(edges.begin() - , edges.end() - , edges_fil.begin() - , idx_u); // number of points labeled from 0 to idx_u-1 - - auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph); - - boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; - for (std::tie(vi, vi_end) = boost::vertices(skel_graph); - vi != vi_end; ++vi) { - boost::put(vertex_prop, *vi, 0.); - } - - return skel_graph; -} diff --git a/src/Simplex_tree/example/CMakeLists.txt b/src/Simplex_tree/example/CMakeLists.txt index 8bc4ad53..b33b2d05 100644 --- a/src/Simplex_tree/example/CMakeLists.txt +++ b/src/Simplex_tree/example/CMakeLists.txt @@ -34,6 +34,16 @@ if(GMP_FOUND AND CGAL_FOUND) "${CMAKE_SOURCE_DIR}/data/points/bunny_5000.off") install(TARGETS Simplex_tree_example_alpha_shapes_3_from_off DESTINATION bin) + + add_executable ( Simplex_tree_example_cech_complex_cgal_mini_sphere_3d cech_complex_cgal_mini_sphere_3d.cpp ) + target_link_libraries(Simplex_tree_example_cech_complex_cgal_mini_sphere_3d ${Boost_PROGRAM_OPTIONS_LIBRARY} ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Simplex_tree_example_cech_complex_cgal_mini_sphere_3d ${TBB_LIBRARIES}) + endif() + add_test(NAME Simplex_tree_example_cech_complex_cgal_mini_sphere_3d COMMAND $<TARGET_FILE:Simplex_tree_example_cech_complex_cgal_mini_sphere_3d> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" -r 0.3 -d 3) + + install(TARGETS Simplex_tree_example_alpha_shapes_3_from_off DESTINATION bin) endif() add_executable ( Simplex_tree_example_graph_expansion_with_blocker graph_expansion_with_blocker.cpp ) diff --git a/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp new file mode 100644 index 00000000..217e251f --- /dev/null +++ b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp @@ -0,0 +1,234 @@ +/* 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): Clément Maria + * + * Copyright (C) 2014 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/graph_simplicial_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Points_off_io.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Min_sphere_of_spheres_d.h> +#include <CGAL/Min_sphere_of_points_d_traits_d.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <utility> // for pair +#include <map> + +// ------------------------------------------------------------------------------- +// cech_complex_cgal_mini_sphere_3d is an example of each step that is required to +// build a Cech over a Simplex_tree. Please refer to cech_persistence to see +// how to do the same thing with the Cech_complex wrapper for less detailed +// steps. +// ------------------------------------------------------------------------------- + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Vertex_handle = Simplex_tree::Vertex_handle; +using Simplex_handle = Simplex_tree::Simplex_handle; +using Filtration_value = Simplex_tree::Filtration_value; +using Siblings = Simplex_tree::Siblings; +using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS +, boost::property < Gudhi::vertex_filtration_t, Filtration_value > +, boost::property < Gudhi::edge_filtration_t, Filtration_value > +>; +using Edge_t = std::pair< Vertex_handle, Vertex_handle >; + +using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >; +using Point = Kernel::Point_d; +using Traits = CGAL::Min_sphere_of_points_d_traits_d<Kernel,Filtration_value,3>; +using Min_sphere = CGAL::Min_sphere_of_spheres_d<Traits>; + +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +class Cech_blocker { + public: + bool operator()(Simplex_handle sh) { + std::vector<Point> points; +#if DEBUG_TRACES + std::cout << "Cech_blocker on ["; +#endif // DEBUG_TRACES + for (auto vertex : simplex_tree_.simplex_vertex_range(sh)) { + points.push_back(point_cloud_[vertex]); +#if DEBUG_TRACES + std::cout << vertex << ", "; +#endif // DEBUG_TRACES + } + Min_sphere ms(points.begin(),points.end()); + Filtration_value radius = ms.radius(); +#if DEBUG_TRACES + std::cout << "] - radius = " << radius << " - returns " << (radius > threshold_) << std::endl; +#endif // DEBUG_TRACES + simplex_tree_.assign_filtration(sh, radius); + return (radius > threshold_); + } + Cech_blocker(Simplex_tree& simplex_tree, Filtration_value threshold, const std::vector<Point>& point_cloud) + : simplex_tree_(simplex_tree), + threshold_(threshold), + point_cloud_(point_cloud) { } + private: + Simplex_tree simplex_tree_; + Filtration_value threshold_; + std::vector<Point> point_cloud_; +}; + +template< typename InputPointRange> +Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold); + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , Filtration_value & threshold + , int & dim_max); + +int main(int argc, char * argv[]) { + std::string off_file_points; + Filtration_value threshold; + int dim_max; + + program_options(argc, argv, off_file_points, threshold, dim_max); + + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + // Compute the proximity graph of the points + Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold); + + //Min_sphere sph1(off_reader.get_point_cloud()[0], off_reader.get_point_cloud()[1], off_reader.get_point_cloud()[2]); + // Construct the Rips complex in a Simplex Tree + Simplex_tree st; + // insert the proximity graph in the simplex tree + st.insert_graph(prox_graph); + // expand the graph until dimension dim_max + st.expansion_with_blockers(dim_max, Cech_blocker(st, threshold, off_reader.get_point_cloud())); + + 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(); + +#if DEBUG_TRACES + std::cout << "********************************************************************\n"; + // Display the Simplex_tree - Can not be done in the middle of 2 inserts + std::cout << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n"; + std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << " " << "[" << st.filtration(f_simplex) << "] "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << static_cast<int>(vertex) << " "; + } + std::cout << std::endl; + } +#endif // DEBUG_TRACES + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , Filtration_value & threshold + , int & dim_max) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a 3d point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Cech complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Cech complex we want to compute."); + + po::positional_options_description pos; + pos.add("input-file", 1); + + 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 << "Construct a Cech complex defined on a set of input points.\n \n"; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} + +/** Output the proximity graph of the points. + * + * If points contains n elements, the proximity graph is the graph + * with n vertices, and an edge [u,v] iff the distance function between + * points u and v is smaller than threshold. + * + * The type PointCloud furnishes .begin() and .end() methods, that return + * iterators with value_type Point. + */ +template< typename InputPointRange> +Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold) { + std::vector< Edge_t > edges; + std::vector< Filtration_value > edges_fil; + + Kernel k; + Vertex_handle idx_u, idx_v; + Filtration_value fil; + idx_u = 0; + for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { + idx_v = idx_u + 1; + for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) { + fil = k.squared_distance_d_object()(*it_u, *it_v); + // For Cech Complex, threshold is a radius (distance /2) + fil = std::sqrt(fil) / 2.; + if (fil <= threshold) { + edges.emplace_back(idx_u, idx_v); + edges_fil.push_back(fil); + } + } + ++idx_u; + } + + Graph_t skel_graph(edges.begin() + , edges.end() + , edges_fil.begin() + , idx_u); // number of points labeled from 0 to idx_u-1 + + auto vertex_prop = boost::get(Gudhi::vertex_filtration_t(), skel_graph); + + boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; + for (std::tie(vi, vi_end) = boost::vertices(skel_graph); + vi != vi_end; ++vi) { + boost::put(vertex_prop, *vi, 0.); + } + + return skel_graph; +} diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 7da767cb..cb6ab309 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -391,13 +391,13 @@ class Simplex_tree { return sh->second.key(); } - /** \brief Returns the simplex associated to a key. + /** \brief Returns the simplex that has index idx in the filtration. * * The filtration must be initialized. * \pre SimplexTreeOptions::store_key */ - Simplex_handle simplex(Simplex_key key) const { - return filtration_vect_[key]; + Simplex_handle simplex(Simplex_key idx) const { + return filtration_vect_[idx]; } /** \brief Returns the filtration value of a simplex. @@ -761,8 +761,12 @@ class Simplex_tree { return &root_; } - /** Set a dimension for the simplicial complex. */ + /** \brief Set a dimension for the simplicial complex. + * \details This function must be used with caution because it disables dimension recomputation when required + * (this recomputation can be triggered by `remove_maximal_simplex()` or `prune_above_filtration()`). + */ void set_dimension(int dimension) { + dimension_to_be_lowered_ = false; dimension_ = dimension; } @@ -1145,7 +1149,7 @@ class Simplex_tree { Siblings * new_sib = new Siblings(siblings, // oncles simplex->first, // parent boost::adaptors::reverse(intersection)); // boost::container::ordered_unique_range_t - std::vector<Simplex_handle> blocked_new_sib_list; + std::vector<Vertex_handle> blocked_new_sib_vertex_list; // As all intersections are inserted, we can call the blocker function on all new_sib members for (auto new_sib_member = new_sib->members().begin(); new_sib_member != new_sib->members().end(); @@ -1153,17 +1157,19 @@ class Simplex_tree { bool blocker_result = block_simplex(new_sib_member); // new_sib member has been blocked by the blocker function // add it to the list to be removed - do not perform it while looping on it - if (blocker_result) - blocked_new_sib_list.push_back(new_sib_member); - } - bool removed = false; - for (auto& blocked_new_sib_member : blocked_new_sib_list){ - removed = removed || remove_maximal_simplex(blocked_new_sib_member); + if (blocker_result) { + blocked_new_sib_vertex_list.push_back(new_sib_member->first); + } } - if (removed) { + if (blocked_new_sib_vertex_list.size() == new_sib->members().size()) { + // Specific case where all have to be deleted + delete new_sib; // ensure the children property simplex->second.assign_children(siblings); } else { + for (auto& blocked_new_sib_member : blocked_new_sib_vertex_list) { + new_sib->members().erase(blocked_new_sib_member); + } // ensure recursive call simplex->second.assign_children(new_sib); siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex); @@ -1338,16 +1344,14 @@ class Simplex_tree { public: /** \brief Remove a maximal simplex. * @param[in] sh Simplex handle on the maximal simplex to remove. - * @return a boolean value that is an implementation detail, and that the user is supposed to ignore * \pre Please check the simplex has no coface before removing it. * \exception std::invalid_argument In debug mode, if sh has children. - * \post Be aware that removing is shifting data in a flat_map (`initialize_filtration()` to be done). + * \post Be aware that removing is shifting data in a flat_map (initialize_filtration to be done). * \post Note that the dimension of the simplicial complex may be lower after calling `remove_maximal_simplex()` * than it was before. However, `upper_bound_dimension()` will return the old value, which remains a valid upper * bound. If you care, you can call `dimension()` to recompute the exact dimension. - * \internal @return true if the leaf's branch has no other leaves (branch's children has been re-assigned), false otherwise. */ - bool remove_maximal_simplex(Simplex_handle sh) { + void remove_maximal_simplex(Simplex_handle sh) { // Guarantee the simplex has no children GUDHI_CHECK(!has_children(sh), std::invalid_argument("Simplex_tree::remove_maximal_simplex - argument has children")); @@ -1365,9 +1369,7 @@ class Simplex_tree { delete child; // dimension may need to be lowered dimension_to_be_lowered_ = true; - return true; } - return false; } private: diff --git a/src/common/include/gudhi/graph_simplicial_complex.h b/src/common/include/gudhi/graph_simplicial_complex.h index 5fe7c826..d84421b2 100644 --- a/src/common/include/gudhi/graph_simplicial_complex.h +++ b/src/common/include/gudhi/graph_simplicial_complex.h @@ -28,6 +28,9 @@ #include <utility> // for pair<> #include <vector> #include <map> +#include <tuple> // for std::tie + +namespace Gudhi { /* Edge tag for Boost PropertyGraph. */ struct edge_filtration_t { @@ -39,4 +42,64 @@ struct vertex_filtration_t { typedef boost::vertex_property_tag kind; }; +template <typename SimplicialComplexForProximityGraph> +using Proximity_graph = typename boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS +, boost::property < vertex_filtration_t, typename SimplicialComplexForProximityGraph::Filtration_value > +, boost::property < edge_filtration_t, typename SimplicialComplexForProximityGraph::Filtration_value >>; + +/** \brief Computes the proximity graph of the points. + * + * If points contains n elements, the proximity graph is the graph with n vertices, and an edge [u,v] iff the + * distance function between points u and v is smaller than threshold. + * + * \tparam ForwardPointRange furnishes `.begin()` and `.end()` methods. + * + * \tparam Distance furnishes `operator()(const Point& p1, const Point& p2)`, where + * `Point` is a point from the `ForwardPointRange`, and that returns a `Filtration_value`. + */ +template< typename SimplicialComplexForProximityGraph + , typename ForwardPointRange + , typename Distance > +Proximity_graph<SimplicialComplexForProximityGraph> compute_proximity_graph( + const ForwardPointRange& points, + typename SimplicialComplexForProximityGraph::Filtration_value threshold, + Distance distance) { + using Vertex_handle = typename SimplicialComplexForProximityGraph::Vertex_handle; + using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; + + std::vector<std::pair< Vertex_handle, Vertex_handle >> edges; + std::vector< Filtration_value > edges_fil; + std::map< Vertex_handle, Filtration_value > vertices; + + Vertex_handle idx_u, idx_v; + Filtration_value fil; + idx_u = 0; + for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { + idx_v = idx_u + 1; + for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) { + fil = distance(*it_u, *it_v); + if (fil <= threshold) { + edges.emplace_back(idx_u, idx_v); + edges_fil.push_back(fil); + } + } + ++idx_u; + } + + // Points are labeled from 0 to idx_u-1 + Proximity_graph<SimplicialComplexForProximityGraph> skel_graph(edges.begin(), edges.end(), edges_fil.begin(), idx_u); + + auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph); + + typename boost::graph_traits<Proximity_graph<SimplicialComplexForProximityGraph>>::vertex_iterator vi, vi_end; + for (std::tie(vi, vi_end) = boost::vertices(skel_graph); + vi != vi_end; ++vi) { + boost::put(vertex_prop, *vi, 0.); + } + + return skel_graph; +} + +} // namespace Gudhi + #endif // GRAPH_SIMPLICIAL_COMPLEX_H_ diff --git a/src/cython/cython/periodic_cubical_complex.pyx b/src/cython/cython/periodic_cubical_complex.pyx index 581c7b69..3025f125 100644 --- a/src/cython/cython/periodic_cubical_complex.pyx +++ b/src/cython/cython/periodic_cubical_complex.pyx @@ -33,7 +33,7 @@ __license__ = "GPL v3" cdef extern from "Cubical_complex_interface.h" namespace "Gudhi": cdef cppclass Periodic_cubical_complex_base_interface "Gudhi::Cubical_complex::Cubical_complex_interface<Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double>>": - Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells) + Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions) Periodic_cubical_complex_base_interface(string perseus_file) int num_simplices() int dimension() @@ -58,7 +58,7 @@ cdef class PeriodicCubicalComplex: # Fake constructor that does nothing but documenting the constructor def __init__(self, dimensions=None, top_dimensional_cells=None, - perseus_file=''): + periodic_dimensions=None, perseus_file=''): """PeriodicCubicalComplex constructor from dimensions and top_dimensional_cells or from a Perseus-style file name. @@ -66,6 +66,8 @@ cdef class PeriodicCubicalComplex: :type dimensions: list of int :param top_dimensional_cells: A list of cells filtration values. :type top_dimensional_cells: list of double + :param periodic_dimensions: A list of top dimensional cells periodicity value. + :type periodic_dimensions: list of boolean Or @@ -75,10 +77,10 @@ cdef class PeriodicCubicalComplex: # The real cython constructor def __cinit__(self, dimensions=None, top_dimensional_cells=None, - perseus_file=''): - if (dimensions is not None) and (top_dimensional_cells is not None) and (perseus_file is ''): - self.thisptr = new Periodic_cubical_complex_base_interface(dimensions, top_dimensional_cells) - elif (dimensions is None) and (top_dimensional_cells is None) and (perseus_file is not ''): + periodic_dimensions=None, perseus_file=''): + if (dimensions is not None) and (top_dimensional_cells is not None) and (periodic_dimensions is not None) and (perseus_file is ''): + self.thisptr = new Periodic_cubical_complex_base_interface(dimensions, top_dimensional_cells, periodic_dimensions) + elif (dimensions is None) and (top_dimensional_cells is None) and (periodic_dimensions is None) and (perseus_file is not ''): if os.path.isfile(perseus_file): self.thisptr = new Periodic_cubical_complex_base_interface(str.encode(perseus_file)) else: diff --git a/src/cython/cython/simplex_tree.pyx b/src/cython/cython/simplex_tree.pyx index 32b91028..8a436619 100644 --- a/src/cython/cython/simplex_tree.pyx +++ b/src/cython/cython/simplex_tree.pyx @@ -37,11 +37,13 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>": Simplex_tree() double simplex_filtration(vector[int] simplex) + void assign_simplex_filtration(vector[int] simplex, double filtration) void initialize_filtration() int num_vertices() int num_simplices() void set_dimension(int dimension) int dimension() + int upper_bound_dimension() bint find_simplex(vector[int] simplex) bint insert_simplex_and_subfaces(vector[int] simplex, double filtration) @@ -50,8 +52,9 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": vector[pair[vector[int], double]] get_star(vector[int] simplex) vector[pair[vector[int], double]] get_cofaces(vector[int] simplex, int dimension) - void remove_maximal_simplex(vector[int] simplex) void expansion(int max_dim) + void remove_maximal_simplex(vector[int] simplex) + bool prune_above_filtration(double filtration) cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>": @@ -113,15 +116,31 @@ cdef class SimplexTree: """ return self.thisptr.simplex_filtration(simplex) + def assign_filtration(self, simplex, filtration): + """This function assigns the simplicial complex filtration value for a + given N-simplex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int. + :param filtration: The simplicial complex filtration value. + :type filtration: float + """ + self.thisptr.assign_simplex_filtration(simplex, filtration) + def initialize_filtration(self): """This function initializes and sorts the simplicial complex filtration vector. .. note:: - This function must be launched before persistence, betti_numbers, - persistent_betti_numbers or get_filtration after inserting or - removing simplices. + This function must be launched before + :func:`persistence()<gudhi.SimplexTree.persistence>`, + :func:`betti_numbers()<gudhi.SimplexTree.betti_numbers>`, + :func:`persistent_betti_numbers()<gudhi.SimplexTree.persistent_betti_numbers>`, + or :func:`get_filtration()<gudhi.SimplexTree.get_filtration>` + after :func:`inserting<gudhi.SimplexTree.insert>` or + :func:`removing<gudhi.SimplexTree.remove_maximal_simplex>` + simplices. """ self.thisptr.initialize_filtration() @@ -148,21 +167,42 @@ cdef class SimplexTree: :returns: the simplicial complex dimension. :rtype: int + + .. note:: + + This function is not constant time because it can recompute + dimension if required (can be triggered by + :func:`remove_maximal_simplex()<gudhi.SimplexTree.remove_maximal_simplex>` + or + :func:`prune_above_filtration()<gudhi.SimplexTree.prune_above_filtration>` + methods). """ return self.thisptr.dimension() - def set_dimension(self, dimension): - """This function sets the dimension of the simplicial complex. + def upper_bound_dimension(self): + """This function returns a valid dimension upper bound of the + simplicial complex. - insert and remove_maximal_simplex functions do not update dimension - value of the `SimplexTree`. + :returns: an upper bound on the dimension of the simplicial complex. + :rtype: int + """ + return self.thisptr.upper_bound_dimension() - `AlphaComplex`, `RipsComplex`, `TangentialComplex` and `WitnessComplex` - automatically sets the correct dimension in their `create_simplex_tree` - functions. + def set_dimension(self, dimension): + """This function sets the dimension of the simplicial complex. :param dimension: The new dimension value. :type dimension: int. + + .. note:: + + This function must be used with caution because it disables + dimension recomputation when required + (this recomputation can be triggered by + :func:`remove_maximal_simplex()<gudhi.SimplexTree.remove_maximal_simplex>` + or + :func:`prune_above_filtration()<gudhi.SimplexTree.prune_above_filtration>` + ). """ self.thisptr.set_dimension(<int>dimension) @@ -286,9 +326,57 @@ cdef class SimplexTree: :param simplex: The N-simplex, represented by a list of vertex. :type simplex: list of int. + + .. note:: + + Be aware that removing is shifting data in a flat_map + (:func:`initialize_filtration()<gudhi.SimplexTree.initialize_filtration>` to be done). + + .. note:: + + The dimension of the simplicial complex may be lower after calling + remove_maximal_simplex than it was before. However, + :func:`upper_bound_dimension()<gudhi.SimplexTree.upper_bound_dimension>` + method will return the old value, which + remains a valid upper bound. If you care, you can call + :func:`dimension()<gudhi.SimplexTree.dimension>` + to recompute the exact dimension. """ self.thisptr.remove_maximal_simplex(simplex) + def prune_above_filtration(self, filtration): + """Prune above filtration value given as parameter. + + :param filtration: Maximum threshold value. + :type filtration: float. + :returns: The filtration modification information. + :rtype: bint + + + .. note:: + + Some simplex tree functions require the filtration to be valid. + prune_above_filtration function is not launching + :func:`initialize_filtration()<gudhi.SimplexTree.initialize_filtration>` + but returns the filtration modification + information. If the complex has changed , please call + :func:`initialize_filtration()<gudhi.SimplexTree.initialize_filtration>` + to recompute it. + + .. note:: + + Note that the dimension of the simplicial complex may be lower + after calling + :func:`prune_above_filtration()<gudhi.SimplexTree.prune_above_filtration>` + than it was before. However, + :func:`upper_bound_dimension()<gudhi.SimplexTree.upper_bound_dimension>` + will return the old value, which remains a + valid upper bound. If you care, you can call + :func:`dimension()<gudhi.SimplexTree.dimension>` + method to recompute the exact dimension. + """ + return self.thisptr.prune_above_filtration(filtration) + def expansion(self, max_dim): """Expands the Simplex_tree containing only its one skeleton until dimension max_dim. @@ -336,8 +424,9 @@ cdef class SimplexTree: :returns: The Betti numbers ([B0, B1, ..., Bn]). :rtype: list of int - :note: betti_numbers function requires persistence function to be - launched first. + :note: betti_numbers function requires + :func:`persistence()<gudhi.SimplexTree.persistence>` + function to be launched first. """ cdef vector[int] bn_result if self.pcohptr != NULL: @@ -361,7 +450,8 @@ cdef class SimplexTree: :returns: The persistent Betti numbers ([B0, B1, ..., Bn]). :rtype: list of int - :note: persistent_betti_numbers function requires persistence + :note: persistent_betti_numbers function requires + :func:`persistence()<gudhi.SimplexTree.persistence>` function to be launched first. """ cdef vector[int] pbn_result @@ -381,8 +471,9 @@ cdef class SimplexTree: :returns: The persistence intervals. :rtype: list of pair of float - :note: intervals_in_dim function requires persistence function to be - launched first. + :note: intervals_in_dim function requires + :func:`persistence()<gudhi.SimplexTree.persistence>` + function to be launched first. """ cdef vector[pair[double,double]] intervals_result if self.pcohptr != NULL: @@ -399,8 +490,9 @@ cdef class SimplexTree: :param persistence_file: The specific dimension. :type persistence_file: string. - :note: intervals_in_dim function requires persistence function to be - launched first. + :note: intervals_in_dim function requires + :func:`persistence()<gudhi.SimplexTree.persistence>` + function to be launched first. """ if self.pcohptr != NULL: if persistence_file != '': diff --git a/src/cython/include/Cubical_complex_interface.h b/src/cython/include/Cubical_complex_interface.h index 7c0148f1..fad92c2c 100644 --- a/src/cython/include/Cubical_complex_interface.h +++ b/src/cython/include/Cubical_complex_interface.h @@ -43,6 +43,12 @@ class Cubical_complex_interface : public Bitmap_cubical_complex<CubicalComplexOp : Bitmap_cubical_complex<CubicalComplexOptions>(dimensions, top_dimensional_cells) { } + Cubical_complex_interface(const std::vector<unsigned>& dimensions, + const std::vector<double>& top_dimensional_cells, + const std::vector<bool>& periodic_dimensions) + : Bitmap_cubical_complex<CubicalComplexOptions>(dimensions, top_dimensional_cells, periodic_dimensions) { + } + Cubical_complex_interface(const std::string& perseus_file) : Bitmap_cubical_complex<CubicalComplexOptions>(perseus_file.c_str()) { } diff --git a/src/cython/include/Simplex_tree_interface.h b/src/cython/include/Simplex_tree_interface.h index 09e7e992..54a4f824 100644 --- a/src/cython/include/Simplex_tree_interface.h +++ b/src/cython/include/Simplex_tree_interface.h @@ -52,6 +52,10 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return (Base::find(vh) != Base::null_simplex()); } + void assign_simplex_filtration(const Simplex& vh, Filtration_value filtration) { + Base::assign_filtration(Base::find(vh), filtration); + } + bool insert(const Simplex& simplex, Filtration_value filtration = 0) { Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration); return (result.second); diff --git a/src/cython/test/test_simplex_tree.py b/src/cython/test/test_simplex_tree.py index 801d52b7..6dec5d94 100755 --- a/src/cython/test/test_simplex_tree.py +++ b/src/cython/test/test_simplex_tree.py @@ -134,3 +134,30 @@ def test_expansion(): ([1, 2], 0.5), ([0, 1, 2], 0.5), ([1, 2, 3], 0.5), ([5], 0.6), ([6], 0.6), ([5, 6], 0.6), ([4], 0.7), ([2, 4], 0.7), ([0, 3], 0.8), ([0, 1, 3], 0.8), ([0, 2, 3], 0.8), ([0, 1, 2, 3], 0.8), ([4, 6], 0.9), ([3, 6], 1.0)] + +def test_automatic_dimension(): + st = SimplexTree() + assert st.__is_defined() == True + assert st.__is_persistence_defined() == False + + # insert test + assert st.insert([0,1,3], filtration=0.5) == True + assert st.insert([0,1,2], filtration=1.) == True + + assert st.num_vertices() == 4 + assert st.num_simplices() == 11 + + assert st.dimension() == 2 + assert st.upper_bound_dimension() == 2 + + assert st.prune_above_filtration(0.6) == True + assert st.dimension() == 2 + assert st.upper_bound_dimension() == 2 + + st.assign_filtration([0, 1, 3], 0.7) + assert st.filtration([0, 1, 3]) == 0.7 + + st.remove_maximal_simplex([0, 1, 3]) + assert st.upper_bound_dimension() == 2 + assert st.dimension() == 1 + assert st.upper_bound_dimension() == 1 |