Merge VR_bitmap

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bitmap@1074 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 4c9ac225ffcfb924371c72fccc44cbd6ecb3d6e3

3 files changed, 1515 insertions, 1669 deletions

diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
index 4dd295e7..1fd36914 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
@@ -20,28 +20,28 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#ifndef BITMAP_CUBICAL_COMPLEX_H_
+#define BITMAP_CUBICAL_COMPLEX_H_
 
-#pragma once
-#include <limits>
-#include <numeric>

-#include "Bitmap_cubical_complex_base.h"
-#include "Bitmap_cubical_complex_periodic_boundary_conditions_base.h"
-
+#include <gudhi/Bitmap_cubical_complex_base.h>
+#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h>
 
+#include <limits>
+#include <utility>  // for pair<>
+#include <algorithm>  // for sort
+#include <vector>
 
-namespace Gudhi
-{
+namespace Gudhi {
 
-namespace Cubical_complex
-{
+namespace Cubical_complex {
 
-//global variable, was used just for debugging.
+// global variable, was used just for debugging.
 const bool globalDbg = false;
 
-template <typename T> class is_before_in_filtration;

-

+template <typename T> class is_before_in_filtration;
+
 /**
-* This is a Bitmap_cubical_complex class. It joints a functionalities of Bitmap_cubical_complex_base and Bitmap_cubical_complex_periodic_boundary_conditions_base classes into

+* This is a Bitmap_cubical_complex class. It joints a functionalities of Bitmap_cubical_complex_base and Bitmap_cubical_complex_periodic_boundary_conditions_base classes into
 * Gudhi persistent homology engine. It is a template class that inherit from its template parameter. The template parameter is supposed to be either Bitmap_cubical_complex_base or Bitmap_cubical_complex_periodic_boundary_conditions_base class.
 **/
 
@@ -51,529 +51,533 @@ template <typename T> class is_before_in_filtration;
  *@ingroup cubical_complex
  */
 template <typename T>
-class Bitmap_cubical_complex : public T
-{
-public:
-//*********************************************//
-//Typedefs and typenames
-//*********************************************//
-    typedef size_t Simplex_key;
-    typedef typename T::filtration_type Filtration_value;
-    typedef Simplex_key Simplex_handle;
-
-
-//*********************************************//
-//Constructors
-//*********************************************//
-    //Over here we need to definie various input types. I am proposing the following ones:
-    //Perseus style
-    //H5 files? TODO
-    //binary files with little endiangs / big endians? TODO
-    //constructor from a vector of elements of a type T. TODO
-
-    /**
-    * Constructor form a Perseus-style file.
-    **/
-    Bitmap_cubical_complex( const char* perseus_style_file ):
-    T(perseus_style_file),key_associated_to_simplex(this->total_number_of_cells+1)
-    {
-        //clock_t begin = clock();
-        if ( globalDbg ){cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n";}
-        for ( size_t i = 0 ; i != this->total_number_of_cells ; ++i )
-        {
-            this->key_associated_to_simplex[i] = i;
-        }
-        //we initialize this only once, in each constructor, when the bitmap is constructed.
-        //If the user decide to change some elements of the bitmap, then this procedure need
-        //to be called again.
-        this->initialize_simplex_associated_to_key();
-        //cerr << "Time of running Bitmap_cubical_complex( const char* perseus_style_file ) constructor : " << double(clock() - begin) / CLOCKS_PER_SEC << endl;
-    }
-
-
-    /**
-    * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells
-    * in the following directions and vector of element of a type T
-    * with filtration on top dimensional cells.
-    **/
-    Bitmap_cubical_complex( const std::vector<unsigned>& dimensions , const std::vector<typename T::filtration_type>& top_dimensional_cells ):
-    T(dimensions,top_dimensional_cells),
-    key_associated_to_simplex(this->total_number_of_cells+1)
-    {
-        for ( size_t i = 0 ; i != this->total_number_of_cells ; ++i )
-        {
-            this->key_associated_to_simplex[i] = i;
-        }
-        //we initialize this only once, in each constructor, when the bitmap is constructed.
-        //If the user decide to change some elements of the bitmap, then this procedure need
-        //to be called again.
-        this->initialize_simplex_associated_to_key();
-    }
-
-    /**
-    * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells
-    * in the following directions and vector of element of a type T::filtration_type
-    * with filtration on top dimensional cells. The last parameter of the constructor is a vector of bools of a length equal to the dimension of cubical complex.
-    * If the position i on this vector is true, then we impose periodic boundary conditions in this direction.
-    **/
-    Bitmap_cubical_complex( const std::vector<unsigned>& dimensions , const std::vector<typename T::filtration_type>& top_dimensional_cells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed  ):
-    T(dimensions,top_dimensional_cells,directions_in_which_periodic_b_cond_are_to_be_imposed),
-    key_associated_to_simplex(this->total_number_of_cells+1)
-    {
-        for ( size_t i = 0 ; i != this->total_number_of_cells ; ++i )
-        {
-            this->key_associated_to_simplex[i] = i;
-        }
-        //we initialize this only once, in each constructor, when the bitmap is constructed.
-        //If the user decide to change some elements of the bitmap, then this procedure need
-        //to be called again.
-        this->initialize_simplex_associated_to_key();
-    }

-

-    /**
-    * Destructor of the Bitmap_cubical_complex class.
-    **/

-    virtual ~Bitmap_cubical_complex(){}
-
-
-//*********************************************//
-//Other 'easy' functions
-//*********************************************//
-    /**
-    * Returns number of all cubes in the complex.
-    **/
-    size_t 	num_simplices()const
-    {
-        return this->total_number_of_cells;
-    }
-
-    /**
-    * Returns a Simplex_handle to a cube that do not exist in this complex.
-    **/
-    static Simplex_handle null_simplex()
-    {
-        if ( globalDbg ){cerr << "Simplex_handle null_simplex()\n";}
-        return std::numeric_limits<Simplex_handle>::max();
-    }
-
-
-    /**
-    * Returns dimension of the complex.
-    **/
-    inline size_t dimension()const
-    {
-        return this->sizes.size();
-    }
-
-    /**
-    * Return dimension of a cell pointed by the Simplex_handle.
-    **/
-    inline unsigned dimension(Simplex_handle sh)const
-    {
-        if ( globalDbg ){cerr << "unsigned dimension(const Simplex_handle& sh)\n";}
-        if ( sh != std::numeric_limits<Simplex_handle>::max() ) return this->get_dimension_of_a_cell( sh );
-        return -1;
-    }
-
-    /**
-    * Return the filtration of a cell pointed by the Simplex_handle.
-    **/
-    typename T::filtration_type filtration(Simplex_handle sh)
-    {
-        if ( globalDbg ){cerr << "T::filtration_type filtration(const Simplex_handle& sh)\n";}
-        //Returns the filtration value of a simplex.
-        if ( sh != std::numeric_limits<Simplex_handle>::max() ) return this->data[sh];
-        return std::numeric_limits<Simplex_handle>::max();
-    }
-
-    /**
-    * Return a key which is not a key of any cube in the considered data structure.
-    **/
-    static Simplex_key null_key()
-    {
-        if ( globalDbg ){cerr << "Simplex_key null_key()\n";}
-        return std::numeric_limits<Simplex_handle>::max();
-    }
-
-    /**
-    * Return the key of a cube pointed by the Simplex_handle.
-    **/
-    Simplex_key key(Simplex_handle sh)const
-    {
-        if ( globalDbg ){cerr << "Simplex_key key(const Simplex_handle& sh)\n";}
-        if ( sh != std::numeric_limits<Simplex_handle>::max() )
-        {
-            return this->key_associated_to_simplex[sh];
-        }
-        return this->null_key();
-    }
-
-    /**
-    * Return the Simplex_handle given the key of the cube.
-    **/
-    Simplex_handle simplex(Simplex_key key)
-    {
-        if ( globalDbg ){cerr << "Simplex_handle simplex(Simplex_key key)\n";}
-        if ( key != std::numeric_limits<Simplex_handle>::max() )
-        {
-            return this->simplex_associated_to_key[ key ];
-        }
-        return null_simplex();
-    }
-
-    /**
-    * Assign key to a cube pointed by the Simplex_handle
-    **/
-    void assign_key(Simplex_handle sh, Simplex_key key)
-    {
-        if ( globalDbg ){cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n";}
-        if ( key == std::numeric_limits<Simplex_handle>::max() ) return;
-        this->key_associated_to_simplex[sh] = key;
-        this->simplex_associated_to_key[key] = sh;
-    }
-
-    /**
-    * Function called from a constructor. It is needed for Filtration_simplex_iterator to work.
-    **/
-    void initialize_simplex_associated_to_key();
-
-
-
-//*********************************************//
-//Iterators
-//*********************************************//
-
-
-
-    /**
-    * Boundary_simplex_range class provides ranges for boundary iterators.
-    **/
-    typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator;
-    typedef typename std::vector< Simplex_handle > Boundary_simplex_range;
-
-
-
-    /**
-    * Filtration_simplex_iterator class provides an iterator though the whole structure in the order of filtration.
-    * Secondary criteria for filtration are:
-    * (1) Dimension of a cube (lower dimensional comes first).
-    * (2) Position in the data structure (the ones that are earlies in the data structure comes first).
-    **/
-    class Filtration_simplex_range;
-    class Filtration_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle >
-    {
-        //Iterator over all simplices of the complex in the order of the indexing scheme.
-        //'value_type' must be 'Simplex_handle'.
-        public:
-            Filtration_simplex_iterator( Bitmap_cubical_complex* b ):b(b),position(0){};
-            Filtration_simplex_iterator():b(NULL){};
-
-            Filtration_simplex_iterator operator++()
-            {
-                if ( globalDbg ){cerr << "Filtration_simplex_iterator operator++\n";}
-                ++this->position;
-                return (*this);
-            }
-            Filtration_simplex_iterator operator++(int)
-            {
-                Filtration_simplex_iterator result = *this;
-                ++(*this);
-                return result;
-            }
-            Filtration_simplex_iterator operator =( const Filtration_simplex_iterator& rhs )
-            {
-                if ( globalDbg ){cerr << "Filtration_simplex_iterator operator =\n";}
-                this->b = rhs.b;
-                this->position = rhs.position;
-            }
-            bool operator == ( const Filtration_simplex_iterator& rhs )const
-            {
-                if ( globalDbg ){cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n";}
-                return ( this->position == rhs.position );
-            }
-
-            bool operator != ( const Filtration_simplex_iterator& rhs )const
-            {
-                if ( globalDbg ){cerr << "bool operator != ( const Filtration_simplex_iterator& rhs )\n";}
-                return !(*this == rhs);
-            }
-            Simplex_handle operator*()
-            {
-                if ( globalDbg ){cerr << "Simplex_handle operator*()\n";}
-                return this->b->simplex_associated_to_key[ this->position ];
-            }
-
-            friend class Filtration_simplex_range;
-        private:
-            Bitmap_cubical_complex<T>* b;
-            size_t position;
-    };
-
-
-    /**
-    * Filtration_simplex_range provides the ranges for Filtration_simplex_iterator.
-    **/
-    class Filtration_simplex_range
-    {
-        //Range over the simplices of the complex in the order of the filtration.
-        //.begin() and .end() return type Filtration_simplex_iterator.
-        public:
-            typedef Filtration_simplex_iterator const_iterator;
-            typedef Filtration_simplex_iterator iterator;
-            Filtration_simplex_range(Bitmap_cubical_complex<T>* b):b(b){};
-            Filtration_simplex_iterator begin()
-            {
-                if ( globalDbg ){cerr << "Filtration_simplex_iterator begin() \n";}
-                return Filtration_simplex_iterator( this->b );
-            }
-            Filtration_simplex_iterator end()
-            {
-                if ( globalDbg ){cerr << "Filtration_simplex_iterator end()\n";}
-                Filtration_simplex_iterator it( this->b );
-                it.position = this->b->simplex_associated_to_key.size();
-                return it;
-            }
-        private:
-            Bitmap_cubical_complex<T>* b;
-    };
-
-
-
-//*********************************************//
-//Methods to access iterators from the container:
-    /**
-    * boundary_simplex_range creates an object of a Boundary_simplex_range class
-    * that provides ranges for the Boundary_simplex_iterator.
-    **/
-    Boundary_simplex_range boundary_simplex_range(Simplex_handle sh)
-    {
-        /*
-        std::vector< size_t > bdry = this->get_boundary_of_a_cell(sh);
-        Boundary_simplex_range result( bdry.size() );
-        for ( size_t i = 0 ; i != bdry.size() ; ++i )
-        {
-            result[i] = this->simplex_associated_to_key[ bdry[i] ];
-        }
-        return result;
-        */
-        return this->get_boundary_of_a_cell(sh);
-    }
-
-    /**
-    * filtration_simplex_range creates an object of a Filtration_simplex_range class
-    * that provides ranges for the Filtration_simplex_iterator.
-    **/
-    Filtration_simplex_range filtration_simplex_range()
-    {
-        if ( globalDbg ){cerr << "Filtration_simplex_range filtration_simplex_range()\n";}
-        //Returns a range over the simplices of the complex in the order of the filtration
-        return Filtration_simplex_range(this);
-    }
-//*********************************************//
-
-
-
-//*********************************************//
-//Elements which are in Gudhi now, but I (and in all the cases I asked also Marc) do not understand why they are there.
-    //TODO -- the file IndexingTag.h in the Gudhi library contains an empty structure, so
-    //I understand that this is something that was planned (for simplicial maps?)
-    //but was never finished. The only idea I have here is to use the same empty structure from
-    //IndexingTag.h file, but only if the compiler needs it. If the compiler
-    //do not need it, then I would rather not add here elements which I do not understand.
-    //typedef Indexing_tag
-    /**
-    * Function needed for compatibility with Gudhi. Not useful for other purposes.
-    **/
-    std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )
-    {
-        std::vector< size_t > bdry = this->get_boundary_of_a_cell( sh    );
-        if ( globalDbg )
-        {
-		cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n";
-              cerr << "bdry.size() : " << bdry.size() << endl;
-	 }
-        //this method returns two first elements from the boundary of sh.
-        if ( bdry.size() < 2 )
-		throw("Error in endpoints in Bitmap_cubical_complex class. The cell have less than two elements in the boundary.");
-        return std::make_pair( bdry[0] , bdry[1] );
-    }
-
-
-    /**
-    * Class needed for compatibility with Gudhi. Not useful for other purposes.
-    **/
-    class Skeleton_simplex_range;
-    class Skeleton_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle >
-    {
-        //Iterator over all simplices of the complex in the order of the indexing scheme.
-        //'value_type' must be 'Simplex_handle'.
-        public:
-            Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , size_t d ):b(b),dimension(d)
-            {
-                if ( globalDbg ){cerr << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , size_t d )\n";}
-                //find the position of the first simplex of a dimension d
-                this->position = 0;
-                while (
-                      (this->position != b->data.size()) &&
-                      ( this->b->get_dimension_of_a_cell( this->position ) != this->dimension )
-                      )
-                {
-                    ++this->position;
-                }
-            };
-            Skeleton_simplex_iterator ():b(NULL),dimension(0){};
-
-            Skeleton_simplex_iterator operator++()
-            {
-                if ( globalDbg ){cerr << "Skeleton_simplex_iterator operator++()\n";}
-                //increment the position as long as you did not get to the next element of the dimension dimension.
-                ++this->position;
-                while (
-			 (this->position != this->b->data.size()) &&
-                      ( this->b->get_dimension_of_a_cell( this->position ) != this->dimension )
-                      )
-                {
-                    ++this->position;
-                }
-                return (*this);
-            }
-            Skeleton_simplex_iterator operator++(int)
-            {
-                Skeleton_simplex_iterator result = *this;
-                ++(*this);
-                return result;
-            }
-            Skeleton_simplex_iterator operator =( const Skeleton_simplex_iterator& rhs )
-            {
-                if ( globalDbg ){cerr << "Skeleton_simplex_iterator operator =\n";}
-                this->b = rhs.b;
-                this->position = rhs.position;
-            }
-            bool operator == ( const Skeleton_simplex_iterator& rhs )const
-            {
-                if ( globalDbg ){cerr << "bool operator ==\n";}
-                return ( this->position == rhs.position );
-            }
-
-            bool operator != ( const Skeleton_simplex_iterator& rhs )const
-            {
-                if ( globalDbg ){cerr << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n";}
-                return !(*this == rhs);
-            }
-            Simplex_handle operator*()
-            {
-                if ( globalDbg ){cerr << "Simplex_handle operator*() \n";}
-                return this->position;
-            }
-
-            friend class Skeleton_simplex_range;
-        private:
-            Bitmap_cubical_complex<T>* b;
-            size_t position;
-            unsigned dimension;
-    };
-    /**
-    * Class needed for compatibility with Gudhi. Not useful for other purposes.
-    **/
-    class Skeleton_simplex_range
-    {
-        //Range over the simplices of the complex in the order of the filtration.
-        //.begin() and .end() return type Filtration_simplex_iterator.
-        public:
-            typedef Skeleton_simplex_iterator const_iterator;
-            typedef Skeleton_simplex_iterator iterator;
-            Skeleton_simplex_range(Bitmap_cubical_complex<T>* b , unsigned dimension):b(b),dimension(dimension){};
-            Skeleton_simplex_iterator begin()
-            {
-                if ( globalDbg ){cerr << "Skeleton_simplex_iterator begin()\n";}
-                return Skeleton_simplex_iterator( this->b , this->dimension );
-            }
-            Skeleton_simplex_iterator end()
-            {
-                if ( globalDbg ){cerr << "Skeleton_simplex_iterator end()\n";}
-                Skeleton_simplex_iterator it( this->b , this->dimension );
-                it.position = this->b->data.size();
-                return it;
-            }
-        private:
-            Bitmap_cubical_complex<T>* b;
-            unsigned dimension;
-    };
-
-    /**
-    * Function needed for compatibility with Gudhi. Not useful for other purposes.
-    **/
-    Skeleton_simplex_range skeleton_simplex_range( unsigned dimension )
-    {
-        if ( globalDbg ){cerr << "Skeleton_simplex_range skeleton_simplex_range( unsigned dimension )\n";}
-        return Skeleton_simplex_range( this , dimension );
-    }
-
-    friend class is_before_in_filtration<T>;
-
-
-protected:
-    std::vector< size_t > key_associated_to_simplex;
-    std::vector< size_t > simplex_associated_to_key;
-};//Bitmap_cubical_complex
+class Bitmap_cubical_complex : public T {
+ public:
+  //*********************************************//
+  // Typedefs and typenames
+  //*********************************************//
+  typedef size_t Simplex_key;
+  typedef typename T::filtration_type Filtration_value;
+  typedef Simplex_key Simplex_handle;
+
+
+  //*********************************************//
+  // Constructors
+  //*********************************************//
+  // Over here we need to define various input types. I am proposing the following ones:
+  // Perseus style
+  // TODO(PD) H5 files?
+  // TODO(PD) binary files with little endiangs / big endians ?
+  // TODO(PD) constructor from a vector of elements of a type T. ?
+
+  /**
+   * Constructor form a Perseus-style file.
+   **/
+  Bitmap_cubical_complex(const char* perseus_style_file) :
+      T(perseus_style_file), key_associated_to_simplex(this->total_number_of_cells + 1) {
+    if (globalDbg) {
+      std::cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n";
+    }
+    for (size_t i = 0; i != this->total_number_of_cells; ++i) {
+      this->key_associated_to_simplex[i] = i;
+    }
+    // we initialize this only once, in each constructor, when the bitmap is constructed.
+    // If the user decide to change some elements of the bitmap, then this procedure need
+    // to be called again.
+    this->initialize_simplex_associated_to_key();
+  }
+
+  /**
+   * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells
+   * in the following directions and vector of element of a type T
+   * with filtration on top dimensional cells.
+   **/
+  Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
+                         const std::vector<typename T::filtration_type>& top_dimensional_cells) :
+      T(dimensions, top_dimensional_cells),
+      key_associated_to_simplex(this->total_number_of_cells + 1) {
+    for (size_t i = 0; i != this->total_number_of_cells; ++i) {
+      this->key_associated_to_simplex[i] = i;
+    }
+    // we initialize this only once, in each constructor, when the bitmap is constructed.
+    // If the user decide to change some elements of the bitmap, then this procedure need
+    // to be called again.
+    this->initialize_simplex_associated_to_key();
+  }
+
+  /**
+   * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells
+   * in the following directions and vector of element of a type T::filtration_type
+   * with filtration on top dimensional cells. The last parameter of the constructor is a vector of bools of a length equal to the dimension of cubical complex.
+   * If the position i on this vector is true, then we impose periodic boundary conditions in this direction.
+   **/
+  Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
+                         const std::vector<typename T::filtration_type>& top_dimensional_cells,
+                         std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed) :
+      T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed),
+      key_associated_to_simplex(this->total_number_of_cells + 1) {
+    for (size_t i = 0; i != this->total_number_of_cells; ++i) {
+      this->key_associated_to_simplex[i] = i;
+    }
+    // we initialize this only once, in each constructor, when the bitmap is constructed.
+    // If the user decide to change some elements of the bitmap, then this procedure need
+    // to be called again.
+    this->initialize_simplex_associated_to_key();
+  }
+
+  /**
+   * Destructor of the Bitmap_cubical_complex class.
+   **/
+  virtual ~Bitmap_cubical_complex(){}
+
+  //*********************************************//
+  // Other 'easy' functions
+  //*********************************************//
+
+  /**
+   * Returns number of all cubes in the complex.
+   **/
+  size_t num_simplices()const {
+    return this->total_number_of_cells;
+  }
+
+  /**
+   * Returns a Simplex_handle to a cube that do not exist in this complex.
+   **/
+  static Simplex_handle null_simplex() {
+    if (globalDbg) {
+      std::cerr << "Simplex_handle null_simplex()\n";
+    }
+    return std::numeric_limits<Simplex_handle>::max();
+  }
+
+  /**
+   * Returns dimension of the complex.
+   **/
+  inline size_t dimension()const {
+    return this->sizes.size();
+  }
+
+  /**
+   * Return dimension of a cell pointed by the Simplex_handle.
+   **/
+  inline unsigned dimension(Simplex_handle sh)const {
+    if (globalDbg) {
+      std::cerr << "unsigned dimension(const Simplex_handle& sh)\n";
+    }
+    if (sh != std::numeric_limits<Simplex_handle>::max()) return this->get_dimension_of_a_cell(sh);
+    return -1;
+  }
+
+  /**
+   * Return the filtration of a cell pointed by the Simplex_handle.
+   **/
+  typename T::filtration_type filtration(Simplex_handle sh) {
+    if (globalDbg) {
+      std::cerr << "T::filtration_type filtration(const Simplex_handle& sh)\n";
+    }
+    // Returns the filtration value of a simplex.
+    if (sh != std::numeric_limits<Simplex_handle>::max()) return this->data[sh];
+    return std::numeric_limits<Simplex_handle>::max();
+  }
+
+  /**
+   * Return a key which is not a key of any cube in the considered data structure.
+   **/
+  static Simplex_key null_key() {
+    if (globalDbg) {
+      std::cerr << "Simplex_key null_key()\n";
+    }
+    return std::numeric_limits<Simplex_handle>::max();
+  }
+
+  /**
+   * Return the key of a cube pointed by the Simplex_handle.
+   **/
+  Simplex_key key(Simplex_handle sh)const {
+    if (globalDbg) {
+      std::cerr << "Simplex_key key(const Simplex_handle& sh)\n";
+    }
+    if (sh != std::numeric_limits<Simplex_handle>::max()) {
+      return this->key_associated_to_simplex[sh];
+    }
+    return this->null_key();
+  }
+
+  /**
+   * Return the Simplex_handle given the key of the cube.
+   **/
+  Simplex_handle simplex(Simplex_key key) {
+    if (globalDbg) {
+      std::cerr << "Simplex_handle simplex(Simplex_key key)\n";
+    }
+    if (key != std::numeric_limits<Simplex_handle>::max()) {
+      return this->simplex_associated_to_key[ key ];
+    }
+    return null_simplex();
+  }
+
+  /**
+   * Assign key to a cube pointed by the Simplex_handle
+   **/
+  void assign_key(Simplex_handle sh, Simplex_key key) {
+    if (globalDbg) {
+      std::cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n";
+    }
+    if (key == std::numeric_limits<Simplex_handle>::max()) return;
+    this->key_associated_to_simplex[sh] = key;
+    this->simplex_associated_to_key[key] = sh;
+  }
+
+  /**
+   * Function called from a constructor. It is needed for Filtration_simplex_iterator to work.
+   **/
+  void initialize_simplex_associated_to_key();
+
+  //*********************************************//
+  // Iterators
+  //*********************************************//
+
+  /**
+   * Boundary_simplex_range class provides ranges for boundary iterators.
+   **/
+  typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator;
+  typedef typename std::vector< Simplex_handle > Boundary_simplex_range;
+
+  /**
+   * Filtration_simplex_iterator class provides an iterator though the whole structure in the order of filtration.
+   * Secondary criteria for filtration are:
+   * (1) Dimension of a cube (lower dimensional comes first).
+   * (2) Position in the data structure (the ones that are earlies in the data structure comes first).
+   **/
+  class Filtration_simplex_range;
+
+  class Filtration_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > {
+    // Iterator over all simplices of the complex in the order of the indexing scheme.
+    // 'value_type' must be 'Simplex_handle'.
+   public:
+    Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) { }
+
+    Filtration_simplex_iterator() : b(NULL) { }
+
+    Filtration_simplex_iterator operator++() {
+      if (globalDbg) {
+        std::cerr << "Filtration_simplex_iterator operator++\n";
+      }
+      ++this->position;
+      return (*this);
+    }
 
-template <typename T>
-void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key()
-{
-    if ( globalDbg )
-    {
-	cerr << "void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration() \n";
-    }
-    this->simplex_associated_to_key = std::vector<size_t>( this->data.size() );
-    std::iota (std::begin(simplex_associated_to_key), std::end(simplex_associated_to_key), 0);
-    std::sort( simplex_associated_to_key.begin() ,
-               simplex_associated_to_key.end() ,
-               is_before_in_filtration<T>(this) );
-
-    //we still need to deal here with a key_associated_to_simplex:
-    for ( size_t i = 0  ; i != simplex_associated_to_key.size() ; ++i )
-    {
-        this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i;
+    Filtration_simplex_iterator operator++(int) {
+      Filtration_simplex_iterator result = *this;
+      ++(*this);
+      return result;
     }
-}
 
+    Filtration_simplex_iterator& operator=(const Filtration_simplex_iterator& rhs) {
+      if (globalDbg) {
+        std::cerr << "Filtration_simplex_iterator operator =\n";
+      }
+      this->b = rhs.b;
+      this->position = rhs.position;
+    }
 
-template <typename T>
-class is_before_in_filtration
-{
-public:
-    explicit is_before_in_filtration(Bitmap_cubical_complex<T> * CC)
-        : CC_(CC) { }
-
-    bool operator()( const typename Bitmap_cubical_complex<T>::Simplex_handle sh1, const typename Bitmap_cubical_complex<T>::Simplex_handle sh2) const
-    {
-      // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
-      typename T::filtration_type fil1 = CC_->data[sh1];
-      typename T::filtration_type fil2 = CC_->data[sh2];
-      if ( fil1 != fil2 )
-      {
-	   return fil1 < fil2;
+    bool operator==(const Filtration_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n";
+      }
+      return ( this->position == rhs.position);
+    }
+
+    bool operator!=(const Filtration_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator != ( const Filtration_simplex_iterator& rhs )\n";
+      }
+      return !(*this == rhs);
+    }
+
+    Simplex_handle operator*() {
+      if (globalDbg) {
+        std::cerr << "Simplex_handle operator*()\n";
+      }
+      return this->b->simplex_associated_to_key[ this->position ];
+    }
+
+    friend class Filtration_simplex_range;
+
+   private:
+    Bitmap_cubical_complex<T>* b;
+    size_t position;
+  };
+
+  /**
+   * Filtration_simplex_range provides the ranges for Filtration_simplex_iterator.
+   **/
+  class Filtration_simplex_range {
+    // Range over the simplices of the complex in the order of the filtration.
+    // .begin() and .end() return type Filtration_simplex_iterator.
+   public:
+    typedef Filtration_simplex_iterator const_iterator;
+    typedef Filtration_simplex_iterator iterator;
+
+    Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) { }
+
+    Filtration_simplex_iterator begin() {
+      if (globalDbg) {
+        std::cerr << "Filtration_simplex_iterator begin() \n";
       }
-      //in this case they are on the same filtration level, so the dimension decide.
-      size_t dim1 = CC_->get_dimension_of_a_cell(sh1);
-      size_t dim2 = CC_->get_dimension_of_a_cell(sh2);
-      if ( dim1 != dim2 )
-      {
-	      return dim1 < dim2;
+      return Filtration_simplex_iterator(this->b);
+    }
+
+    Filtration_simplex_iterator end() {
+      if (globalDbg) {
+        std::cerr << "Filtration_simplex_iterator end()\n";
       }
-      //in this case both filtration and dimensions of the considered cubes are the same. To have stable sort, we simply compare their positions in the bitmap:
-      return sh1 < sh2;
-     }
-protected:
-    Bitmap_cubical_complex<T>* CC_;
+      Filtration_simplex_iterator it(this->b);
+      it.position = this->b->simplex_associated_to_key.size();
+      return it;
+    }
+
+   private:
+    Bitmap_cubical_complex<T>* b;
   };
 
 
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-//****************************************************************************************************************//
 
+  //*********************************************//
+  // Methods to access iterators from the container:
 
+  /**
+   * boundary_simplex_range creates an object of a Boundary_simplex_range class
+   * that provides ranges for the Boundary_simplex_iterator.
+   **/
+  Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) {
+    return this->get_boundary_of_a_cell(sh);
+  }
+
+  /**
+   * filtration_simplex_range creates an object of a Filtration_simplex_range class
+   * that provides ranges for the Filtration_simplex_iterator.
+   **/
+  Filtration_simplex_range filtration_simplex_range() {
+    if (globalDbg) {
+      std::cerr << "Filtration_simplex_range filtration_simplex_range()\n";
+    }
+    // Returns a range over the simplices of the complex in the order of the filtration
+    return Filtration_simplex_range(this);
+  }
+  //*********************************************//
+
+
+
+  //*********************************************//
+  // Elements which are in Gudhi now, but I (and in all the cases I asked also Marc) do not understand why they are
+  // there.
+  // TODO(PD) the file IndexingTag.h in the Gudhi library contains an empty structure, so
+  // I understand that this is something that was planned (for simplicial maps?)
+  // but was never finished. The only idea I have here is to use the same empty structure from
+  // IndexingTag.h file, but only if the compiler needs it. If the compiler
+  // do not need it, then I would rather not add here elements which I do not understand.
+  // typedef Indexing_tag
+
+  /**
+   * Function needed for compatibility with Gudhi. Not useful for other purposes.
+   **/
+  std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) {
+    std::vector< size_t > bdry = this->get_boundary_of_a_cell(sh);
+    if (globalDbg) {
+      std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n";
+      std::cerr << "bdry.size() : " << bdry.size() << std::endl;
+    }
+    // this method returns two first elements from the boundary of sh.
+    if (bdry.size() < 2)
+      throw("Error in endpoints in Bitmap_cubical_complex class. The cell have less than two elements in the boundary.");
+    return std::make_pair(bdry[0], bdry[1]);
+  }
+
+
+  /**
+   * Class needed for compatibility with Gudhi. Not useful for other purposes.
+   **/
+  class Skeleton_simplex_range;
+
+  class Skeleton_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > {
+    // Iterator over all simplices of the complex in the order of the indexing scheme.
+    // 'value_type' must be 'Simplex_handle'.
+   public:
+    Skeleton_simplex_iterator(Bitmap_cubical_complex* b, size_t d) : b(b), dimension(d) {
+      if (globalDbg) {
+        std::cerr << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , size_t d )\n";
+      }
+      // find the position of the first simplex of a dimension d
+      this->position = 0;
+      while (
+             (this->position != b->data.size()) &&
+             (this->b->get_dimension_of_a_cell(this->position) != this->dimension)
+             ) {
+        ++this->position;
+      }
+    }
+
+    Skeleton_simplex_iterator() : b(NULL), position(0), dimension(0) { }
+
+    Skeleton_simplex_iterator operator++() {
+      if (globalDbg) {
+        std::cerr << "Skeleton_simplex_iterator operator++()\n";
+      }
+      // increment the position as long as you did not get to the next element of the dimension dimension.
+      ++this->position;
+      while (
+             (this->position != this->b->data.size()) &&
+             (this->b->get_dimension_of_a_cell(this->position) != this->dimension)
+             ) {
+        ++this->position;
+      }
+      return (*this);
+    }
+
+    Skeleton_simplex_iterator operator++(int) {
+      Skeleton_simplex_iterator result = *this;
+      ++(*this);
+      return result;
+    }
+
+    Skeleton_simplex_iterator& operator=(const Skeleton_simplex_iterator& rhs) {
+      if (globalDbg) {
+        std::cerr << "Skeleton_simplex_iterator operator =\n";
+      }
+      this->b = rhs.b;
+      this->position = rhs.position;
+      this->dimension = rhs.dimension;
+    }
+
+    bool operator==(const Skeleton_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator ==\n";
+      }
+      return ( this->position == rhs.position);
+    }
+
+    bool operator!=(const Skeleton_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n";
+      }
+      return !(*this == rhs);
+    }
+
+    Simplex_handle operator*() {
+      if (globalDbg) {
+        std::cerr << "Simplex_handle operator*() \n";
+      }
+      return this->position;
+    }
+
+    friend class Skeleton_simplex_range;
+   private:
+    Bitmap_cubical_complex<T>* b;
+    size_t position;
+    unsigned dimension;
+  };
+
+  /**
+   * Class needed for compatibility with Gudhi. Not useful for other purposes.
+   **/
+  class Skeleton_simplex_range {
+    // Range over the simplices of the complex in the order of the filtration.
+    // .begin() and .end() return type Filtration_simplex_iterator.
+   public:
+    typedef Skeleton_simplex_iterator const_iterator;
+    typedef Skeleton_simplex_iterator iterator;
+
+    Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, unsigned dimension) : b(b), dimension(dimension) { }
+
+    Skeleton_simplex_iterator begin() {
+      if (globalDbg) {
+        std::cerr << "Skeleton_simplex_iterator begin()\n";
+      }
+      return Skeleton_simplex_iterator(this->b, this->dimension);
+    }
+
+    Skeleton_simplex_iterator end() {
+      if (globalDbg) {
+        std::cerr << "Skeleton_simplex_iterator end()\n";
+      }
+      Skeleton_simplex_iterator it(this->b, this->dimension);
+      it.position = this->b->data.size();
+      return it;
+    }
+
+   private:
+    Bitmap_cubical_complex<T>* b;
+    unsigned dimension;
+  };
+
+  /**
+   * Function needed for compatibility with Gudhi. Not useful for other purposes.
+   **/
+  Skeleton_simplex_range skeleton_simplex_range(unsigned dimension) {
+    if (globalDbg) {
+      std::cerr << "Skeleton_simplex_range skeleton_simplex_range( unsigned dimension )\n";
+    }
+    return Skeleton_simplex_range(this, dimension);
+  }
+
+  friend class is_before_in_filtration<T>;
+
+ protected:
+  std::vector< size_t > key_associated_to_simplex;
+  std::vector< size_t > simplex_associated_to_key;
+};  // Bitmap_cubical_complex
+
+template <typename T>
+void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key() {
+  if (globalDbg) {
+    std::cerr << "void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration() \n";
+  }
+  this->simplex_associated_to_key = std::vector<size_t>(this->data.size());
+  std::iota(std::begin(simplex_associated_to_key), std::end(simplex_associated_to_key), 0);
+  std::sort(simplex_associated_to_key.begin(),
+            simplex_associated_to_key.end(),
+            is_before_in_filtration<T>(this));
+
+  // we still need to deal here with a key_associated_to_simplex:
+  for ( size_t i = 0  ; i != simplex_associated_to_key.size() ; ++i )
+  {
+    this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i;
+  }
 }
