Loading...
Searching...
No Matches
Bitmap_cubical_complex_base.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Pawel Dlotko
4 *
5 * Copyright (C) 2015 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
11#ifndef BITMAP_CUBICAL_COMPLEX_BASE_H_
12#define BITMAP_CUBICAL_COMPLEX_BASE_H_
13
14#include <gudhi/Debug_utils.h>
15
16#include <boost/config.hpp>
17#include <boost/iterator/counting_iterator.hpp>
18#include <boost/range/iterator_range.hpp>
19
20#include <iostream>
21#include <vector>
22#include <string>
23#include <fstream>
24#include <algorithm>
25#include <iterator>
26#include <limits>
27#include <utility>
28#include <stdexcept>
29#include <cstddef>
30#include <numeric>
31#include <functional>
32
33namespace Gudhi {
34
35namespace cubical_complex {
36
56template <typename T>
58 public:
59 typedef T filtration_type;
60
64 Bitmap_cubical_complex_base() : total_number_of_cells(0) {}
71 explicit Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes);
78 explicit Bitmap_cubical_complex_base(const char* perseus_style_file);
83 Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& cells, bool input_top_cells = true);
84
89
101 virtual inline std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const;
116 virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;
117
126 inline size_t get_top_dimensional_coface_of_a_cell(size_t splx);
127
136 inline size_t get_vertex_of_a_cell(size_t splx);
137
157 virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) const {
158 // first get the counters for coface and face:
159 std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface);
160 std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face);
161
162 // coface_counter and face_counter should agree at all positions except from one:
163 int number_of_position_in_which_counters_do_not_agree = -1;
164 std::size_t number_of_full_faces_that_comes_before = 0;
165 for (std::size_t i = 0; i != coface_counter.size(); ++i) {
166 if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) {
167 ++number_of_full_faces_that_comes_before;
168 }
169 if (coface_counter[i] != face_counter[i]) {
170 if (number_of_position_in_which_counters_do_not_agree != -1) {
171 std::cerr << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
172 throw std::logic_error(
173 "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.");
174 }
175 number_of_position_in_which_counters_do_not_agree = i;
176 }
177 }
178
179 int incidence = 1;
180 if (number_of_full_faces_that_comes_before % 2) incidence = -1;
181 // if the face cell is on the right from coface cell:
182 if (coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 ==
183 face_counter[number_of_position_in_which_counters_do_not_agree]) {
184 incidence *= -1;
185 }
186
187 return incidence;
188 }
189
198 inline unsigned get_dimension_of_a_cell(std::size_t cell) const;
199
206 inline T& get_cell_data(std::size_t cell);
207
216 void impose_lower_star_filtration(); // assume that top dimensional cells are already set.
217
222 void impose_lower_star_filtration_from_vertices(); // assume that vertices are already set.
223
227 inline unsigned dimension() const { return sizes.size(); }
228
232 inline std::size_t size() const { return this->data.size(); }
233
242 void put_data_to_bins(std::size_t number_of_bins);
243
254 void put_data_to_bins(T diameter_of_bin);
255
259 std::pair<T, T> min_max_filtration();
260
261 // ITERATORS
262
267 typedef boost::counting_iterator<std::size_t> All_cells_iterator;
268
273
278
282 typedef boost::iterator_range<All_cells_iterator> All_cells_range;
283
286
290 typedef typename std::vector<std::size_t>::const_iterator Boundary_iterator;
291 typedef typename std::vector<std::size_t> Boundary_range;
292
297 Boundary_range boundary_range(std::size_t sh) { return this->get_boundary_of_a_cell(sh); }
298
302 typedef typename std::vector<std::size_t>::const_iterator Coboundary_iterator;
303 typedef typename std::vector<std::size_t> Coboundary_range;
304
309 Coboundary_range coboundary_range(std::size_t sh) { return this->get_coboundary_of_a_cell(sh); }
310
316 public:
317 typedef std::input_iterator_tag iterator_category;
318 typedef std::size_t value_type;
319 typedef std::ptrdiff_t difference_type;
320 typedef value_type* pointer;
321 typedef value_type reference;
322
324
325 Top_dimensional_cells_iterator operator++() {
326 // first find first element of the counter that can be increased:
327 std::size_t dim = 0;
328 while ((dim != this->b->dimension()) && (this->counter[dim] == this->b->sizes[dim] - 1)) ++dim;
329
330 if (dim != this->b->dimension()) {
331 ++this->counter[dim];
332 for (std::size_t i = 0; i != dim; ++i) {
333 this->counter[i] = 0;
334 }
335 } else {
336 ++this->counter[0];
337 }
338 return *this;
339 }
340
341 Top_dimensional_cells_iterator operator++(int) {
342 Top_dimensional_cells_iterator result = *this;
343 ++(*this);
344 return result;
345 }
346
347 bool operator==(const Top_dimensional_cells_iterator& rhs) const {
348 if (this->b != rhs.b) return false;
349 if (this->counter.size() != rhs.counter.size()) return false;
350 for (std::size_t i = 0; i != this->counter.size(); ++i) {
351 if (this->counter[i] != rhs.counter[i]) return false;
352 }
353 return true;
354 }
355
356 bool operator!=(const Top_dimensional_cells_iterator& rhs) const { return !(*this == rhs); }
357
358 /*
359 * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as
360 * an argument of the following functions:
361 * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell
362 * boundary and coboundary and dimension
363 * and in function get_cell_data to get a filtration of a cell.
364 */
365 std::size_t operator*() { return this->compute_index_in_bitmap(); }
366
367 std::size_t compute_index_in_bitmap() const {
368 std::size_t index = 0;
369 for (std::size_t i = 0; i != this->counter.size(); ++i) {
370 index += (2 * this->counter[i] + 1) * this->b->multipliers[i];
371 }
372 return index;
373 }
374
375 void print_counter() const {
376 for (std::size_t i = 0; i != this->counter.size(); ++i) {
377 std::clog << this->counter[i] << " ";
378 }
379 }
380 friend class Bitmap_cubical_complex_base;
381
382 protected:
383 std::vector<std::size_t> counter;
385 };
386
392 return a;
393 }
394
400 for (std::size_t i = 0; i != this->dimension(); ++i) {
401 a.counter[i] = this->sizes[i] - 1;
402 }
403 a.counter[0]++;
404 return a;
405 }
406
411 public:
413
415
417
418 private:
420 };
421
424
425 //****************************************************************************************************************//
426 //****************************************************************************************************************//
427 //****************************************************************************************************************//
428 //****************************************************************************************************************//
429 class Vertices_iterator {
430 public:
431 typedef std::input_iterator_tag iterator_category;
432 typedef std::size_t value_type;
433 typedef std::ptrdiff_t difference_type;
434 typedef value_type* pointer;
435 typedef value_type reference;
436
437 Vertices_iterator(Bitmap_cubical_complex_base* b) : counter(b->dimension()), b(b) {}
438
439 Vertices_iterator operator++() {
440 // first find first element of the counter that can be increased:
441 std::size_t dim = 0;
442 while ((dim != this->b->dimension()) && (this->counter[dim] == this->b->sizes[dim])) ++dim;
443
444 if (dim != this->b->dimension()) {
445 ++this->counter[dim];
446 for (std::size_t i = 0; i != dim; ++i) {
447 this->counter[i] = 0;
448 }
449 } else {
450 ++this->counter[0];
451 }
452 return *this;
453 }
454
455 Vertices_iterator operator++(int) {
456 Vertices_iterator result = *this;
457 ++(*this);
458 return result;
459 }
460
461 bool operator==(const Vertices_iterator& rhs) const {
462 if (this->b != rhs.b) return false;
463 GUDHI_CHECK(this->counter.size() == rhs.counter.size(), "impossible");
464 for (std::size_t i = 0; i != this->counter.size(); ++i) {
465 if (this->counter[i] != rhs.counter[i]) return false;
466 }
467 return true;
468 }
469
470 bool operator!=(const Vertices_iterator& rhs) const { return !(*this == rhs); }
471
472 /*
473 * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as
474 * an argument of the following functions:
475 * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell
476 * boundary and coboundary and dimension
477 * and in function get_cell_data to get a filtration of a cell.
478 */
479 std::size_t operator*() const { return this->compute_index_in_bitmap(); }
480
481 std::size_t compute_index_in_bitmap() const {
482 std::size_t index = 0;
483 for (std::size_t i = 0; i != this->counter.size(); ++i) {
484 index += 2 * this->counter[i] * this->b->multipliers[i];
485 }
486 return index;
487 }
488
489 void print_counter() const {
490 for (std::size_t i = 0; i != this->counter.size(); ++i) {
491 std::clog << this->counter[i] << " ";
492 }
493 }
494 friend class Bitmap_cubical_complex_base;
495
496 protected:
497 std::vector<std::size_t> counter;
498 Bitmap_cubical_complex_base* b;
499 };
500
501 /*
502 * Function returning a Vertices_iterator to the first vertex of the bitmap.
503 */
504 Vertices_iterator vertices_iterator_begin() {
505 Vertices_iterator a(this);
506 return a;
507 }
508
509 /*
510 * Function returning a Vertices_iterator to the last vertex of the bitmap.
511 */
512 Vertices_iterator vertices_iterator_end() {
513 Vertices_iterator a(this);
514 for (std::size_t i = 0; i != this->dimension(); ++i) {
515 a.counter[i] = this->sizes[i];
516 }
517 a.counter[0]++;
518 return a;
519 }
520
521 /*
522 * @brief Vertices_iterator_range class provides ranges for Vertices_iterator_range
523 */
524 class Vertices_range {
525 public:
526 Vertices_range(Bitmap_cubical_complex_base* b) : b(b) {}
527
528 Vertices_iterator begin() { return b->vertices_iterator_begin(); }
529
530 Vertices_iterator end() { return b->vertices_iterator_end(); }
531
532 private:
534 };
535
536 /* Returns a range over all vertices. */
537 Vertices_range vertices_range() { return Vertices_range(this); }
538
539 //****************************************************************************************************************//
540 //****************************************************************************************************************//
541 //****************************************************************************************************************//
542 //****************************************************************************************************************//
543
544 inline std::size_t number_cells() const { return this->total_number_of_cells; }
545
546 //****************************************************************************************************************//
547 //****************************************************************************************************************//
548 //****************************************************************************************************************//
549 //****************************************************************************************************************//
550
552 template <class F> void for_each_vertex(F&&f) {
553 for_each_vertex_rec(f, 0, multipliers.size()-1);
554 }
555
556 protected:
557 std::vector<unsigned> sizes;
558 std::vector<unsigned> multipliers;
559 std::vector<T> data;
560 std::size_t total_number_of_cells;
561
562 template <class F> void for_each_vertex_rec(F&&f, std::size_t base, int dim);
563 void propagate_from_vertices_rec(int special_dim, int current_dim, std::size_t base);
564
565 void set_up_containers(const std::vector<unsigned>& sizes, bool is_pos_inf) {
566 // The fact that multipliers[0]=1 is relied on by optimizations in other functions
567 unsigned multiplier = 1;
568 for (std::size_t i = 0; i != sizes.size(); ++i) {
569 this->sizes.push_back(sizes[i]);
570 this->multipliers.push_back(multiplier);
571 multiplier *= 2 * sizes[i] + 1;
572 }
573 if(is_pos_inf)
574 this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity());
575 else
576 this->data = std::vector<T>(multiplier, -std::numeric_limits<T>::infinity());
577 this->total_number_of_cells = multiplier;
578 }
579
580 std::size_t compute_position_in_bitmap(const std::vector<unsigned>& counter) {
581 std::size_t position = 0;
582 for (std::size_t i = 0; i != this->multipliers.size(); ++i) {
583 position += this->multipliers[i] * counter[i];
584 }
585 return position;
586 }
587
588 std::vector<unsigned> compute_counter_for_given_cell(std::size_t cell) const {
589 std::vector<unsigned> counter;
590 counter.reserve(this->sizes.size());
591 for (std::size_t dim = this->sizes.size(); dim > 1; --dim) {
592 std::size_t quot = cell / this->multipliers[dim - 1];
593 cell = cell % this->multipliers[dim - 1];
594 counter.push_back(quot);
595 }
596 // Split out the last iteration to save a costly division by multipliers[0]=1
597 counter.push_back(cell);
598 std::reverse(counter.begin(), counter.end());
599 return counter;
600 }
601 void read_perseus_style_file(const char* perseus_style_file);
602 void setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions,
603 const std::vector<T>& top_dimensional_cells);
604 void setup_bitmap_based_on_vertices(const std::vector<unsigned>& sizes_in_following_directions,
605 const std::vector<T>& vertices);
606 Bitmap_cubical_complex_base(const char* perseus_style_file, const std::vector<bool>& directions);
607 Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, const std::vector<bool>& directions);
608 Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& cells,
609 const std::vector<bool>& directions, bool input_top_cells);
610};
611
612template <typename T>
613void Bitmap_cubical_complex_base<T>::put_data_to_bins(std::size_t number_of_bins) {
614
615 std::pair<T, T> min_max = this->min_max_filtration();
616 T dx = (min_max.second - min_max.first) / (T)number_of_bins;
617
618 // now put the data into the appropriate bins:
619 for (std::size_t i = 0; i != this->data.size(); ++i) {
620#ifdef DEBUG_TRACES
621 std::clog << "Before binning : " << this->data[i] << std::endl;
622#endif
623 this->data[i] = min_max.first + dx * (this->data[i] - min_max.first) / number_of_bins;
624#ifdef DEBUG_TRACES
625 std::clog << "After binning : " << this->data[i] << std::endl;
626#endif
627 }
628}
629
630template <typename T>
632 std::pair<T, T> min_max = this->min_max_filtration();
633
634 std::size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin;
635 // now put the data into the appropriate bins:
636 for (std::size_t i = 0; i != this->data.size(); ++i) {
637#ifdef DEBUG_TRACES
638 std::clog << "Before binning : " << this->data[i] << std::endl;
639#endif
640 this->data[i] = min_max.first + diameter_of_bin * (this->data[i] - min_max.first) / number_of_bins;
641#ifdef DEBUG_TRACES
642 std::clog << "After binning : " << this->data[i] << std::endl;
643#endif
644 }
645}
646
647template <typename T>
649 std::pair<T, T> min_max(std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity());
650 for (std::size_t i = 0; i != this->data.size(); ++i) {
651 if (this->data[i] < min_max.first) min_max.first = this->data[i];
652 if (this->data[i] > min_max.second) min_max.second = this->data[i];
653 }
654 return min_max;
655}
656
657template <typename T>
659 this->set_up_containers(sizes, true);
660}
661
662template <typename T>
664 const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) {
665 this->set_up_containers(sizes_in_following_directions, true);
666 std::size_t number_of_top_dimensional_elements = std::accumulate(std::begin(sizes_in_following_directions),
667 std::end(sizes_in_following_directions), std::size_t(1),
668 std::multiplies<std::size_t>());
669 if (number_of_top_dimensional_elements != top_dimensional_cells.size()) {
670 std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<unsigned> "
671 << "sizes_in_following_directions, std::vector<T> top_dimensional_cells ). Number of top dimensional "
672 << "elements that follow from sizes_in_following_directions vector is different from the size of "
673 << "top_dimensional_cells vector."
674 << std::endl;
675 throw std::invalid_argument(
676 "Error in constructor Bitmap_cubical_complex_base( std::vector<unsigned> sizes_in_following_directions,"
677 "std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from "
678 "sizes_in_following_directions vector is different from the size of top_dimensional_cells vector.");
679 }
680
681 std::size_t index = 0;
682 for (auto it = this->top_dimensional_cells_iterator_begin();
683 it != this->top_dimensional_cells_iterator_end(); ++it) {
684 this->get_cell_data(*it) = top_dimensional_cells[index];
685 ++index;
686 }
687 this->impose_lower_star_filtration();
688}
689
690template <typename T>
691void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_vertices(const std::vector<unsigned>& sizes_in_following_directions,
692 const std::vector<T>& vertices) {
693 std::vector<unsigned> top_cells_sizes;
694 std::transform (sizes_in_following_directions.begin(), sizes_in_following_directions.end(), std::back_inserter(top_cells_sizes),
695 [](int i){ return i-1;});
696 this->set_up_containers(top_cells_sizes, false);
697 std::size_t number_of_vertices = std::accumulate(std::begin(sizes_in_following_directions),
698 std::end(sizes_in_following_directions), (std::size_t)1,
699 std::multiplies<std::size_t>());
700 if (number_of_vertices != vertices.size()) {
701 std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<unsigned> "
702 << "sizes_in_following_directions, std::vector<T> vertices ). Number of vertices "
703 << "elements that follow from sizes_in_following_directions vector is different from the size of "
704 << "vertices vector."
705 << std::endl;
706 throw std::invalid_argument(
707 "Error in constructor Bitmap_cubical_complex_base( std::vector<unsigned> sizes_in_following_directions,"
708 "std::vector<T> vertices ). Number of vertices elements that follow from "
709 "sizes_in_following_directions vector is different from the size of vertices vector.");
710 }
711
712 for_each_vertex([this, &vertices, index=(std::size_t)0] (auto cell) mutable { get_cell_data(cell) = vertices[index++]; });
713 this->impose_lower_star_filtration_from_vertices();
714}
715
716template <typename T>
718 if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;}
719 else {
720 for (auto v : this->get_coboundary_of_a_cell(splx)){
721 if(this->get_cell_data(v) == this->get_cell_data(splx)){
722 return this->get_top_dimensional_coface_of_a_cell(v);
723 }
724 }
725 }
726 BOOST_UNREACHABLE_RETURN(-2);
727}
728
729template <typename T>
731 if (this->get_dimension_of_a_cell(splx) == 0){return splx;}
732 else {
733 for (auto v : this->get_boundary_of_a_cell(splx)){
734 if(this->get_cell_data(v) == this->get_cell_data(splx)){
735 return this->get_vertex_of_a_cell(v);
736 }
737 }
738 }
739 BOOST_UNREACHABLE_RETURN(-2);
740}
741
742template <typename T>
743Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes_in_following_directions,
744 const std::vector<T>& cells, bool input_top_cells) {
745 if (input_top_cells) {
746 this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, cells);
747 } else {
748 this->setup_bitmap_based_on_vertices(sizes_in_following_directions, cells);
749 }
750}
751
752template <typename T>
753void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus_style_file) {
754 std::ifstream inFiltration(perseus_style_file);
755 if(!inFiltration) throw std::ios_base::failure(std::string("Could not open the file ") + perseus_style_file);
756 unsigned dimensionOfData;
757 inFiltration >> dimensionOfData;
758
759#ifdef DEBUG_TRACES
760 std::clog << "dimensionOfData : " << dimensionOfData << std::endl;
761#endif
762
763 std::vector<unsigned> sizes;
764 sizes.reserve(dimensionOfData);
765 // all dimensions multiplied
766 std::size_t dimensions = 1;
767 for (std::size_t i = 0; i != dimensionOfData; ++i) {
768 unsigned size_in_this_dimension;
769 inFiltration >> size_in_this_dimension;
770 sizes.push_back(size_in_this_dimension);
771 dimensions *= size_in_this_dimension;
772#ifdef DEBUG_TRACES
773 std::clog << "size_in_this_dimension : " << size_in_this_dimension << std::endl;
774#endif
775 }
776 this->set_up_containers(sizes, true);
777
778 Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it = this->top_dimensional_cells_iterator_begin();
779
780 T filtrationLevel = 0.;
781 std::size_t filtration_counter = 0;
782 while (!inFiltration.eof()) {
783 std::string line;
784 getline(inFiltration, line);
785 if (line.length() != 0) {
786 int n = sscanf(line.c_str(), "%lf", &filtrationLevel);
787 if (n != 1) {
788 std::string perseus_error("Bad Perseus file format. This line is incorrect : " + line);
789 throw std::ios_base::failure(perseus_error.c_str());
790 }
791
792#ifdef DEBUG_TRACES
793 std::clog << "Cell of an index : " << it.compute_index_in_bitmap()
794 << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap())
795 << " get the value : " << filtrationLevel << std::endl;
796#endif
797 this->get_cell_data(*it) = filtrationLevel;
798 ++it;
799 ++filtration_counter;
800 }
801 }
802
803 if (filtration_counter != dimensions) {
804 std::string perseus_error("Bad Perseus file format. Read " + std::to_string(filtration_counter) + " expected " + \
805 std::to_string(dimensions) + " values");
806 throw std::ios_base::failure(perseus_error);
807 }
808
809 inFiltration.close();
810 this->impose_lower_star_filtration();
811}
812
813template <typename T>
815 const std::vector<bool>& directions) {
816 // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
817 // conditions.
818 // It ignores the last parameter of the function.
819 this->read_perseus_style_file(perseus_style_file);
820}
821
822template <typename T>
823Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes,
824 const std::vector<bool>& directions) {
825 // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
826 // conditions.
827 // It ignores the last parameter of the function.
828 this->set_up_containers(sizes, true);
829}
830
831template <typename T>
832Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions,
833 const std::vector<T>& cells,
834 const std::vector<bool>& directions,
835 bool input_top_cells) {
836 // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
837 // conditions.
838 // It ignores the last parameter of the function.
839 if (input_top_cells) {
840 this->setup_bitmap_based_on_top_dimensional_cells_list(dimensions, cells);
841 } else {
842 this->setup_bitmap_based_on_vertices(dimensions, cells);
843 }
844}
845
846template <typename T>
848 this->read_perseus_style_file(perseus_style_file);
849}
850
851template <typename T>
852std::vector<std::size_t> Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(std::size_t cell) const {
853 std::vector<std::size_t> boundary_elements;
854
855 // Speed traded of for memory. Check if it is better in practice.
856 boundary_elements.reserve(this->dimension() * 2);
857
858 std::size_t sum_of_dimensions = 0;
859 std::size_t cell1 = cell;
860 for (std::size_t i = this->multipliers.size(); i > 1; --i) {
861 unsigned position = cell1 / this->multipliers[i - 1];
862 cell1 = cell1 % this->multipliers[i - 1];
863 if (position % 2 == 1) {
864 if (sum_of_dimensions % 2) {
865 boundary_elements.push_back(cell + this->multipliers[i - 1]);
866 boundary_elements.push_back(cell - this->multipliers[i - 1]);
867 } else {
868 boundary_elements.push_back(cell - this->multipliers[i - 1]);
869 boundary_elements.push_back(cell + this->multipliers[i - 1]);
870 }
871 ++sum_of_dimensions;
872 }
873 }
874 // Split out the last iteration to save a costly division by multipliers[0]=1
875 if (cell1 % 2 == 1) {
876 if (sum_of_dimensions % 2) {
877 boundary_elements.push_back(cell + 1);
878 boundary_elements.push_back(cell - 1);
879 } else {
880 boundary_elements.push_back(cell - 1);
881 boundary_elements.push_back(cell + 1);
882 }
883 ++sum_of_dimensions;
884 }
885
886 return boundary_elements;
887}
888
889template <typename T>
890std::vector<std::size_t> Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(std::size_t cell) const {
891 std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
892 std::vector<std::size_t> coboundary_elements;
893 std::size_t cell1 = cell;
894 for (std::size_t i = this->multipliers.size(); i > 1; --i) {
895 // It is a bit sad to recompute those divisions when we just did them in compute_counter_for_given_cell.
896 unsigned position = cell1 / this->multipliers[i - 1];
897 cell1 = cell1 % this->multipliers[i - 1];
898 if (position % 2 == 0) {
899 if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) {
900 coboundary_elements.push_back(cell - this->multipliers[i - 1]);
901 }
902 if ((cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) {
903 coboundary_elements.push_back(cell + this->multipliers[i - 1]);
904 }
905 }
906 }
907 if (cell1 % 2 == 0) {
908 if ((cell > 1) && (counter[0] != 0)) {
909 coboundary_elements.push_back(cell - 1);
910 }
911 if ((cell + 1 < this->data.size()) && (counter[0] != 2 * this->sizes[0])) {
912 coboundary_elements.push_back(cell + 1);
913 }
914 }
915 return coboundary_elements;
916}
917
918template <typename T>
920#ifdef DEBUG_TRACES
921 std::clog << "\n\n\n Computing position o a cell of an index : " << cell << std::endl;
922#endif
923 unsigned dimension = 0;
924 for (std::size_t i = this->multipliers.size(); i > 1; --i) {
925 unsigned position = cell / this->multipliers[i - 1];
926 std::size_t newcell = cell % this->multipliers[i - 1];
927
928#ifdef DEBUG_TRACES
929 std::clog << "i-1 :" << i - 1 << std::endl;
930 std::clog << "cell : " << cell << std::endl;
931 std::clog << "position : " << position << std::endl;
932 std::clog << "multipliers[" << i - 1 << "] = " << this->multipliers[i - 1] << std::endl;
933#endif
934
935 if (position % 2 == 1) {
936#ifdef DEBUG_TRACES
937 std::clog << "Nonzero length in this direction \n";
938#endif
939 dimension++;
940 }
941 cell = newcell;
942 }
943 // Split out the last iteration to save a costly division by multipliers[0]=1
944#ifdef DEBUG_TRACES
945 std::clog << "i-1 :" << 0 << std::endl;
946 std::clog << "cell : " << cell << std::endl;
947 std::clog << "position : " << cell << std::endl;
948 std::clog << "multipliers[" << 0 << "] = " << 1 << std::endl;
949#endif
950
951 if (cell % 2 == 1) {
952#ifdef DEBUG_TRACES
953 std::clog << "Nonzero length in this direction \n";
954#endif
955 dimension++;
956 }
957
958 return dimension;
959}
960
961template <typename T>
963 return this->data[cell];
964}
965
966template <typename T>
968 // this vector will be used to check which elements have already been taken care of in imposing lower star filtration
969 std::vector<bool> is_this_cell_considered(this->data.size(), false);
970
971 std::vector<std::size_t> indices_to_consider;
972 // we assume here that we already have a filtration on the top dimensional cells and
973 // we have to extend it to lower ones.
974 for (auto it = this->top_dimensional_cells_iterator_begin();
975 it != this->top_dimensional_cells_iterator_end(); ++it) {
976 indices_to_consider.push_back(it.compute_index_in_bitmap());
977 }
978
979 while (indices_to_consider.size()) {
980#ifdef DEBUG_TRACES
981 std::clog << "indices_to_consider in this iteration \n";
982 for (auto index : indices_to_consider) {
983 std::clog << index << " ";
984 }
985#endif
986 std::vector<std::size_t> new_indices_to_consider;
987 for (auto index : indices_to_consider) {
988 std::vector<std::size_t> bd = this->get_boundary_of_a_cell(index);
989 for (auto boundary : bd) {
990#ifdef DEBUG_TRACES
991 std::clog << "filtration of a cell : " << boundary << " is : " << this->data[boundary]
992 << " while of a cell: " << index << " is: " << this->data[index]
993 << std::endl;
994#endif
995 if (this->data[boundary] > this->data[index]) {
996 this->data[boundary] = this->data[index];
997#ifdef DEBUG_TRACES
998 std::clog << "Setting the value of a cell : " << boundary
999 << " to : " << this->data[index] << std::endl;
1000#endif
1001 }
1002 if (is_this_cell_considered[boundary] == false) {
1003 new_indices_to_consider.push_back(boundary);
1004 is_this_cell_considered[boundary] = true;
1005 }
1006 }
1007 }
1008 indices_to_consider.swap(new_indices_to_consider);
1009 }
1010}
1011
1012template <typename T>
1013template <class F>
1014void Bitmap_cubical_complex_base<T>::for_each_vertex_rec(F&&f, std::size_t base, int dim) {
1015 if (dim > 0)
1016 for(std::size_t i = 0; i < sizes[dim] + 1; ++i)
1017 for_each_vertex_rec(f, base + multipliers[dim] * 2 * i, dim - 1);
1018 else
1019 for(std::size_t i = 0; i < sizes[0] + 1; ++i)
1020 f(base + 2 * i);
1021}
1022
1023template <typename T>
1025 int max_dim = multipliers.size()-1;
1026 for (int dim = max_dim; dim >= 0; --dim)
1027 propagate_from_vertices_rec(dim, max_dim, 0);
1028}
1029
1030template <typename T>
1031void Bitmap_cubical_complex_base<T>::propagate_from_vertices_rec (int special_dim, int current_dim, std::size_t base) {
1032 if (special_dim == current_dim) {
1033 propagate_from_vertices_rec(special_dim, current_dim - 1, base);
1034 return;
1035 }
1036 if (current_dim < 0) {
1037 std::size_t step = multipliers[special_dim];
1038 for(std::size_t i = 0; i < sizes[special_dim]; ++i) {
1039 std::size_t ref = base + step * 2 * i;
1040 data[ref + step] = std::max(data[ref], data[ref + 2 * step]);
1041 }
1042 return;
1043 }
1044 if (current_dim < special_dim)
1045 for(std::size_t i = 0; i < sizes[current_dim] + 1; ++i)
1046 propagate_from_vertices_rec(special_dim, current_dim - 1, base + multipliers[current_dim] * 2 * i);
1047 else
1048 for(std::size_t i = 0; i < 2 * sizes[current_dim] + 1; ++i)
1049 propagate_from_vertices_rec(special_dim, current_dim - 1, base + multipliers[current_dim] * i);
1050}
1051
1052template <typename T>
1053bool compareFirstElementsOfTuples(const std::pair<std::pair<T, std::size_t>, char>& first,
1054 const std::pair<std::pair<T, std::size_t>, char>& second) {
1055 if (first.first.first < second.first.first) {
1056 return true;
1057 } else {
1058 if (first.first.first > second.first.first) {
1059 return false;
1060 }
1061 // in this case first.first.first == second.first.first, so we need to compare dimensions
1062 return first.second < second.second;
1063 }
1064}
1065
1066} // namespace cubical_complex
1067
1068namespace Cubical_complex = cubical_complex;
1069
1070} // namespace Gudhi
1071
1072#endif // BITMAP_CUBICAL_COMPLEX_BASE_H_
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
Range corresponding to Top_dimensional_cells_iterator.
Definition: Bitmap_cubical_complex_base.h:410
Cubical complex represented as a bitmap, class with basic implementation.
Definition: Bitmap_cubical_complex_base.h:57
std::vector< std::size_t >::const_iterator Coboundary_iterator
Definition: Bitmap_cubical_complex_base.h:302
Bitmap_cubical_complex_base()
Definition: Bitmap_cubical_complex_base.h:64
void put_data_to_bins(std::size_t number_of_bins)
Definition: Bitmap_cubical_complex_base.h:613
void impose_lower_star_filtration()
Definition: Bitmap_cubical_complex_base.h:967
boost::counting_iterator< std::size_t > All_cells_iterator
Iterator through all cells in the complex (in order they appear in the structure – i....
Definition: Bitmap_cubical_complex_base.h:267
std::size_t size() const
Definition: Bitmap_cubical_complex_base.h:232
boost::iterator_range< All_cells_iterator > All_cells_range
Range corresponding to All_cells_iterator.
Definition: Bitmap_cubical_complex_base.h:282
virtual std::vector< std::size_t > get_coboundary_of_a_cell(std::size_t cell) const
Definition: Bitmap_cubical_complex_base.h:890
unsigned get_dimension_of_a_cell(std::size_t cell) const
Definition: Bitmap_cubical_complex_base.h:919
unsigned dimension() const
Definition: Bitmap_cubical_complex_base.h:227
Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin()
Definition: Bitmap_cubical_complex_base.h:390
Top_dimensional_cells_range top_dimensional_cells_range()
Definition: Bitmap_cubical_complex_base.h:423
virtual std::vector< std::size_t > get_boundary_of_a_cell(std::size_t cell) const
Definition: Bitmap_cubical_complex_base.h:852
Coboundary_range coboundary_range(std::size_t sh)
Definition: Bitmap_cubical_complex_base.h:309
All_cells_range all_cells_range() const
Definition: Bitmap_cubical_complex_base.h:285
std::pair< T, T > min_max_filtration()
Definition: Bitmap_cubical_complex_base.h:648
Top_dimensional_cells_iterator top_dimensional_cells_iterator_end()
Definition: Bitmap_cubical_complex_base.h:398
std::vector< std::size_t >::const_iterator Boundary_iterator
Definition: Bitmap_cubical_complex_base.h:290
size_t get_top_dimensional_coface_of_a_cell(size_t splx)
Definition: Bitmap_cubical_complex_base.h:717
All_cells_iterator all_cells_iterator_end() const
Definition: Bitmap_cubical_complex_base.h:277
T & get_cell_data(std::size_t cell)
Definition: Bitmap_cubical_complex_base.h:962
void impose_lower_star_filtration_from_vertices()
Definition: Bitmap_cubical_complex_base.h:1024
Boundary_range boundary_range(std::size_t sh)
Definition: Bitmap_cubical_complex_base.h:297
virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) const
Definition: Bitmap_cubical_complex_base.h:157
size_t get_vertex_of_a_cell(size_t splx)
Definition: Bitmap_cubical_complex_base.h:730
All_cells_iterator all_cells_iterator_begin() const
Definition: Bitmap_cubical_complex_base.h:272
virtual ~Bitmap_cubical_complex_base()
Definition: Bitmap_cubical_complex_base.h:88