11#ifndef BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
12#define BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
14#include <gudhi/Bitmap_cubical_complex_base.h>
15#include <gudhi/Debug_utils.h>
25namespace cubical_complex {
57 const std::vector<unsigned>& sizes,
58 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
72 const std::vector<unsigned>& dimensions,
const std::vector<T>& cells,
73 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed,
74 bool input_top_cells);
121 std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface);
122 std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face);
125 int number_of_position_in_which_counters_do_not_agree = -1;
126 std::size_t number_of_full_faces_that_comes_before = 0;
127 for (std::size_t i = 0; i != coface_counter.size(); ++i) {
128 if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) {
129 ++number_of_full_faces_that_comes_before;
131 if (coface_counter[i] != face_counter[i]) {
132 if (number_of_position_in_which_counters_do_not_agree != -1) {
133 std::cerr <<
"Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
134 throw std::logic_error(
135 "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.");
137 number_of_position_in_which_counters_do_not_agree = i;
142 if (number_of_full_faces_that_comes_before % 2) incidence = -1;
144 if ((coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 ==
145 face_counter[number_of_position_in_which_counters_do_not_agree]) ||
146 ((coface_counter[number_of_position_in_which_counters_do_not_agree] != 1) &&
147 (face_counter[number_of_position_in_which_counters_do_not_agree] == 0))) {
155 class Vertices_iterator {
157 typedef std::input_iterator_tag iterator_category;
158 typedef std::size_t value_type;
159 typedef std::ptrdiff_t difference_type;
160 typedef value_type* pointer;
161 typedef value_type reference;
165 Vertices_iterator operator++() {
168 while ((dim != this->b->dimension()) && (this->counter[dim] == this->b->sizes[dim] - this->b->directions_in_which_periodic_b_cond_are_to_be_imposed[dim])) ++dim;
170 if (dim != this->b->dimension()) {
171 ++this->counter[dim];
172 for (std::size_t i = 0; i != dim; ++i) {
173 this->counter[i] = 0;
181 Vertices_iterator operator++(
int) {
182 Vertices_iterator result = *
this;
187 bool operator==(
const Vertices_iterator& rhs)
const {
188 if (this->b != rhs.b)
return false;
189 GUDHI_CHECK(this->counter.size() == rhs.counter.size(),
"impossible");
190 for (std::size_t i = 0; i != this->counter.size(); ++i) {
191 if (this->counter[i] != rhs.counter[i])
return false;
196 bool operator!=(
const Vertices_iterator& rhs)
const {
return !(*
this == rhs); }
205 std::size_t operator*()
const {
return this->compute_index_in_bitmap(); }
207 std::size_t compute_index_in_bitmap()
const {
208 std::size_t index = 0;
209 for (std::size_t i = 0; i != this->counter.size(); ++i) {
210 index += 2 * this->counter[i] * this->b->multipliers[i];
215 void print_counter()
const {
216 for (std::size_t i = 0; i != this->counter.size(); ++i) {
217 std::clog << this->counter[i] <<
" ";
220 friend class Bitmap_cubical_complex_periodic_boundary_conditions_base;
223 std::vector<std::size_t> counter;
224 Bitmap_cubical_complex_periodic_boundary_conditions_base* b;
230 Vertices_iterator vertices_iterator_begin() {
231 Vertices_iterator a(
this);
238 Vertices_iterator vertices_iterator_end() {
239 Vertices_iterator a(
this);
240 for (std::size_t i = 0; i != this->
dimension(); ++i) {
241 a.counter[i] = this->sizes[i] - this->directions_in_which_periodic_b_cond_are_to_be_imposed[i];
250 class Vertices_range {
254 Vertices_iterator begin() {
return b->vertices_iterator_begin(); }
256 Vertices_iterator end() {
return b->vertices_iterator_end(); }
263 Vertices_range vertices_range() {
return Vertices_range(
this); }
272 std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed;
274 void set_up_containers(
const std::vector<unsigned>& sizes,
bool is_pos_inf) {
276 unsigned multiplier = 1;
277 for (std::size_t i = 0; i != sizes.size(); ++i) {
278 this->sizes.push_back(sizes[i]);
279 this->multipliers.push_back(multiplier);
281 if (directions_in_which_periodic_b_cond_are_to_be_imposed[i]) {
282 multiplier *= 2 * sizes[i];
284 multiplier *= 2 * sizes[i] + 1;
289 this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity());
291 this->data = std::vector<T>(multiplier, -std::numeric_limits<T>::infinity());
292 this->total_number_of_cells = multiplier;
296 const std::vector<T>& cells,
297 bool input_top_cells);
302 void construct_complex_based_on_top_dimensional_cells(
303 const std::vector<unsigned>& dimensions,
const std::vector<T>& topDimensionalCells,
304 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
305 void construct_complex_based_on_vertices(
const std::vector<unsigned>& dimensions,
const std::vector<T>& vertices,
306 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
310void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells(
311 const std::vector<unsigned>& dimensions,
const std::vector<T>& topDimensionalCells,
312 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
313 this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed;
314 this->set_up_containers(dimensions,
true);
317 for (
auto it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
318 this->get_cell_data(*it) = topDimensionalCells[i];
321 this->impose_lower_star_filtration();
325void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_vertices(
326 const std::vector<unsigned>& dimensions,
const std::vector<T>& vertices,
327 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
328 this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed;
330 std::vector<unsigned> top_cells_sizes;
331 std::transform (dimensions.begin(), dimensions.end(), directions_in_which_periodic_b_cond_are_to_be_imposed.begin(),
332 std::back_inserter(top_cells_sizes), [](
unsigned i,
bool b){ return i - !b;});
333 this->set_up_containers(top_cells_sizes,
false);
336 for (
auto it = this->vertices_iterator_begin(); it != this->vertices_iterator_end(); ++it) {
337 this->get_cell_data(*it) = vertices[i];
340 this->impose_lower_star_filtration_from_vertices();
345 const std::vector<unsigned>& sizes,
346 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
347 this->directions_in_which_periodic_b_cond_are_to_be_imposed(directions_in_which_periodic_b_cond_are_to_be_imposed);
348 this->set_up_containers(sizes,
true);
353 const char* perseus_style_file) {
356 std::ifstream inFiltration(perseus_style_file);
357 if(!inFiltration)
throw std::ios_base::failure(std::string(
"Could not open the file ") + perseus_style_file);
358 unsigned dimensionOfData;
359 inFiltration >> dimensionOfData;
361 this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData,
false);
363 std::vector<unsigned> sizes;
364 sizes.reserve(dimensionOfData);
365 for (std::size_t i = 0; i != dimensionOfData; ++i) {
366 int size_in_this_dimension;
367 inFiltration >> size_in_this_dimension;
368 if (size_in_this_dimension < 0) {
369 this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] =
true;
371 sizes.push_back(abs(size_in_this_dimension));
373 this->set_up_containers(sizes,
true);
376 it = this->top_dimensional_cells_iterator_begin();
378 while (!inFiltration.eof()) {
379 double filtrationLevel;
380 inFiltration >> filtrationLevel;
381 if (inFiltration.eof())
break;
384 std::clog <<
"Cell of an index : " << it.compute_index_in_bitmap()
385 <<
" and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap())
386 <<
" get the value : " << filtrationLevel << std::endl;
388 this->get_cell_data(*it) = filtrationLevel;
391 inFiltration.close();
392 this->impose_lower_star_filtration();
397 const std::vector<unsigned>& sizes) {
398 this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(),
false);
399 this->set_up_containers(sizes,
true);
404 const std::vector<unsigned>& dimensions,
const std::vector<T>& cells,
bool input_top_cells) {
405 std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(),
false);
406 if (input_top_cells) {
407 this->construct_complex_based_on_top_dimensional_cells(dimensions, cells,
408 directions_in_which_periodic_b_cond_are_to_be_imposed);
410 this->construct_complex_based_on_vertices(dimensions, cells,
411 directions_in_which_periodic_b_cond_are_to_be_imposed);
417 const std::vector<unsigned>& dimensions,
const std::vector<T>& cells,
418 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed,
419 bool input_top_cells) {
420 if(input_top_cells) {
421 this->construct_complex_based_on_top_dimensional_cells(dimensions, cells,
422 directions_in_which_periodic_b_cond_are_to_be_imposed);
424 this->construct_complex_based_on_vertices(dimensions, cells,
425 directions_in_which_periodic_b_cond_are_to_be_imposed);
433 std::size_t cell)
const {
435 std::clog <<
"Computations of boundary of a cell : " << cell << std::endl;
438 std::vector<std::size_t> boundary_elements;
439 boundary_elements.reserve(this->dimension() * 2);
440 std::size_t cell1 = cell;
441 std::size_t sum_of_dimensions = 0;
443 for (std::size_t i = this->multipliers.size(); i != 0; --i) {
444 unsigned position = cell1 / this->multipliers[i - 1];
445 cell1 = cell1 % this->multipliers[i - 1];
447 if (position % 2 == 1) {
449 if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
450 if (sum_of_dimensions % 2) {
451 boundary_elements.push_back(cell - this->multipliers[i - 1]);
452 boundary_elements.push_back(cell + this->multipliers[i - 1]);
454 boundary_elements.push_back(cell + this->multipliers[i - 1]);
455 boundary_elements.push_back(cell - this->multipliers[i - 1]);
458 std::clog << cell - this->multipliers[i - 1] <<
" " << cell + this->multipliers[i - 1] <<
" ";
462 if (position != 2 * this->sizes[i - 1] - 1) {
463 if (sum_of_dimensions % 2) {
464 boundary_elements.push_back(cell - this->multipliers[i - 1]);
465 boundary_elements.push_back(cell + this->multipliers[i - 1]);
467 boundary_elements.push_back(cell + this->multipliers[i - 1]);
468 boundary_elements.push_back(cell - this->multipliers[i - 1]);
471 std::clog << cell - this->multipliers[i - 1] <<
" " << cell + this->multipliers[i - 1] <<
" ";
474 if (sum_of_dimensions % 2) {
475 boundary_elements.push_back(cell - this->multipliers[i - 1]);
476 boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
478 boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
479 boundary_elements.push_back(cell - this->multipliers[i - 1]);
482 std::clog << cell - this->multipliers[i - 1] <<
" "
483 << cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1] <<
" ";
490 return boundary_elements;
495 std::size_t cell)
const {
496 std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
497 std::vector<std::size_t> coboundary_elements;
498 std::size_t cell1 = cell;
499 for (std::size_t i = this->multipliers.size(); i != 0; --i) {
500 unsigned position = cell1 / this->multipliers[i - 1];
501 cell1 = cell1 % this->multipliers[i - 1];
503 if (position % 2 == 0) {
504 if (!this->directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
506 if ((counter[i - 1] != 0) && (cell > this->multipliers[i - 1])) {
507 coboundary_elements.push_back(cell - this->multipliers[i - 1]);
509 if ((counter[i - 1] != 2 * this->sizes[i - 1]) && (cell + this->multipliers[i - 1] < this->data.size())) {
510 coboundary_elements.push_back(cell + this->multipliers[i - 1]);
514 if (counter[i - 1] != 0) {
515 coboundary_elements.push_back(cell - this->multipliers[i - 1]);
516 coboundary_elements.push_back(cell + this->multipliers[i - 1]);
519 coboundary_elements.push_back(cell + this->multipliers[i - 1]);
520 coboundary_elements.push_back(cell + (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
525 return coboundary_elements;
532 std::vector<bool> is_this_cell_considered(this->data.size(),
false);
534 std::vector<std::size_t> indices_to_consider;
537 for (
auto it = this->vertices_iterator_begin();
538 it != this->vertices_iterator_end(); ++it) {
539 indices_to_consider.push_back(it.compute_index_in_bitmap());
542 while (indices_to_consider.size()) {
544 std::clog <<
"indices_to_consider in this iteration \n";
545 for (
auto index : indices_to_consider) {
546 std::clog << index <<
" ";
549 std::vector<std::size_t> new_indices_to_consider;
550 for (
auto index : indices_to_consider) {
551 std::vector<std::size_t> cbd = this->get_coboundary_of_a_cell(index);
552 for (
auto coboundary : cbd) {
554 std::clog <<
"filtration of a cell : " << coboundary <<
" is : " << this->data[coboundary]
555 <<
" while of a cell: " << index <<
" is: " << this->data[index]
558 if (this->data[coboundary] < this->data[index]) {
559 this->data[coboundary] = this->data[index];
561 std::clog <<
"Setting the value of a cell : " << coboundary
562 <<
" to : " << this->data[index] << std::endl;
565 if (is_this_cell_considered[coboundary] ==
false) {
566 new_indices_to_consider.push_back(coboundary);
567 is_this_cell_considered[coboundary] =
true;
571 indices_to_consider.swap(new_indices_to_consider);
576namespace Cubical_complex = cubical_complex;
Iterator through top dimensional cells of the complex. The cells appear in order they are stored in t...
Definition: Bitmap_cubical_complex_base.h:315
Cubical complex represented as a bitmap, class with basic implementation.
Definition: Bitmap_cubical_complex_base.h:57
unsigned dimension() const
Definition: Bitmap_cubical_complex_base.h:227
Cubical complex with periodic boundary conditions represented as a bitmap.
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:41
void impose_lower_star_filtration_from_vertices()
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:530
Bitmap_cubical_complex_periodic_boundary_conditions_base()
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:48
virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base()
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:79
virtual std::vector< std::size_t > get_coboundary_of_a_cell(std::size_t cell) const override
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:494
virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) const override
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:119
virtual std::vector< std::size_t > get_boundary_of_a_cell(std::size_t cell) const override
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:432