-}
+
+template <typename T>
+class is_before_in_filtration {
+ public:
+  explicit is_before_in_filtration(Bitmap_cubical_complex<T> * CC)
+      : CC_(CC) { }
+
+  bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1,
+                  const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const {
+    // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
+    typename T::filtration_type fil1 = CC_->data[sh1];
+    typename T::filtration_type fil2 = CC_->data[sh2];
+    if (fil1 != fil2) {
+      return fil1 < fil2;
+    }
+    // in this case they are on the same filtration level, so the dimension decide.
+    size_t dim1 = CC_->get_dimension_of_a_cell(sh1);
+    size_t dim2 = CC_->get_dimension_of_a_cell(sh2);
+    if (dim1 != dim2) {
+      return dim1 < dim2;
+    }
+    // in this case both filtration and dimensions of the considered cubes are the same. To have stable sort, we simply
+    // compare their positions in the bitmap:
+    return sh1 < sh2;
+  }
+
+ protected:
+  Bitmap_cubical_complex<T>* CC_;
+};
+
+}  // namespace Cubical_complex
+
+}  // namespace Gudhi
+
+#endif  // BITMAP_CUBICAL_COMPLEX_H_
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
index f0517a86..62776019 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -20,35 +20,30 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#pragma once
+#ifndef BITMAP_CUBICAL_COMPLEX_BASE_H_
+#define BITMAP_CUBICAL_COMPLEX_BASE_H_
 
