diff options
Diffstat (limited to 'src/Bitmap_cubical_complex/include')
4 files changed, 1867 insertions, 0 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 new file mode 100644 index 00000000..67e1fed3 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -0,0 +1,596 @@ +/* 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): Pawel Dlotko + * + * Copyright (C) 2015 INRIA Sophia-Saclay (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 BITMAP_CUBICAL_COMPLEX_H_ +#define BITMAP_CUBICAL_COMPLEX_H_ + +#include <gudhi/Bitmap_cubical_complex_base.h> +#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h> + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_sort.h> +#endif + +#include <limits> +#include <utility> // for pair<> +#include <algorithm> // for sort +#include <vector> +#include <numeric> // for iota + +namespace Gudhi { + +namespace Cubical_complex { + +// global variable, was used just for debugging. +const bool globalDbg = false; + +template <typename T> class is_before_in_filtration; + +/** + * This is a Bitmap_cubical_complex class. It joints a functionalities of Bitmap_cubical_complex_base and + * Bitmap_cubical_complex_periodic_boundary_conditions_base classes into + * Gudhi persistent homology engine. It is a template class that inherit from its template parameter. The template + * parameter is supposed to be either Bitmap_cubical_complex_base or + * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + **/ + +/** + *@class Bitmap_cubical_complex + *@brief Cubical complex represented as a bitmap. + *@ingroup cubical_complex + */ +template <typename T> +class Bitmap_cubical_complex : public T { + public: + //*********************************************// + // Typedefs and typenames + //*********************************************// + typedef size_t Simplex_key; + typedef typename T::filtration_type Filtration_value; + typedef Simplex_key Simplex_handle; + + + //*********************************************// + // Constructors + //*********************************************// + // Over here we need to define various input types. I am proposing the following ones: + // Perseus style + // TODO(PD) H5 files? + // TODO(PD) binary files with little endiangs / big endians ? + // TODO(PD) constructor from a vector of elements of a type T. ? + + /** + * Constructor form a Perseus-style file. + **/ + Bitmap_cubical_complex(const char* perseus_style_file) : + T(perseus_style_file), key_associated_to_simplex(this->total_number_of_cells + 1) { + 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) { + this->key_associated_to_simplex[i] = i; + } + // we initialize this only once, in each constructor, when the bitmap is constructed. + // If the user decide to change some elements of the bitmap, then this procedure need + // to be called again. + this->initialize_simplex_associated_to_key(); + } + + /** + * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells + * in the following directions and vector of element of a type T + * with filtration on top dimensional cells. + **/ + Bitmap_cubical_complex(const std::vector<unsigned>& dimensions, + const std::vector<typename T::filtration_type>& 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) { + this->key_associated_to_simplex[i] = i; + } + // we initialize this only once, in each constructor, when the bitmap is constructed. + // If the user decide to change some elements of the bitmap, then this procedure need + // to be called again. + this->initialize_simplex_associated_to_key(); + } + + /** + * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells + * in the following directions and vector of element of a type T::filtration_type + * with filtration on top dimensional cells. The last parameter of the constructor is a vector of boolean of a length + * equal to the dimension of cubical complex. + * If the position i on this vector is true, then we impose periodic boundary conditions in this direction. + **/ + Bitmap_cubical_complex(const std::vector<unsigned>& dimensions, + const std::vector<typename T::filtration_type>& top_dimensional_cells, + 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) { + this->key_associated_to_simplex[i] = i; + } + // we initialize this only once, in each constructor, when the bitmap is constructed. + // If the user decide to change some elements of the bitmap, then this procedure need + // to be called again. + this->initialize_simplex_associated_to_key(); + } + + /** + * Destructor of the Bitmap_cubical_complex class. + **/ + virtual ~Bitmap_cubical_complex() {} + + //*********************************************// + // Other 'easy' functions + //*********************************************// + + /** + * Returns number of all cubes in the complex. + **/ + 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. + **/ + static Simplex_handle null_simplex() { + if (globalDbg) { + std::cerr << "Simplex_handle null_simplex()\n"; + } + return std::numeric_limits<Simplex_handle>::max(); + } + + /** + * Returns dimension of the complex. + **/ + inline size_t dimension()const { + return this->sizes.size(); + } + + /** + * Return dimension of a cell pointed by the Simplex_handle. + **/ + inline unsigned dimension(Simplex_handle sh)const { + if (globalDbg) { + std::cerr << "unsigned dimension(const Simplex_handle& sh)\n"; + } + if (sh != std::numeric_limits<Simplex_handle>::max()) return this->get_dimension_of_a_cell(sh); + return -1; + } + + /** + * Return the filtration of a cell pointed by the Simplex_handle. + **/ + typename T::filtration_type filtration(Simplex_handle sh) { + if (globalDbg) { + std::cerr << "T::filtration_type filtration(const Simplex_handle& sh)\n"; + } + // Returns the filtration value of a simplex. + if (sh != std::numeric_limits<Simplex_handle>::max()) return this->data[sh]; + return std::numeric_limits<Simplex_handle>::max(); + } + + /** + * Return a key which is not a key of any cube in the considered data structure. + **/ + static Simplex_key null_key() { + if (globalDbg) { + std::cerr << "Simplex_key null_key()\n"; + } + return std::numeric_limits<Simplex_handle>::max(); + } + + /** + * Return the key of a cube pointed by the Simplex_handle. + **/ + Simplex_key key(Simplex_handle sh)const { + if (globalDbg) { + std::cerr << "Simplex_key key(const Simplex_handle& sh)\n"; + } + if (sh != std::numeric_limits<Simplex_handle>::max()) { + return this->key_associated_to_simplex[sh]; + } + return this->null_key(); + } + + /** + * Return the Simplex_handle given the key of the cube. + **/ + Simplex_handle simplex(Simplex_key key) { + if (globalDbg) { + std::cerr << "Simplex_handle simplex(Simplex_key key)\n"; + } + if (key != std::numeric_limits<Simplex_handle>::max()) { + return this->simplex_associated_to_key[ key ]; + } + return null_simplex(); + } + + /** + * Assign key to a cube pointed by the Simplex_handle + **/ + void assign_key(Simplex_handle sh, Simplex_key key) { + if (globalDbg) { + std::cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n"; + } + if (key == std::numeric_limits<Simplex_handle>::max()) return; + this->key_associated_to_simplex[sh] = key; + this->simplex_associated_to_key[key] = sh; + } + + /** + * Function called from a constructor. It is needed for Filtration_simplex_iterator to work. + **/ + void initialize_simplex_associated_to_key(); + + //*********************************************// + // Iterators + //*********************************************// + + /** + * Boundary_simplex_range class provides ranges for boundary iterators. + **/ + typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator; + typedef typename std::vector< Simplex_handle > Boundary_simplex_range; + + /** + * Filtration_simplex_iterator class provides an iterator though the whole structure in the order of filtration. + * Secondary criteria for filtration are: + * (1) Dimension of a cube (lower dimensional comes first). + * (2) Position in the data structure (the ones that are earlies in the data structure comes first). + **/ + class Filtration_simplex_range; + + class Filtration_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > { + // Iterator over all simplices of the complex in the order of the indexing scheme. + // 'value_type' must be 'Simplex_handle'. + public: + Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) { } + + Filtration_simplex_iterator() : b(NULL), position(0) { } + + Filtration_simplex_iterator operator++() { + if (globalDbg) { + std::cerr << "Filtration_simplex_iterator operator++\n"; + } + ++this->position; + return (*this); + } + + Filtration_simplex_iterator operator++(int) { + Filtration_simplex_iterator result = *this; + ++(*this); + return result; + } + + Filtration_simplex_iterator& operator=(const Filtration_simplex_iterator& rhs) { + if (globalDbg) { + std::cerr << "Filtration_simplex_iterator operator =\n"; + } + this->b = rhs.b; + this->position = rhs.position; + return (*this); + } + + bool operator==(const Filtration_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n"; + } + return ( this->position == rhs.position); + } + + bool operator!=(const Filtration_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator != ( const Filtration_simplex_iterator& rhs )\n"; + } + return !(*this == rhs); + } + + Simplex_handle operator*() { + if (globalDbg) { + std::cerr << "Simplex_handle operator*()\n"; + } + return this->b->simplex_associated_to_key[ this->position ]; + } + + friend class Filtration_simplex_range; + + private: + Bitmap_cubical_complex<T>* b; + size_t position; + }; + + /** + * Filtration_simplex_range provides the ranges for Filtration_simplex_iterator. + **/ + class Filtration_simplex_range { + // Range over the simplices of the complex in the order of the filtration. + // .begin() and .end() return type Filtration_simplex_iterator. + public: + typedef Filtration_simplex_iterator const_iterator; + typedef Filtration_simplex_iterator iterator; + + Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) { } + + Filtration_simplex_iterator begin() { + if (globalDbg) { + std::cerr << "Filtration_simplex_iterator begin() \n"; + } + return Filtration_simplex_iterator(this->b); + } + + Filtration_simplex_iterator end() { + if (globalDbg) { + std::cerr << "Filtration_simplex_iterator end()\n"; + } + Filtration_simplex_iterator it(this->b); + it.position = this->b->simplex_associated_to_key.size(); + return it; + } + + private: + Bitmap_cubical_complex<T>* b; + }; + + + + //*********************************************// + // Methods to access iterators from the container: + + /** + * boundary_simplex_range creates an object of a Boundary_simplex_range class + * that provides ranges for the Boundary_simplex_iterator. + **/ + Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { + return this->get_boundary_of_a_cell(sh); + } + + /** + * filtration_simplex_range creates an object of a Filtration_simplex_range class + * that provides ranges for the Filtration_simplex_iterator. + **/ + Filtration_simplex_range filtration_simplex_range() { + if (globalDbg) { + std::cerr << "Filtration_simplex_range filtration_simplex_range()\n"; + } + // Returns a range over the simplices of the complex in the order of the filtration + return Filtration_simplex_range(this); + } + //*********************************************// + + + + //*********************************************// + // Elements which are in Gudhi now, but I (and in all the cases I asked also Marc) do not understand why they are + // there. + // TODO(PD) the file IndexingTag.h in the Gudhi library contains an empty structure, so + // I understand that this is something that was planned (for simplicial maps?) + // but was never finished. The only idea I have here is to use the same empty structure from + // IndexingTag.h file, but only if the compiler needs it. If the compiler + // do not need it, then I would rather not add here elements which I do not understand. + // typedef Indexing_tag + + /** + * 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); + if (globalDbg) { + std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n"; + std::cerr << "bdry.size() : " << bdry.size() << std::endl; + } + // this method returns two first elements from the boundary of sh. + if (bdry.size() < 2) + throw("Error in endpoints in Bitmap_cubical_complex class. The cell have less than two elements in the " + "boundary."); + return std::make_pair(bdry[0], bdry[1]); + } + + + /** + * Class needed for compatibility with Gudhi. Not useful for other purposes. + **/ + class Skeleton_simplex_range; + + class Skeleton_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > { + // 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) { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , size_t d )\n"; + } + // find the position of the first simplex of a dimension d + this->position = 0; + while ( + (this->position != b->data.size()) && + (this->b->get_dimension_of_a_cell(this->position) != this->dimension) + ) { + ++this->position; + } + } + + Skeleton_simplex_iterator() : b(NULL), position(0), dimension(0) { } + + Skeleton_simplex_iterator operator++() { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator operator++()\n"; + } + // increment the position as long as you did not get to the next element of the dimension dimension. + ++this->position; + while ( + (this->position != this->b->data.size()) && + (this->b->get_dimension_of_a_cell(this->position) != this->dimension) + ) { + ++this->position; + } + return (*this); + } + + Skeleton_simplex_iterator operator++(int) { + Skeleton_simplex_iterator result = *this; + ++(*this); + return result; + } + + Skeleton_simplex_iterator& operator=(const Skeleton_simplex_iterator& rhs) { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator operator =\n"; + } + this->b = rhs.b; + this->position = rhs.position; + this->dimension = rhs.dimension; + return (*this); + } + + bool operator==(const Skeleton_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator ==\n"; + } + return ( this->position == rhs.position); + } + + bool operator!=(const Skeleton_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n"; + } + return !(*this == rhs); + } + + Simplex_handle operator*() { + if (globalDbg) { + std::cerr << "Simplex_handle operator*() \n"; + } + return this->position; + } + + friend class Skeleton_simplex_range; + private: + Bitmap_cubical_complex<T>* b; + size_t position; + unsigned dimension; + }; + + /** + * Class needed for compatibility with Gudhi. Not useful for other purposes. + **/ + class Skeleton_simplex_range { + // Range over the simplices of the complex in the order of the filtration. + // .begin() and .end() return type Filtration_simplex_iterator. + public: + typedef Skeleton_simplex_iterator const_iterator; + typedef Skeleton_simplex_iterator iterator; + + Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, unsigned dimension) : b(b), dimension(dimension) { } + + Skeleton_simplex_iterator begin() { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator begin()\n"; + } + return Skeleton_simplex_iterator(this->b, this->dimension); + } + + Skeleton_simplex_iterator end() { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator end()\n"; + } + Skeleton_simplex_iterator it(this->b, this->dimension); + it.position = this->b->data.size(); + return it; + } + + private: + Bitmap_cubical_complex<T>* b; + unsigned dimension; + }; + + /** + * Function needed for compatibility with Gudhi. Not useful for other purposes. + **/ + Skeleton_simplex_range skeleton_simplex_range(unsigned dimension) { + if (globalDbg) { + std::cerr << "Skeleton_simplex_range skeleton_simplex_range( unsigned dimension )\n"; + } + return Skeleton_simplex_range(this, dimension); + } + + friend class is_before_in_filtration<T>; + + protected: + std::vector< size_t > key_associated_to_simplex; + std::vector< size_t > simplex_associated_to_key; +}; // Bitmap_cubical_complex + +template <typename T> +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()); + 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, is_before_in_filtration<T>(this)); +#else + std::sort(simplex_associated_to_key.begin(), simplex_associated_to_key.end(), is_before_in_filtration<T>(this)); +#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 ) { + this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i; + } +} + +template <typename T> +class is_before_in_filtration { + public: + explicit is_before_in_filtration(Bitmap_cubical_complex<T> * CC) + : CC_(CC) { } + + bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1, + const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const { + // Not using st_->filtration(sh1) because it uselessly tests for null_simplex. + typename T::filtration_type fil1 = CC_->data[sh1]; + typename T::filtration_type fil2 = CC_->data[sh2]; + if (fil1 != fil2) { + 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); + if (dim1 != dim2) { + return dim1 < dim2; + } + // in this case both filtration and dimensions of the considered cubes are the same. To have stable sort, we simply + // compare their positions in the bitmap: + return sh1 < sh2; + } + + protected: + Bitmap_cubical_complex<T>* CC_; +}; + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_H_ 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 new file mode 100644 index 00000000..266ce051 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h @@ -0,0 +1,143 @@ +/* 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): Pawel Dlotko + * + * Copyright (C) 2015 INRIA Sophia-Saclay (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 BITMAP_CUBICAL_COMPLEX_COUNTER_H_ +#define BITMAP_CUBICAL_COMPLEX_COUNTER_H_ + +#include <iostream> +#include <vector> + +namespace Gudhi { + +namespace Cubical_complex { + +/** + * This is an implementation of a counter being a vector of integers. + * The constructor of the class takes as an input two vectors W and V. + * It assumes that W < V coordinatewise. + * If the initial counter W is not specified, it is assumed to be vector of zeros. + * The class allows to iterate between W and V by using increment() function. + * The increment() function returns a bool value. + * The current counter reach the end counter V if the value returned by the increment function is FALSE. + * This class is needed for the implementation of a bitmapCubicalComplex. + **/ + +class counter { + public: + /** + * Constructor of a counter class. It takes only the parameter which is the end value of the counter. + * The default beginning value is a vector of the same length as the endd, filled-in with zeros. + **/ + counter(const std::vector<unsigned>& endd) : begin(endd.size(), 0), end(endd), current(endd.size(), 0) { } + + /** + * Constructor of a counter class. It takes as the input beginn and end vector. + * It assumes that begin vector is lexicographically below the end vector. + **/ + counter(const std::vector< unsigned >& beginn, const std::vector< unsigned >& endd) : begin(beginn), end(endd), current(endd.size(), 0) { + if (beginn.size() != endd.size()) + throw "In constructor of a counter, begin and end vectors do not have the same size. Program terminate"; + } + + /** + * Function to increment the counter. If the value returned by the function is true, + * then the incrementation process was successful. + * If the value of the function is false, that means, that the counter have reached its end-value. + **/ + bool increment() { + 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) { + this->current[j] = this->begin[j]; + } + return true; + } + + /** + * Function to check if we are at the end of counter. + **/ + bool isFinal() { + for (size_t i = 0; i != this->current.size(); ++i) { + if (this->current[i] == this->end[i])return true; + } + return false; + } + + /** + * Function required in the implementation of bitmapCubicalComplexWPeriodicBoundaryCondition. + * Its aim is to find an counter corresponding to the element the following + * boundary element is identified with when periodic boundary conditions are imposed. + **/ + std::vector< unsigned > find_opposite(const std::vector< bool >& directionsForPeriodicBCond) { + std::vector< unsigned > result; + for (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 { + result.push_back(this->current[i]); + } + } + return result; + } + + /** + * Function checking at which positions the current value of a counter is the final value of the counter. + **/ + std::vector< bool > directions_of_finals() { + std::vector< bool > result; + for (size_t i = 0; i != this->current.size(); ++i) { + if (this->current[i] == this->end[i]) { + result.push_back(true); + } else { + result.push_back(false); + } + } + return result; + } + + /** + * Function to write counter to the stream. + **/ + 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) { + out << c.current[i] << " "; + } + return out; + } + + private: + std::vector< unsigned > begin; + std::vector< unsigned > end; + std::vector< unsigned > current; +}; + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_COUNTER_H_ 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 new file mode 100644 index 00000000..7294da98 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -0,0 +1,819 @@ +/* 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): Pawel Dlotko + * + * Copyright (C) 2015 INRIA Sophia-Saclay (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 BITMAP_CUBICAL_COMPLEX_BASE_H_ +#define BITMAP_CUBICAL_COMPLEX_BASE_H_ + +#include <gudhi/Bitmap_cubical_complex/counter.h> + +#include <iostream> +#include <vector> +#include <string> +#include <fstream> +#include <algorithm> +#include <iterator> +#include <limits> +#include <utility> // for pair<> + +namespace Gudhi { + +namespace Cubical_complex { + +/** + * @class Bitmap_cubical_complex_base + * @brief Cubical complex represented as a bitmap, class with basic implementation. + * @ingroup cubical_complex + */ + +/** + * This is a class implementing a basic bitmap data structure to store cubical complexes. + * It implements only the most basic subroutines. + * The idea of the bitmap is the following. Our aim is to have a memory efficient + * data structure to store d-dimensional cubical complex + * C being a cubical decomposition + * of a rectangular region of a space. This is achieved by storing C as a + * vector of bits (this is where the name 'bitmap' came from). + * Each cell is represented by a single + * bit (in case of black and white bitmaps, or by a single element of a type T + * (here T is a filtration type of a bitmap, typically a double). + * All the informations needed for homology and + * persistent homology computations (like dimension of a cell, boundary and + * coboundary elements of a cell, are then obtained from the + * position of the element in C. + * The default filtration used in this implementation is the lower star filtration. + */ +template <typename T> +class Bitmap_cubical_complex_base { + public: + typedef T filtration_type; + + /** + *Default constructor + **/ + Bitmap_cubical_complex_base() : + total_number_of_cells(0) { } + /** + * There are a few constructors of a Bitmap_cubical_complex_base class. + * First one, that takes vector<unsigned>, creates an empty bitmap of a dimension equal + * the number of elements in the + * input vector and size in the i-th dimension equal the number in the position i-of the input vector. + */ + Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes); + /** + * The second constructor takes as a input a Perseus style file. For more details, + * please consult the documentations of + * Perseus software as well as examples attached to this + * implementation. + **/ + Bitmap_cubical_complex_base(const char* perseus_style_file); + /** + * The last constructor of a Bitmap_cubical_complex_base class accepts vector of dimensions (as the first one) + * together with vector of filtration values of top dimensional cells. + **/ + Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells); + + /** + * Destructor of the Bitmap_cubical_complex_base class. + **/ + virtual ~Bitmap_cubical_complex_base() { } + + /** + * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell + * and get_cell_data are the basic + * functions that compute boundary / coboundary / dimension and the filtration + * value form a position of a cell in the structure of a bitmap. The input parameter of all of those function is a + * non-negative integer, indicating a position of a cube in the data structure. + * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to + * the positions of (co)boundary element of the input cell. + */ + virtual inline std::vector< size_t > get_boundary_of_a_cell(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 + * functions that compute boundary / coboundary / dimension and the filtration + * value form a position of a cell in the structure of a bitmap. + * The input parameter of all of those function is a non-negative integer, + * indicating a position of a cube in the data structure. + * In the case of functions that compute (co)boundary, the output is a vector if + * non-negative integers pointing to the + * positions of (co)boundary element of the input cell. + **/ + virtual inline std::vector< size_t > get_coboundary_of_a_cell(size_t cell)const; + /** + * In the case of get_dimension_of_a_cell function, the output is a non-negative integer + * indicating the dimension of a cell. + **/ + inline unsigned get_dimension_of_a_cell(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. + * This allows reading and changing the value of filtration. Note that if the value of a filtration is changed, the + * 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); + + + /** + * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells. + * Then, there are a few ways one can pick the filtration of lower dimensional + * cells. The most typical one is by so called lower star filtration. This function is always called by any + * constructor which takes the top dimensional cells. If you use such a constructor, + * then there is no need to call this function. Call it only if you are putting the filtration + * of the cells by your own (for instance by using Top_dimensional_cells_iterator). + **/ + void impose_lower_star_filtration(); // assume that top dimensional cells are already set. + + /** + * Returns dimension of a complex. + **/ + inline unsigned dimension()const { + return sizes.size(); + } + + /** + * Returns number of all cubes in the data structure. + **/ + inline unsigned size()const { + return this->data.size(); + } + + /** + * Writing to stream operator. By using it we get the values T of cells in order in which they are stored in the + * structure. This procedure is used for debugging purposes. + **/ + template <typename K> + friend std::ostream& operator<<(std::ostream & os, const Bitmap_cubical_complex_base<K>& b); + + /** + * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values + * 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. + * 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); + + /** + * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values + * 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( T diameter_of_bin ) is designed for that purpose. + * The parameter of it is the diameter of each bin. Note that the bottleneck distance between the persistence + * diagram of the cubical complex before and after using such a function will be bounded by the parameter + * diameter_of_bin. + **/ + void put_data_to_bins(T diameter_of_bin); + + /** + * Functions to find min and max values of filtration. + **/ + std::pair< T, T > min_max_filtration(); + + // ITERATORS + + /** + * Iterator through all cells in the complex (in order they appear in the structure -- i.e. + * in lexicographical order). + **/ + class All_cells_iterator : std::iterator< std::input_iterator_tag, T > { + public: + All_cells_iterator() { + this->counter = 0; + } + + All_cells_iterator operator++() { + // first find first element of the counter that can be increased: + ++this->counter; + return *this; + } + + All_cells_iterator operator++(int) { + All_cells_iterator result = *this; + ++(*this); + return result; + } + + All_cells_iterator& operator=(const All_cells_iterator& rhs) { + this->counter = rhs.counter; + return *this; + } + + bool operator==(const All_cells_iterator& rhs)const { + if (this->counter != rhs.counter)return false; + return true; + } + + bool operator!=(const All_cells_iterator& rhs)const { + return !(*this == rhs); + } + + /* + * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as + * an argument of the following functions: + * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell + * boundary and coboundary and dimension + * and in function get_cell_data to get a filtration of a cell. + */ + size_t operator*() { + return this->counter; + } + friend class Bitmap_cubical_complex_base; + protected: + size_t counter; + }; + + /** + * Function returning a All_cells_iterator to the first cell of the bitmap. + **/ + All_cells_iterator all_cells_iterator_begin() { + All_cells_iterator a; + return a; + } + + /** + * Function returning a All_cells_iterator to the last cell of the bitmap. + **/ + All_cells_iterator all_cells_iterator_end() { + All_cells_iterator a; + a.counter = this->data.size(); + return a; + } + + /** + * All_cells_range class provides ranges for All_cells_iterator + **/ + class All_cells_range { + public: + All_cells_range(Bitmap_cubical_complex_base* b) : b(b) { } + + All_cells_iterator begin() { + return b->all_cells_iterator_begin(); + } + + All_cells_iterator end() { + return b->all_cells_iterator_end(); + } + private: + Bitmap_cubical_complex_base<T>* b; + }; + + All_cells_range all_cells_range() { + return All_cells_range(this); + } + + + /** + * 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; + + /** + * 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); + } + + /** + * 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; + + /** + * 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); + } + + /** + * Iterator through top dimensional cells of the complex. The cells appear in order they are stored + * in the structure (i.e. in lexicographical order) + **/ + 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()); + // 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; + 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) { + this->counter[i] = 0; + } + } else { + ++this->counter[0]; + } + return *this; + } + + Top_dimensional_cells_iterator operator++(int) { + Top_dimensional_cells_iterator result = *this; + ++(*this); + return result; + } + + Top_dimensional_cells_iterator& operator=(const Top_dimensional_cells_iterator& rhs) { + this->counter = rhs.counter; + this->b = rhs.b; + return *this; + } + + 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) { + if (this->counter[i] != rhs.counter[i])return false; + } + return true; + } + + bool operator!=(const Top_dimensional_cells_iterator& rhs)const { + return !(*this == rhs); + } + + /* + * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as + * an argument of the following functions: + * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell + * 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(); + } + + size_t compute_index_in_bitmap()const { + size_t index = 0; + for (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) { + std::cout << this->counter[i] << " "; + } + } + friend class Bitmap_cubical_complex_base; + protected: + std::vector< size_t > counter; + Bitmap_cubical_complex_base& b; + }; + + /** + * Function returning a Top_dimensional_cells_iterator to the first top dimensional cell of the bitmap. + **/ + Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin() { + Top_dimensional_cells_iterator a(*this); + return a; + } + + /** + * Function returning a Top_dimensional_cells_iterator to the last top dimensional cell of the bitmap. + **/ + Top_dimensional_cells_iterator top_dimensional_cells_iterator_end() { + Top_dimensional_cells_iterator a(*this); + for (size_t i = 0; i != this->dimension(); ++i) { + a.counter[i] = this->sizes[i] - 1; + } + a.counter[0]++; + return a; + } + + /** + * Top_dimensional_cells_iterator_range class provides ranges for Top_dimensional_cells_iterator_range + **/ + class Top_dimensional_cells_range { + public: + Top_dimensional_cells_range(Bitmap_cubical_complex_base* b) : b(b) { } + + Top_dimensional_cells_iterator begin() { + return b->top_dimensional_cells_iterator_begin(); + } + + Top_dimensional_cells_iterator end() { + return b->top_dimensional_cells_iterator_end(); + } + private: + Bitmap_cubical_complex_base<T>* b; + }; + + Top_dimensional_cells_range top_dimensional_cells_range() { + return Top_dimensional_cells_range(this); + } + + + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + + inline size_t number_cells()const { + return this->total_number_of_cells; + } + + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + + protected: + std::vector<unsigned> sizes; + std::vector<unsigned> multipliers; + std::vector<T> data; + 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) { + this->sizes.push_back(sizes[i]); + this->multipliers.push_back(multiplier); + multiplier *= 2 * sizes[i] + 1; + } + this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max()); + 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) { + position += this->multipliers[i] * counter[i]; + } + return position; + } + + std::vector<unsigned> compute_counter_for_given_cell(size_t cell)const { + std::vector<unsigned> counter; + counter.reserve(this->sizes.size()); + for (size_t dim = this->sizes.size(); dim != 0; --dim) { + counter.push_back(cell / this->multipliers[dim - 1]); + cell = cell % this->multipliers[dim - 1]; + } + std::reverse(counter.begin(), counter.end()); + return counter; + } + void read_perseus_style_file(const char* perseus_style_file); + void setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions, + const std::vector<T>& top_dimensional_cells); + Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector<bool> directions); + Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, std::vector<bool> directions); + Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& top_dimensional_cells, + std::vector<bool> directions); +}; + +template <typename T> +void Bitmap_cubical_complex_base<T>::put_data_to_bins(size_t number_of_bins) { + bool bdg = 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) { + 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) { + 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; + std::pair< T, T > min_max = this->min_max_filtration(); + + 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) { + 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) { + std::cerr << "After binning : " << this->data[i] << std::endl; + getchar(); + } + } +} + +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) { + 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]; + } + return min_max; +} + +template <typename K> +std::ostream& operator<<(std::ostream & out, const Bitmap_cubical_complex_base<K>& b) { + for (typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator + it = b.all_cells_const_begin(); it != b.all_cells_const_end(); ++it) { + out << *it << " "; + } + return out; +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base +(const std::vector<unsigned>& sizes) { + this->set_up_containers(sizes); +} + +template <typename T> +void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list(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) { + 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::endl; + throw("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."); + } + + Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + 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; + } + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base +(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) { + this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells); +} + +template <typename T> +void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus_style_file) { + bool dbg = false; + std::ifstream inFiltration; + inFiltration.open(perseus_style_file); + unsigned dimensionOfData; + inFiltration >> dimensionOfData; + + if (dbg) { + std::cerr << "dimensionOfData : " << dimensionOfData << std::endl; + getchar(); + } + + std::vector<unsigned> sizes; + sizes.reserve(dimensionOfData); + for (size_t i = 0; i != dimensionOfData; ++i) { + unsigned size_in_this_dimension; + inFiltration >> size_in_this_dimension; + sizes.push_back(size_in_this_dimension); + if (dbg) { + std::cerr << "size_in_this_dimension : " << size_in_this_dimension << std::endl; + } + } + this->set_up_containers(sizes); + + 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; + 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()) + << " get the value : " << filtrationLevel << std::endl; + } + this->get_cell_data(*it) = filtrationLevel; + ++it; + } + inFiltration.close(); + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file, + std::vector<bool> directions) { + // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary + // conditions. + // It ignores the last parameter of the function. + this->read_perseus_style_file(perseus_style_file); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, + std::vector<bool> directions) { + // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary + // conditions. + // It ignores the last parameter of the function. + this->set_up_containers(sizes); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& top_dimensional_cells, + std::vector<bool> directions) { + // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary + // conditions. + // It ignores the last parameter of the function. + this->setup_bitmap_based_on_top_dimensional_cells_list(dimensions, top_dimensional_cells); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file) { + this->read_perseus_style_file(perseus_style_file); +} + +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; + + // Speed traded of for memory. Check if it is better in practice. + boundary_elements.reserve(this->dimension()*2); + + size_t cell1 = cell; + for (size_t i = this->multipliers.size(); i != 0; --i) { + unsigned position = cell1 / this->multipliers[i - 1]; + if (position % 2 == 1) { + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + } + cell1 = cell1 % this->multipliers[i - 1]; + } + return boundary_elements; +} + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(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) { + unsigned position = cell1 / this->multipliers[i - 1]; + if (position % 2 == 0) { + if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) { + coboundary_elements.push_back(cell - this->multipliers[i - 1]); + } + if ( + (cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) { + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + } + } + cell1 = cell1 % this->multipliers[i - 1]; + } + return coboundary_elements; +} + +template <typename T> +unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(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) { + unsigned position = cell / this->multipliers[i - 1]; + + if (dbg) { + std::cerr << "i-1 :" << i - 1 << std::endl; + 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) { + if (dbg) std::cerr << "Nonzero length in this direction \n"; + dimension++; + } + cell = cell % this->multipliers[i - 1]; + } + return dimension; +} + +template <typename T> +inline T& Bitmap_cubical_complex_base<T>::get_cell_data(size_t cell) { + return this->data[cell]; +} + +template <typename T> +void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { + bool dbg = false; + + // 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::vector<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. + typename Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) { + indices_to_consider.push_back(it.compute_index_in_bitmap()); + } + + 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) { + 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) { + 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) { + new_indices_to_consider.push_back(bd[boundaryIt]); + is_this_cell_considered[ bd[boundaryIt] ] = true; + } + } + } + indices_to_consider.swap(new_indices_to_consider); + } +} + +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) { + if (first.first.first < second.first.first) { + return true; + } else { + if (first.first.first > second.first.first) { + return false; + } + // in this case first.first.first == second.first.first, so we need to compare dimensions + return first.second < second.second; + } +} + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_BASE_H_ 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 new file mode 100644 index 00000000..a446c0e8 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h @@ -0,0 +1,309 @@ +/* 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): Pawel Dlotko + * + * Copyright (C) 2015 INRIA Sophia-Saclay (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 BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_ +#define BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_ + +#include <gudhi/Bitmap_cubical_complex_base.h> + +#include <cmath> +#include <limits> // for numeric_limits<> +#include <vector> + +namespace Gudhi { + +namespace Cubical_complex { + +// in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic +// boundary conditions). But, we set up the iterators and the procedures to compute boundary and coboundary in the way +// that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all +// the cells on the right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base + +/** + * @class Bitmap_cubical_complex_periodic_boundary_conditions_base + * @brief Cubical complex with periodic boundary conditions represented as a bitmap. + * @ingroup cubical_complex + */ +/** + * This is a class implementing a bitmap data structure with periodic boundary conditions. Most of the functions are + * identical to the functions from Bitmap_cubical_complex_base. + * The ones that needed to be updated are the constructors and get_boundary_of_a_cell and get_coboundary_of_a_cell. + */ +template <typename T> +class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T> { + public: + // constructors that take an extra parameter: + + /** + * Default constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base() { } + /** + * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following + * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of booleans. If + * at i-th position of this vector there is true value, that means that periodic boundary conditions are to be + * imposed in this direction. In case of false, the periodic boundary conditions will not be imposed in the direction + * i. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); + /** + * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the name of Perseus + * style file as an input. Please consult the documentation about the specification of the file. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseusStyleFile); + /** + * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following + * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of top dimensional + * cells (ordered lexicographically) and (3) vector of booleans. If at i-th position of this vector there is true + * value, that means that periodic boundary conditions are to be imposed in this direction. In case of false, the + * periodic boundary conditions will not be imposed in the direction i. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector< bool >& directions_in_which_periodic_b_cond_are_to_be_imposed); + + /** + * Destructor of the Bitmap_cubical_complex_periodic_boundary_conditions_base class. + **/ + virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base() {} + + // overwritten methods co compute boundary and coboundary + /** + * A version of a function that return boundary of a given cell for an object of + * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + */ + virtual std::vector< size_t > get_boundary_of_a_cell(size_t cell) const; + + /** + * A version of a function that return coboundary of a given cell for an object of + * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + */ + virtual std::vector< size_t > get_coboundary_of_a_cell(size_t cell) const; + + protected: + std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed; + + void set_up_containers(const std::vector<unsigned>& sizes) { + unsigned multiplier = 1; + for (size_t i = 0; i != sizes.size(); ++i) { + this->sizes.push_back(sizes[i]); + this->multipliers.push_back(multiplier); + + if (directions_in_which_periodic_b_cond_are_to_be_imposed[i]) { + multiplier *= 2 * sizes[i]; + } else { + multiplier *= 2 * sizes[i] + 1; + } + } + // std::reverse( this->sizes.begin() , this->sizes.end() ); + this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max()); + this->total_number_of_cells = multiplier; + } + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes); + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells); + void construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); +}; + +template <typename T> +void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { + 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; + 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; + } + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { + 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(sizes); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseus_style_file) { + // for Perseus style files: + bool dbg = false; + + std::ifstream inFiltration; + inFiltration.open(perseus_style_file); + unsigned dimensionOfData; + inFiltration >> dimensionOfData; + + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData, false); + + std::vector<unsigned> sizes; + sizes.reserve(dimensionOfData); + for (size_t i = 0; i != dimensionOfData; ++i) { + int size_in_this_dimension; + inFiltration >> size_in_this_dimension; + if (size_in_this_dimension < 0) { + this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] = true; + } + sizes.push_back(abs(size_in_this_dimension)); + } + this->set_up_containers(sizes); + + typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator it(*this); + it = this->top_dimensional_cells_iterator_begin(); + + while (!inFiltration.eof()) { + double filtrationLevel; + inFiltration >> filtrationLevel; + if (inFiltration.eof())break; + + 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()) + << " get the value : " << filtrationLevel << std::endl; + } + this->get_cell_data(*it) = filtrationLevel; + ++it; + } + inFiltration.close(); + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes) { + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(), false); + this->set_up_containers(sizes); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells) { + std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(), false); + this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells, + directions_in_which_periodic_b_cond_are_to_be_imposed); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>:: +Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { + this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells, + directions_in_which_periodic_b_cond_are_to_be_imposed); +} + +// ***********************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 { + bool dbg = false; + if (dbg) { + std::cerr << "Computations of boundary of a cell : " << cell << std::endl; + } + + std::vector< size_t > boundary_elements; + size_t cell1 = cell; + for (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) { + // if there are no periodic boundary conditions in this direction, we do not have to do anything. + if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) { + // std::cerr << "A\n"; + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + if (dbg) { + std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " "; + } + } else { + // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end. + if (position != 2 * this->sizes[ i - 1 ] - 1) { + // std::cerr << "B\n"; + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + if (dbg) { + std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " "; + } + } else { + // std::cerr << "C\n"; + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ]); + if (dbg) { + std::cerr << cell - this->multipliers[ i - 1 ] << " " << + cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ] << " "; + } + } + } + } + cell1 = cell1 % this->multipliers[i - 1]; + } + return boundary_elements; +} + +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<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) { + 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) { + if (!this->directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) { + // no periodic boundary conditions in this direction + if ((counter[i - 1] != 0) && (cell > this->multipliers[i - 1])) { + coboundary_elements.push_back(cell - this->multipliers[i - 1]); + } + if ((counter[i - 1] != 2 * this->sizes[i - 1]) && (cell + this->multipliers[i - 1] < this->data.size())) { + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + } + } else { + // we want to have periodic boundary conditions in this direction + if (counter[i - 1] != 0) { + coboundary_elements.push_back(cell - this->multipliers[i - 1]); + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + } else { + // in this case counter[i-1] == 0. + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + coboundary_elements.push_back(cell + (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[i - 1]); + } + } + } + + cell1 = cell1 % this->multipliers[i - 1]; + } + return coboundary_elements; +} + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_ |