Bitmap_cubical_complex_periodic_boundary_conditions_base.h

/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Pawel Dlotko
 *
 *    Copyright (C) 2015 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */

#ifndef BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
#define BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_

#include <gudhi/Bitmap_cubical_complex_base.h>

#include <cmath>
#include <limits>  // for numeric_limits<>
#include <vector>
#include <stdexcept>
#include <cstddef>

namespace Gudhi {

namespace cubical_complex {

// in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic
// boundary conditions). But, we set up the iterators and the procedures to compute boundary and coboundary in the way
// that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all
// the cells on the right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base

template <typename T>
class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T> {
 public:
  // constructors that take an extra parameter:

  Bitmap_cubical_complex_periodic_boundary_conditions_base() {}
  Bitmap_cubical_complex_periodic_boundary_conditions_base(
      const std::vector<unsigned>& sizes,
      const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
  Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseusStyleFile);
  Bitmap_cubical_complex_periodic_boundary_conditions_base(
      const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
      const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);

  virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base() {}

  // overwritten methods co compute boundary and coboundary
  virtual std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const;

  virtual std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;

  virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) {
    // first get the counters for coface and face:
    std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface);
    std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face);

    // coface_counter and face_counter should agree at all positions except from one:
    int number_of_position_in_which_counters_do_not_agree = -1;
    std::size_t number_of_full_faces_that_comes_before = 0;
    for (std::size_t i = 0; i != coface_counter.size(); ++i) {
      if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) {
        ++number_of_full_faces_that_comes_before;
      }
      if (coface_counter[i] != face_counter[i]) {
        if (number_of_position_in_which_counters_do_not_agree != -1) {
          std::cout << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
          throw std::logic_error(
              "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.");
        }
        number_of_position_in_which_counters_do_not_agree = i;
      }
    }

    int incidence = 1;
    if (number_of_full_faces_that_comes_before % 2) incidence = -1;
    // if the face cell is on the right from coface cell:
    if ((coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 ==
         face_counter[number_of_position_in_which_counters_do_not_agree]) ||
        ((coface_counter[number_of_position_in_which_counters_do_not_agree] != 1) &&
         (face_counter[number_of_position_in_which_counters_do_not_agree] == 0))) {
      incidence *= -1;
    }

    return incidence;
  }

 protected:
  std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed;

  void set_up_containers(const std::vector<unsigned>& sizes) {
    unsigned multiplier = 1;
    for (std::size_t i = 0; i != sizes.size(); ++i) {
      this->sizes.push_back(sizes[i]);
      this->multipliers.push_back(multiplier);

      if (directions_in_which_periodic_b_cond_are_to_be_imposed[i]) {
        multiplier *= 2 * sizes[i];
      } else {
        multiplier *= 2 * sizes[i] + 1;
      }
    }
    // std::reverse( this->sizes.begin() , this->sizes.end() );
    this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity());
    this->total_number_of_cells = multiplier;
  }
  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes);
  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
                                                           const std::vector<T>& topDimensionalCells);

  void construct_complex_based_on_top_dimensional_cells(
      const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
      const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
};

template <typename T>
void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells(
    const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
    const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
  this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed;
  this->set_up_containers(dimensions);

  std::size_t i = 0;
  for (auto it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
    this->get_cell_data(*it) = topDimensionalCells[i];
    ++i;
  }
  this->impose_lower_star_filtration();
}

template <typename T>
Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(
    const std::vector<unsigned>& sizes,
    const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
  this->directions_in_which_periodic_b_cond_are_to_be_imposed(directions_in_which_periodic_b_cond_are_to_be_imposed);
  this->set_up_containers(sizes);
}

template <typename T>
Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(
    const char* perseus_style_file) {
  // for Perseus style files:
  bool dbg = false;

  std::ifstream inFiltration;
  inFiltration.open(perseus_style_file);
  unsigned dimensionOfData;
  inFiltration >> dimensionOfData;

  this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData, false);

  std::vector<unsigned> sizes;
  sizes.reserve(dimensionOfData);
  for (std::size_t i = 0; i != dimensionOfData; ++i) {
    int size_in_this_dimension;
    inFiltration >> size_in_this_dimension;
    if (size_in_this_dimension < 0) {
      this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] = true;
    }
    sizes.push_back(abs(size_in_this_dimension));
  }
  this->set_up_containers(sizes);

  typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator it(*this);
  it = this->top_dimensional_cells_iterator_begin();

  while (!inFiltration.eof()) {
    double filtrationLevel;
    inFiltration >> filtrationLevel;
    if (inFiltration.eof()) break;

    if (dbg) {
      std::cerr << "Cell of an index : " << it.compute_index_in_bitmap()
                << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap())
                << " get the value : " << filtrationLevel << std::endl;
    }
    this->get_cell_data(*it) = filtrationLevel;
    ++it;
  }
  inFiltration.close();
  this->impose_lower_star_filtration();
}