-#include <iostream>

-#include <numeric>
-#include <string>
+#include <gudhi/Bitmap_cubical_complex/counter.h>
+
+#include <iostream>
 #include <vector>
 #include <string>
 #include <fstream>
 #include <algorithm>
 #include <iterator>
-#include <limits>

-#include <ctime>
-#include "Bitmap_cubical_complex/counter.h"
-
-
-using namespace std;
+#include <limits>
+#include <utility>  // for pair<>
 
-namespace Gudhi
-{
-
-namespace Cubical_complex
-{
+namespace Gudhi {
 
+namespace Cubical_complex {
 
 /**
- *@class Bitmap_cubical_complex_base
- *@brief Cubical complex represented as a bitmap, class with basic implementation.
- *@ingroup cubical_complex
+ * @class Bitmap_cubical_complex_base
+ * @brief Cubical complex represented as a bitmap, class with basic implementation.
+ * @ingroup cubical_complex
  */
+
 /**
  * This is a class implementing a basic bitmap data structure to store cubical complexes.
  * It implements only the most basic subroutines.
@@ -67,844 +62,750 @@ namespace Cubical_complex
  * The default filtration used in this implementation is the lower star filtration.
  */
 template <typename T>
-class Bitmap_cubical_complex_base
-{
-public:

-    typedef T filtration_type;

-    /**

-    *Default constructor

-    **/

-    Bitmap_cubical_complex_base()

-    {

-    }
-    /**
-    * There are a few constructors of a Bitmap_cubical_complex_base class.
-    * First one, that takes vector<unsigned>, creates an empty bitmap of a dimension equal
-    * the number of elements in the
-    * input vector and size in the i-th dimension equal the number in the position i-of the input vector.
-    */
-    Bitmap_cubical_complex_base( const std::vector<unsigned>& sizes );
-    /**
-    * The second constructor takes as a input a Perseus style file. For more details,
-    * please consult the documentations of
-    * Perseus software as well as examples attached to this
-    * implementation.
-    **/
-    Bitmap_cubical_complex_base( const char* perseus_style_file );
-    /**
-    * The last constructor of a Bitmap_cubical_complex_base class accepts vector of dimensions (as the first one)
-    * together with vector of filtration values of top dimensional cells.
-    **/
-    Bitmap_cubical_complex_base( const std::vector<unsigned>& dimensions , const std::vector<T>& top_dimensional_cells );
-

-    /**
-    * Destructor of the Bitmap_cubical_complex_base class.
-    **/

-    virtual ~Bitmap_cubical_complex_base(){}

-
-    /**
-    * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell
-    * and get_cell_data are the basic
-     * functions that compute boundary / coboundary / dimension and the filtration
-     * value form a position of a cell in the structure of a bitmap. The input parameter of all of those function is a
-     * non-negative integer, indicating a position of a cube in the data structure.
-     * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to
-     * the positions of (co)boundary element of the input cell.
-    */
-    virtual inline std::vector< size_t > get_boundary_of_a_cell( size_t cell )const;
-    /**
-    * The functions get_coboundary_of_a_cell, get_coboundary_of_a_cell,
-    * get_dimension_of_a_cell and get_cell_data are the basic
-    * functions that compute boundary / coboundary / dimension and the filtration
-    * value form a position of a cell in the structure of a bitmap.
-    * The input parameter of all of those function is a non-negative integer,
-    * indicating a position of a cube in the data structure.
-    * In the case of functions that compute (co)boundary, the output is a vector if
-    * non-negative integers pointing to the
-    * positions of (co)boundary element of the input cell.
-    **/
-    virtual inline std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const;
-    /**
-    * In the case of get_dimension_of_a_cell function, the output is a non-negative integer
-    * indicating the dimension of a cell.
-    **/
-    inline unsigned get_dimension_of_a_cell( size_t cell )const;
-    /**
-    * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position. This allows reading and changing the value of filtration.

-    * Note that if the value of a filtration is changed, the code do not check if we have a filtration or not. i.e. it do not check if the value of a filtration of a cell is

-    * not smaller than the value of a filtration of its boundary and not greater than the value of its coboundary.
-    **/
-    inline T& get_cell_data( size_t cell );
-
-
-    /**
-    * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells.
-    * Then, there are a few ways one can pick the filtration of lower dimensional
-    * cells. The most typical one is by so called lower star filtration. This function is always called by any
-    * constructor which takes the top dimensional cells. If you use such a constructor,
-    * then there is no need to call this function. Call it only if you are putting the filtration
-    * of the cells by your own (for instance by using Top_dimensional_cells_iterator).
-    **/
-    void impose_lower_star_filtration();//assume that top dimensional cells are already set.
-
-    /**
-    * Returns dimension of a complex.
-    **/
-    inline unsigned dimension()const{ return sizes.size(); }
-
-    /**
-    * Returns number of all cubes in the data structure.
-    **/
-    inline unsigned size()const
-    {
-        return this->data.size();
-    }
-
-    /**
-    * Writing to stream operator. By using it we get the values T of cells in order in which they are stored in the structure. This procedure is used for debugging purposes.
-    **/
-    template <typename K>
-    friend ostream& operator << ( ostream & os , const Bitmap_cubical_complex_base<K>& b );

-

-

-    /**
-    * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values equally distributed in the range of data.

-    * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute persistence gets

-    * worst. When dealing with this type of data, one may want to put different values on cells to some number of bins. The function put_data_to_bins( size_t number_of_bins )

-    * ais designed for that purpose. The parameter of the function is the number of bins (distinct values) we want to have in the cubical complex.
-    **/
-    void put_data_to_bins( size_t number_of_bins );

-

-    /**
-    * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values equally distributed in the range of data.

-    * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute persistence gets

-    * worst. When dealing with this type of data, one may want to put different values on cells to some number of bins. The function put_data_to_bins( T diameter_of_bin ) is

-    * designed for that purpose. The parameter of it is the diameter of each bin. Note that the bottleneck distance between the persistence diagram of the cubical complex

-    * before and after using such a function will be bounded by the parameter diameter_of_bin.
-    **/
-    void put_data_to_bins( T diameter_of_bin );

-

-    /**
-    * Functions to find min and max values of filtration.
-    **/

-    std::pair< T ,T > min_max_filtration();
-
-    //ITERATORS

-    /**
-    * Iterator through all cells in the complex (in order they appear in the structure -- i.e.
-    * in lexicographical order).
-    **/

-    class All_cells_iterator : std::iterator< std::input_iterator_tag, T >
-    {
-        public:
-            All_cells_iterator()
-            {
-               this->counter = 0;
-             }
-            All_cells_iterator operator++()
-            {
-                //first find first element of the counter that can be increased:
-               ++this->counter;
-                return *this;
-            }
-            All_cells_iterator operator++(int)
-            {
-                All_cells_iterator result = *this;
-                ++(*this);
-                return result;
-            }
-            All_cells_iterator operator =( const All_cells_iterator& rhs )
-            {
-                this->counter = rhs.counter;
-                return *this;
-            }
-            bool operator == ( const All_cells_iterator& rhs )const
-            {
-                if ( this->counter != rhs.counter )return false;
-                return true;
-            }
-            bool operator != ( const All_cells_iterator& rhs )const
-            {
-                return !(*this == rhs);
-            }
-

-            /*

-            * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as an argument of the following functions:

-            * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell boundary and coboundary and dimension

-            * and in function get_cell_data to get a filtration of a cell.

-            */
-            size_t operator*()
-            {
-                return this->counter;
-            }
-            friend class Bitmap_cubical_complex_base;
-        protected:
-            size_t counter;
-    };

-

-    /**
-    * Function returning a All_cells_iterator to the first cell of the bitmap.
-    **/
-    All_cells_iterator all_cells_iterator_begin()
-    {
-        All_cells_iterator a;
-        return a;
-    }

-

-    /**
-    * Function returning a All_cells_iterator to the last cell of the bitmap.
-    **/
-    All_cells_iterator all_cells_iterator_end()
-    {
-        All_cells_iterator a;
-        a.counter = this->data.size();
-        return a;
-    }

-

-

-    /**
-    * All_cells_range class provides ranges for All_cells_iterator
-    **/

-    class All_cells_range
-    {
-        public:
-            All_cells_range(Bitmap_cubical_complex_base* b):b(b){};
-            All_cells_iterator begin()
-            {
-                return b->all_cells_iterator_begin();
-            }
-            All_cells_iterator end()
-            {
-                return b->all_cells_iterator_end();
-            }
-        private:
-            Bitmap_cubical_complex_base<T>* b;
-    };

-

-    All_cells_range all_cells_range()

-    {

-        return All_cells_range(this);

-    }
-

-

-    /**
-    * Boundary_range class provides ranges for boundary iterators.
-    **/
-    typedef typename std::vector< size_t >::const_iterator Boundary_iterator;
-    typedef typename std::vector< size_t > Boundary_range;

-

-    /**
-    * boundary_simplex_range creates an object of a Boundary_simplex_range class
-    * that provides ranges for the Boundary_simplex_iterator.
-    **/
-    Boundary_range boundary_range(size_t sh)
-    {
-        return this->get_boundary_of_a_cell(sh);
-    }

-

-    /**
-    * Coboundary_range class provides ranges for boundary iterators.
-    **/
-    typedef typename std::vector< size_t >::const_iterator Coboundary_iterator;
-    typedef typename std::vector< size_t > Coboundary_range;

-

-    /**
-    * boundary_simplex_range creates an object of a Boundary_simplex_range class
-    * that provides ranges for the Boundary_simplex_iterator.
-    **/
-    Coboundary_range coboundary_range(size_t sh)
-    {
-        return this->get_coboundary_of_a_cell(sh);
-    }

-

-

-

-

-

-

-
-
-    /**
-    * Iterator through top dimensional cells of the complex. The cells appear in order they are stored
-    * in the structure (i.e. in lexicographical order)
-    **/
-    class Top_dimensional_cells_iterator : std::iterator< std::input_iterator_tag, T >
-    {
-        public:
-            Top_dimensional_cells_iterator( Bitmap_cubical_complex_base& b ):b(b)
-            {
-               this->counter = std::vector<size_t>(b.dimension());
-               //std::fill( this->counter.begin() , this->counter.end() , 0 );
-             }
-            Top_dimensional_cells_iterator operator++()
-            {
-                //first find first element of the counter that can be increased:
-                size_t dim = 0;
-                while ( ( dim != this->b.dimension() ) && ( this->counter[dim] == this->b.sizes[dim]-1 ) )++dim;
-
-                if ( dim != this->b.dimension() )
-                {
-                    ++this->counter[dim];
-                    for ( size_t i = 0 ; i != dim ; ++i )
-                    {
-                        this->counter[i] = 0;
-                    }
-                }
-                else
-                {
-                    ++this->counter[0];
-                }
-                return *this;
-            }
-            Top_dimensional_cells_iterator operator++(int)
-            {
-                Top_dimensional_cells_iterator result = *this;
-                ++(*this);
-                return result;
-            }
-            Top_dimensional_cells_iterator operator =( const Top_dimensional_cells_iterator& rhs )
-            {
-                this->counter = rhs.counter;
-                this->b = rhs.b;
-                return *this;
-            }
-            bool operator == ( const Top_dimensional_cells_iterator& rhs )const
-            {
-                if ( &this->b != &rhs.b )return false;
-                if ( this->counter.size() != rhs.counter.size() )return false;
-                for ( size_t i = 0 ; i != this->counter.size() ; ++i )
-                {
-                    if ( this->counter[i] != rhs.counter[i] )return false;
-                }
-                return true;
-            }
-            bool operator != ( const Top_dimensional_cells_iterator& rhs )const
-            {
-                return !(*this == rhs);
-            }
-
-            //T& operator*()
-            //{
-            //    //given the counter, compute the index in the array and return this element.
-            //    unsigned index = 0;
-            //    for ( size_t i = 0 ; i != this->counter.size() ; ++i )
-            //    {
-            //        index += (2*this->counter[i]+1)*this->b.multipliers[i];
-            //    }
-            //    return this->b.data[index];
-            //}

-

-            /*

-            * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as an argument of the following functions:

-            * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell boundary and coboundary and dimension

-            * and in function get_cell_data to get a filtration of a cell.

-            */

-            size_t operator*()
-            {

-                return this->compute_index_in_bitmap();

-            }
-
-            size_t compute_index_in_bitmap()const
-            {
-                size_t index = 0;
-                for ( size_t i = 0 ; i != this->counter.size() ; ++i )
-                {
-                    index += (2*this->counter[i]+1)*this->b.multipliers[i];
-                }
-                return index;
-            }

-
-
-            void print_counter()const
-            {
-                for ( size_t i = 0 ; i != this->counter.size() ; ++i )
-                {
-                    cout << this->counter[i] << " ";
-                }
-            }
-            friend class Bitmap_cubical_complex_base;
-        protected:
-            std::vector< size_t > counter;
-            Bitmap_cubical_complex_base& b;
-    };

-

-    /**
-    * Function returning a Top_dimensional_cells_iterator to the first top dimensional cell of the bitmap.
-    **/
-    Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin()
-    {
-        Top_dimensional_cells_iterator a(*this);
-        return a;
-    }

-

-    /**
-    * Function returning a Top_dimensional_cells_iterator to the last top dimensional cell of the bitmap.
-    **/
-    Top_dimensional_cells_iterator top_dimensional_cells_iterator_end()
-    {
-        Top_dimensional_cells_iterator a(*this);
-        for ( size_t i = 0 ; i != this->dimension() ; ++i )
-        {
-            a.counter[i] = this->sizes[i]-1;
-        }
-        a.counter[0]++;
-        return a;
-    }

-

-    /**
-    * Top_dimensional_cells_iterator_range class provides ranges for Top_dimensional_cells_iterator_range
-    **/

-    class Top_dimensional_cells_range
-    {
-        public:
-            Top_dimensional_cells_range(Bitmap_cubical_complex_base* b):b(b){};
-            Top_dimensional_cells_iterator begin()
-            {
-                return b->top_dimensional_cells_iterator_begin();
-            }
-            Top_dimensional_cells_iterator end()
-            {
-                return b->top_dimensional_cells_iterator_end();
-            }
-        private:
-            Bitmap_cubical_complex_base<T>* b;
-    };

-

-    Top_dimensional_cells_range top_dimensional_cells_range()

-    {

-        return Top_dimensional_cells_range(this);

-    }
-
-
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-

-

-inline size_t number_cells()const

-{

-    return this->total_number_of_cells;

-}
+class Bitmap_cubical_complex_base {
+ public:
+  typedef T filtration_type;
+
+  /**
+   *Default constructor
+   **/
+  Bitmap_cubical_complex_base() :
+      total_number_of_cells(0) { }
+  /**
+   * There are a few constructors of a Bitmap_cubical_complex_base class.
+   * First one, that takes vector<unsigned>, creates an empty bitmap of a dimension equal
+   * the number of elements in the
+   * input vector and size in the i-th dimension equal the number in the position i-of the input vector.
+   */
+  Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes);
+  /**
+   * The second constructor takes as a input a Perseus style file. For more details,
+   * please consult the documentations of
+   * Perseus software as well as examples attached to this
+   * implementation.
+   **/
+  Bitmap_cubical_complex_base(const char* perseus_style_file);
+  /**
+   * The last constructor of a Bitmap_cubical_complex_base class accepts vector of dimensions (as the first one)
+   * together with vector of filtration values of top dimensional cells.
+   **/
+  Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells);
+
+  /**
+   * Destructor of the Bitmap_cubical_complex_base class.
+   **/
+  virtual ~Bitmap_cubical_complex_base() { }
+
+  /**
+   * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell
+   * and get_cell_data are the basic
+   * functions that compute boundary / coboundary / dimension and the filtration
+   * value form a position of a cell in the structure of a bitmap. The input parameter of all of those function is a
+   * non-negative integer, indicating a position of a cube in the data structure.
+   * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to
+   * the positions of (co)boundary element of the input cell.
+   */
+  virtual inline std::vector< size_t > get_boundary_of_a_cell(size_t cell)const;
+  /**
+   * The functions get_coboundary_of_a_cell, get_coboundary_of_a_cell,
+   * get_dimension_of_a_cell and get_cell_data are the basic
+   * functions that compute boundary / coboundary / dimension and the filtration
+   * value form a position of a cell in the structure of a bitmap.
+   * The input parameter of all of those function is a non-negative integer,
+   * indicating a position of a cube in the data structure.
+   * In the case of functions that compute (co)boundary, the output is a vector if
+   * non-negative integers pointing to the
+   * positions of (co)boundary element of the input cell.
+   **/
+  virtual inline std::vector< size_t > get_coboundary_of_a_cell(size_t cell)const;
+  /**
+   * In the case of get_dimension_of_a_cell function, the output is a non-negative integer
+   * indicating the dimension of a cell.
+   **/
+  inline unsigned get_dimension_of_a_cell(size_t cell)const;
+  /**
+   * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position.
+   * This allows reading and changing the value of filtration. Note that if the value of a filtration is changed, the
+   * code do not check if we have a filtration or not. i.e. it do not check if the value of a filtration of a cell is
+   * not smaller than the value of a filtration of its boundary and not greater than the value of its coboundary.
+   **/
+  inline T& get_cell_data(size_t cell);
+
+
+  /**
+   * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells.
+   * Then, there are a few ways one can pick the filtration of lower dimensional
+   * cells. The most typical one is by so called lower star filtration. This function is always called by any
+   * constructor which takes the top dimensional cells. If you use such a constructor,
+   * then there is no need to call this function. Call it only if you are putting the filtration
+   * of the cells by your own (for instance by using Top_dimensional_cells_iterator).
+   **/
+  void impose_lower_star_filtration(); // assume that top dimensional cells are already set.
+
+  /**
+   * Returns dimension of a complex.
+   **/
+  inline unsigned dimension()const {
+    return sizes.size();
+  }
+
+  /**
+   * Returns number of all cubes in the data structure.
+   **/
+  inline unsigned size()const {
+    return this->data.size();
+  }
+
+  /**
+   * Writing to stream operator. By using it we get the values T of cells in order in which they are stored in the
+   * structure. This procedure is used for debugging purposes.
+   **/
+  template <typename K>
+  friend std::ostream& operator<<(std::ostream & os, const Bitmap_cubical_complex_base<K>& b);
+
+  /**
+   * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values
+   * equally distributed in the range of data.
+   * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute
+   * persistence gets worst. When dealing with this type of data, one may want to put different values on cells to
+   * some number of bins. The function put_data_to_bins( size_t number_of_bins ) is designed for that purpose.
+   * The parameter of the function is the number of bins (distinct values) we want to have in the cubical complex.
+   **/
+  void put_data_to_bins(size_t number_of_bins);
+
+  /**
+   * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values
+   * equally distributed in the range of data.
+   * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute
+   * persistence gets worst. When dealing with this type of data, one may want to put different values on cells to
+   * some number of bins. The function put_data_to_bins( T diameter_of_bin ) is designed for that purpose.
+   * The parameter of it is the diameter of each bin. Note that the bottleneck distance between the persistence
+   * diagram of the cubical complex before and after using such a function will be bounded by the parameter
+   * diameter_of_bin.
+   **/
+  void put_data_to_bins(T diameter_of_bin);
+
+  /**
+   * Functions to find min and max values of filtration.
+   **/
+  std::pair< T, T > min_max_filtration();
+
+  // ITERATORS
+
+  /**
+   * Iterator through all cells in the complex (in order they appear in the structure -- i.e.
+   * in lexicographical order).
+   **/
+  class All_cells_iterator : std::iterator< std::input_iterator_tag, T > {
+   public:
+
+    All_cells_iterator() {
+      this->counter = 0;
+    }
 
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-//****************************************************************************************************************//
-
-protected:
-    std::vector<unsigned> sizes;
-    std::vector<unsigned> multipliers;
-    std::vector<T> data;
-    size_t total_number_of_cells;
-    void set_up_containers( const std::vector<unsigned>& sizes )
-    {
-        unsigned multiplier = 1;
-        for ( size_t i = 0 ; i != sizes.size() ; ++i )
-        {
-            this->sizes.push_back(sizes[i]);
-            this->multipliers.push_back(multiplier);
-            multiplier *= 2*sizes[i]+1;
-        }
-        this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max());
-        this->total_number_of_cells = multiplier;
+    All_cells_iterator operator++() {
+      //first find first element of the counter that can be increased:
+      ++this->counter;
+      return *this;
     }
 
-    size_t compute_position_in_bitmap( const std::vector< unsigned >& counter )
-    {
-        size_t position = 0;
-        for ( size_t i = 0 ; i != this->multipliers.size() ; ++i )
-        {
-            position += this->multipliers[i]*counter[i];
-        }
-        return position;
+    All_cells_iterator operator++(int) {
+      All_cells_iterator result = *this;
+      ++(*this);
+      return result;
+    }
+
+    All_cells_iterator operator=(const All_cells_iterator& rhs) {
+      this->counter = rhs.counter;
+      return *this;
     }
 
-     std::vector<unsigned> compute_counter_for_given_cell( size_t cell )const
-    {
-        std::vector<unsigned> counter;

-        counter.reserve( this->sizes.size() );
-        for ( size_t dim = this->sizes.size() ; dim != 0 ; --dim )
-        {
-            counter.push_back(cell/this->multipliers[dim-1]);
-            cell = cell%this->multipliers[dim-1];
+    bool operator==(const All_cells_iterator& rhs)const {
+      if (this->counter != rhs.counter)return false;
+      return true;
+    }
+
+    bool operator!=(const All_cells_iterator& rhs)const {
+      return !(*this == rhs);
+    }
+
+    /*
+     * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as an argument of the following functions:
+     * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell boundary and coboundary and dimension
+     * and in function get_cell_data to get a filtration of a cell.
+     */
+    size_t operator*() {
+      return this->counter;
+    }
+    friend class Bitmap_cubical_complex_base;
+   protected:
+    size_t counter;
+  };
+
+  /**
+   * Function returning a All_cells_iterator to the first cell of the bitmap.
+   **/
+  All_cells_iterator all_cells_iterator_begin() {
+    All_cells_iterator a;
+    return a;
+  }
+
+  /**
+   * Function returning a All_cells_iterator to the last cell of the bitmap.
+   **/
+  All_cells_iterator all_cells_iterator_end() {
+    All_cells_iterator a;
+    a.counter = this->data.size();
+    return a;
+  }
+
+  /**
+   * All_cells_range class provides ranges for All_cells_iterator
+   **/
+  class All_cells_range {
+   public:
+
+    All_cells_range(Bitmap_cubical_complex_base* b) : b(b) { };
+
+    All_cells_iterator begin() {
+      return b->all_cells_iterator_begin();
+    }
+
+    All_cells_iterator end() {
+      return b->all_cells_iterator_end();
+    }
+   private:
+    Bitmap_cubical_complex_base<T>* b;
+  };
+
+  All_cells_range all_cells_range() {
+    return All_cells_range(this);
+  }
+
+
+  /**
+   * Boundary_range class provides ranges for boundary iterators.
+   **/
+  typedef typename std::vector< size_t >::const_iterator Boundary_iterator;
+  typedef typename std::vector< size_t > Boundary_range;
+
+  /**
+   * boundary_simplex_range creates an object of a Boundary_simplex_range class
+   * that provides ranges for the Boundary_simplex_iterator.
+   **/
+  Boundary_range boundary_range(size_t sh) {
+    return this->get_boundary_of_a_cell(sh);
+  }
+
+  /**
+   * Coboundary_range class provides ranges for boundary iterators.
+   **/
+  typedef typename std::vector< size_t >::const_iterator Coboundary_iterator;
+  typedef typename std::vector< size_t > Coboundary_range;
+
+  /**
+   * boundary_simplex_range creates an object of a Boundary_simplex_range class
+   * that provides ranges for the Boundary_simplex_iterator.
+   **/
+  Coboundary_range coboundary_range(size_t sh) {
+    return this->get_coboundary_of_a_cell(sh);
+  }
+
+  /**
+   * Iterator through top dimensional cells of the complex. The cells appear in order they are stored
+   * in the structure (i.e. in lexicographical order)
+   **/
+  class Top_dimensional_cells_iterator : std::iterator< std::input_iterator_tag, T > {
+   public:
+
+    Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) {
+      this->counter = std::vector<size_t>(b.dimension());
+      //std::fill( this->counter.begin() , this->counter.end() , 0 );
+    }
+
+    Top_dimensional_cells_iterator operator++() {
+      //first find first element of the counter that can be increased:
+      size_t dim = 0;
+      while ((dim != this->b.dimension()) && (this->counter[dim] == this->b.sizes[dim] - 1))++dim;
+
+      if (dim != this->b.dimension()) {
+        ++this->counter[dim];
+        for (size_t i = 0; i != dim; ++i) {
+          this->counter[i] = 0;
         }
-        std::reverse( counter.begin() , counter.end() );
-        return counter;
-    }

-    void read_perseus_style_file( const char* perseus_style_file );

-    void setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions , const std::vector<T>& top_dimensional_cells);

-    Bitmap_cubical_complex_base( const char* perseus_style_file , std::vector<bool> directions );

-    Bitmap_cubical_complex_base( const std::vector<unsigned>& sizes , std::vector<bool> directions );

-    Bitmap_cubical_complex_base( const std::vector<unsigned>& dimensions , const std::vector<T>& top_dimensional_cells , std::vector<bool> directions );
+      } else {
+        ++this->counter[0];
+      }
+      return *this;
+    }
+
+    Top_dimensional_cells_iterator operator++(int) {
+      Top_dimensional_cells_iterator result = *this;
+      ++(*this);
+      return result;
+    }
+
+    Top_dimensional_cells_iterator operator=(const Top_dimensional_cells_iterator& rhs) {
+      this->counter = rhs.counter;
+      this->b = rhs.b;
+      return *this;
+    }
+
+    bool operator==(const Top_dimensional_cells_iterator& rhs)const {
+      if (&this->b != &rhs.b)return false;
+      if (this->counter.size() != rhs.counter.size())return false;
+      for (size_t i = 0; i != this->counter.size(); ++i) {
+        if (this->counter[i] != rhs.counter[i])return false;
+      }
+      return true;
+    }
+
+    bool operator!=(const Top_dimensional_cells_iterator& rhs)const {
+      return !(*this == rhs);
+    }
+
+    /*
+     * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as an argument of the following functions:
+     * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell boundary and coboundary and dimension
+     * and in function get_cell_data to get a filtration of a cell.
+     */
+    size_t operator*() {
+      return this->compute_index_in_bitmap();
+    }
+
+    size_t compute_index_in_bitmap()const {
+      size_t index = 0;
+      for (size_t i = 0; i != this->counter.size(); ++i) {
+        index += (2 * this->counter[i] + 1) * this->b.multipliers[i];
+      }
+      return index;
+    }
+
+    void print_counter()const {
+      for (size_t i = 0; i != this->counter.size(); ++i) {
+        std::cout << this->counter[i] << " ";
+      }
+    }
+    friend class Bitmap_cubical_complex_base;
+   protected:
+    std::vector< size_t > counter;
+    Bitmap_cubical_complex_base& b;
+  };
+
+  /**
+   * Function returning a Top_dimensional_cells_iterator to the first top dimensional cell of the bitmap.
+   **/
+  Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin() {
+    Top_dimensional_cells_iterator a(*this);
+    return a;
+  }
+
+  /**
+   * Function returning a Top_dimensional_cells_iterator to the last top dimensional cell of the bitmap.
+   **/
+  Top_dimensional_cells_iterator top_dimensional_cells_iterator_end() {
+    Top_dimensional_cells_iterator a(*this);
+    for (size_t i = 0; i != this->dimension(); ++i) {
+      a.counter[i] = this->sizes[i] - 1;
+    }
+    a.counter[0]++;
+    return a;
+  }
+
+  /**
+   * Top_dimensional_cells_iterator_range class provides ranges for Top_dimensional_cells_iterator_range
+   **/
+  class Top_dimensional_cells_range {
+   public:
+
+    Top_dimensional_cells_range(Bitmap_cubical_complex_base* b) : b(b) { };
+
+    Top_dimensional_cells_iterator begin() {
+      return b->top_dimensional_cells_iterator_begin();
+    }
+
+    Top_dimensional_cells_iterator end() {
+      return b->top_dimensional_cells_iterator_end();
+    }
+   private:
+    Bitmap_cubical_complex_base<T>* b;
+  };
+
+  Top_dimensional_cells_range top_dimensional_cells_range() {
+    return Top_dimensional_cells_range(this);
+  }
+
+
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+
+  inline size_t number_cells()const {
+    return this->total_number_of_cells;
+  }
+
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+
+ protected:
+  std::vector<unsigned> sizes;
+  std::vector<unsigned> multipliers;
+  std::vector<T> data;
+  size_t total_number_of_cells;
+
+  void set_up_containers(const std::vector<unsigned>& sizes) {
+    unsigned multiplier = 1;
+    for (size_t i = 0; i != sizes.size(); ++i) {
+      this->sizes.push_back(sizes[i]);
+      this->multipliers.push_back(multiplier);
+      multiplier *= 2 * sizes[i] + 1;
+    }
+    this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max());
+    this->total_number_of_cells = multiplier;
+  }
+
+  size_t compute_position_in_bitmap(const std::vector< unsigned >& counter) {
+    size_t position = 0;
+    for (size_t i = 0; i != this->multipliers.size(); ++i) {
+      position += this->multipliers[i] * counter[i];
+    }
+    return position;
+  }
+
+  std::vector<unsigned> compute_counter_for_given_cell(size_t cell)const {
+    std::vector<unsigned> counter;
+    counter.reserve(this->sizes.size());
+    for (size_t dim = this->sizes.size(); dim != 0; --dim) {
+      counter.push_back(cell / this->multipliers[dim - 1]);
+      cell = cell % this->multipliers[dim - 1];
+    }
+    std::reverse(counter.begin(), counter.end());
+    return counter;
+  }
+  void read_perseus_style_file(const char* perseus_style_file);
+  void setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells);
+  Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector<bool> directions);
+  Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, std::vector<bool> directions);
+  Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells, std::vector<bool> directions);
 };
 
