Bitmap_cubical_complex_periodic_boundary_conditions_base.h
1 /* This file is part of the Gudhi Library. The Gudhi library
2  * (Geometric Understanding in Higher Dimensions) is a generic C++
3  * library for computational topology.
4  *
5  * Author(s): Pawel Dlotko
6  *
7  * Copyright (C) 2015 Inria
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <http://www.gnu.org/licenses/>.
21  */
22 
23 #ifndef BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
24 #define BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
25 
26 #include <gudhi/Bitmap_cubical_complex_base.h>
27 
28 #include <cmath>
29 #include <limits> // for numeric_limits<>
30 #include <vector>
31 #include <stdexcept>
32 #include <cstddef>
33 
34 namespace Gudhi {
35 
36 namespace cubical_complex {
37 
38 // in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic
39 // boundary conditions). But, we set up the iterators and the procedures to compute boundary and coboundary in the way
40 // that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all
41 // the cells on the right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base
42 
51 template <typename T>
53  public:
54  // constructors that take an extra parameter:
55 
68  const std::vector<unsigned>& sizes,
69  const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
83  const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
84  const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
85 
90 
91  // overwritten methods co compute boundary and coboundary
98  virtual std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const;
99 
108  virtual std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;
109 
129  virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) {
130  // first get the counters for coface and face:
131  std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface);
132  std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face);
133 
134  // coface_counter and face_counter should agree at all positions except from one:
135  int number_of_position_in_which_counters_do_not_agree = -1;
136  std::size_t number_of_full_faces_that_comes_before = 0;
137  for (std::size_t i = 0; i != coface_counter.size(); ++i) {
138  if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) {
139  ++number_of_full_faces_that_comes_before;
140  }
141  if (coface_counter[i] != face_counter[i]) {
142  if (number_of_position_in_which_counters_do_not_agree != -1) {
143  std::cout << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
144  throw std::logic_error(
145  "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.");
146  }
147  number_of_position_in_which_counters_do_not_agree = i;
148  }
149  }
150 
151  int incidence = 1;
152  if (number_of_full_faces_that_comes_before % 2) incidence = -1;
153  // if the face cell is on the right from coface cell:
154  if ((coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 ==
155  face_counter[number_of_position_in_which_counters_do_not_agree]) ||
156  ((coface_counter[number_of_position_in_which_counters_do_not_agree] != 1) &&
157  (face_counter[number_of_position_in_which_counters_do_not_agree] == 0))) {
158  incidence *= -1;
159  }
160 
161  return incidence;
162  }
163 
164  protected:
165  std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed;
166 
167  void set_up_containers(const std::vector<unsigned>& sizes) {
168  unsigned multiplier = 1;
169  for (std::size_t i = 0; i != sizes.size(); ++i) {
170  this->sizes.push_back(sizes[i]);
171  this->multipliers.push_back(multiplier);
172 
173  if (directions_in_which_periodic_b_cond_are_to_be_imposed[i]) {
174  multiplier *= 2 * sizes[i];
175  } else {
176  multiplier *= 2 * sizes[i] + 1;
177  }
178  }
179  // std::reverse( this->sizes.begin() , this->sizes.end() );
180  this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity());
181  this->total_number_of_cells = multiplier;
182  }
183  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes);
184  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
185  const std::vector<T>& topDimensionalCells);
186 
190  void construct_complex_based_on_top_dimensional_cells(
191  const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
192  const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
193 };
194 
195 template <typename T>
197  const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
198  const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
199  this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed;
200  this->set_up_containers(dimensions);
201 
202  std::size_t i = 0;
203  for (auto it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
204  this->get_cell_data(*it) = topDimensionalCells[i];
205  ++i;
206  }
208 }
209 
210 template <typename T>
212  const std::vector<unsigned>& sizes,
213  const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
214  this->directions_in_which_periodic_b_cond_are_to_be_imposed(directions_in_which_periodic_b_cond_are_to_be_imposed);
215  this->set_up_containers(sizes);
216 }
217 
218 template <typename T>
220  const char* perseus_style_file) {
221  // for Perseus style files:
222  bool dbg = false;
223 
224  std::ifstream inFiltration;
225  inFiltration.open(perseus_style_file);
226  unsigned dimensionOfData;
227  inFiltration >> dimensionOfData;
228 
229  this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData, false);
230 
231  std::vector<unsigned> sizes;
232  sizes.reserve(dimensionOfData);
233  for (std::size_t i = 0; i != dimensionOfData; ++i) {
234  int size_in_this_dimension;
235  inFiltration >> size_in_this_dimension;
236  if (size_in_this_dimension < 0) {
237  this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] = true;
238  }
239  sizes.push_back(abs(size_in_this_dimension));
240  }
241  this->set_up_containers(sizes);
242 
245 
246  while (!inFiltration.eof()) {
247  double filtrationLevel;
248  inFiltration >> filtrationLevel;
249  if (inFiltration.eof()) break;
250 
251  if (dbg) {
252  std::cerr << "Cell of an index : " << it.compute_index_in_bitmap()
253  << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap())
254  << " get the value : " << filtrationLevel << std::endl;
255  }
256  this->get_cell_data(*it) = filtrationLevel;
257  ++it;
258  }
259  inFiltration.close();
261 }
262 
263 template <typename T>
265  const std::vector<unsigned>& sizes) {
266  this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(), false);
267  this->set_up_containers(sizes);
268 }
269 
270 template <typename T>
272  const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells) {
273  std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(), false);
274  this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells,
275  directions_in_which_periodic_b_cond_are_to_be_imposed);
276 }
277 
278 template <typename T>
280  const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
281  const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
282  this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells,
283  directions_in_which_periodic_b_cond_are_to_be_imposed);
284 }
285 
286 // ***********************Methods************************ //
287 
288 template <typename T>
290  std::size_t cell) const {
291  bool dbg = false;
292  if (dbg) {
293  std::cerr << "Computations of boundary of a cell : " << cell << std::endl;
294  }
295 
296  std::vector<std::size_t> boundary_elements;
297  boundary_elements.reserve(this->dimension() * 2);
298  std::size_t cell1 = cell;
299  std::size_t sum_of_dimensions = 0;
300 
301  for (std::size_t i = this->multipliers.size(); i != 0; --i) {
302  unsigned position = cell1 / this->multipliers[i - 1];
303  // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction.
304  if (position % 2 == 1) {
305  // if there are no periodic boundary conditions in this direction, we do not have to do anything.
306  if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
307  // std::cerr << "A\n";
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::cerr << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
317  }
318  } else {
319  // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end.
320  if (position != 2 * this->sizes[i - 1] - 1) {
321  // std::cerr << "B\n";
322  if (sum_of_dimensions % 2) {
323  boundary_elements.push_back(cell - this->multipliers[i - 1]);
324  boundary_elements.push_back(cell + this->multipliers[i - 1]);
325  } else {
326  boundary_elements.push_back(cell + this->multipliers[i - 1]);
327  boundary_elements.push_back(cell - this->multipliers[i - 1]);
328  }
329  if (dbg) {
330  std::cerr << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
331  }
332  } else {
333  // std::cerr << "C\n";
334  if (sum_of_dimensions % 2) {
335  boundary_elements.push_back(cell - this->multipliers[i - 1]);
336  boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
337  } else {
338  boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
339  boundary_elements.push_back(cell - this->multipliers[i - 1]);
340  }
341  if (dbg) {
342  std::cerr << cell - this->multipliers[i - 1] << " "
343  << cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1] << " ";
344  }
345  }
346  }
347  ++sum_of_dimensions;
348  }
349  cell1 = cell1 % this->multipliers[i - 1];
350  }
351  return boundary_elements;
352 }
353 
354 template <typename T>
356  std::size_t cell) const {
357  std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
358  std::vector<std::size_t> coboundary_elements;
359  std::size_t cell1 = cell;
360  for (std::size_t i = this->multipliers.size(); i != 0; --i) {
361  unsigned position = cell1 / this->multipliers[i - 1];
362  // if the cell has zero length in this direction, then it will have cbd in this direction.
363  if (position % 2 == 0) {
364  if (!this->directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
365  // no periodic boundary conditions in this direction
366  if ((counter[i - 1] != 0) && (cell > this->multipliers[i - 1])) {
367  coboundary_elements.push_back(cell - this->multipliers[i - 1]);
368  }
369  if ((counter[i - 1] != 2 * this->sizes[i - 1]) && (cell + this->multipliers[i - 1] < this->data.size())) {
370  coboundary_elements.push_back(cell + this->multipliers[i - 1]);
371  }
372  } else {
373  // we want to have periodic boundary conditions in this direction
374  if (counter[i - 1] != 0) {
375  coboundary_elements.push_back(cell - this->multipliers[i - 1]);
376  coboundary_elements.push_back(cell + this->multipliers[i - 1]);
377  } else {
378  // in this case counter[i-1] == 0.
379  coboundary_elements.push_back(cell + this->multipliers[i - 1]);
380  coboundary_elements.push_back(cell + (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
381  }
382  }
383  }
384 
385  cell1 = cell1 % this->multipliers[i - 1];
386  }
387  return coboundary_elements;
388 }
389 
390 } // namespace cubical_complex
391 
392 namespace Cubical_complex = cubical_complex;
393 
394 } // namespace Gudhi
395 
396 #endif // BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
void impose_lower_star_filtration()
Definition: Bitmap_cubical_complex_base.h:781
unsigned dimension() const
Definition: Bitmap_cubical_complex_base.h:207
T & get_cell_data(std::size_t cell)
Definition: Bitmap_cubical_complex_base.h:776
Cubical complex with periodic boundary conditions represented as a bitmap.
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:52
Cubical complex represented as a bitmap, class with basic implementation.
Definition: Bitmap_cubical_complex_base.h:63
unsigned get_dimension_of_a_cell(std::size_t cell) const
Definition: Bitmap_cubical_complex_base.h:752
Definition: SimplicialComplexForAlpha.h:26
virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face)
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:129
This is an implementation of a counter being a vector of integers.
Definition: counter.h:44
Bitmap_cubical_complex_periodic_boundary_conditions_base()
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:59
Top_dimensional_cells_iterator top_dimensional_cells_iterator_end()
Definition: Bitmap_cubical_complex_base.h:444
Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin()
Definition: Bitmap_cubical_complex_base.h:436
virtual std::vector< std::size_t > get_boundary_of_a_cell(std::size_t cell) const
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:289
virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base()
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:89
virtual std::vector< std::size_t > get_coboundary_of_a_cell(std::size_t cell) const
Definition: Bitmap_cubical_complex_periodic_boundary_conditions_base.h:355
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
GUDHI  Version 2.3.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : GPL v3 Generated on Tue Sep 4 2018 14:32:59 for GUDHI by Doxygen 1.8.13