template <typename T>
Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(
    const std::vector<unsigned>& sizes) {
  this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(), false);
  this->set_up_containers(sizes);
}

template <typename T>
Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(
    const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells) {
  std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(), false);
  this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells,
                                                         directions_in_which_periodic_b_cond_are_to_be_imposed);
}

template <typename T>
Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(
    const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells,
    const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) {
  this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells,
                                                         directions_in_which_periodic_b_cond_are_to_be_imposed);
}

// ***********************Methods************************ //

template <typename T>
std::vector<std::size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell(
    std::size_t cell) const {
  bool dbg = false;
  if (dbg) {
    std::cerr << "Computations of boundary of a cell : " << cell << std::endl;
  }

  std::vector<std::size_t> boundary_elements;
  boundary_elements.reserve(this->dimension() * 2);
  std::size_t cell1 = cell;
  std::size_t sum_of_dimensions = 0;

  for (std::size_t i = this->multipliers.size(); i != 0; --i) {
    unsigned position = cell1 / this->multipliers[i - 1];
    // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction.
    if (position % 2 == 1) {
      // if there are no periodic boundary conditions in this direction, we do not have to do anything.
      if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
        // std::cerr << "A\n";
        if (sum_of_dimensions % 2) {
          boundary_elements.push_back(cell - this->multipliers[i - 1]);
          boundary_elements.push_back(cell + this->multipliers[i - 1]);
        } else {
          boundary_elements.push_back(cell + this->multipliers[i - 1]);
          boundary_elements.push_back(cell - this->multipliers[i - 1]);
        }
        if (dbg) {
          std::cerr << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
        }
      } else {
        // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end.
        if (position != 2 * this->sizes[i - 1] - 1) {
          // std::cerr << "B\n";
          if (sum_of_dimensions % 2) {
            boundary_elements.push_back(cell - this->multipliers[i - 1]);
            boundary_elements.push_back(cell + this->multipliers[i - 1]);
          } else {
            boundary_elements.push_back(cell + this->multipliers[i - 1]);
            boundary_elements.push_back(cell - this->multipliers[i - 1]);
          }
          if (dbg) {
            std::cerr << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " ";
          }
        } else {
          // std::cerr << "C\n";
          if (sum_of_dimensions % 2) {
            boundary_elements.push_back(cell - this->multipliers[i - 1]);
            boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
          } else {
            boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
            boundary_elements.push_back(cell - this->multipliers[i - 1]);
          }
          if (dbg) {
            std::cerr << cell - this->multipliers[i - 1] << " "
                      << cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1] << " ";
          }
        }
      }
      ++sum_of_dimensions;
    }
    cell1 = cell1 % this->multipliers[i - 1];
  }
  return boundary_elements;
}

template <typename T>
std::vector<std::size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell(
    std::size_t cell) const {
  std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
  std::vector<std::size_t> coboundary_elements;
  std::size_t cell1 = cell;
  for (std::size_t i = this->multipliers.size(); i != 0; --i) {
    unsigned position = cell1 / this->multipliers[i - 1];
    // if the cell has zero length in this direction, then it will have cbd in this direction.
    if (position % 2 == 0) {
      if (!this->directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
        // no periodic boundary conditions in this direction
        if ((counter[i - 1] != 0) && (cell > this->multipliers[i - 1])) {
          coboundary_elements.push_back(cell - this->multipliers[i - 1]);
        }
        if ((counter[i - 1] != 2 * this->sizes[i - 1]) && (cell + this->multipliers[i - 1] < this->data.size())) {
          coboundary_elements.push_back(cell + this->multipliers[i - 1]);
        }
      } else {
        // we want to have periodic boundary conditions in this direction
        if (counter[i - 1] != 0) {
          coboundary_elements.push_back(cell - this->multipliers[i - 1]);
          coboundary_elements.push_back(cell + this->multipliers[i - 1]);
        } else {
          // in this case counter[i-1] == 0.
          coboundary_elements.push_back(cell + this->multipliers[i - 1]);
          coboundary_elements.push_back(cell + (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]);
        }
      }
    }

    cell1 = cell1 % this->multipliers[i - 1];
  }
  return coboundary_elements;
}

}  // namespace cubical_complex

namespace Cubical_complex = cubical_complex;

}  // namespace Gudhi

#endif  // BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_