+template <typename T>
+void Bitmap_cubical_complex_base<T>::put_data_to_bins(size_t number_of_bins) {
+  bool bdg = false;
+
+  std::pair< T, T > min_max = this->min_max_filtration();
+  T dx = (min_max.second - min_max.first) / (T) number_of_bins;
+
+  //now put the data into the appropriate bins:
+  for (size_t i = 0; i != this->data.size(); ++i) {
+    if (bdg) {
+      std::cerr << "Before binning : " << this->data[i] << std::endl;
+    }
+    this->data[i] = min_max.first + dx * (this->data[i] - min_max.first) / number_of_bins;
+    if (bdg) {
+      std::cerr << "After binning : " << this->data[i] << std::endl;
+      getchar();
+    }
+  }
+}
 
-template <typename T>

-void Bitmap_cubical_complex_base<T>::put_data_to_bins( size_t number_of_bins )

-{

-    bool bdg = false;

-

-    std::pair< T ,T > min_max = this->min_max_filtration();

-    T dx = (min_max.second-min_max.first)/(T)number_of_bins;

-

-    //now put the data into the appropriate bins:

-    for ( size_t i = 0 ; i != this->data.size() ; ++i )

-    {

-        if ( bdg ){cerr << "Before binning : " << this->data[i] << endl;}

-        this->data[i] = min_max.first + dx*(this->data[i]-min_max.first)/number_of_bins;

-        if ( bdg ){cerr << "After binning : " << this->data[i] << endl;getchar();}

-    }

-}

