/* 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
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*/
#ifndef BITMAP_CUBICAL_COMPLEX_BASE_H_
#define BITMAP_CUBICAL_COMPLEX_BASE_H_
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
namespace Gudhi {
namespace cubical_complex {
/**
* @brief Cubical complex represented as a bitmap, class with basic implementation.
* @ingroup cubical_complex
* @details 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
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, 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& 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& dimensions, const std::vector& 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.
* The boundary elements are guaranteed to be returned so that the
* incidence coefficients of boundary elements are alternating.
*/
virtual inline std::vector 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
* 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.
* Note that unlike in the case of boundary, over here the elements are
* not guaranteed to be returned with alternating incidence numbers.
*
**/
virtual inline std::vector get_coboundary_of_a_cell(std::size_t cell) const;
/**
* This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of
* dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence
* between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$.
* Note that first parameter is a cube of dimension n,
* and the second parameter is an adjusted cube in dimension n-1.
* Given \f$A = [b_1,e_1] \times \ldots \ [b_{j-1},e_{j-1}] \times [b_{j},e_{j}] \times [b_{j+1},e_{j+1}] \times \ldots
*\times [b_{n},e_{n}] \f$
* such that \f$ b_{j} \neq e_{j} \f$
* and \f$B = [b_1,e_1] \times \ldots \ [b_{j-1},e_{j-1}] \times [a,a] \times [b_{j+1},e_{j+1}] \times \ldots \times
*[b_{n},e_{n}] \f$
* where \f$ a = b_{j}\f$ or \f$ a = e_{j}\f$, the incidence between \f$A\f$ and \f$B\f$
* computed by this procedure is given by formula:
* \f$ c\ (-1)^{\sum_{i=1}^{j-1} dim [b_{i},e_{i}]} \f$
* Where \f$ dim [b_{i},e_{i}] = 0 \f$ if \f$ b_{i}=e_{i} \f$ and 1 in other case.
* c is -1 if \f$ a = b_{j}\f$ and 1 if \f$ a = e_{j}\f$.
* @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(std::size_t coface, std::size_t face) const {
// first get the counters for coface and face:
std::vector coface_counter = this->compute_counter_for_given_cell(coface);
std::vector 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;
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;
}
if (coface_counter[i] != face_counter[i]) {
if (number_of_position_in_which_counters_do_not_agree != -1) {
std::cout << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
throw std::logic_error(
"Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.");
}
number_of_position_in_which_counters_do_not_agree = i;
}
}
int incidence = 1;
if (number_of_full_faces_that_comes_before % 2) incidence = -1;
// if the face cell is on the right from coface cell:
if (coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 ==
face_counter[number_of_position_in_which_counters_do_not_agree]) {
incidence *= -1;
}
return incidence;
}
/**
* In the case of get_dimension_of_a_cell function, the output is a non-negative integer
* indicating the dimension of a cell.
* Note that unlike in the case of boundary, over here the elements are
* not guaranteed to be returned with alternating incidence numbers.
* To compute incidence between cells use compute_incidence_between_cells
* procedure
**/
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.
* 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(std::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
friend std::ostream& operator<<(std::ostream& os, const Bitmap_cubical_complex_base& 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( 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(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
* 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 min_max_filtration();
// ITERATORS
/**
* @brief 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 {
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.
*/
std::size_t operator*() { return this->counter; }
friend class Bitmap_cubical_complex_base;
protected:
std::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;
}
/**
* @brief 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* b;
};
All_cells_range all_cells_range() { return All_cells_range(this); }
/**
* Boundary_range class provides ranges for boundary iterators.
**/
typedef typename std::vector::const_iterator Boundary_iterator;
typedef typename std::vector 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(std::size_t sh) { return this->get_boundary_of_a_cell(sh); }
/**
* Coboundary_range class provides ranges for boundary iterators.
**/
typedef typename std::vector::const_iterator Coboundary_iterator;
typedef typename std::vector 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(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
* in the structure (i.e. in lexicographical order)
**/
class Top_dimensional_cells_iterator : std::iterator {
public:
Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) {
this->counter = std::vector(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:
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 (std::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 (std::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.
*/
std::size_t operator*() { return this->compute_index_in_bitmap(); }
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 (std::size_t i = 0; i != this->counter.size(); ++i) {
std::cout << this->counter[i] << " ";
}
}
friend class Bitmap_cubical_complex_base;
protected:
std::vector 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 (std::size_t i = 0; i != this->dimension(); ++i) {
a.counter[i] = this->sizes[i] - 1;
}
a.counter[0]++;
return a;
}
/**
* @brief 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* b;
};
Top_dimensional_cells_range top_dimensional_cells_range() { return Top_dimensional_cells_range(this); }
//****************************************************************************************************************//
//****************************************************************************************************************//
//****************************************************************************************************************//
//****************************************************************************************************************//
inline std::size_t number_cells() const { return this->total_number_of_cells; }
//****************************************************************************************************************//
//****************************************************************************************************************//
//****************************************************************************************************************//
//****************************************************************************************************************//
protected:
std::vector sizes;
std::vector multipliers;
std::vector data;
std::size_t total_number_of_cells;
void set_up_containers(const std::vector& sizes) {
unsigned multiplier = 1;
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;
}
this->data = std::vector(multiplier, std::numeric_limits::infinity());
this->total_number_of_cells = multiplier;
}
std::size_t compute_position_in_bitmap(const std::vector& 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 compute_counter_for_given_cell(std::size_t cell) const {
std::vector counter;
counter.reserve(this->sizes.size());
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];
}
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& sizes_in_following_directions,
const std::vector& top_dimensional_cells);
Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector directions);
Bitmap_cubical_complex_base(const std::vector& sizes, std::vector directions);
Bitmap_cubical_complex_base(const std::vector& dimensions, const std::vector& top_dimensional_cells,
std::vector directions);
};
template
void Bitmap_cubical_complex_base::put_data_to_bins(std::size_t number_of_bins) {
bool dbg = false;
std::pair 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 (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 (dbg) {
std::cerr << "After binning : " << this->data[i] << std::endl;
}
}
}
template
void Bitmap_cubical_complex_base::put_data_to_bins(T diameter_of_bin) {
bool dbg = false;
std::pair min_max = this->min_max_filtration();
std::size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin;
// now put the data into the appropriate bins:
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 (dbg) {
std::cerr << "After binning : " << this->data[i] << std::endl;
}
}
}
template
std::pair Bitmap_cubical_complex_base::min_max_filtration() {
std::pair min_max(std::numeric_limits::infinity(), -std::numeric_limits::infinity());
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];
}
return min_max;
}
template
std::ostream& operator<<(std::ostream& out, const Bitmap_cubical_complex_base& b) {
for (typename Bitmap_cubical_complex_base::all_cells_const_iterator it = b.all_cells_const_begin();
it != b.all_cells_const_end(); ++it) {
out << *it << " ";
}
return out;
}
template
Bitmap_cubical_complex_base::Bitmap_cubical_complex_base(const std::vector& sizes) {
this->set_up_containers(sizes);
}
template
void Bitmap_cubical_complex_base::setup_bitmap_based_on_top_dimensional_cells_list(
const std::vector& sizes_in_following_directions, const std::vector& top_dimensional_cells) {
this->set_up_containers(sizes_in_following_directions);
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 "
<< "sizes_in_following_directions, std::vector 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 sizes_in_following_directions,"
"std::vector 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::Top_dimensional_cells_iterator it(*this);
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;
}
this->impose_lower_star_filtration();
}
template
Bitmap_cubical_complex_base::Bitmap_cubical_complex_base(const std::vector& sizes_in_following_directions,
const std::vector& top_dimensional_cells) {
this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells);
}
template
void Bitmap_cubical_complex_base::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;
}
std::vector sizes;
sizes.reserve(dimensionOfData);
// 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;
}
}
this->set_up_containers(sizes);
Bitmap_cubical_complex_base::Top_dimensional_cells_iterator it(*this);
it = this->top_dimensional_cells_iterator_begin();
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())
<< " get the value : " << filtrationLevel << std::endl;
}
this->get_cell_data(*it) = filtrationLevel;
++it;
}
inFiltration.close();
this->impose_lower_star_filtration();
}
template
Bitmap_cubical_complex_base::Bitmap_cubical_complex_base(const char* perseus_style_file,
std::vector 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
Bitmap_cubical_complex_base::Bitmap_cubical_complex_base(const std::vector& sizes,
std::vector 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
Bitmap_cubical_complex_base::Bitmap_cubical_complex_base(const std::vector& dimensions,
const std::vector& top_dimensional_cells,
std::vector 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
Bitmap_cubical_complex_base::Bitmap_cubical_complex_base(const char* perseus_style_file) {
this->read_perseus_style_file(perseus_style_file);
}
template
std::vector Bitmap_cubical_complex_base::get_boundary_of_a_cell(std::size_t cell) const {
std::vector boundary_elements;
// Speed traded of for memory. Check if it is better in practice.
boundary_elements.reserve(this->dimension() * 2);
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) {
boundary_elements.push_back(cell + this->multipliers[i - 1]);
boundary_elements.push_back(cell - this->multipliers[i - 1]);
} else {
boundary_elements.push_back(cell - this->multipliers[i - 1]);
boundary_elements.push_back(cell + this->multipliers[i - 1]);
}
++sum_of_dimensions;
}
cell1 = cell1 % this->multipliers[i - 1];
}
return boundary_elements;
}
template
std::vector Bitmap_cubical_complex_base::get_coboundary_of_a_cell(std::size_t cell) const {
std::vector counter = this->compute_counter_for_given_cell(cell);
std::vector 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)) {
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
unsigned Bitmap_cubical_complex_base::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 (std::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;
}
if (position % 2 == 1) {
if (dbg) std::cerr << "Nonzero length in this direction \n";
dimension++;
}
cell = cell % this->multipliers[i - 1];
}
return dimension;
}
template
inline T& Bitmap_cubical_complex_base::get_cell_data(std::size_t cell) {
return this->data[cell];
}
template
void Bitmap_cubical_complex_base::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 is_this_cell_considered(this->data.size(), false);
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 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::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 (std::size_t i = 0; i != indices_to_consider.size(); ++i) {
std::cout << indices_to_consider[i] << " ";
}
}
std::vector new_indices_to_consider;
for (std::size_t i = 0; i != indices_to_consider.size(); ++i) {
std::vector 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;
}
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;
}
}
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
bool compareFirstElementsOfTuples(const std::pair, char>& first,
const std::pair, 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 Cubical_complex = cubical_complex;
} // namespace Gudhi
#endif // BITMAP_CUBICAL_COMPLEX_BASE_H_