Loading...
Searching...
No Matches
Bitmap_cubical_complex_periodic_boundary_conditions_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_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
12#define BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
13
14#include <gudhi/Bitmap_cubical_complex_base.h>
15#include <gudhi/Debug_utils.h>
16
17#include <cmath>
18#include <limits> // for numeric_limits<>
19#include <vector>
20#include <stdexcept>
21#include <cstddef>
22
23namespace Gudhi {
24
25namespace cubical_complex {
26
27// in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic
28// boundary conditions). But, we set up the iterators and the procedures to compute boundary and coboundary in the way
29// that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all
30// the cells on the right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base
31
40template <typename T>
42 public:
43 // constructors that take an extra parameter:
44
57 const std::vector<unsigned>& sizes,
58 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
63 explicit Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseusStyleFile);
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);
75
80
81 // overwritten methods co compute boundary and coboundary
88 virtual std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const override;
89
98 virtual std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const override;
99
119 virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) const override {
120 // first get the counters for coface and face:
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);
123
124 // coface_counter and face_counter should agree at all positions except from one:
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;
130 }
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.");
136 }
137 number_of_position_in_which_counters_do_not_agree = i;
138 }
139 }
140
141 int incidence = 1;
142 if (number_of_full_faces_that_comes_before % 2) incidence = -1;
143 // if the face cell is on the right from coface cell:
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))) {
148 incidence *= -1;
149 }
150
151 return incidence;
152 }
153
154 // The non-periodic code works for top dimensional cells, but not vertices.
155 class Vertices_iterator {
156 public:
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;
162
163 Vertices_iterator(Bitmap_cubical_complex_periodic_boundary_conditions_base* b) : counter(b->dimension()), b(b) {}
164
165 Vertices_iterator operator++() {
166 // first find first element of the counter that can be increased:
167 std::size_t dim = 0;
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;
169
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;
174 }
175 } else {
176 ++this->counter[0];
177 }
178 return *this;
179 }
180
181 Vertices_iterator operator++(int) {
182 Vertices_iterator result = *this;
183 ++(*this);
184 return result;
185 }
186
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;
192 }
193 return true;
194 }
195
196 bool operator!=(const Vertices_iterator& rhs) const { return !(*this == rhs); }
197
198 /*
199 * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as
200 * an argument of the following functions:
201 * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell
202 * boundary and coboundary and dimension
203 * and in function get_cell_data to get a filtration of a cell.
204 */
205 std::size_t operator*() const { return this->compute_index_in_bitmap(); }
206
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];
211 }
212 return index;
213 }
214
215 void print_counter() const {
216 for (std::size_t i = 0; i != this->counter.size(); ++i) {
217 std::clog << this->counter[i] << " ";
218 }
219 }
220 friend class Bitmap_cubical_complex_periodic_boundary_conditions_base;
221
222 protected:
223 std::vector<std::size_t> counter;
224 Bitmap_cubical_complex_periodic_boundary_conditions_base* b;
225 };
226
227 /*
228 * Function returning a Vertices_iterator to the first vertex of the bitmap.
229 */
230 Vertices_iterator vertices_iterator_begin() {
231 Vertices_iterator a(this);
232 return a;
233 }
234
235 /*
236 * Function returning a Vertices_iterator to the last vertex of the bitmap.
237 */
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];
242 }
243 a.counter[0]++;
244 return a;
245 }
246
247 /*
248 * @brief Vertices_iterator_range class provides ranges for Vertices_iterator_range
249 */
250 class Vertices_range {
251 public:
253
254 Vertices_iterator begin() { return b->vertices_iterator_begin(); }
255
256 Vertices_iterator end() { return b->vertices_iterator_end(); }
257
258 private:
260 };
261
262 /* Returns a range over all vertices. */
263 Vertices_range vertices_range() { return Vertices_range(this); }
264
270
271 protected:
272 std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed;
273
274 void set_up_containers(const std::vector<unsigned>& sizes, bool is_pos_inf) {
275 // The fact that multipliers[0]=1 is relied on by optimizations in other functions
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);
280
281 if (directions_in_which_periodic_b_cond_are_to_be_imposed[i]) {
282 multiplier *= 2 * sizes[i];
283 } else {
284 multiplier *= 2 * sizes[i] + 1;
285 }
286 }
287 // std::reverse( this->sizes.begin() , this->sizes.end() );
288 if (is_pos_inf)
289 this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity());
290 else
291 this->data = std::vector<T>(multiplier, -std::numeric_limits<T>::infinity());
292 this->total_number_of_cells = multiplier;
293 }
294 Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes);
295 Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
296 const std::vector<T>& cells,
297 bool input_top_cells);
298
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);
307};
308
309template <typename T>
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);
315
316 std::size_t i = 0;
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];
319 ++i;
320 }
321 this->impose_lower_star_filtration();
322}
323
324template <typename T>
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;
329
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);
334
335 std::size_t i = 0;
336 for (auto it = this->vertices_iterator_begin(); it != this->vertices_iterator_end(); ++it) {
337 this->get_cell_data(*it) = vertices[i];
338 ++i;
339 }
340 this->impose_lower_star_filtration_from_vertices();
341}
342
343template <typename T>
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);
349}
350
351template <typename T>
353 const char* perseus_style_file) {
354 // for Perseus style files:
355
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;
360
361 this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData, false);
362
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;
370 }
371 sizes.push_back(abs(size_in_this_dimension));
372 }
373 this->set_up_containers(sizes, true);
374
376 it = this->top_dimensional_cells_iterator_begin();
377
378 while (!inFiltration.eof()) {
379 double filtrationLevel;
380 inFiltration >> filtrationLevel;
381 if (inFiltration.eof()) break;
382
383#ifdef DEBUG_TRACES
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;
387#endif
388 this->get_cell_data(*it) = filtrationLevel;
389 ++it;
390 }
391 inFiltration.close();
392 this->impose_lower_star_filtration();
393}
394
395template <typename T>
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);
400}
401
402template <typename T>
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);
409 } else {
410 this->construct_complex_based_on_vertices(dimensions, cells,
411 directions_in_which_periodic_b_cond_are_to_be_imposed);
412 }
413}
414
415template <typename T>
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);
423 } else {
424 this->construct_complex_based_on_vertices(dimensions, cells,
425 directions_in_which_periodic_b_cond_are_to_be_imposed);
426 }
427}
428
429// ***********************Methods************************ //
430
431template <typename T>
433 std::size_t cell) const {
434#ifdef DEBUG_TRACES
435 std::clog << "Computations of boundary of a cell : " << cell << std::endl;
436#endif
437
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;
442
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];
446 // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction.
447 if (position % 2 == 1) {
448 // if there are no periodic boundary conditions in this direction, we do not have to do anything.
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]);
453 } else {
454 boundary_elements.push_back(cell + this->multipliers[i - 1]);
455 boundary_elements.push_back(cell - this->multipliers[i - 1]);
456 }
457#ifdef DEBUG_TRACES
458 std::clog << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
459#endif
460 } else {
461 // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end.
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]);
466 } else {
467 boundary_elements.push_back(cell + this->multipliers[i - 1]);
468 boundary_elements.push_back(cell - this->multipliers[i - 1]);
469 }
470#ifdef DEBUG_TRACES
471 std::clog << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
472#endif
473 } else {
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]);
477 } else {
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]);
480 }
481#ifdef DEBUG_TRACES
482 std::clog << cell - this->multipliers[i - 1] << " "
483 << cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1] << " ";
484#endif
485 }
486 }
487 ++sum_of_dimensions;
488 }
489 }
490 return boundary_elements;
491}
492
493template <typename T>
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];
502 // if the cell has zero length in this direction, then it will have cbd in this direction.
503 if (position % 2 == 0) {
504 if (!this->directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
505 // no periodic boundary conditions in this direction
506 if ((counter[i - 1] != 0) && (cell > this->multipliers[i - 1])) {
507 coboundary_elements.push_back(cell - this->multipliers[i - 1]);
508 }
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]);
511 }
512 } else {
513 // we want to have periodic boundary conditions in this direction
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]);
517 } else {
518 // in this case counter[i-1] == 0.
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]);
521 }
522 }
523 }
524 }
525 return coboundary_elements;
526}
527
528// Exact same code as the non-periodic version, duplicated for now to ensure it calls the right version of vertices_iterator_begin.
529template <typename T>
531 // this vector will be used to check which elements have already been taken care of in imposing lower star filtration
532 std::vector<bool> is_this_cell_considered(this->data.size(), false);
533
534 std::vector<std::size_t> indices_to_consider;
535 // we assume here that we already have a filtration on the vertices and
536 // we have to extend it to higher ones.
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());
540 }
541
542 while (indices_to_consider.size()) { // Iteration on the dimension
543#ifdef DEBUG_TRACES
544 std::clog << "indices_to_consider in this iteration \n";
545 for (auto index : indices_to_consider) {
546 std::clog << index << " ";
547 }
548#endif
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) {
553#ifdef DEBUG_TRACES
554 std::clog << "filtration of a cell : " << coboundary << " is : " << this->data[coboundary]
555 << " while of a cell: " << index << " is: " << this->data[index]
556 << std::endl;
557#endif
558 if (this->data[coboundary] < this->data[index]) {
559 this->data[coboundary] = this->data[index];
560#ifdef DEBUG_TRACES
561 std::clog << "Setting the value of a cell : " << coboundary
562 << " to : " << this->data[index] << std::endl;
563#endif
564 }
565 if (is_this_cell_considered[coboundary] == false) {
566 new_indices_to_consider.push_back(coboundary);
567 is_this_cell_considered[coboundary] = true;
568 }
569 }
570 }
571 indices_to_consider.swap(new_indices_to_consider);
572 }
573}
574} // namespace cubical_complex
575
576namespace Cubical_complex = cubical_complex;
577
578} // namespace Gudhi
579
580#endif // BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_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
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