-

 template <typename T>
-void Bitmap_cubical_complex_base<T>::put_data_to_bins( T diameter_of_bin )

-{

-    bool bdg = false;

-    std::pair< T ,T > min_max = this->min_max_filtration();

-

-    size_t number_of_bins = (min_max.second - min_max.first)/diameter_of_bin;

-    //now put the data into the appropriate bins:

-    for ( size_t i = 0 ; i != this->data.size() ; ++i )

-    {

-        if ( bdg ){cerr << "Before binning : " << this->data[i] << endl;}

-        this->data[i] = min_max.first + diameter_of_bin*(this->data[i]-min_max.first)/number_of_bins;

-        if ( bdg ){cerr << "After binning : " << this->data[i] << endl;getchar();}

-    }

-}

-

-template <typename T>

-std::pair< T ,T > Bitmap_cubical_complex_base<T>::min_max_filtration()

-{

-    std::pair< T ,T > min_max( std::numeric_limits<T>::max() , std::numeric_limits<T>::min() );

-    for ( size_t i = 0 ; i != this->data.size() ; ++i )

-    {

-        if ( this->data[i] < min_max.first )min_max.first = this->data[i];

-        if ( this->data[i] > min_max.second )min_max.second = this->data[i];

-    }

-    return min_max;

-}

