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
16#include <cmath>
17#include <limits> // for numeric_limits<>
18#include <vector>
19#include <stdexcept>
20#include <cstddef>
21
22namespace Gudhi {
23
24namespace cubical_complex {
25
26// in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic
27// boundary conditions). But, we set up the iterators and the procedures to compute boundary and coboundary in the way
28// that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all
29// the cells on the right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base
30
39template <typename T>
41 public:
42 // constructors that take an extra parameter:
43
56 const std::vector<unsigned>& sizes,
57 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
71 const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
72 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
73
78
79 // overwritten methods co compute boundary and coboundary
86 virtual std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const override;
87
96 virtual std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const override;
97
117 virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) const override {
118 // first get the counters for coface and face:
119 std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface);
120 std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face);
121
122 // coface_counter and face_counter should agree at all positions except from one:
123 int number_of_position_in_which_counters_do_not_agree = -1;
124 std::size_t number_of_full_faces_that_comes_before = 0;
125 for (std::size_t i = 0; i != coface_counter.size(); ++i) {
126 if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) {
127 ++number_of_full_faces_that_comes_before;
128 }
129 if (coface_counter[i] != face_counter[i]) {
130 if (number_of_position_in_which_counters_do_not_agree != -1) {
131 std::cerr << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
132 throw std::logic_error(
133 "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.");
134 }
135 number_of_position_in_which_counters_do_not_agree = i;
136 }
137 }
138
139 int incidence = 1;
140 if (number_of_full_faces_that_comes_before % 2) incidence = -1;
141 // if the face cell is on the right from coface cell:
142 if ((coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 ==
143 face_counter[number_of_position_in_which_counters_do_not_agree]) ||
144 ((coface_counter[number_of_position_in_which_counters_do_not_agree] != 1) &&
145 (face_counter[number_of_position_in_which_counters_do_not_agree] == 0))) {
146 incidence *= -1;
147 }
148
149 return incidence;
150 }
151
152 protected:
153 std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed;
154
155 void set_up_containers(const std::vector<unsigned>& sizes) {
156 unsigned multiplier = 1;
157 for (std::size_t i = 0; i != sizes.size(); ++i) {
158 this->sizes.push_back(sizes[i]);
159 this->multipliers.push_back(multiplier);
160
161 if (directions_in_which_periodic_b_cond_are_to_be_imposed[i]) {
162 multiplier *= 2 * sizes[i];
163 } else {
164 multiplier *= 2 * sizes[i] + 1;
165 }
166 }
167 // std::reverse( this->sizes.begin() , this->sizes.end() );
168 this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity());
169 this->total_number_of_cells = multiplier;
170 }
171 Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes);
172 Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
173 const std::vector<T>& topDimensionalCells);
174
178 void construct_complex_based_on_top_dimensional_cells(
179 const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
180 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
181};
182
183template <typename T>
184void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells(
185 const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
186 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
187 this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed;
188 this->set_up_containers(dimensions);
189
190 std::size_t i = 0;
191 for (auto it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
192 this->get_cell_data(*it) = topDimensionalCells[i];
193 ++i;
194 }
195 this->impose_lower_star_filtration();
196}
197
198template <typename T>
200 const std::vector<unsigned>& sizes,
201 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
202 this->directions_in_which_periodic_b_cond_are_to_be_imposed(directions_in_which_periodic_b_cond_are_to_be_imposed);
203 this->set_up_containers(sizes);
204}
205
206template <typename T>
208 const char* perseus_style_file) {
209 // for Perseus style files:
210 bool dbg = false;
211
212 std::ifstream inFiltration;
213 inFiltration.open(perseus_style_file);
214 unsigned dimensionOfData;
215 inFiltration >> dimensionOfData;
216
217 this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData, false);
218
219 std::vector<unsigned> sizes;
220 sizes.reserve(dimensionOfData);
221 for (std::size_t i = 0; i != dimensionOfData; ++i) {
222 int size_in_this_dimension;
223 inFiltration >> size_in_this_dimension;
224 if (size_in_this_dimension < 0) {
225 this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] = true;
226 }
227 sizes.push_back(abs(size_in_this_dimension));
228 }
229 this->set_up_containers(sizes);
230
232 it = this->top_dimensional_cells_iterator_begin();
233
234 while (!inFiltration.eof()) {
235 double filtrationLevel;
236 inFiltration >> filtrationLevel;
237 if (inFiltration.eof()) break;
238
239 if (dbg) {
240 std::clog << "Cell of an index : " << it.compute_index_in_bitmap()
241 << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap())
242 << " get the value : " << filtrationLevel << std::endl;
243 }
244 this->get_cell_data(*it) = filtrationLevel;
245 ++it;
246 }
247 inFiltration.close();
248 this->impose_lower_star_filtration();
249}
250
251template <typename T>
253 const std::vector<unsigned>& sizes) {
254 this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(), false);
255 this->set_up_containers(sizes);
256}
257
258template <typename T>
260 const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells) {
261 std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(), false);
262 this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells,
263 directions_in_which_periodic_b_cond_are_to_be_imposed);
264}
265
266template <typename T>
268 const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
269 const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
270 this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells,
271 directions_in_which_periodic_b_cond_are_to_be_imposed);
272}
273
274// ***********************Methods************************ //
275
276template <typename T>
278 std::size_t cell) const {
279 bool dbg = false;
280 if (dbg) {
281 std::clog << "Computations of boundary of a cell : " << cell << std::endl;
282 }
283
284 std::vector<std::size_t> boundary_elements;
285 boundary_elements.reserve(this->dimension() * 2);
286 std::size_t cell1 = cell;
287 std::size_t sum_of_dimensions = 0;
288
289 for (std::size_t i = this->multipliers.size(); i != 0; --i) {
290 unsigned position = cell1 / this->multipliers[i - 1];
291 // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction.
292 if (position % 2 == 1) {
293 // if there are no periodic boundary conditions in this direction, we do not have to do anything.
294 if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
295 if (sum_of_dimensions % 2) {
296 boundary_elements.push_back(cell - this->multipliers[i - 1]);
297 boundary_elements.push_back(cell + this->multipliers[i - 1]);
298 } else {
299 boundary_elements.push_back(cell + this->multipliers[i - 1]);
300 boundary_elements.push_back(cell - this->multipliers[i - 1]);
301 }
302 if (dbg) {
303 std::clog << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
304 }
305 } else {
306 // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end.
307 if (position != 2 * this->sizes[i - 1] - 1) {
308 if (sum_of_dimensions % 2) {
309 boundary_elements.push_back(cell - this->multipliers[i - 1]);
310 boundary_elements.push_back(cell + this->multipliers[i - 1]);
311 } else {
312 boundary_elements.push_back(cell + this->multipliers[i - 1]);
313 boundary_elements.push_back(cell - this->multipliers[i - 1]);
314 }
315 if (dbg) {
316 std::clog << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
317 }
318 } else {
319 if (sum_of_dimensions % 2) {
320 boundary_elements.push_back(cell - this->multipliers[i - 1]);
321 boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
322 } else {
323 boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
324 boundary_elements.push_back(cell - this->multipliers[i - 1]);
325 }
326 if (dbg) {
327 std::clog << cell - this->multipliers[i - 1] << " "
328 << cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1] << " ";
329 }
330 }
331 }
332 ++sum_of_dimensions;
333 }
334 cell1 = cell1 % this->multipliers[i - 1];
335 }
336 return boundary_elements;
337}
338
339template <typename T>
341 std::size_t cell) const {
342 std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
343 std::vector<std::size_t> coboundary_elements;
344 std::size_t cell1 = cell;
345 for (std::size_t i = this->multipliers.size(); i != 0; --i) {
346 unsigned position = cell1 / this->multipliers[i - 1];
347 // if the cell has zero length in this direction, then it will have cbd in this direction.
348 if (position % 2 == 0) {
349 if (!this->directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
350 // no periodic boundary conditions in this direction
351 if ((counter[i - 1] != 0) && (cell > this->multipliers[i - 1])) {
352 coboundary_elements.push_back(cell - this->multipliers[i - 1]);
353 }
354 if ((counter[i - 1] != 2 * this->sizes[i - 1]) && (cell + this->multipliers[i - 1] < this->data.size())) {
355 coboundary_elements.push_back(cell + this->multipliers[i - 1]);
356 }
357 } else {
358 // we want to have periodic boundary conditions in this direction
359 if (counter[i - 1] != 0) {
360 coboundary_elements.push_back(cell - this->multipliers[i - 1]);
361 coboundary_elements.push_back(cell + this->multipliers[i - 1]);
362 } else {
363 // in this case counter[i-1] == 0.
364 coboundary_elements.push_back(cell + this->multipliers[i - 1]);
365 coboundary_elements.push_back(cell + (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
366 }
367 }
368 }
369
370 cell1 = cell1 % this->multipliers[i - 1];
371 }
372 return coboundary_elements;
373}
374
375} // namespace cubical_complex
376
377namespace Cubical_complex = cubical_complex;
378
379} // namespace Gudhi
380
381#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:358
Cubical complex represented as a bitmap, class with basic implementation.
Definition: Bitmap_cubical_complex_base.h:53
Cubical complex with periodic boundary conditions represented as a bitmap.
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:40
Bitmap_cubical_complex_periodic_boundary_conditions_base()
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:47
virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base()
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:77
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:340
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:117
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:277
This is an implementation of a counter being a vector of integers.
Definition: counter.h:32