+void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) {
+  bool bdg = false;
+  std::pair< T, T > min_max = this->min_max_filtration();
+
+  size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin;
+  //now put the data into the appropriate bins:
+  for (size_t i = 0; i != this->data.size(); ++i) {
+    if (bdg) {
+      std::cerr << "Before binning : " << this->data[i] << std::endl;
+    }
+    this->data[i] = min_max.first + diameter_of_bin * (this->data[i] - min_max.first) / number_of_bins;
+    if (bdg) {
+      std::cerr << "After binning : " << this->data[i] << std::endl;
+      getchar();
+    }
+  }
+}
 
+template <typename T>
+std::pair< T, T > Bitmap_cubical_complex_base<T>::min_max_filtration() {
+  std::pair< T, T > min_max(std::numeric_limits<T>::max(), std::numeric_limits<T>::min());
+  for (size_t i = 0; i != this->data.size(); ++i) {
+    if (this->data[i] < min_max.first)min_max.first = this->data[i];
+    if (this->data[i] > min_max.second)min_max.second = this->data[i];
+  }
+  return min_max;
+}
 
 template <typename K>
-ostream& operator << ( ostream & out , const Bitmap_cubical_complex_base<K>& b )
-{
-    for ( typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator
-          it = b.all_cells_const_begin() ; it != b.all_cells_const_end() ; ++it )
-    {
-        out << *it << " ";
-    }
-    return out;
+ostream& operator<<(ostream & out, const Bitmap_cubical_complex_base<K>& b) {
+  for (typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator
+       it = b.all_cells_const_begin(); it != b.all_cells_const_end(); ++it) {
+    out << *it << " ";
+  }
+  return out;
 }
 
-
 template <typename T>
 Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base
-( const std::vector<unsigned>& sizes )
-{
-    this->set_up_containers( sizes );
+(const std::vector<unsigned>& sizes) {
+  this->set_up_containers(sizes);
 }
-

-template <typename T>

-void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions , const std::vector<T>& top_dimensional_cells)

-{

-    this->set_up_containers( sizes_in_following_directions );
-
-    size_t number_of_top_dimensional_elements = 1;
-    for ( size_t i = 0 ; i != sizes_in_following_directions.size() ; ++i )
-    {
-        number_of_top_dimensional_elements *= sizes_in_following_directions[i];
-    }
-    if ( number_of_top_dimensional_elements != top_dimensional_cells.size() )
-    {
-        cerr <<
-       "Error in constructor\
+
+template <typename T>
+void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) {
+  this->set_up_containers(sizes_in_following_directions);
+
+  size_t number_of_top_dimensional_elements = 1;
+  for (size_t i = 0; i != sizes_in_following_directions.size(); ++i) {
+    number_of_top_dimensional_elements *= sizes_in_following_directions[i];
+  }
+  if (number_of_top_dimensional_elements != top_dimensional_cells.size()) {
+    std::cerr <<
+        "Error in constructor\
         Bitmap_cubical_complex_base\
        ( std::vector<size_t> sizes_in_following_directions , std::vector<T> top_dimensional_cells ).\
        Number of top dimensional elements that follow from sizes_in_following_directions vector is different\
-       than the size of top_dimensional_cells vector." << endl;
-        throw("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions,\
+       than the size of top_dimensional_cells vector." << std::endl;
+    throw ("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions,\
         std::vector<T> top_dimensional_cells )\
       . Number of top dimensional elements that follow from sizes_in_following_directions vector is different than the\
        size of top_dimensional_cells vector.");
-    }
-
-    Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
-    size_t index = 0;
-    for ( it = this->top_dimensional_cells_iterator_begin() ; it != this->top_dimensional_cells_iterator_end() ; ++it )
-    {
-        this->get_cell_data(*it) = top_dimensional_cells[index];
-        ++index;
-    }
-    this->impose_lower_star_filtration();

-}

+  }
+
+  Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
+  size_t index = 0;
+  for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
+    this->get_cell_data(*it) = top_dimensional_cells[index];
+    ++index;
+  }
+  this->impose_lower_star_filtration();
+}
 
 template <typename T>
 Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base
-( const std::vector<unsigned>& sizes_in_following_directions , const std::vector<T>& top_dimensional_cells )
-{
-    this->setup_bitmap_based_on_top_dimensional_cells_list( sizes_in_following_directions , top_dimensional_cells );
-}

-

-template <typename T>

-void Bitmap_cubical_complex_base<T>::read_perseus_style_file( const char* perseus_style_file )

-{

-    bool dbg = false;
-    ifstream inFiltration, inIds;
-    inFiltration.open( perseus_style_file );
-    unsigned dimensionOfData;
-    inFiltration >> dimensionOfData;
-
-    if (dbg){cerr << "dimensionOfData : " << dimensionOfData << endl;getchar();}
-
-    std::vector<unsigned> sizes;

-    sizes.reserve( dimensionOfData );
-    for ( size_t i = 0 ; i != dimensionOfData ; ++i )
-    {
-        unsigned size_in_this_dimension;
-        inFiltration >> size_in_this_dimension;
-        size_in_this_dimension = size_in_this_dimension;
-        sizes.push_back( size_in_this_dimension );
-        if (dbg){cerr << "size_in_this_dimension : " << size_in_this_dimension << endl;}
-    }
-    this->set_up_containers( sizes );
-
-    Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
-    it = this->top_dimensional_cells_iterator_begin();
-
-    while ( !inFiltration.eof() )
-    {
-        T filtrationLevel;
-        inFiltration >> filtrationLevel;
-        if ( dbg )
-        {
-            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 << endl;
-        }
-        this->get_cell_data(*it) = filtrationLevel;
-        ++it;
-    }
-    inFiltration.close();
-    this->impose_lower_star_filtration();

-}

-

-template <typename T>

-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const char* perseus_style_file , std::vector<bool> directions )

-{

-    //this constructor is here just for compatibility with a class that creates cubical complexes with periodic bundary conditions.

-    //It ignores the last parameter of the function.

-    this->read_perseus_style_file( perseus_style_file );

-}

-

-template <typename T>

-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const std::vector<unsigned>& sizes , std::vector<bool> directions )

-{

-    //this constructor is here just for compatibility with a class that creates cubical complexes with periodic bundary conditions.

-    //It ignores the last parameter of the function.

-    this->set_up_containers( sizes );

-}

-

-template <typename T>

-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const std::vector<unsigned>& dimensions , const std::vector<T>& top_dimensional_cells , std::vector<bool> directions )
-{

-    //this constructor is here just for compatibility with a class that creates cubical complexes with periodic bundary conditions.

-    //It ignores the last parameter of the function.

-    this->setup_bitmap_based_on_top_dimensional_cells_list( dimensions , top_dimensional_cells );

+(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) {
+  this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells);
 }
 
 template <typename T>
-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const char* perseus_style_file )
-{
-    this->read_perseus_style_file( perseus_style_file );
+void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus_style_file) {
+  bool dbg = false;
+  std::ifstream inFiltration, inIds;
+  inFiltration.open(perseus_style_file);
+  unsigned dimensionOfData;
+  inFiltration >> dimensionOfData;
+
+  if (dbg) {
+    std::cerr << "dimensionOfData : " << dimensionOfData << std::endl;
+    getchar();
+  }
+
+  std::vector<unsigned> sizes;
+  sizes.reserve(dimensionOfData);
+  for (size_t i = 0; i != dimensionOfData; ++i) {
+    unsigned size_in_this_dimension;
+    inFiltration >> size_in_this_dimension;
+    size_in_this_dimension = size_in_this_dimension;
+    sizes.push_back(size_in_this_dimension);
+    if (dbg) {
+      std::cerr << "size_in_this_dimension : " << size_in_this_dimension << std::endl;
+    }
+  }
+  this->set_up_containers(sizes);
+
+  Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
+  it = this->top_dimensional_cells_iterator_begin();
+
+  while (!inFiltration.eof()) {
+    T filtrationLevel;
+    inFiltration >> filtrationLevel;
+    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>
-std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell( size_t cell )const
-{
-    std::vector< size_t > boundary_elements;

-

-    //Speed traded of for memory. Check if it is better in practice.

-    boundary_elements.reserve( this->dimension()*2 );

-
-    size_t cell1 = cell;
-    for ( size_t i = this->multipliers.size() ; i != 0 ; --i )
-    {
-        unsigned position = cell1/this->multipliers[i-1];
-        if ( position%2 == 1 )
-        {

-            boundary_elements.push_back( cell - this->multipliers[ i-1 ] );

-            boundary_elements.push_back( cell + this->multipliers[ i-1 ] );
-        }
-        cell1 = cell1%this->multipliers[i-1];
-    }
-    return boundary_elements;
-}

-
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector<bool> directions) {
+  //this constructor is here just for compatibility with a class that creates cubical complexes with periodic bundary conditions.
+  //It ignores the last parameter of the function.
+  this->read_perseus_style_file(perseus_style_file);
+}
 
-

-

-

 template <typename T>
-std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell( size_t cell )const
-{
-    std::vector<unsigned> counter = this->compute_counter_for_given_cell( cell );
-    std::vector< size_t > coboundary_elements;
-    size_t cell1 = cell;
-    for ( size_t i = this->multipliers.size() ; i != 0 ; --i )
-    {
-        unsigned position = cell1/this->multipliers[i-1];
-        if ( position%2 == 0 )
-        {
-            if ( (cell > this->multipliers[i-1]) && (counter[i-1] != 0) )
-            {
-                coboundary_elements.push_back( cell - this->multipliers[i-1] );
-            }
-            if (
-               (cell + this->multipliers[i-1] < this->data.size()) && (counter[i-1] != 2*this->sizes[i-1]) )
-            {
-                coboundary_elements.push_back( cell + this->multipliers[i-1] );
-            }
-        }
-        cell1 = cell1%this->multipliers[i-1];
-    }
-    return coboundary_elements;
-}

-

-
-
-
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, std::vector<bool> directions) {
+  //this constructor is here just for compatibility with a class that creates cubical complexes with periodic bundary conditions.
+  //It ignores the last parameter of the function.
+  this->set_up_containers(sizes);
+}
 
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells, std::vector<bool> directions) {
+  //this constructor is here just for compatibility with a class that creates cubical complexes with periodic bundary conditions.
+  //It ignores the last parameter of the function.
+  this->setup_bitmap_based_on_top_dimensional_cells_list(dimensions, top_dimensional_cells);
+}
 
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file) {
+  this->read_perseus_style_file(perseus_style_file);
+}
 
+template <typename T>
+std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell)const {
+  std::vector< size_t > boundary_elements;
+
+  //Speed traded of for memory. Check if it is better in practice.
+  boundary_elements.reserve(this->dimension()*2);
+
+  size_t cell1 = cell;
+  for (size_t i = this->multipliers.size(); i != 0; --i) {
+    unsigned position = cell1 / this->multipliers[i - 1];
+    if (position % 2 == 1) {
+      boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+      boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+    }
+    cell1 = cell1 % this->multipliers[i - 1];
+  }
+  return boundary_elements;
+}
 
 template <typename T>
-unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell( size_t cell )const
-{
-    bool dbg = false;
-    if (dbg)cerr << "\n\n\n Computing position o a cell of an index : " << cell << endl;
-    unsigned dimension = 0;
-    for ( size_t i = this->multipliers.size() ; i != 0 ; --i )
-    {
-        unsigned position = cell/this->multipliers[i-1];
-
-        if (dbg)cerr << "i-1 :" << i-1 << endl;
-        if (dbg)cerr << "cell : " << cell << endl;
-        if (dbg)cerr << "position : " << position << endl;
-        if (dbg)cerr << "multipliers["<<i-1<<"] = " << this->multipliers[i-1] << endl;
-        if (dbg)getchar();
-
-        if ( position%2 == 1 )
-        {
-            if (dbg)cerr << "Nonzero length in this direction \n";
-            dimension++;
-        }
-        cell = cell%this->multipliers[i-1];
+std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell)const {
+  std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
+  std::vector< size_t > coboundary_elements;
+  size_t cell1 = cell;
+  for (size_t i = this->multipliers.size(); i != 0; --i) {
+    unsigned position = cell1 / this->multipliers[i - 1];
+    if (position % 2 == 0) {
+      if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) {
+        coboundary_elements.push_back(cell - this->multipliers[i - 1]);
+      }
+      if (
+          (cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) {
+        coboundary_elements.push_back(cell + this->multipliers[i - 1]);
+      }
     }
-    return dimension;
+    cell1 = cell1 % this->multipliers[i - 1];
+  }
+  return coboundary_elements;
 }
 
 template <typename T>
-inline T& Bitmap_cubical_complex_base<T>::get_cell_data( size_t cell )
-{
-    return this->data[cell];
+unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell)const {
+  bool dbg = false;
+  if (dbg) std::cerr << "\n\n\n Computing position o a cell of an index : " << cell << std::endl;
+  unsigned dimension = 0;
+  for (size_t i = this->multipliers.size(); i != 0; --i) {
+    unsigned position = cell / this->multipliers[i - 1];
+
+    if (dbg) {
+      std::cerr << "i-1 :" << i - 1 << std::endl;
+      std::cerr << "cell : " << cell << std::endl;
+      std::cerr << "position : " << position << std::endl;
+      std::cerr << "multipliers[" << i - 1 << "] = " << this->multipliers[i - 1] << std::endl;
+      getchar();
+    }
+
+    if (position % 2 == 1) {
+      if (dbg) std::cerr << "Nonzero length in this direction \n";
+      dimension++;
+    }
+    cell = cell % this->multipliers[i - 1];
+  }
+  return dimension;
 }
 
+template <typename T>
+inline T& Bitmap_cubical_complex_base<T>::get_cell_data(size_t cell) {
+  return this->data[cell];
+}
 
 template <typename T>
-void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration()
-{
-    bool dbg = false;
-
-    //this vector will be used to check which elements have already been taken care of
-    //in imposing lower star filtration:
-    std::vector<bool> is_this_cell_considered( this->data.size() , false );
-

-    size_t size_to_reserve = 1;

-    for ( size_t i = 0 ; i != this->multipliers.size() ; ++i )

-    {

-        size_to_reserve *= (size_t)((this->multipliers[i]-1)/2);

-    }

-
-    std::vector<size_t> indices_to_consider;

-    indices_to_consider.reserve( size_to_reserve );
-    //we assume here that we already have a filtration on the top dimensional cells and
-    //we have to extend it to lower ones.
-    typename Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
-    for ( it = this->top_dimensional_cells_iterator_begin() ; it != this->top_dimensional_cells_iterator_end() ; ++it )
-    {
-        indices_to_consider.push_back( it.compute_index_in_bitmap() );
-    }
-
-    while ( indices_to_consider.size() )
-    {
-        if ( dbg )
-        {
-            cerr << "indices_to_consider in this iteration \n";
-            for ( size_t i = 0 ; i != indices_to_consider.size() ; ++i )
-            {
-                cout << indices_to_consider[i] << "  ";
-            }
+void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() {
+  bool dbg = false;
+
+  //this vector will be used to check which elements have already been taken care of
+  //in imposing lower star filtration:
+  std::vector<bool> is_this_cell_considered(this->data.size(), false);
+
+  size_t size_to_reserve = 1;
+  for (size_t i = 0; i != this->multipliers.size(); ++i) {
+    size_to_reserve *= (size_t) ((this->multipliers[i] - 1) / 2);
+  }
+
+  std::vector<size_t> indices_to_consider;
+  indices_to_consider.reserve(size_to_reserve);
+  //we assume here that we already have a filtration on the top dimensional cells and
+  //we have to extend it to lower ones.
+  typename Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
+  for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
+    indices_to_consider.push_back(it.compute_index_in_bitmap());
+  }
+
+  while (indices_to_consider.size()) {
+    if (dbg) {
+      std::cerr << "indices_to_consider in this iteration \n";
+      for (size_t i = 0; i != indices_to_consider.size(); ++i) {
+        std::cout << indices_to_consider[i] << "  ";
+      }
+      getchar();
+    }
+    std::vector<size_t> new_indices_to_consider;
+    for (size_t i = 0; i != indices_to_consider.size(); ++i) {
+      std::vector<size_t> bd = this->get_boundary_of_a_cell(indices_to_consider[i]);
+      for (size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) {
+        if (dbg) {
+          std::cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[ bd[boundaryIt] ] << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[ indices_to_consider[i] ] << std::endl;
+          getchar();
+
+        }
+        if (this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ]) {
+          this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ];
+          if (dbg) {
+            std::cerr << "Setting the value of a cell : " << bd[boundaryIt] << " to : " << this->data[ indices_to_consider[i] ] << std::endl;
             getchar();
+          }
         }
-        std::vector<size_t> new_indices_to_consider;
-        for ( size_t i = 0 ; i != indices_to_consider.size() ; ++i )
-        {
-            std::vector<size_t> bd = this->get_boundary_of_a_cell( indices_to_consider[i] );
-            for ( size_t boundaryIt = 0 ; boundaryIt != bd.size() ; ++boundaryIt )
-            {

-                if ( dbg )

-                {

-                    cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[ bd[boundaryIt] ] << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[ indices_to_consider[i] ] << endl;

-                    getchar();

-

-                }
-                if ( this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ] )
-                {
-                    this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ];

-                    if ( dbg )

-                    {

-                        cerr << "Setting the value of a cell : " << bd[boundaryIt] << " to : " << this->data[ indices_to_consider[i] ] << endl;

-                        getchar();

-                    }
-                }
-                if ( is_this_cell_considered[ bd[boundaryIt] ] == false )
-                {
-                    new_indices_to_consider.push_back( bd[boundaryIt] );
-                    is_this_cell_considered[ bd[boundaryIt] ] = true;
-                }
-            }
+        if (is_this_cell_considered[ bd[boundaryIt] ] == false) {
+          new_indices_to_consider.push_back(bd[boundaryIt]);
+          is_this_cell_considered[ bd[boundaryIt] ] = true;
         }
-        indices_to_consider.swap(new_indices_to_consider);
+      }
     }
+    indices_to_consider.swap(new_indices_to_consider);
+  }
 }
 
-
 template <typename T>
-bool compareFirstElementsOfTuples( const std::pair< std::pair< T , size_t > , char >& first ,
-                                   const std::pair< std::pair< T , size_t > , char >& second )
-{
-    if ( first.first.first < second.first.first )
-    {
-        return true;
-    }
-    else
-    {
-        if ( first.first.first > second.first.first )
-        {
-            return false;
-        }
-        //in this case first.first.first == second.first.first, so we need to compare dimensions
-        return first.second < second.second;
+bool compareFirstElementsOfTuples(const std::pair< std::pair< T, size_t >, char >& first,
+                                  const std::pair< std::pair< T, size_t >, char >& second) {
+  if (first.first.first < second.first.first) {
+    return true;
+  } else {
+    if (first.first.first > second.first.first) {
+      return false;
     }
+    //in this case first.first.first == second.first.first, so we need to compare dimensions
+    return first.second < second.second;
+  }
 }
 
+} // namespace Cubical_complex
 
+} // namespace Gudhi
 
-
-}
-
-}

+#endif  // BITMAP_CUBICAL_COMPLEX_BASE_H_
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
index 956e74a7..2c0d77fe 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
@@ -18,348 +18,289 @@
  *
  *    You should have received a copy of the GNU General Public License
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */

-

-#pragma once

-#include <cmath>

-#include "Bitmap_cubical_complex_base.h"

-

-using namespace std;

-

-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

-

-

+ */
+
+#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>
+
+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
+
 /**
- *@class Bitmap_cubical_complex_periodic_boundary_conditions_base
- *@brief Cubical complex with periodic boundary conditions represented as a bitmap.
- *@ingroup cubical_complex
- */

+ * @class Bitmap_cubical_complex_periodic_boundary_conditions_base
+ * @brief Cubical complex with periodic boundary conditions represented as a bitmap.
+ * @ingroup cubical_complex
+ */
 /**
-* This is a class implementing a bitmap data structure with periodic boundary conditions. Most of the functions are identical to the functions from Bitmap_cubical_complex_base.

+* This is a class implementing a bitmap data structure with periodic boundary conditions. Most of the functions are identical to the functions from Bitmap_cubical_complex_base.
 * The ones that needed to be updated are the constructors and get_boundary_of_a_cell and get_coboundary_of_a_cell.
-*/

-template <typename T>

-class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T>

-{

-public:

-    //constructors that take an extra parameter:

-    /**
-    * Default constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class.
-    */

-    Bitmap_cubical_complex_periodic_boundary_conditions_base(){};

-

-    /**
-    * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of booleans. If at i-th position of this vector there is true value, that means that periodic boundary conditions are to be imposed in this direction. In case of false, the periodic boundary conditions will not be imposed in the direction i.
-    */

-    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 );

-    /**
-    * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the name of Perseus style file as an input. Please consult the documentation about the specification of the file.
-    */

-    Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseusStyleFile );

-    /**
-    * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of top dimensional cells (ordered lexicographically) and (3) vector of booleans. If at i-th position of this vector there is true value, that means that periodic boundary conditions are to be imposed in this direction. In case of false, the periodic boundary conditions will not be imposed in the direction i.
-    */

-    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 );

-

-    /**
-    * Destructor of the Bitmap_cubical_complex_periodic_boundary_conditions_base class.
-    **/

-    virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base(){}

-

-    //overwritten methods co compute boundary and coboundary

-    /**
-    * A version of a function that return boundary of a given cell for an object of Bitmap_cubical_complex_periodic_boundary_conditions_base class.
-    */

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

-    /**
-    * A version of a function that return coboundary of a given cell for an object of Bitmap_cubical_complex_periodic_boundary_conditions_base class.
-    */

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

-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 ( 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>::max());

-        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 );

-

-    size_t i = 0;

-    for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator 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;

-    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 ( 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 )

-        {

-            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 << endl;

-        }

-        this->get_cell_data(*it) = filtrationLevel;

-        ++it;

-    }

-    inFiltration.close();

-    this->impose_lower_star_filtration();

-

-/*

-     char* filename = (char*)perseus_style_file;

-     //char* filename = "combustionWithPeriodicBoundaryConditions/v0/tV0_000000.float";

-     ifstream file( filename , ios::binary | ios::ate );

-     unsigned realSizeOfFile = file.tellg();

-     file.close();

-     realSizeOfFile = realSizeOfFile/sizeof(T);

-

-     unsigned w, h, d;

-

-     w = h = d = ceil(pow( realSizeOfFile , (double)(1/(double)3) ));

-

-     T* slice = new T[w*h*d];

-     if (slice == NULL)

-     {

-        cerr << "Allocation error, cannot allocate " << w*h*d*sizeof(T) << " bytes to store the data from the file. The program will now terminate \n";

-        exit(EXIT_FAILURE);

-     }

-

-     FILE* fp;

-     if ((fp=fopen( filename, "rb" )) == NULL )

-     {

-         cerr << "Cannot open the file: " << filename << ". The program will now terminate \n";

-         exit(1);

-     }

-

-     clock_t read_begin = clock();

-     fread( slice,4,w*h*d,fp );

-     fclose(fp);

-     cerr << "Time of reading the file : " << double(clock() - read_begin) / CLOCKS_PER_SEC << endl;

-

-

-     clock_t begin_creation_bitap = clock();

-     std::vector<T> data(slice,slice+w*h*d);

-     delete[] slice;

-     std::vector< unsigned > sizes;

-     sizes.push_back(w);

-     sizes.push_back(w);

-     sizes.push_back(w);

-

-    this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true );

-    this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true );

-    this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true );

-    this->set_up_containers( sizes );

-

-    size_t i = 0;

-    for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator it = this->top_dimensional_cells_iterator_begin() ; it != this->top_dimensional_cells_iterator_end() ; ++it )

-    {

-        *it = data[i];

-        ++i;

-    }

-    this->impose_lower_star_filtration();

-    cerr << "Time of creation of a bitmap : " << double(clock() - begin_creation_bitap ) / CLOCKS_PER_SEC << endl;

-    */

-}

-

-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< size_t >  Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell( size_t cell )const

-{

-    bool dbg = false;

-    if ( dbg ){cerr << "Computations of boundary of a cell : " << cell << endl;}

-

-    std::vector< size_t > boundary_elements;

-    size_t cell1 = cell;

-    for ( 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] )

-            {

-                //cerr << "A\n";

-                boundary_elements.push_back( cell - this->multipliers[ i-1 ] );

-                boundary_elements.push_back( cell + this->multipliers[ i-1 ] );

-                if (dbg){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 )

-                {

-                    //cerr << "B\n";

-                    boundary_elements.push_back( cell - this->multipliers[ i-1 ] );

-                    boundary_elements.push_back( cell + this->multipliers[ i-1 ] );

-                    if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";}

-                }

-                else

-                {

-                    //cerr << "C\n";

-                    boundary_elements.push_back( cell - this->multipliers[ i-1 ] );

-                    boundary_elements.push_back( cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] );

-                    if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] << " ";}

-                }

-            }

-        }

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

-    }

-    return boundary_elements;

-}

-

-template <typename T>

-std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell( size_t cell )const

-{

-    std::vector<unsigned> counter = this->compute_counter_for_given_cell( cell );

-    std::vector< size_t > coboundary_elements;

-    size_t cell1 = cell;

-    for ( 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;

-}

-

-

-

-}//Cubical_complex

-}//namespace Gudhi

+*/
+template <typename T>
+class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T> {
+ public:
+  // constructors that take an extra parameter:
+
+  /**
+   * Default constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   */
+  Bitmap_cubical_complex_periodic_boundary_conditions_base() { }
+  /**
+   * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following
+   * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of booleans. If
+   * at i-th position of this vector there is true value, that means that periodic boundary conditions are to be
+   * imposed in this direction. In case of false, the periodic boundary conditions will not be imposed in the direction
+   * i.
+   */
+  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);
+  /**
+   * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the name of Perseus
+   * style file as an input. Please consult the documentation about the specification of the file.
+   */
+  Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseusStyleFile);
+  /**
+   * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following
+   * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of top dimensional
+   * cells (ordered lexicographically) and (3) vector of booleans. If at i-th position of this vector there is true
+   * value, that means that periodic boundary conditions are to be imposed in this direction. In case of false, the
+   * periodic boundary conditions will not be imposed in the direction i.
+   */
+  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);
+
+  /**
+   * Destructor of the Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   **/
+  virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base(){}
+
+  // overwritten methods co compute boundary and coboundary
+  /**
+   * A version of a function that return boundary of a given cell for an object of Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   */
+  virtual std::vector< size_t > get_boundary_of_a_cell(size_t cell)const;
+
+  /**
+   * A version of a function that return coboundary of a given cell for an object of Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   */
+  virtual std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const;
+
+ 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 (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>::max());
+    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);
+
+  size_t i = 0;
+  for (typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator 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 (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< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell(size_t cell) const {
+  bool dbg = false;
+  if (dbg) {
+    std::cerr << "Computations of boundary of a cell : " << cell << std::endl;
+  }
+
+  std::vector< size_t > boundary_elements;
+  size_t cell1 = cell;
+  for (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";
+        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";
+          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";
+          boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+          boundary_elements.push_back(cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ]);
+          if (dbg) {
+            std::cerr << cell - this->multipliers[ i - 1 ] << " " <<
+                cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ] << " ";
+          }
+        }
+      }
+    }
+    cell1 = cell1 % this->multipliers[i - 1];
+  }
+  return boundary_elements;
+}
+
+template <typename T>
+std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell(size_t cell) const {
+  std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
+  std::vector< size_t > coboundary_elements;
+  size_t cell1 = cell;
+  for (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 Gudhi
+
+#endif  // BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_