Adding changes to the cubcial complex class according to Marc and Vincent's comments.

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

Former-commit-id: cb9134a5c18fae3e2e63bbaf8329dd168d08d3b3

7 files changed, 2463 insertions, 2232 deletions

diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
index 37c16618..31da3609 100644
--- a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
+++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
@@ -1,71 +1,76 @@
-/*    This file is part of the Gudhi Library. The Gudhi library

- *    (Geometric Understanding in Higher Dimensions) is a generic C++

- *    library for computational topology.

- *

- *    Author(s):       Pawel Dlotko

- *

- *    Copyright (C) 2015  INRIA Sophia-Saclay (France)

- *

- *    This program is free software: you can redistribute it and/or modify

- *    it under the terms of the GNU General Public License as published by

- *    the Free Software Foundation, either version 3 of the License, or

- *    (at your option) any later version.

- *

- *    This program is distributed in the hope that it will be useful,

- *    but WITHOUT ANY WARRANTY; without even the implied warranty of

- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

- *    GNU General Public License for more details.

- *

- *    You should have received a copy of the GNU General Public License

- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.

- */

-

-// for persistence algorithm

-#include <gudhi/reader_utils.h>

-#include <gudhi/Bitmap_cubical_complex.h>

-#include <gudhi/Persistent_cohomology.h>

-

-#include <boost/program_options.hpp>

-

-// standard stuff

-#include <iostream>

-#include <sstream>

-

-using namespace Gudhi;

-using namespace Gudhi::persistent_cohomology;

-using namespace std;

-

-int main(int argc, char** argv) {

-  cout << "This program computes persistent homology, by using Bitmap_cubical_complex class, of cubical complexes "

-      "provided in text files in Perseus style (the only numbed in the first line is a dimension D of a cubical "

-      "complex. In the lines I between 2 and D+1 there are numbers of top dimensional cells in the direction I. Let N "

-      "denote product of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are filtrations of "

-      "top dimensional cells. We assume that the cells are in the lexicographical order. See CubicalOneSphere.txt or "

-      "CubicalTwoSphere.txt for example." << endl;

-

-  int p = 2;

-  double min_persistence = 0;

-

-  if (argc != 2) {

-    cout << "Wrong number of parameters. Please provide the name of a file with a Perseus style cubical complex at the "

-        "input. The program will now terminate.\n";

-    return 1;

-  }

-

-  Bitmap_cubical_complex<double> b(argv[1]);

-

-

-  // Compute the persistence diagram of the complex

-  persistent_cohomology::Persistent_cohomology< Bitmap_cubical_complex<double>, Field_Zp > pcoh(b);

-  pcoh.init_coefficients(p);  // initializes the coefficient field for homology

-  pcoh.compute_persistent_cohomology(min_persistence);

-

-

-  stringstream ss;

-  ss << argv[1] << "_persistence";

-  std::ofstream out((char*) ss.str().c_str());

-  pcoh.output_diagram(out);

-  out.close();

-

-  return 0;

-}

+ /*    This file is part of the Gudhi Library. The Gudhi library
+  *    (Geometric Understanding in Higher Dimensions) is a generic C++
+  *    library for computational topology.
+  *
+  *    Author(s):       Pawel Dlotko
+  *
+  *    Copyright (C) 2015  INRIA Saclay (France)
+  *
+  *    This program is free software: you can redistribute it and/or modify
+  *    it under the terms of the GNU General Public License as published by
+  *    the Free Software Foundation, either version 3 of the License, or
+  *    (at your option) any later version.
+  *
+  *    This program is distributed in the hope that it will be useful,
+  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  *    GNU General Public License for more details.
+  *
+  *    You should have received a copy of the GNU General Public License
+  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+  */
+
+
+//for persistence algorithm
+#include <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+
+using namespace Gudhi;
+using namespace Gudhi::Cubical_complex;
+using namespace Gudhi::persistent_cohomology;
+
+//standard stuff
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+int main( int argc , char** argv )
+{
+    cout << "This program computes persistent homology, by using bitmap_cubical_complex class,\
+ of cubical complexes provided in text files in Perseus style (the only numbered in \
+the first line is a dimension D of a bitmap. In the lines I between 2 and D+1 there are\
+ numbers of top dimensional cells in the direction I. Let N denote product \
+of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are\
+ filtrations of top dimensional cells. We assume that the cells are in the \
+lexicographical order. See CubicalOneSphere.txt or CubicalTwoSphere.txt for example." << endl;
+
+    int p = 2;
+    double min_persistence = 0;
+
+    if ( argc != 2 )
+    {
+        cout << "Wrong number of parameters. Please provide the name of a file with a\
+                 Perseus style bitmap at the input. The program will now terminate.\n";
+        return 1;
+    }
+
+    Bitmap_cubical_complex<double> b( argv[1] );
+
+
+    // Compute the persistence diagram of the complex
+    persistent_cohomology::Persistent_cohomology< Bitmap_cubical_complex<double>, Field_Zp > pcoh(b);
+    pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
+    pcoh.compute_persistent_cohomology( min_persistence );
+
+
+    stringstream ss;
+    ss << argv[1] << "_persistence";
+    std::ofstream out((char*)ss.str().c_str());
+    pcoh.output_diagram(out);
+    out.close();
+
+    return 0;
+}
diff --git a/src/Bitmap_cubical_complex/example/CMakeLists.txt b/src/Bitmap_cubical_complex/example/CMakeLists.txt
index 05ef1319..dd252a79 100644
--- a/src/Bitmap_cubical_complex/example/CMakeLists.txt
+++ b/src/Bitmap_cubical_complex/example/CMakeLists.txt
@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 2.6)
-project(GUDHISimplexTreeFromFile)
+project(GUDHIBitmap)
 
 add_executable ( Bitmap_cubical_complex Bitmap_cubical_complex.cpp )
 target_link_libraries(Bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY})	
diff --git a/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp
index ac7557ce..60cfc113 100644
--- a/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp
+++ b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp
@@ -1,81 +1,93 @@
-/*    This file is part of the Gudhi Library. The Gudhi library

- *    (Geometric Understanding in Higher Dimensions) is a generic C++

- *    library for computational topology.

- *

- *    Author(s):       Pawel Dlotko

- *

- *    Copyright (C) 2015  INRIA Saclay (France)

- *

- *    This program is free software: you can redistribute it and/or modify

- *    it under the terms of the GNU General Public License as published by

- *    the Free Software Foundation, either version 3 of the License, or

- *    (at your option) any later version.

- *

- *    This program is distributed in the hope that it will be useful,

- *    but WITHOUT ANY WARRANTY; without even the implied warranty of

- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

- *    GNU General Public License for more details.

- *

- *    You should have received a copy of the GNU General Public License

- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.

- */

-

-

-// for persistence algorithm

-#include <gudhi/reader_utils.h>

-#include <gudhi/Bitmap_cubical_complex.h>

-#include <gudhi/Persistent_cohomology.h>

-

-#include <boost/program_options.hpp>

-

-// standard stuff

-#include <iostream>

-#include <sstream>

-#include <vector>

-#include <cstdlib>

-#include <ctime>

-

-using namespace Gudhi;

-using namespace Gudhi::persistent_cohomology;

-using namespace std;

-

-int main(int argc, char** argv) {

-  srand(time(0));

-

-  cout << "This program computes persistent homology, by using Bitmap_cubical_complex class, of cubical complexes. "

-      "The first parameter of the program is the dimension D of the cubical complex. The next D parameters are number "

-      "of top dimensional cubes in each dimension of the cubical complex. The program will create random cubical "

-      "complex of that sizes and compute persistent homology of it." << endl;

-

-  int p = 2;

-  double min_persistence = 0;

-

-  size_t dimensionOfBitmap = (size_t) atoi(argv[1]);

-  std::vector< unsigned > sizes;

-  size_t multipliers = 1;

-  for (size_t dim = 0; dim != dimensionOfBitmap; ++dim) {

-    unsigned sizeInThisDimension = (unsigned) atoi(argv[2 + dim]);

-    sizes.push_back(sizeInThisDimension);

-    multipliers *= sizeInThisDimension;

-  }

-

-  std::vector< double > data;

-  for (size_t i = 0; i != multipliers; ++i) {

-    data.push_back(rand() / (double) RAND_MAX);

-  }

-

-  Bitmap_cubical_complex<double> b(sizes, data);

-

-  // Compute the persistence diagram of the complex

-  persistent_cohomology::Persistent_cohomology< Bitmap_cubical_complex<double>, Field_Zp > pcoh(b);

-  pcoh.init_coefficients(p);  // initializes the coefficient field for homology

-  pcoh.compute_persistent_cohomology(min_persistence);

-

-  stringstream ss;

-  ss << "randomComplex_persistence";

-  std::ofstream out((char*) ss.str().c_str());

-  pcoh.output_diagram(out);

-  out.close();

-

-  return 0;

-}

+ /*    This file is part of the Gudhi Library. The Gudhi library
+  *    (Geometric Understanding in Higher Dimensions) is a generic C++
+  *    library for computational topology.
+  *
+  *    Author(s):       Pawel Dlotko
+  *
+  *    Copyright (C) 2015  INRIA Saclay (France)
+  *
+  *    This program is free software: you can redistribute it and/or modify
+  *    it under the terms of the GNU General Public License as published by
+  *    the Free Software Foundation, either version 3 of the License, or
+  *    (at your option) any later version.
+  *
+  *    This program is distributed in the hope that it will be useful,
+  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  *    GNU General Public License for more details.
+  *
+  *    You should have received a copy of the GNU General Public License
+  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+  */
+
+
+//for persistence algorithm
+#include <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+
+using namespace Gudhi;
+using namespace Gudhi::Cubical_complex;
+using namespace Gudhi::persistent_cohomology;
+
+//standard stuff
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+using namespace std;
+
+int main( int argc , char** argv )
+{
+    srand( time(0) );
+
+    cout << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical complexes. \
+The first parameter of the program is the dimension D of the bitmap. \
+The next D parameters are number of top dimensional cubes in each dimension of the bitmap.\
+The program will create random cubical complex of that sizes and compute persistent homology of it." << endl;
+
+    int p = 2;
+    double min_persistence = 0;
+
+    if ( argc < 3 )
+    {
+        cerr << "Wrong number of parameters, the program will now terminate\n";
+        return 1;
+    }
+
+    size_t dimensionOfBitmap = (size_t)atoi( argv[1] );
+    std::vector< unsigned > sizes;
+    size_t multipliers = 1;
+    for ( size_t dim = 0 ; dim != dimensionOfBitmap ; ++dim )
+    {
+        unsigned sizeInThisDimension = (unsigned)atoi( argv[2+dim] );
+        sizes.push_back( sizeInThisDimension );
+        multipliers *= sizeInThisDimension;
+    }
+
+    std::vector< double > data;
+    for ( size_t i = 0 ; i != multipliers ; ++i )
+    {
+        data.push_back( rand()/(double)RAND_MAX );
+    }
+
+
+
+    Bitmap_cubical_complex<double> b( sizes , data );
+
+
+    // Compute the persistence diagram of the complex
+    persistent_cohomology::Persistent_cohomology< Bitmap_cubical_complex<double>, Field_Zp > pcoh(b);
+    pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
+    pcoh.compute_persistent_cohomology( min_persistence );
+
+
+    stringstream ss;
+    ss << "randomComplex_persistence";
+    std::ofstream out((char*)ss.str().c_str());
+    pcoh.output_diagram(out);
+    out.close();
+
+    return 0;
+}
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 2f8cb0a3..f2c753d9 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
@@ -1,724 +1,719 @@
-/*    This file is part of the Gudhi Library. The Gudhi library

- *    (Geometric Understanding in Higher Dimensions) is a generic C++

- *    library for computational topology.

- *

- *    Author(s):       Pawel Dlotko

- *

- *    Copyright (C) 2015  INRIA Sophia-Saclay (France)

- *

- *    This program is free software: you can redistribute it and/or modify

- *    it under the terms of the GNU General Public License as published by

- *    the Free Software Foundation, either version 3 of the License, or

- *    (at your option) any later version.

- *

- *    This program is distributed in the hope that it will be useful,

- *    but WITHOUT ANY WARRANTY; without even the implied warranty of

- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

- *    GNU General Public License for more details.

- *

- *    You should have received a copy of the GNU General Public License

- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.

- */

-

-#ifndef BITMAP_CUBICAL_COMPLEX_H_

-#define BITMAP_CUBICAL_COMPLEX_H_

-

-#include <gudhi/Bitmap_cubical_complex_base.h>

-

-#include <limits>

-#include <utility>  // for pair

-#include <algorithm>  // for sort

-#include <vector>  // for vector

-

-// global variable, was used just for debugging.

-bool globalDbg = false;

-

-template <typename T = double>

-class Bitmap_cubical_complex : public Bitmap_cubical_complex_base<T> {

- public:

-  //******************************************************************************************************************//

-  // Typedefs and typenames

-  //******************************************************************************************************************//

-  friend class Simplex_handle;

-  typedef size_t Simplex_key;

-  typedef T Filtration_value;

-

-

-  //******************************************************************************************************************//

-  // Simplex handle class

-  //******************************************************************************************************************//

-

-  /**

-   * Handle of a cell, required for compatibility with the function to compute persistence in Gudhi. Elements of this

-   * class are: the pointer to the bitmap B in which the considered cell is together with a position of this cell in B.

-   * Given this data, one can get all the information about the considered cell.

-   **/

-  class Simplex_handle {

-   public:

-    Simplex_handle() {

-      if (globalDbg) {

-        std::cerr << "Simplex_handle()\n";

-      }

-      this->b = 0;

-      this->position = 0;

-    }

-

-    Simplex_handle(Bitmap_cubical_complex<T>* b) {

-      if (globalDbg) {

-        std::cerr << "Simplex_handle(Bitmap_cubical_complex<T>* b)\n";

-      }

-      this->b = b;

-      this->position = 0;

-    }

-

-    Simplex_handle(const Simplex_handle& org) : b(org.b) {

-      if (globalDbg) {

-        std::cerr << "Simplex_handle( const Simplex_handle& org )\n";

-      }

-      this->position = org.position;

-    }

-

-    Simplex_handle& operator=(const Simplex_handle& rhs) {

-      if (globalDbg) {

-        std::cerr << "Simplex_handle operator =  \n";

-      }

-      this->position = rhs.position;

-      this->b = rhs.b;

-      return *this;

-    }

-

-    Simplex_handle(Bitmap_cubical_complex<T>* b, Simplex_key position) {

-      if (globalDbg) {

-        std::cerr << "Simplex_handle(Bitmap_cubical_complex<T>* b , Simplex_key position)\n";

-        std::cerr << "Position : " << position << std::endl;

-      }

-      this->b = b;

-      this->position = position;

-    }

-    friend class Bitmap_cubical_complex<T>;

-   private:

-    Bitmap_cubical_complex<T>* b;

-    Simplex_key position;

-    // Assumption -- this field always keep the REAL position of simplex in the bitmap, no matter what keys have been.

-    // to deal with the keys, the class Bitmap_cubical_complex have extra vectors: keyAssociatedToSimplex and

-    // simplexAssociatedToKey that allow to move between actual cell and the key assigned to it.

-  };

-

-

-  //******************************************************************************************************************//

-  // Constructors

-  //******************************************************************************************************************//

-  // Over here we need to definie various input types. I am proposing the following ones:

-  // Perseus style

-  // TODO(Pawel Dlotko): H5 files?

-  // TODO(Pawel Dlotko): binary files with little endiangs / big endians?

-  // TODO(Pawel Dlotko): constructor from a vector of elements of a type T.

-

-  /**

-   * Constructor form a Perseus-style file.

-   **/

-  Bitmap_cubical_complex(char* perseusStyleFile) : Bitmap_cubical_complex_base<T>(perseusStyleFile) {

-    if (globalDbg) {

-      std::cerr << "Bitmap_cubical_complex( char* perseusStyleFile )\n";

-    }

-    std::vector< size_t > keyAssociatedToSimplex(this->totalNumberOfCells + 1);

-    std::vector< size_t > simplexAssociatedToKey(this->totalNumberOfCells + 1);

-

-    for (size_t i = 0; i != this->totalNumberOfCells; ++i) {

-      keyAssociatedToSimplex[i] = simplexAssociatedToKey[i] = i;

-    }

-    this->keyAssociatedToSimplex = keyAssociatedToSimplex;

-    this->simplexAssociatedToKey = simplexAssociatedToKey;

-    // 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->initializeElementsOrderedAccordingToFiltration();

-  }

-

-  /**

-   * 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(std::vector<unsigned> dimensions, std::vector<T> topDimensionalCells)

-  : Bitmap_cubical_complex_base<T>(dimensions, topDimensionalCells) {

-    std::vector< size_t > keyAssociatedToSimplex(this->totalNumberOfCells + 1);

-    std::vector< size_t > simplexAssociatedToKey(this->totalNumberOfCells + 1);

-

-    for (size_t i = 0; i != this->totalNumberOfCells; ++i) {

-      keyAssociatedToSimplex[i] = simplexAssociatedToKey[i] = i;

-    }

-    this->keyAssociatedToSimplex = keyAssociatedToSimplex;

-    this->simplexAssociatedToKey = simplexAssociatedToKey;

-    // 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->initializeElementsOrderedAccordingToFiltration();

-  }

-

-  //******************************************************************************************************************//

-  // Other 'easy' functions

-  //******************************************************************************************************************//

-

-  /**

-   * Returns number of all cubes in the complex.

-   **/

-  size_t num_simplices()const {

-    return this->totalNumberOfCells;

-  }

-

-  /**

-   * Returns a Simplex_handle to a cube that do not exist in this complex.

-   **/

-  Simplex_handle null_simplex() {

-    return Simplex_handle(this, this->data.size());

-  }

-

-  /**

-   * Returns dimension of the complex.

-   **/

-  size_t dimension() {

-    return this->sizes.size();

-  }

-

-  /**

-   * Return dimension of a cell pointed by the Simplex_handle.

-   **/

-  size_t dimension(const Simplex_handle& sh) {

-    if (globalDbg) {

-      std::cerr << "int dimension(const Simplex_handle& sh)\n";

-    }

-    if (sh.position != this->data.size()) return sh.b->get_dimension_of_a_cell(sh.position);

-    return std::numeric_limits<size_t>::max();

-  }

-

-  /**

-   * Return the filtration of a cell pointed by the Simplex_handle.

-   **/

-  T filtration(const Simplex_handle& sh) {

-    if (globalDbg) {

-      std::cerr << "T filtration(const Simplex_handle& sh)\n";

-    }

-    // Returns the filtration value of a simplex.

-    if (sh.position != this->data.size()) return sh.b->data[ sh.position ];

-    return INT_MAX;

-  }

-

-  /**

-   * Return a key which is not a key of any cube in the considered data structure.

-   **/

-  Simplex_key null_key() {

-    if (globalDbg) {

-      std::cerr << "Simplex_key null_key()\n";

-    }

-    return this->data.size();

-  }

-

-  /**

-   * Return the key of a cube pointed by the Simplex_handle.

-   **/

-  Simplex_key key(const Simplex_handle& sh) {

-    if (globalDbg) {

-      std::cerr << "Simplex_key key(const Simplex_handle& sh)\n";

-    }

-    return sh.b->keyAssociatedToSimplex[ sh.position ];

-  }

-

-  /**

-   * 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";

-    }

-    return Simplex_handle(this, this->simplexAssociatedToKey[ key ]);

-  }

-

-  /**

-   * 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";

-    }

-    this->keyAssociatedToSimplex[sh.position] = key;

-    this->simplexAssociatedToKey[key] = sh.position;

-  }

-

-  /**

-   * Function called from a constructor. It is needed for Filtration_simplex_iterator to work.

-   **/

-  void initializeElementsOrderedAccordingToFiltration();

-

-

-

-  //******************************************************************************************************************//

-  // Iterators

-  //******************************************************************************************************************//

-

-  /**

-   * Boundary_simplex_iterator class allows iteration on boundary of each cube.

-   **/

-  class Boundary_simplex_range;

-

-  class Boundary_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > {

-  // Iterator on the simplices belonging to the boundary of a simplex.

-  // value_type must be 'Simplex_handle'.

-   public:

-    Boundary_simplex_iterator(Simplex_handle& sh) : sh(sh) {

-      if (globalDbg) {

-        std::cerr << "Boundary_simplex_iterator( Simplex_handle& sh )\n";

-      }

-      this->position = 0;

-      this->boundaryElements = this->sh.b->get_boundary_of_a_cell(this->sh.position);

-    }

-

-    Boundary_simplex_iterator operator++() {

-      if (globalDbg) {

-        std::cerr << "Boundary_simplex_iterator operator++()\n";

-      }

-      ++this->position;

-      return *this;

-    }

-

-    Boundary_simplex_iterator operator++(int) {

-      Boundary_simplex_iterator result = *this;

-      ++(*this);

-      return result;

-    }

-

-    Boundary_simplex_iterator operator=(const Boundary_simplex_iterator& rhs) {

-      if (globalDbg) {

-        std::cerr << "Boundary_simplex_iterator operator =\n";

-      }

-      this->sh = rhs.sh;

-      this->boundaryElements.clear();

-      this->boundaryElementsinsert(this->boundaryElements.end(),

-                                   rhs.boundaryElements.begin(), rhs.boundaryElements.end());

-    }

-

-    bool operator==(const Boundary_simplex_iterator& rhs) {

-      if (globalDbg) {

-        std::cerr << "bool operator ==\n";

-      }

-      if (this->position == rhs.position) {

-        if (this->boundaryElements.size() != rhs.boundaryElements.size())return false;

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

-          if (this->boundaryElements[i] != rhs.boundaryElements[i])return false;

-        }

-        return true;

-      }

-      return false;

-    }

-

-    bool operator!=(const Boundary_simplex_iterator& rhs) {

-      if (globalDbg) {

-        std::cerr << "bool operator != \n";

-      }

-      return !(*this == rhs);

-    }

-

-    Simplex_handle operator*() {

-      if (globalDbg) {

-        std::cerr << "Simplex_handle operator*\n";

-      }

-      return Simplex_handle(this->sh.b, this->boundaryElements[this->position]);

-    }

-

-    friend class Boundary_simplex_range;

-   private:

-    Simplex_handle sh;

-    std::vector< size_t > boundaryElements;

-    size_t position;

-  };

-

-  /**

-   * Boundary_simplex_range class provides ranges for boundary iterators.

-   **/

-  class Boundary_simplex_range {

-  // Range giving access to the simplices in the boundary of a simplex.

-  // .begin() and .end() return type Boundary_simplex_iterator.

-   public:

-    Boundary_simplex_range(const Simplex_handle& sh) : sh(sh) { }

-

-    Boundary_simplex_iterator begin() {

-      if (globalDbg) {

-        std::cerr << "Boundary_simplex_iterator begin\n";

-      }

-      Boundary_simplex_iterator it(this->sh);

-      return it;

-    }

-

-    Boundary_simplex_iterator end() {

-      if (globalDbg) {

-        std::cerr << "Boundary_simplex_iterator end()\n";

-      }

-      Boundary_simplex_iterator it(this->sh);

-      it.position = it.boundaryElements.size();

-      return it;

-    }

-

-   private:

-    Simplex_handle sh;

-  };

-

-

-  /**

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

-    }

-

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

-    }

-

-    bool operator==(const Filtration_simplex_iterator& rhs) {

-      if (globalDbg) {

-        std::cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n";

-      }

-      if (this->position == rhs.position) {

-        return true;

-      }

-      return false;

-    }

-

-    bool operator!=(const Filtration_simplex_iterator& rhs) {

-      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 Simplex_handle(this->b, this->b->elementsOrderedAccordingToFiltration[ 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:

-    Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) { }

-

-    Filtration_simplex_iterator begin() {

-      if (globalDbg) {

-        std::cerr << "Filtration_simplex_iterator begin() \n";

-      }

-      return Filtration_simplex_iterator(this->b);

-    }

-

-    Filtration_simplex_iterator end() {

-      if (globalDbg) {

-        std::cerr << "Filtration_simplex_iterator end()\n";

-      }

-      Filtration_simplex_iterator it(this->b);

-      it.position = this->b->elementsOrderedAccordingToFiltration.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) {

-    if (globalDbg) {

-      std::cerr << "Boundary_simplex_range boundary_simplex_range(Simplex_handle& sh)\n";

-    }

-    // Returns a range giving access to all simplices of the boundary of a simplex, i.e. the set of

-    // codimension 1 subsimplices of the Simplex.

-    return Boundary_simplex_range(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(Pawel Dlotko): 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.position);

-    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 for which this method was called have less than two elements in the boundary.");

-    return std::make_pair(Simplex_handle(this, bdry[0]), Simplex_handle(this, 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), 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;

-    }

-

-    bool operator==(const Skeleton_simplex_iterator& rhs) {

-      if (globalDbg) {

-        std::cerr << "bool operator ==\n";

-      }

-      if (this->position == rhs.position) {

-        return true;

-      }

-      return false;

-    }

-

-    bool operator!=(const Skeleton_simplex_iterator& rhs) {

-      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 Simplex_handle(this->b, this->position);

-    }

-

-    friend class Skeleton_simplex_range;

-   private:

-    Bitmap_cubical_complex<T>* b;

-    size_t position;

-    int 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:

-    Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, int 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;

-    int dimension;

-  };

-

-  /**

-   * Function needed for compatibility with Gudhi. Not useful for other purposes.

-   **/

-  Skeleton_simplex_range skeleton_simplex_range(int dimension) {

-    if (globalDbg) {

-      std::cerr << "Skeleton_simplex_range skeleton_simplex_range( int dimension )\n";

-    }

-    return Skeleton_simplex_range(this, dimension);

-  }

-

-

-

-  //******************************************************************************************************************//

-  // functions used for debugging:

-

-  /**

-   * Function used for debugging purposes.

-   **/

-  void printKeyAssociatedToSimplex() {

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

-      std::cerr << i << " -> " << this->simplexAssociatedToKey[i] << std::endl;

-    }

-  }

-

-  /**

-   * Function used for debugging purposes.

-   **/

-  size_t printRealPosition(const Simplex_handle& sh) {

-    return sh.position;

-  }

-

- private:

-  std::vector< size_t > keyAssociatedToSimplex;

-  std::vector< size_t > simplexAssociatedToKey;

-  // needed by Filtration_simplex_iterator. If this iterator is not used, this field is not initialized.

-  std::vector< size_t > elementsOrderedAccordingToFiltration;

-};

-

-template <typename T>

-bool compareElementsForElementsOrderedAccordingToFiltration(const std::pair< size_t,

-                                                            std::pair< T, char > >& f,

-                                                            const std::pair< size_t,

-                                                            std::pair< T, char > >& s) {

-  if (globalDbg) {

-    std::cerr << "ompareElementsForElementsOrderedAccordingToFiltration\n";

-  }

-  if (f.second.first < s.second.first) {

-    return true;

-  } else {

-    if (f.second.first > s.second.first) {

-      return false;

-    } else {

-      // in this case f.second.first == s.second.first, and we use dimension to compare:

-      if (f.second.second < s.second.second) {

-        return true;

-      } else {

-        if (f.second.second > s.second.second) {

-          return false;

-        } else {

-          // in this case, both the filtration value and the dimensions for those cells are the same.

-          // Since it may be nice to have a stable sorting procedure, in this case, we compare positions in the bitmap:

-          return ( f.first < s.first);

-        }

-      }

-    }

-  }

-}

-

-template <typename T>

-void Bitmap_cubical_complex<T>::initializeElementsOrderedAccordingToFiltration() {

-  if (globalDbg) {

-    std::cerr << "void Bitmap_cubical_complex<T>::initializeElementsOrderedAccordingToFiltration() \n";

-  }

-  // ( position , (filtration , dimension) )

-  std::vector< std::pair< size_t, std::pair< T, char > > > dataOfElementsFromBitmap(this->data.size());

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

-    // TODO(Pawel Dlotko): This can be optimized by having a counter here. We do not need to re-compute the dimension

-    // for every cell from scratch

-    dataOfElementsFromBitmap[i] = std::make_pair(i, std::make_pair(this->data[i], this->get_dimension_of_a_cell(i)));

-  }

-  std::sort(dataOfElementsFromBitmap.begin(), dataOfElementsFromBitmap.end(),

-            compareElementsForElementsOrderedAccordingToFiltration<T>);

-

-  // Elements of bitmap ordered according to filtration then according to dimension then according to position in bitmap

-  std::vector< size_t > elements_of_bitmap_ordered(this->data.size());

-  for (size_t i = 0; i != dataOfElementsFromBitmap.size(); ++i) {

-    elements_of_bitmap_ordered[i] = dataOfElementsFromBitmap[i].first;

-  }

-  this->elementsOrderedAccordingToFiltration = elements_of_bitmap_ordered;

-}

-

-

-//****************************************************************************************************************//

-//****************************************************************************************************************//

-//****************************************************************************************************************//

-//****************************************************************************************************************//

-

-#endif  // BITMAP_CUBICAL_COMPLEX_H_

+	/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    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 <limits>
+#include "Bitmap_cubical_complex_base.h"
+
+
+
+namespace Gudhi
+{
+
+namespace Cubical_complex
+{
+
+//global variable, was used just for debugging.
+const bool globalDbg = false;
+
+template <typename T = double>
+class Bitmap_cubical_complex : public Bitmap_cubical_complex_base<T>
+{
+public:
+//*********************************************//
+//Typedefs and typenames
+//*********************************************//
+    friend class Simplex_handle;
+    typedef size_t Simplex_key;
+    typedef T Filtration_value;
+
+
+//*********************************************//
+//Simplex handle class
+//*********************************************//
+    /**
+    * Handle of a cell, required for compatibility with the function to compute persistence in Gudhi. 
+    * Elements of this class are: the pointer to the bitmap B in which the considered cell is
+    * together with a position of this cell in B. Given this data, 
+    * one can get all the information about the considered cell.
+    **/
+    class Simplex_handle
+    {
+    public:
+        Simplex_handle()
+        {
+            if ( globalDbg ){cerr << "Simplex_handle()\n";}
+            this->b = 0;
+            this->position = 0;
+        }
+
+        Simplex_handle(Bitmap_cubical_complex<T>* b)
+        {
+            if ( globalDbg )
+	     {
+		 cerr << "Simplex_handle(Bitmap_cubical_complex<T>* b)\n";
+	     }
+            this->b = b;
+            this->position = 0;
+        }
+
+        //Simplex_handle( const Simplex_handle& org ):b(org.b)
+        //{
+        //    if ( globalDbg ){cerr << "Simplex_handle( const Simplex_handle& org )\n";}
+        //    this->position = org.position;
+        //}
+
+        Simplex_handle operator = ( const Simplex_handle& rhs )
+        {
+            if ( globalDbg ){cerr << "Simplex_handle operator =  \n";}
+            this->position = rhs.position;
+            this->b = rhs.b;
+	     return *this;
+        }
+
+        Simplex_handle(Bitmap_cubical_complex<T>* b , Simplex_key position)
+        {
+            if ( globalDbg )
+            {
+                 cerr << "Simplex_handle(Bitmap_cubical_complex<T>* b , Simplex_key position)\n";
+                 cerr << "Position : " << position << endl;
+            }
+            this->b = b;
+            this->position = position;
+        }
+        friend class Bitmap_cubical_complex<T>;
+    private:
+        Bitmap_cubical_complex<T>* b;
+        Simplex_key position;
+        //Assumption -- field above always keep the REAL position of simplex in the bitmap, 
+        //no matter what keys have been.
+	 //to deal with the keys, the class Bitmap_cubical_complex have extra vectors: key_associated_to_simplex and 
+        //simplex_associated_to_key that allow to move between actual cell and the key assigned to it.
+    };
+
+
+//*********************************************//
+//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 ):
+    Bitmap_cubical_complex_base<T>(perseus_style_file),key_associated_to_simplex(this->total_number_of_cells+1),
+    simplex_associated_to_key(this->total_number_of_cells+1)
+    {
+        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] = this->simplex_associated_to_key[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_elements_ordered_according_to_filtration();
+    }
+
+
+    /**
+    * 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( std::vector<unsigned>& dimensions , std::vector<T>& top_dimensional_cells ):
+    Bitmap_cubical_complex_base<T>(dimensions,top_dimensional_cells),
+    key_associated_to_simplex(this->total_number_of_cells+1),
+    simplex_associated_to_key(this->total_number_of_cells+1)
+    {
+        for ( size_t i = 0 ; i != this->total_number_of_cells ; ++i )
+        {
+            this->key_associated_to_simplex[i] = this->simplex_associated_to_key[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_elements_ordered_according_to_filtration();
+    }
+
+//*********************************************//
+//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.
+    **/
+    Simplex_handle null_simplex()
+    {
+        return Simplex_handle(this,this->data.size());
+    }
+
+    /**
+    * Returns dimension of the complex.
+    **/
+    size_t dimension()
+    {
+        return this->sizes.size();
+    }
+
+    /**
+    * Return dimension of a cell pointed by the Simplex_handle.
+    **/
+    unsigned dimension(const Simplex_handle& sh)
+    {
+        if ( globalDbg ){cerr << "unsigned dimension(const Simplex_handle& sh)\n";}
+        if ( sh.position != this->data.size() ) return sh.b->get_dimension_of_a_cell( sh.position );
+        return -1;
+    }
+
+    /**
+    * Return the filtration of a cell pointed by the Simplex_handle.
+    **/
+    T filtration(const Simplex_handle& sh)
+    {
+        if ( globalDbg ){cerr << "T filtration(const Simplex_handle& sh)\n";}
+        //Returns the filtration value of a simplex.
+        if ( sh.position != this->data.size() ) return sh.b->data[ sh.position ];
+        return std::numeric_limits<int>::max();
+    }
+
+    /**
+    * Return a key which is not a key of any cube in the considered data structure.
+    **/
+    Simplex_key null_key()
+    {
+        if ( globalDbg ){cerr << "Simplex_key null_key()\n";}
+        return this->data.size();
+    }
+
+    /**
+    * Return the key of a cube pointed by the Simplex_handle.
+    **/
+    Simplex_key key(const Simplex_handle& sh)
+    {
+        if ( globalDbg ){cerr << "Simplex_key key(const Simplex_handle& sh)\n";}
+        return sh.b->key_associated_to_simplex[ sh.position ];
+    }
+
+    /**
+    * 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";}
+        return Simplex_handle( this , this->simplex_associated_to_key[ key ] );
+    }
+
+    /**
+    * 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";}
+        this->key_associated_to_simplex[sh.position] = key;
+        this->simplex_associated_to_key[key] = sh.position;
+    }
+
+    /**
+    * Function called from a constructor. It is needed for Filtration_simplex_iterator to work.
+    **/
+    void initialize_elements_ordered_according_to_filtration();
+
+
+
+//*********************************************//
+//Iterators
+//*********************************************//
+
+    /**
+    * Boundary_simplex_iterator class allows iteration on boundary of each cube.
+    **/
+    class Boundary_simplex_range;
+    class Boundary_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle >
+    {
+        //Iterator on the simplices belonging to the boundary of a simplex.
+        //value_type must be 'Simplex_handle'.
+        public:
+            Boundary_simplex_iterator( Simplex_handle& sh ):sh(sh)
+            {
+                if ( globalDbg ){cerr << "Boundary_simplex_iterator( Simplex_handle& sh )\n";}
+                this->position = 0;
+                this->boundary_elements = this->sh.b->get_boundary_of_a_cell( this->sh.position );
+            }
+            Boundary_simplex_iterator operator++()
+            {
+                if ( globalDbg ){cerr << "Boundary_simplex_iterator operator++()\n";}
+                ++this->position;
+                return *this;
+            }
+            Boundary_simplex_iterator operator++(int)
+            {
+                Boundary_simplex_iterator result = *this;
+                ++(*this);
+                return result;
+            }
+            Boundary_simplex_iterator operator =( const Boundary_simplex_iterator& rhs )
+            {
+                if ( globalDbg ){cerr << "Boundary_simplex_iterator operator =\n";}
+                this->sh = rhs.sh;
+                this->boundary_elements.clear();
+                this->boundary_elementsinsert
+		  (this->boundary_elements.end(), rhs.boundary_elements.begin(), rhs.boundary_elements.end());
+            }
+            bool operator == ( const Boundary_simplex_iterator& rhs )
+            {
+                if ( globalDbg ){cerr << "bool operator ==\n";}
+                if ( this->position == rhs.position )
+                {
+                    if ( this->boundary_elements.size() != rhs.boundary_elements.size() )return false;
+                    for ( size_t i = 0 ; i != this->boundary_elements.size() ; ++i )
+                    {
+                        if ( this->boundary_elements[i] != rhs.boundary_elements[i] )return false;
+                    }
+                    return true;
+                }
+                return false;
+            }
+
+            bool operator != ( const Boundary_simplex_iterator& rhs )
+            {
+                if ( globalDbg ){cerr << "bool operator != \n";}
+                return !(*this == rhs);
+            }
+            Simplex_handle operator*()
+            {
+                if ( globalDbg ){cerr << "Simplex_handle operator*\n";}
+                return Simplex_handle( this->sh.b , this->boundary_elements[this->position]  );
+            }
+
+            friend class Boundary_simplex_range;
+    private:
+            Simplex_handle sh;
+            std::vector< size_t > boundary_elements;
+            size_t position;
+    };
+
+
+    /**
+    * Boundary_simplex_range class provides ranges for boundary iterators.
+    **/
+    class Boundary_simplex_range
+    {
+        //Range giving access to the simplices in the boundary of a simplex.
+        //.begin() and .end() return type Boundary_simplex_iterator.
+        public:
+            Boundary_simplex_range(const Simplex_handle& sh):sh(sh){};
+            Boundary_simplex_iterator begin()
+            {
+                if ( globalDbg ){cerr << "Boundary_simplex_iterator begin\n";}
+                Boundary_simplex_iterator it( this->sh );
+                return it;
+            }
+            Boundary_simplex_iterator end()
+            {
+                if ( globalDbg ){cerr << "Boundary_simplex_iterator end()\n";}
+                Boundary_simplex_iterator it( this->sh );
+                it.position = it.boundary_elements.size();
+                return it;
+            }
+        private:
+            Simplex_handle sh;
+    };
+
+
+    /**
+    * 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 )
+            {
+                if ( globalDbg ){cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n";}
+                if ( this->position == rhs.position )
+                {
+                    return true;
+                }
+                return false;
+            }
+
+            bool operator != ( const Filtration_simplex_iterator& rhs )
+            {
+                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 Simplex_handle( this->b , this->b->elements_ordered_according_to_filtration[ 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:
+            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->elements_ordered_according_to_filtration.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)
+    {
+        if ( globalDbg ){cerr << "Boundary_simplex_range boundary_simplex_range(Simplex_handle& sh)\n";}
+        //Returns a range giving access to all simplices of the boundary of a simplex, 
+        //i.e. the set of codimension 1 subsimplices of the Simplex.
+        return Boundary_simplex_range(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.position );
+        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( Simplex_handle(this,bdry[0]) , Simplex_handle(this,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 )
+            {
+                if ( globalDbg ){cerr << "bool operator ==\n";}
+                if ( this->position == rhs.position )
+                {
+                    return true;
+                }
+                return false;
+            }
+
+            bool operator != ( const Skeleton_simplex_iterator& rhs )
+            {
+                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 Simplex_handle( this->b , 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:
+            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 );
+    }
+
+
+
+//*********************************************//
+//functions used for debugging:
+    /**
+    * Function used for debugging purposes.
+    **/
+    //void printkey_associated_to_simplex()
+    //{
+    //    for ( size_t i = 0 ; i != this->data.size() ; ++i )
+    //    {
+    //        cerr << i << " -> " << this->simplex_associated_to_key[i] << endl;
+    //    }
+    //}
+
+    /**
+    * Function used for debugging purposes.
+    **/
+    size_t printRealPosition( const Simplex_handle& sh )
+    {
+        return sh.position;
+    }
+
+private:
+    std::vector< size_t > key_associated_to_simplex;
+    std::vector< size_t > simplex_associated_to_key;
+    std::vector< size_t > elements_ordered_according_to_filtration;
+    //filed above is needed by Filtration_simplex_iterator. If this iterator is not used, this field is not initialized.
+};//Bitmap_cubical_complex
+
+template <typename T>
+bool compare_elements_for_elements_ordered_according_to_filtration
+( const std::pair< size_t , std::pair< T , char > >& f , const std::pair< size_t , std::pair< T , char > >& s )
+{
+    if ( globalDbg ){cerr << "compare_elements_for_elements_ordered_according_to_filtration\n";}
+    if ( f.second.first < s.second.first )
+    {
+        return true;
+    }
+    else
+    {
+        if ( f.second.first > s.second.first )
+        {
+            return false;
+        }
+        else
+        {
+            //in this case f.second.first == s.second.first, and we use dimension to compare:
+            if ( f.second.second < s.second.second )
+            {
+                return true;
+            }
+            else
+            {
+                if ( f.second.second > s.second.second )
+                {
+                    return false;
+                }
+                else
+                {
+                    //in this case, both the filtration value and the dimensions for those cells are the same. 
+                    //Since it may be nice to have a stable sorting procedure, in this case, 
+                    //we compare positions in the bitmap:
+                    return ( f.first < s.first );
+                }
+            }
+        }
+    }
+}
+
+template <typename T>
+void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration()
+{
+    if ( globalDbg )
+    {
+	cerr << "void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration() \n";
+    }
+    //( position , (filtration , dimension) )
+    std::vector< std::pair< size_t , std::pair< T , char > > > data_of_elements_from_bitmap( this->data.size() );
+    for ( size_t i = 0 ; i != this->data.size() ; ++i )
+    {
+        //TODO -- this can be optimized by having a counter here. 
+        //We do not need to re-compute the dimension for every cell from scratch
+        data_of_elements_from_bitmap[i] = 
+        std::make_pair( i , std::make_pair( this->data[i] , this->get_dimension_of_a_cell(i) ) );
+    }
+    std::sort( data_of_elements_from_bitmap.begin() , 
+               data_of_elements_from_bitmap.end() , 
+               compare_elements_for_elements_ordered_according_to_filtration<T> );
+
+    std::vector< size_t > 
+    elements_ordered_according_to_filtration_then_to_dimension_then_to_position
+    ( this->data.size() );
+    for ( size_t i = 0 ; i != data_of_elements_from_bitmap.size() ; ++i )
+    {
+        elements_ordered_according_to_filtration_then_to_dimension_then_to_position[i] 
+        = data_of_elements_from_bitmap[i].first;
+    }
+    this->elements_ordered_according_to_filtration = 
+    elements_ordered_according_to_filtration_then_to_dimension_then_to_position;
+}
+
+
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+
+
+}
+}
\ No newline at end of file
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 d9c91832..2c2bd481 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
@@ -1,593 +1,759 @@
-/*    This file is part of the Gudhi Library. The Gudhi library

- *    (Geometric Understanding in Higher Dimensions) is a generic C++

- *    library for computational topology.

- *

- *    Author(s):       Pawel Dlotko

- *

- *    Copyright (C) 2015  INRIA Sophia-Saclay (France)

- *

- *    This program is free software: you can redistribute it and/or modify

- *    it under the terms of the GNU General Public License as published by

- *    the Free Software Foundation, either version 3 of the License, or

- *    (at your option) any later version.

- *

- *    This program is distributed in the hope that it will be useful,

- *    but WITHOUT ANY WARRANTY; without even the implied warranty of

- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

- *    GNU General Public License for more details.

- *

- *    You should have received a copy of the GNU General Public License

- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.

- */

-

-#ifndef BITMAP_CUBICAL_COMPLEX_BASE_H_

-#define BITMAP_CUBICAL_COMPLEX_BASE_H_

-

-#include <gudhi/counter.h>

-

-#include <iostream>

-#include <string>

-#include <vector>

-#include <cstdlib>

-#include <climits>

-#include <fstream>

-#include <algorithm>

-#include <iterator>

-

-/**

- * This is a class implementing a basic bitmap data structure to store cubical complexes. It implements only the most basic subroutines.

- * The idea of the bitmap is the following. Our aim is to have a memory efficient data structure to store d-dimensional cubical complex C being a cubical decomposition

- * of a rectangular region of a space. This is achieved by storing C as a vector of bits (this is where the name 'bitmap' came from). Each cell is represented by a single

- * bit (in case of black and white bitmaps, or by a single element of a type T (here T is a filtration type of a bitmap, typically a double). All the informations needed for homology and

- * persistent homology computations (like dimension of a cell, boundary and coboundary elements of a cell, are then obtained from the position of the element in C.

- */

-template <typename T>

-class Bitmap_cubical_complex_base {

- public:

-  /**

-   * 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 to 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(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(char* perseusStyleFile_);

-  /**

-   * 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(std::vector<unsigned> dimensions_, std::vector<T> topDimensionalCells_);

-

-  /**

-   * The functions get_boundary_of_a_cell, get_coboundary_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.

-   */

-  inline std::vector< size_t > get_boundary_of_a_cell(size_t cell_);

-  /**

-   * 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.

-   **/

-  inline std::vector< size_t > get_coboundary_of_a_cell(size_t cell_);

-  /**

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

-  /**

-   * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position.

-   **/

-  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 topDimensionalCellsIterator).

-   **/

-  void impose_lower_star_filtration();  // assume that top dimensional cells are already set.

-

-  /**

-   * Returns dimension of a complex.

-   **/

-  inline unsigned dimension() {

-    return sizes.size();

-  }

-

-  /**

-   * Returns number of all cubes in the data structure.

-   **/

-  inline unsigned size_of_bitmap() {

-    return this->data.size();

-  }

-

-  /**

-   * Writing to stream operator.

-   **/

-  template <typename K>

-  friend std::ostream& operator<<(std::ostream & os_, const Bitmap_cubical_complex_base<K>& b_);

-

-  // ITERATORS

-

-  /**

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

-   **/

-  typedef typename std::vector< T >::iterator all_cells_iterator;

-

-  all_cells_iterator all_cells_begin()const {

-    return this->data.begin();

-  }

-

-  all_cells_iterator all_cells_end()const {

-    return this->data.end();

-  }

-

-

-  typedef typename std::vector< T >::const_iterator all_cells_const_iterator;

-

-  all_cells_const_iterator all_cells_const_begin()const {

-    return this->data.begin();

-  }

-

-  all_cells_const_iterator all_cells_const_end()const {

-    return this->data.end();

-  }

-

-  /**

-   * 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, double > {

-   public:

-    Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b_) : b(b_) {

-      for (size_t i = 0; i != b_.dimension(); ++i) {

-        this->counter.push_back(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_) {

-      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_) {

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

-    }

-

-    size_t computeIndexInBitmap() {

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

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

-        std::cout << this->counter[i] << " ";

-      }

-    }

-    friend class Bitmap_cubical_complex_base;

-   protected:

-    std::vector< unsigned > counter;

-    Bitmap_cubical_complex_base& b;

-  };

-

-  Top_dimensional_cells_iterator top_dimensional_cells_begin() {

-    Top_dimensional_cells_iterator a(*this);

-    return a;

-  }

-

-  Top_dimensional_cells_iterator top_dimensional_cells_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;

-  }

-

-

-  //****************************************************************************************************************//

-  //****************************************************************************************************************//

-  //****************************************************************************************************************//

-  //****************************************************************************************************************//

-

-

-  //****************************************************************************************************************//

-  //****************************************************************************************************************//

-  //****************************************************************************************************************//

-  //****************************************************************************************************************//

-

- protected:

-  std::vector<unsigned> sizes;

-  std::vector<unsigned> multipliers;

-  std::vector<T> data;

-  size_t totalNumberOfCells;

-

-  void set_up_containers(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)+1;

-      multiplier *= 2 * sizes_[i] + 1;

-    }

-    // std::reverse( this->sizes.begin() , this->sizes.end() );

-    std::vector<T> data(multiplier);

-    std::fill(data.begin(), data.end(), INT_MAX);

-    this->totalNumberOfCells = multiplier;

-    this->data = data;

-  }

-

-  size_t compute_position_in_bitmap(std::vector< int > 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_) {

-    std::vector<unsigned> counter;

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

-  }

-

-  std::vector< size_t > generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions(std::vector< bool > directionsForPeriodicBCond_);

-};

-

-template <typename K>

-std::ostream& operator<<(std::ostream & out_, const Bitmap_cubical_complex_base<K>& b_) {

-  // for ( typename bitmap<K>::all_cells_const_iterator it = b.all_cells_const_begin() ;

-  // it != b.all_cells_const_end() ; ++it )

-  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(std::vector<unsigned> sizes_) {

-  this->set_up_containers(sizes_);

-}

-

-template <typename T>

-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(std::vector<unsigned> sizesInFollowingDirections_,

-                                                            std::vector<T> topDimensionalCells_) {

-  this->set_up_containers(sizesInFollowingDirections_);

-

-  size_t numberOfTopDimensionalElements = 1;

-  for (size_t i = 0; i != sizesInFollowingDirections_.size(); ++i) {

-    numberOfTopDimensionalElements *= sizesInFollowingDirections_[i];

-  }

-  if (numberOfTopDimensionalElements != topDimensionalCells_.size()) {

-    std::cerr << "Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizesInFollowingDirections_ , "

-        "std::vector<float> topDimensionalCells_ ). Number of top dimensional elements that follow from "

-        "sizesInFollowingDirections vector is different than the size of topDimensionalCells vector." << std::endl;

-    throw("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizesInFollowingDirections_ , "

-           "std::vector<float> topDimensionalCells_ ). Number of top dimensional elements that follow from "

-           "sizesInFollowingDirections vector is different than the size of topDimensionalCells vector.");

-  }

-

-  Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);

-  size_t index = 0;

-  for (it = this->top_dimensional_cells_begin(); it != this->top_dimensional_cells_end(); ++it) {

-    (*it) = topDimensionalCells_[index];

-    ++index;

-  }

-  this->impose_lower_star_filtration();

-}

-

-template <typename T>

-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(char* perseusStyleFile_) {

-  bool dbg = false;

-  std::ifstream inFiltration, inIds;

-  inFiltration.open(perseusStyleFile_);

-  unsigned dimensionOfData;

-  inFiltration >> dimensionOfData;

-

-  if (dbg) {

-    std::cerr << "dimensionOfData : " << dimensionOfData << std::endl;

-  }

-

-  std::vector<unsigned> sizes;

-  for (size_t i = 0; i != dimensionOfData; ++i) {

-    int sizeInThisDimension;

-    inFiltration >> sizeInThisDimension;

-    sizeInThisDimension = abs(sizeInThisDimension);

-    sizes.push_back(sizeInThisDimension);

-    if (dbg) {

-      std::cerr << "sizeInThisDimension : " << sizeInThisDimension << std::endl;

-    }

-  }

-  this->set_up_containers(sizes);

-

-  Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);

-  it = this->top_dimensional_cells_begin();

-

-  // TODO(Pawel Dlotko): Over here we also need to read id's of cell and put them to bitmapElement structure!

-  while (!inFiltration.eof()) {

-    double filtrationLevel;

-    inFiltration >> filtrationLevel;

-    if (dbg) {

-      std::cerr << "Cell of an index : " << it.computeIndexInBitmap() << " and dimension: " <<

-          this->get_dimension_of_a_cell(it.computeIndexInBitmap()) << " get the value : " <<

-          filtrationLevel << std::endl;

-    }

-    *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_) {

-  bool bdg = false;

-  // First of all, we need to take the list of coordinates in which the cell has nonzero length.

-  // We do it by using modified version to compute dimension of a cell:

-  std::vector< unsigned > dimensionsInWhichCellHasNonzeroLength;

-  unsigned dimension = 0;

-  size_t cell1 = cell_;

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

-    unsigned position = cell1 / multipliers[i - 1];

-    if (position % 2 == 1) {

-      dimensionsInWhichCellHasNonzeroLength.push_back(i - 1);

-      dimension++;

-    }

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

-  }

-

-  if (bdg) {

-    std::cerr << "dimensionsInWhichCellHasNonzeroLength : \n";

-    for (size_t i = 0; i != dimensionsInWhichCellHasNonzeroLength.size(); ++i) {

-      std::cerr << dimensionsInWhichCellHasNonzeroLength[i] << std::endl;

-    }

-    getchar();

-  }

-

-  std::vector< size_t > boundaryElements;

-  if (dimensionsInWhichCellHasNonzeroLength.size() == 0)return boundaryElements;

-  for (size_t i = 0; i != dimensionsInWhichCellHasNonzeroLength.size(); ++i) {

-    boundaryElements.push_back(cell_ - multipliers[ dimensionsInWhichCellHasNonzeroLength[i] ]);

-    boundaryElements.push_back(cell_ + multipliers[ dimensionsInWhichCellHasNonzeroLength[i] ]);

-

-    if (bdg) std::cerr << "multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " <<

-        multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << std::endl;

-    if (bdg) std::cerr << "cell_ - multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " <<

-        cell_ - multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << std::endl;

-    if (bdg) std::cerr << "cell_ + multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " <<

-        cell_ + multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << std::endl;

-  }

-  return boundaryElements;

-}

-

-template <typename T>

-std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell_) {

-  bool bdg = false;

-  // First of all, we need to take the list of coordinates in which the cell has nonzero length.

-  // We do it by using modified version to compute dimension of a cell:

-  std::vector< unsigned > dimensionsInWhichCellHasZeroLength;

-  unsigned dimension = 0;

-  size_t cell1 = cell_;

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

-    unsigned position = cell1 / multipliers[i - 1];

-    if (position % 2 == 0) {

-      dimensionsInWhichCellHasZeroLength.push_back(i - 1);

-      dimension++;

-    }

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

-  }

-

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

-  // reverse(counter.begin() , counter.end());

-

-  if (bdg) {

-    std::cerr << "dimensionsInWhichCellHasZeroLength : \n";

-    for (size_t i = 0; i != dimensionsInWhichCellHasZeroLength.size(); ++i) {

-      std::cerr << dimensionsInWhichCellHasZeroLength[i] << std::endl;

-    }

-    std::cerr << "\n counter : " << std::endl;

-    for (size_t i = 0; i != counter.size(); ++i) {

-      std::cerr << counter[i] << std::endl;

-    }

-    getchar();

-  }

-

-  std::vector< size_t > coboundaryElements;

-  if (dimensionsInWhichCellHasZeroLength.size() == 0)return coboundaryElements;

-  for (size_t i = 0; i != dimensionsInWhichCellHasZeroLength.size(); ++i) {

-    if (bdg) {

-      std::cerr << "Dimension : " << i << std::endl;

-      if (counter[dimensionsInWhichCellHasZeroLength[i]] == 0) {

-        std::cerr << "In dimension : " << i <<

-            " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n";

-      }

-      if (counter[dimensionsInWhichCellHasZeroLength[i]] == 2 * this->sizes[dimensionsInWhichCellHasZeroLength[i]]) {

-        std::cerr << "In dimension : " << i <<

-            " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n";

-      }

-    }

-

-

-    if ((cell_ > multipliers[dimensionsInWhichCellHasZeroLength[i]]) &&

-        (counter[dimensionsInWhichCellHasZeroLength[i]] != 0)) {

-      // if ( counter[dimensionsInWhichCellHasZeroLength[i]] != 0 )

-      if (bdg)std::cerr << "Subtracting : " << cell_ - multipliers[dimensionsInWhichCellHasZeroLength[i]] << std::endl;

-      coboundaryElements.push_back(cell_ - multipliers[dimensionsInWhichCellHasZeroLength[i]]);

-    }

-    if ((cell_ + multipliers[dimensionsInWhichCellHasZeroLength[i]] < this->data.size()) &&

-        (counter[dimensionsInWhichCellHasZeroLength[i]] != 2 * this->sizes[dimensionsInWhichCellHasZeroLength[i]])) {

-      // if ( counter[dimensionsInWhichCellHasZeroLength[i]] != 2*this->sizes[dimensionsInWhichCellHasZeroLength[i]] )

-      coboundaryElements.push_back(cell_ + multipliers[dimensionsInWhichCellHasZeroLength[i]]);

-      if (bdg)std::cerr << "Adding : " << cell_ + multipliers[dimensionsInWhichCellHasZeroLength[i]] << std::endl;

-    }

-  }

-  return coboundaryElements;

-}

-

-template <typename T>

-unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell_) {

-  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_ / multipliers[i - 1];

-

-    if (dbg)std::cerr << "i-1 :" << i - 1 << std::endl;

-    if (dbg)std::cerr << "cell_ : " << cell_ << std::endl;

-    if (dbg)std::cerr << "position : " << position << std::endl;

-    if (dbg)std::cerr << "multipliers[" << i - 1 << "] = " << multipliers[i - 1] << std::endl;

-    if (dbg)getchar();

-

-    if (position % 2 == 1) {

-      if (dbg)std::cerr << "Nonzero length in this direction \n";

-      dimension++;

-    }

-    cell_ = cell_ % multipliers[i - 1];

-  }

-  return dimension;

-}

-

-template <typename T>

-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> isThisCellConsidered(this->data.size(), false);

-

-  std::vector<size_t> indicesToConsider;

-  // 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_begin(); it != this->top_dimensional_cells_end(); ++it) {

-    indicesToConsider.push_back(it.computeIndexInBitmap());

-  }

-

-  while (indicesToConsider.size()) {

-    if (dbg) {

-      std::cerr << "indicesToConsider in this iteration \n";

-      for (size_t i = 0; i != indicesToConsider.size(); ++i) {

-        std::cout << indicesToConsider[i] << "  ";

-      }

-      getchar();

-    }

-    std::vector<size_t> newIndicesToConsider;

-    for (size_t i = 0; i != indicesToConsider.size(); ++i) {

-      std::vector<size_t> bd = this->get_boundary_of_a_cell(indicesToConsider[i]);

-      for (size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) {

-        if (this->data[ bd[boundaryIt] ] > this->data[ indicesToConsider[i] ]) {

-          this->data[ bd[boundaryIt] ] = this->data[ indicesToConsider[i] ];

-        }

-        if (isThisCellConsidered[ bd[boundaryIt] ] == false) {

-          newIndicesToConsider.push_back(bd[boundaryIt]);

-          isThisCellConsidered[ bd[boundaryIt] ] = true;

-        }

-      }

-    }

-    indicesToConsider.swap(newIndicesToConsider);

-  }

-}

-

-template <typename T>

-std::vector< size_t >

-Bitmap_cubical_complex_base<T>::generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions(std::vector< bool > directionsForPeriodicBCond_) {

-  bool dbg = false;

-  if (this->sizes.size() != directionsForPeriodicBCond_.size())

-    throw ("directionsForPeriodicBCond_ vector size is different from the size of the bitmap. "

-           "The program will now terminate \n");

-

-  std::vector<int> sizes(this->sizes.size());

-  for (size_t i = 0; i != this->sizes.size(); ++i)sizes[i] = 2 * this->sizes[i];

-

-  counter c(sizes);

-

-  std::vector< size_t > result;

-

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

-    size_t position;

-    if (!c.isFinal()) {

-      position = i;

-      // result.push_back( i );

-    } else {

-      std::vector< bool > finals = c.directionsOfFinals();

-      bool jumpInPosition = false;

-      for (size_t dir = 0; dir != finals.size(); ++dir) {

-        if (finals[dir] == false)continue;

-        if (directionsForPeriodicBCond_[dir]) {

-          jumpInPosition = true;

-        }

-      }

-      if (jumpInPosition == true) {

-        // in this case this guy is final, so we need to find 'the opposite one'

-        position = compute_position_in_bitmap(c.findOpposite(directionsForPeriodicBCond_));

-      } else {

-        position = i;

-      }

-    }

-    result.push_back(position);

-    if (dbg) {

-      std::cerr << " position : " << position << std::endl;

-      std::cerr << c << std::endl;

-      getchar();

-    }

-

-    c.increment();

-  }

-

-  return result;

-}

-

-#endif  // BITMAP_CUBICAL_COMPLEX_BASE_H_

+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    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 <iostream>
+#include <string>
+#include <vector>
+#include <string>
+#include <fstream>
+#include <algorithm>
+#include <iterator>
+#include <limits>
+#include "counter.h"
+
+
+using namespace std;
+
+namespace Gudhi
+{
+
+namespace Cubical_complex
+{
+
+
+
+/**
+ * This is a class implementing a basic bitmap data structure to store cubical complexes. 
+ * It implements only the most basic subroutines.
+ * The idea of the bitmap is the following. Our aim is to have a memory efficient 
+ * data structure to store d-dimensional cubical complex 
+ * C being a cubical decomposition
+ * of a rectangular region of a space. This is achieved by storing C as a 
+ * vector of bits (this is where the name 'bitmap' came from). 
+ * Each cell is represented by a single
+ * bit (in case of black and white bitmaps, or by a single element of a type T 
+ * (here T is a filtration type of a bitmap, typically a double). 
+ * All the informations needed for homology and
+ * persistent homology computations (like dimension of a cell, boundary and 
+ * coboundary elements of a cell, are then obtained from the 
+ * position of the element in C.
+ * The default filtration used in this implementation is the lower star filtration.
+ */
+template <typename T>
+class Bitmap_cubical_complex_base
+{
+public:
+    /**
+    * 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( 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( std::vector<unsigned>& dimensions , const std::vector<T>& top_dimensional_cells );
+
+    /**
+    * 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.
+    */
+    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.
+    **/
+    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.
+    **/
+    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_of_bitmap()const
+    {
+        return this->data.size();
+    }
+
+    /**
+    * Writing to stream operator.
+    **/
+    template <typename K>
+    friend ostream& operator << ( ostream & os , const Bitmap_cubical_complex_base<K>& b );
+
+    //ITERATORS
+
+    /**
+    * Iterator through all cells in the complex (in order they appear in the structure -- i.e. 
+    * in lexicographical order).
+    **/
+    typedef typename std::vector< T >::iterator all_cells_iterator;
+    all_cells_iterator all_cells_begin()const
+    {
+        return this->data.begin();
+    }
+    all_cells_iterator all_cells_end()const
+    {
+        return this->data.end();
+    }
+
+
+    typedef typename std::vector< T >::const_iterator all_cells_const_iterator;
+    all_cells_const_iterator all_cells_const_begin()const
+    {
+        return this->data.begin();
+    }
+    all_cells_const_iterator all_cells_const_end()const
+    {
+        return this->data.end();
+    }
+
+    /**
+    * 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, double >
+    {
+        public:
+            Top_dimensional_cells_iterator( Bitmap_cubical_complex_base& b ):b(b)
+            {
+                for ( size_t i = 0 ; i != b.dimension() ; ++i )
+                {
+                    this->counter.push_back(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];
+            }
+
+            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< unsigned > counter;
+            Bitmap_cubical_complex_base& b;
+    };
+    Top_dimensional_cells_iterator top_dimensional_cells_begin()
+    {
+        Top_dimensional_cells_iterator a(*this);
+        return a;
+    }
+    Top_dimensional_cells_iterator top_dimensional_cells_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;
+    }
+
+
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+
+
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+
+protected:
+    std::vector<unsigned> sizes;
+    std::vector<unsigned> multipliers;
+    std::vector<T> data;
+    size_t total_number_of_cells;
+    void set_up_containers( 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)+1;
+            multiplier *= 2*sizes[i]+1;
+        }
+        //std::reverse( this->sizes.begin() , this->sizes.end() );
+        std::vector<T> data(multiplier);
+        std::fill( data.begin() , data.end() , std::numeric_limits<int>::max() );
+        this->total_number_of_cells = multiplier;
+        this->data = data;
+    }
+
+    size_t compute_position_in_bitmap( 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;
+        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;
+    }
+
+    std::vector< size_t > 
+    generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions
+    ( std::vector< bool >& directions_for_periodic_b_cond );
+};
+
+
+
+
+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;
+}
+
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base
+( std::vector<unsigned>& sizes )
+{
+    this->set_up_containers( sizes );
+}
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base
+( 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\
+        Bitmap_cubical_complex_base\
+       ( std::vector<size_t> sizes_in_following_directions , std::vector<float> 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,\
+        std::vector<float> 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_begin() ; it != this->top_dimensional_cells_end() ; ++it )
+    {
+        (*it) = top_dimensional_cells[index];
+        ++index;
+    }
+    this->impose_lower_star_filtration();
+}
+
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( 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;}
+
+    std::vector<unsigned> sizes;
+    for ( size_t i = 0 ; i != dimensionOfData ; ++i )
+    {
+        unsigned size_in_this_dimension;
+        inFiltration >> size_in_this_dimension;
+        size_in_this_dimension = abs(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_begin();
+
+    //TODO -- over here we also need to read id's of cell and put them to bitmapElement structure!
+    while ( !inFiltration.eof() )
+    {
+        double 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;
+        }
+        *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
+{
+    bool bdg = false;
+    //first of all, we need to take the list of coordinates in which the cell has nonzero length. 
+    //We do it by using modified version to compute dimension of a cell:
+    std::vector< unsigned > dimensions_in_which_cell_has_nonzero_length;
+    unsigned dimension = 0;
+    size_t cell1 = cell;
+    for ( size_t i = this->multipliers.size() ; i != 0 ; --i )
+    {
+        unsigned position = cell1/multipliers[i-1];
+        if ( position%2 == 1 )
+        {
+            dimensions_in_which_cell_has_nonzero_length.push_back(i-1);
+            dimension++;
+        }
+        cell1 = cell1%multipliers[i-1];
+    }
+
+    if (bdg)
+    {
+        cerr << "dimensions_in_which_cell_has_nonzero_length : \n";
+        for ( size_t i = 0 ; i != dimensions_in_which_cell_has_nonzero_length.size() ; ++i )
+        {
+            cerr << dimensions_in_which_cell_has_nonzero_length[i] << endl;
+        }
+        getchar();
+    }
+
+    std::vector< size_t > boundary_elements;
+    if ( dimensions_in_which_cell_has_nonzero_length.size() == 0 )return boundary_elements;
+    for ( size_t i = 0 ; i != dimensions_in_which_cell_has_nonzero_length.size() ; ++i )
+    {
+        boundary_elements.push_back( cell - multipliers[ dimensions_in_which_cell_has_nonzero_length[i] ] );
+        boundary_elements.push_back( cell + multipliers[ dimensions_in_which_cell_has_nonzero_length[i] ] );
+
+        if (bdg) cerr << "multipliers[dimensions_in_which_cell_has_nonzero_length[i]] : " 
+        << multipliers[dimensions_in_which_cell_has_nonzero_length[i]] << endl;
+        if (bdg) cerr << "cell - multipliers[dimensions_in_which_cell_has_nonzero_length[i]] : " 
+        << cell - multipliers[dimensions_in_which_cell_has_nonzero_length[i]] << endl;
+        if (bdg) cerr << "cell + multipliers[dimensions_in_which_cell_has_nonzero_length[i]] : " 
+        << cell + multipliers[dimensions_in_which_cell_has_nonzero_length[i]] << endl;
+    }
+    return boundary_elements;
+}
+
+
+
+
+template <typename T>
+std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell( size_t cell )const
+{
+    bool bdg = false;
+    //first of all, we need to take the list of coordinates in which the cell has nonzero length. 
+    //We do it by using modified version to compute dimension of a cell:
+    std::vector< unsigned > dimensions_in_which_cell_has_zero_length;
+    unsigned dimension = 0;
+    size_t cell1 = cell;
+    for ( size_t i = this->multipliers.size() ; i != 0 ; --i )
+    {
+        unsigned position = cell1/multipliers[i-1];
+        if ( position%2 == 0 )
+        {
+            dimensions_in_which_cell_has_zero_length.push_back(i-1);
+            dimension++;
+        }
+        cell1 = cell1%multipliers[i-1];
+    }
+
+    std::vector<unsigned> counter = this->compute_counter_for_given_cell( cell );
+    //reverse(counter.begin() , counter.end());
+
+    if (bdg)
+    {
+        cerr << "dimensions_in_which_cell_has_zero_length : \n";
+        for ( size_t i = 0 ; i != dimensions_in_which_cell_has_zero_length.size() ; ++i )
+        {
+            cerr << dimensions_in_which_cell_has_zero_length[i] << endl;
+        }
+        cerr << "\n counter : " << endl;
+        for ( size_t i = 0 ; i != counter.size() ; ++i )
+        {
+            cerr << counter[i] << endl;
+        }
+        getchar();
+    }
+
+    std::vector< size_t > coboundary_elements;
+    if ( dimensions_in_which_cell_has_zero_length.size() == 0 )return coboundary_elements;
+    for ( size_t i = 0 ; i != dimensions_in_which_cell_has_zero_length.size() ; ++i )
+    {
+       if ( bdg )
+        {
+            cerr << "Dimension : " << i << endl;
+            if (counter[dimensions_in_which_cell_has_zero_length[i]] == 0)
+            {
+                cerr << "In dimension : " << i 
+                << " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n";
+            }
+            if ( counter[dimensions_in_which_cell_has_zero_length[i]] 
+                 == 
+                 2*this->sizes[dimensions_in_which_cell_has_zero_length[i]] )
+            {
+                cerr << "In dimension : " << i 
+                     << " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n";
+            }
+        }
+
+
+        if ( (cell > multipliers[dimensions_in_which_cell_has_zero_length[i]]) 
+              && (counter[dimensions_in_which_cell_has_zero_length[i]] != 0) )
+        //if ( counter[dimensions_in_which_cell_has_zero_length[i]] != 0 )
+        {
+            if ( bdg )
+            {
+                 cerr << "Subtracting : " << cell - multipliers[dimensions_in_which_cell_has_zero_length[i]] << endl;
+            }
+            coboundary_elements.push_back( cell - multipliers[dimensions_in_which_cell_has_zero_length[i]] );
+        }
+        if ( 
+           (cell + multipliers[dimensions_in_which_cell_has_zero_length[i]] < this->data.size()) && 
+           (counter[dimensions_in_which_cell_has_zero_length[i]] 
+           != 
+           2*this->sizes[dimensions_in_which_cell_has_zero_length[i]]) 
+           )
+        //if ( counter[dimensions_in_which_cell_has_zero_length[i]] != 
+        //2*this->sizes[dimensions_in_which_cell_has_zero_length[i]] )
+        {
+            coboundary_elements.push_back( cell + multipliers[dimensions_in_which_cell_has_zero_length[i]] );
+            if ( bdg )cerr << "Adding : " << cell + multipliers[dimensions_in_which_cell_has_zero_length[i]] << endl;
+        }
+    }
+    return coboundary_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/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<<"] = " << multipliers[i-1] << endl;
+        if (dbg)getchar();
+
+        if ( position%2 == 1 )
+        {
+            if (dbg)cerr << "Nonzero length in this direction \n";
+            dimension++;
+        }
+        cell = cell%multipliers[i-1];
+    }
+    return dimension;
+}
+
+template <typename T>
+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 );
+
+    std::vector<size_t> indices_to_consider;
+    //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_begin() ; it != this->top_dimensional_cells_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] << "  ";
+            }
+            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 ( this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ] )
+                {
+                    this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ];
+                }
+                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);
+    }
+}
+
+
+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;
+    }
+}
+
+
+
+template <typename T>
+std::vector< size_t > Bitmap_cubical_complex_base<T>::
+generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions
+( std::vector< bool >& directions_for_periodic_b_cond )
+{
+    bool dbg = false;
+    if ( this->sizes.size() != directions_for_periodic_b_cond.size() )
+    throw "directions_for_periodic_b_cond vector size is different from the size of the bitmap. Program terminate \n";
+
+    std::vector<unsigned> sizes( this->sizes.size() );
+    for ( size_t i = 0 ; i != this->sizes.size() ; ++i )sizes[i] = 2*this->sizes[i];
+
+    counter c( sizes );
+
+    std::vector< size_t > result;
+
+    for ( size_t i = 0 ; i != this->data.size() ; ++i )
+    {
+        size_t position;
+        if ( !c.isFinal() )
+        {
+            position = i;
+            //result.push_back( i );
+        }
+        else
+        {
+            std::vector< bool > finals = c.directions_of_finals();
+            bool jump_in_position = false;
+            for ( size_t dir = 0 ; dir != finals.size() ; ++dir )
+            {
+                if ( finals[dir] == false )continue;
+                if ( directions_for_periodic_b_cond[dir] )
+                {
+                    jump_in_position = true;
+                }
+            }
+            if ( jump_in_position == true )
+            {
+                //in this case this guy is final, so we need to find 'the opposite one'
+                position = compute_position_in_bitmap(  c.find_opposite( directions_for_periodic_b_cond ) );
+            }
+            else
+            {
+                position = i;
+            }
+        }
+        result.push_back( position );
+        if ( dbg )
+        {
+            cerr << " position : " << position << endl;
+            cerr << c << endl;
+            getchar();
+        }
+
+        c.increment();
+    }
+
+    return result;
+}
+
+}
+
+}
\ No newline at end of file
diff --git a/src/Bitmap_cubical_complex/include/gudhi/counter.h b/src/Bitmap_cubical_complex/include/gudhi/counter.h
index 9445a422..a5fda36f 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/counter.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/counter.h
@@ -1,139 +1,177 @@
-/*    This file is part of the Gudhi Library. The Gudhi library

- *    (Geometric Understanding in Higher Dimensions) is a generic C++

- *    library for computational topology.

- *

- *    Author(s):       Pawel Dlotko

- *

- *    Copyright (C) 2015  INRIA Sophia-Saclay (France)

- *

- *    This program is free software: you can redistribute it and/or modify

- *    it under the terms of the GNU General Public License as published by

- *    the Free Software Foundation, either version 3 of the License, or

- *    (at your option) any later version.

- *

- *    This program is distributed in the hope that it will be useful,

- *    but WITHOUT ANY WARRANTY; without even the implied warranty of

- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

- *    GNU General Public License for more details.

- *

- *    You should have received a copy of the GNU General Public License

- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.

- */

-

-#ifndef COUNTER_H_

-#define COUNTER_H_

-

-#include <iostream>

-#include <vector>

-

-/**

- * This is an implementation of a simple counter. It is needed for the implementation of a bitmapCubicalComplex.

- **/

-

-class counter {

- public:

-  /**

-   * Constructor of a counter class. It takes only the parameter which is the end value of the counter.

-   * The default beginning value is a vector of the same length as the end, filled-in with zeros.

-   **/

-  counter(std::vector< int > endd) {

-    for (size_t i = 0; i != endd.size(); ++i) {

-      this->current.push_back(0);

-      this->begin.push_back(0);

-      this->end.push_back(endd[i]);

-    }

-  }

-

-  /**

-   * Constructor of a counter class. It takes as the input beginn and endd vector. It assumes that begin vector is

-   * lexicographically below the end vector.

-   **/

-  counter(std::vector< int > beginn, std::vector< int > endd) {

-    if (beginn.size() != endd.size())

-      throw("In constructor of a counter, begin and end vectors do not have the same size. Program terminate");

-    for (size_t i = 0; i != endd.size(); ++i) {

-      this->current.push_back(0);

-      this->begin.push_back(0);

-      this->end.push_back(endd[i]);

-    }

-  }

-

-  /**

-   * Function to increment the counter. If the value returned by the function is true, then the incrementation process

-   * was successful.

-   * If the value of the function is false, that means, that the counter have reached its end-value.

-   **/

-  bool increment() {

-    size_t i = 0;

-    while ((i != this->end.size()) && (this->current[i] == this->end[i])) {

-      ++i;

-    }

-

-    if (i == this->end.size())return false;

-    ++this->current[i];

-    for (size_t j = 0; j != i; ++j) {

-      this->current[j] = this->begin[j];

-    }

-    return true;

-  }

-

-  /**

-   * Function to check if we are at the end of counter.

-   **/

-  bool isFinal() {

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

-      if (this->current[i] == this->end[i])return true;

-    }

-    return false;

-  }

-

-  /**

-   * Function required in the implementation of bitmapCubicalComplexWPeriodicBoundaryCondition.

-   * Its aim is to find an counter corresponding to the element the following boundary element is identified with

-   * when periodic boundary conditions are imposed.

-   **/

-  std::vector< int > findOpposite(std::vector< bool > directionsForPeriodicBCond) {

-    std::vector< int > result;

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

-      if ((this->current[i] == this->end[i]) && (directionsForPeriodicBCond[i] == true)) {

-        result.push_back(this->begin[i]);

-      } else {

-        result.push_back(this->current[i]);

-      }

-    }

-    return result;

-  }

-

-  /**

-   * Function checking at which positions the current value of a counter is the final value of the counter.

-   **/

-  std::vector< bool > directionsOfFinals() {

-    std::vector< bool > result;

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

-      if (this->current[i] == this->end[i]) {

-        result.push_back(true);

-      } else {

-        result.push_back(false);

-      }

-    }

-    return result;

-  }

-

-  /**

-   * Function to write counter to the stream.

-   **/

-  friend std::ostream& operator<<(std::ostream& out, const counter& c) {

-    // std::cerr << "c.current.size() : " << c.current.size() << std::endl;

-    for (size_t i = 0; i != c.current.size(); ++i) {

-      out << c.current[i] << " ";

-    }

-    return out;

-  }

-

- private:

-  std::vector< int > begin;

-  std::vector< int > end;

-  std::vector< int > current;

-};

-

-#endif  // COUNTER_H_

+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    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 <iostream>
+#include <vector>
+
+using namespace std;
+
+namespace Gudhi
+{
+
+namespace Cubical_complex
+{
+
+/**
+* This is an implementation of a counter being a vector of integers. 
+* The constructor of the class takes as an input two vectors W and V. 
+* It assumes that W < V coordinatewise.
+* If the initial counter W is not specified, it is assumed to be vector of zeros. 
+* The class allows to iterate between W and V by using increment() function.
+* The increment() function returns a bool value. 
+* The current counter reach the end counter V if the value returned by the increment function is FALSE.
+* This class is needed for the implementation of a bitmapCubicalComplex.
+**/
+
+class counter
+{
+public:
+    /**
+    * Constructor of a counter class. It takes only the parameter which is the end value of the counter. 
+    * The default beginning value is a vector of the same length as the endd, filled-in with zeros.
+    **/
+    counter(std::vector<unsigned> const& endd): begin(endd.size(),0), end(endd), current(endd.size(),0){}
+	//counter(std::vector< int >& endd)
+	//{
+	//    for ( size_t i = 0 ; i != endd.size() ; ++i )
+    //    {
+    //        this->current.push_back(0);
+    //        this->begin.push_back(0);
+    //        this->end.push_back( endd[i] );
+    //    }
+	//}
+
+
+    /**
+    * Constructor of a counter class. It takes as the input beginn and end vector. 
+    * It assumes that begin vector is lexicographically below the end vector.
+    **/
+	counter(std::vector< unsigned >& beginn , std::vector< unsigned >& endd)
+	{
+	    if ( beginn.size() != endd.size() )
+		throw "In constructor of a counter, begin and end vectors do not have the same size. Program terminate";
+	    for ( size_t i = 0 ; i != endd.size() ; ++i )
+        {
+            this->current.push_back(0);
+            this->begin.push_back(0);
+            this->end.push_back( endd[i] );
+        }
+	}
+
+    /**
+    * Function to increment the counter. If the value returned by the function is true, 
+    * then the incrementation process was successful.
+    * If the value of the function is false, that means, that the counter have reached its end-value.
+    **/
+	bool increment()
+	{
+	    size_t i = 0;
+	    while( (i != this->end.size()) && (this->current[i] == this->end[i]) )
+        {
+            ++i;
+        }
+
+        if ( i == this->end.size() )return false;
+        ++this->current[i];
+        for ( size_t j = 0 ; j != i ; ++j )
+        {
+            this->current[j] = this->begin[j];
+        }
+        return true;
+	}
+
+    /**
+    * Function to check if we are at the end of counter.
+    **/
+	bool isFinal()
+	{
+	    for ( size_t i = 0 ; i != this->current.size() ; ++i )
+        {
+            if ( this->current[i] == this->end[i] )return true;
+        }
+        return false;
+	}
+
+    /**
+    * Function required in the implementation of bitmapCubicalComplexWPeriodicBoundaryCondition. 
+    * Its aim is to find an counter corresponding to the element the following
+    * boundary element is identified with when periodic boundary conditions are imposed.
+    **/
+    std::vector< unsigned > find_opposite( std::vector< bool >& directionsForPeriodicBCond )
+    {
+        std::vector< unsigned > result;
+        for ( size_t i = 0 ; i != this->current.size() ; ++i )
+        {
+            if ( (this->current[i] == this->end[i]) && (directionsForPeriodicBCond[i] == true) )
+            {
+                result.push_back( this->begin[i] );
+            }
+            else
+            {
+                result.push_back( this->current[i] );
+            }
+        }
+        return result;
+    }
+
+    /**
+    * Function checking at which positions the current value of a counter is the final value of the counter.
+    **/
+    std::vector< bool > directions_of_finals()
+    {
+        std::vector< bool > result;
+        for ( size_t i = 0 ; i != this->current.size() ; ++i )
+        {
+            if ( this->current[i] == this->end[i] )
+            {
+                result.push_back( true );
+            }
+            else
+            {
+                result.push_back( false );
+            }
+        }
+        return result;
+    }
+
+    /**
+    * Function to write counter to the stream.
+    **/
+    friend std::ostream& operator<<(std::ostream& out , const counter& c )
+    {
+        //cerr << "c.current.size() : " << c.current.size() << endl;
+        for ( size_t i = 0 ; i != c.current.size() ; ++i )
+        {
+            out << c.current[i] << " ";
+        }
+        return out;
+    }
+private:
+	std::vector< unsigned > begin;
+	std::vector< unsigned > end;
+	std::vector< unsigned > current;
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
index 1c204bae..183b856a 100644
--- a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
+++ b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
@@ -1,623 +1,638 @@
-#include "gudhi/Bitmap_cubical_complex.h"

-

-#define BOOST_TEST_DYN_LINK

-#define BOOST_TEST_MODULE "cubical_complex"

-#include <boost/test/unit_test.hpp>

-

-using namespace std;

-

-BOOST_AUTO_TEST_CASE(check_dimension) {

-  std::vector< double > increasingFiltrationOfTopDimensionalCells;

-  increasingFiltrationOfTopDimensionalCells.push_back(1);

-  increasingFiltrationOfTopDimensionalCells.push_back(2);

-  increasingFiltrationOfTopDimensionalCells.push_back(3);

-  increasingFiltrationOfTopDimensionalCells.push_back(4);

-  increasingFiltrationOfTopDimensionalCells.push_back(5);

-  increasingFiltrationOfTopDimensionalCells.push_back(6);

-  increasingFiltrationOfTopDimensionalCells.push_back(7);

-  increasingFiltrationOfTopDimensionalCells.push_back(8);

-  increasingFiltrationOfTopDimensionalCells.push_back(9);

-

-  std::vector<unsigned> dimensions;

-  dimensions.push_back(3);

-  dimensions.push_back(3);

-

-  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);

-  BOOST_CHECK(increasing.dimension() == 2);

-}

-

-BOOST_AUTO_TEST_CASE(topDimensionalCellsIterator_test) {

-  std::vector< double > expectedFiltrationValues1;

-  expectedFiltrationValues1.push_back(0);

-  expectedFiltrationValues1.push_back(0);

-  expectedFiltrationValues1.push_back(0);

-  expectedFiltrationValues1.push_back(0);

-  expectedFiltrationValues1.push_back(100);

-  expectedFiltrationValues1.push_back(0);

-  expectedFiltrationValues1.push_back(0);

-  expectedFiltrationValues1.push_back(0);

-  expectedFiltrationValues1.push_back(0);

-

-  std::vector< double > expectedFiltrationValues2;

-  expectedFiltrationValues2.push_back(1);

-  expectedFiltrationValues2.push_back(2);

-  expectedFiltrationValues2.push_back(3);

-  expectedFiltrationValues2.push_back(4);

-  expectedFiltrationValues2.push_back(5);

-  expectedFiltrationValues2.push_back(6);

-  expectedFiltrationValues2.push_back(7);

-  expectedFiltrationValues2.push_back(8);

-  expectedFiltrationValues2.push_back(9);

-

-  std::vector< double > increasingFiltrationOfTopDimensionalCells;

-  increasingFiltrationOfTopDimensionalCells.push_back(1);

-  increasingFiltrationOfTopDimensionalCells.push_back(2);

-  increasingFiltrationOfTopDimensionalCells.push_back(3);

-  increasingFiltrationOfTopDimensionalCells.push_back(4);

-  increasingFiltrationOfTopDimensionalCells.push_back(5);

-  increasingFiltrationOfTopDimensionalCells.push_back(6);

-  increasingFiltrationOfTopDimensionalCells.push_back(7);

-  increasingFiltrationOfTopDimensionalCells.push_back(8);

-  increasingFiltrationOfTopDimensionalCells.push_back(9);

-

-  std::vector< double > oneDimensionalCycle;

-  oneDimensionalCycle.push_back(0);

-  oneDimensionalCycle.push_back(0);

-  oneDimensionalCycle.push_back(0);

-  oneDimensionalCycle.push_back(0);

-  oneDimensionalCycle.push_back(100);

-  oneDimensionalCycle.push_back(0);

-  oneDimensionalCycle.push_back(0);

-  oneDimensionalCycle.push_back(0);

-  oneDimensionalCycle.push_back(0);

-

-  std::vector<unsigned> dimensions;

-  dimensions.push_back(3);

-  dimensions.push_back(3);

-

-  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);

-  Bitmap_cubical_complex<double> hole(dimensions, oneDimensionalCycle);

-

-

-  int i = 0;

-  for (Bitmap_cubical_complex<double>::Top_dimensional_cells_iterator it = increasing.top_dimensional_cells_begin(); it != increasing.top_dimensional_cells_end(); ++it) {

-    BOOST_CHECK(*it == expectedFiltrationValues2[i]);

-    ++i;

-  }

-  i = 0;

-  for (Bitmap_cubical_complex<double>::Top_dimensional_cells_iterator it = hole.top_dimensional_cells_begin(); it != hole.top_dimensional_cells_end(); ++it) {

-    BOOST_CHECK(*it == expectedFiltrationValues1[i]);

-    ++i;

-  }

-}

-

-BOOST_AUTO_TEST_CASE(compute_boundary_test_1) {

-

-  std::vector<double> boundary0;

-  std::vector<double> boundary1;

-  boundary1.push_back(0);

-  boundary1.push_back(2);

-  std::vector<double> boundary2;

-  std::vector<double> boundary3;

-  boundary3.push_back(2);

-  boundary3.push_back(4);

-  std::vector<double> boundary4;

-  std::vector<double> boundary5;

-  boundary5.push_back(4);

-  boundary5.push_back(6);

-  std::vector<double> boundary6;

-  std::vector<double> boundary7;

-  boundary7.push_back(0);

-  boundary7.push_back(14);

-  std::vector<double> boundary8;

-  boundary8.push_back(1);

-  boundary8.push_back(15);

-  boundary8.push_back(7);

-  boundary8.push_back(9);

-  std::vector<double> boundary9;

-  boundary9.push_back(2);

-  boundary9.push_back(16);

-  std::vector<double> boundary10;

-  boundary10.push_back(3);

-  boundary10.push_back(17);

-  boundary10.push_back(9);

-  boundary10.push_back(11);

-  std::vector<double> boundary11;

-  boundary11.push_back(4);

-  boundary11.push_back(18);

-  std::vector<double> boundary12;

-  boundary12.push_back(5);

-  boundary12.push_back(19);

-  boundary12.push_back(11);

-  boundary12.push_back(13);

-  std::vector<double> boundary13;

-  boundary13.push_back(6);

-  boundary13.push_back(20);

-  std::vector<double> boundary14;

-  std::vector<double> boundary15;

-  boundary15.push_back(14);

-  boundary15.push_back(16);

-  std::vector<double> boundary16;

-  std::vector<double> boundary17;

-  boundary17.push_back(16);

-  boundary17.push_back(18);

-  std::vector<double> boundary18;

-  std::vector<double> boundary19;

-  boundary19.push_back(18);

-  boundary19.push_back(20);

-  std::vector<double> boundary20;

-  std::vector<double> boundary21;

-  boundary21.push_back(14);

-  boundary21.push_back(28);

-  std::vector<double> boundary22;

-  boundary22.push_back(15);

-  boundary22.push_back(29);

-  boundary22.push_back(21);

-  boundary22.push_back(23);

-  std::vector<double> boundary23;

-  boundary23.push_back(16);

-  boundary23.push_back(30);

-  std::vector<double> boundary24;

-  boundary24.push_back(17);

-  boundary24.push_back(31);

-  boundary24.push_back(23);

-  boundary24.push_back(25);

-  std::vector<double> boundary25;

-  boundary25.push_back(18);

-  boundary25.push_back(32);

-  std::vector<double> boundary26;

-  boundary26.push_back(19);

-  boundary26.push_back(33);

-  boundary26.push_back(25);

-  boundary26.push_back(27);

-  std::vector<double> boundary27;

-  boundary27.push_back(20);

-  boundary27.push_back(34);

-  std::vector<double> boundary28;

-  std::vector<double> boundary29;

-  boundary29.push_back(28);

-  boundary29.push_back(30);

-  std::vector<double> boundary30;

-  std::vector<double> boundary31;

-  boundary31.push_back(30);

-  boundary31.push_back(32);

-  std::vector<double> boundary32;

-  std::vector<double> boundary33;

-  boundary33.push_back(32);

-  boundary33.push_back(34);

-  std::vector<double> boundary34;

-  std::vector<double> boundary35;

-  boundary35.push_back(28);

-  boundary35.push_back(42);

-  std::vector<double> boundary36;

-  boundary36.push_back(29);

-  boundary36.push_back(43);

-  boundary36.push_back(35);

-  boundary36.push_back(37);

-  std::vector<double> boundary37;

-  boundary37.push_back(30);

-  boundary37.push_back(44);

-  std::vector<double> boundary38;

-  boundary38.push_back(31);

-  boundary38.push_back(45);

-  boundary38.push_back(37);

-  boundary38.push_back(39);

-  std::vector<double> boundary39;

-  boundary39.push_back(32);

-  boundary39.push_back(46);

-  std::vector<double> boundary40;

-  boundary40.push_back(33);

-  boundary40.push_back(47);

-  boundary40.push_back(39);

-  boundary40.push_back(41);

-  std::vector<double> boundary41;

-  boundary41.push_back(34);

-  boundary41.push_back(48);

-  std::vector<double> boundary42;

-  std::vector<double> boundary43;

-  boundary43.push_back(42);

-  boundary43.push_back(44);

-  std::vector<double> boundary44;

-  std::vector<double> boundary45;

-  boundary45.push_back(44);

-  boundary45.push_back(46);

-  std::vector<double> boundary46;

-  std::vector<double> boundary47;

-  boundary47.push_back(46);

-  boundary47.push_back(48);

-  std::vector<double> boundary48;

-  std::vector< std::vector<double> > boundaries;

-  boundaries.push_back(boundary0);

-  boundaries.push_back(boundary1);

-  boundaries.push_back(boundary2);

-  boundaries.push_back(boundary3);

-  boundaries.push_back(boundary4);

-  boundaries.push_back(boundary5);

-  boundaries.push_back(boundary6);

-  boundaries.push_back(boundary7);

-  boundaries.push_back(boundary8);

-  boundaries.push_back(boundary9);

-  boundaries.push_back(boundary10);

-  boundaries.push_back(boundary11);

-  boundaries.push_back(boundary12);

-  boundaries.push_back(boundary13);

-  boundaries.push_back(boundary14);

-  boundaries.push_back(boundary15);

-  boundaries.push_back(boundary16);

-  boundaries.push_back(boundary17);

-  boundaries.push_back(boundary18);

-  boundaries.push_back(boundary19);

-  boundaries.push_back(boundary20);

-  boundaries.push_back(boundary21);

-  boundaries.push_back(boundary22);

-  boundaries.push_back(boundary23);

-  boundaries.push_back(boundary24);

-  boundaries.push_back(boundary25);

-  boundaries.push_back(boundary26);

-  boundaries.push_back(boundary27);

-  boundaries.push_back(boundary28);

-  boundaries.push_back(boundary29);

-  boundaries.push_back(boundary30);

-  boundaries.push_back(boundary31);

-  boundaries.push_back(boundary32);

-  boundaries.push_back(boundary33);

-  boundaries.push_back(boundary34);

-  boundaries.push_back(boundary35);

-  boundaries.push_back(boundary36);

-  boundaries.push_back(boundary37);

-  boundaries.push_back(boundary38);

-  boundaries.push_back(boundary39);

-  boundaries.push_back(boundary40);

-  boundaries.push_back(boundary41);

-  boundaries.push_back(boundary42);

-  boundaries.push_back(boundary43);

-  boundaries.push_back(boundary44);

-  boundaries.push_back(boundary45);

-  boundaries.push_back(boundary46);

-  boundaries.push_back(boundary47);

-  boundaries.push_back(boundary48);

-

-

-

-  std::vector< double > increasingFiltrationOfTopDimensionalCells;

-  increasingFiltrationOfTopDimensionalCells.push_back(1);

-  increasingFiltrationOfTopDimensionalCells.push_back(2);

-  increasingFiltrationOfTopDimensionalCells.push_back(3);

-  increasingFiltrationOfTopDimensionalCells.push_back(4);

-  increasingFiltrationOfTopDimensionalCells.push_back(5);

-  increasingFiltrationOfTopDimensionalCells.push_back(6);

-  increasingFiltrationOfTopDimensionalCells.push_back(7);

-  increasingFiltrationOfTopDimensionalCells.push_back(8);

-  increasingFiltrationOfTopDimensionalCells.push_back(9);

-

-  std::vector<unsigned> dimensions;

-  dimensions.push_back(3);

-  dimensions.push_back(3);

-

-  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);

-  for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) {

-    std::vector< size_t > bd = increasing.get_boundary_of_a_cell(i);

-    for (size_t j = 0; j != bd.size(); ++j) {

-      BOOST_CHECK(boundaries[i][j] == bd[j]);

-    }

-  }

-}

-

-BOOST_AUTO_TEST_CASE(compute_boundary_test_2) {

-  std::vector< double > increasingFiltrationOfTopDimensionalCells;

-  increasingFiltrationOfTopDimensionalCells.push_back(1);

-  increasingFiltrationOfTopDimensionalCells.push_back(2);

-  increasingFiltrationOfTopDimensionalCells.push_back(3);

-  increasingFiltrationOfTopDimensionalCells.push_back(4);

-  increasingFiltrationOfTopDimensionalCells.push_back(5);

-  increasingFiltrationOfTopDimensionalCells.push_back(6);

-  increasingFiltrationOfTopDimensionalCells.push_back(7);

-  increasingFiltrationOfTopDimensionalCells.push_back(8);

-  increasingFiltrationOfTopDimensionalCells.push_back(9);

-

-  std::vector<unsigned> dimensions;

-  dimensions.push_back(3);

-  dimensions.push_back(3);

-

-  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);

-

-

-  std::vector<double> coboundaryElements;

-  coboundaryElements.push_back(7);

-  coboundaryElements.push_back(1);

-  coboundaryElements.push_back(8);

-  coboundaryElements.push_back(9);

-  coboundaryElements.push_back(1);

-  coboundaryElements.push_back(3);

-  coboundaryElements.push_back(10);

-  coboundaryElements.push_back(11);

-  coboundaryElements.push_back(3);

-  coboundaryElements.push_back(5);

-  coboundaryElements.push_back(12);

-  coboundaryElements.push_back(13);

-  coboundaryElements.push_back(5);

-  coboundaryElements.push_back(8);

-  coboundaryElements.push_back(8);

-  coboundaryElements.push_back(10);

-  coboundaryElements.push_back(10);

-  coboundaryElements.push_back(12);

-  coboundaryElements.push_back(12);

-  coboundaryElements.push_back(7);

-  coboundaryElements.push_back(21);

-  coboundaryElements.push_back(15);

-  coboundaryElements.push_back(8);

-  coboundaryElements.push_back(22);

-  coboundaryElements.push_back(9);

-  coboundaryElements.push_back(23);

-  coboundaryElements.push_back(15);

-  coboundaryElements.push_back(17);

-  coboundaryElements.push_back(10);

-  coboundaryElements.push_back(24);

-  coboundaryElements.push_back(11);

-  coboundaryElements.push_back(25);

-  coboundaryElements.push_back(17);

-  coboundaryElements.push_back(19);

-  coboundaryElements.push_back(12);

-  coboundaryElements.push_back(26);

-  coboundaryElements.push_back(13);

-  coboundaryElements.push_back(27);

-  coboundaryElements.push_back(19);

-  coboundaryElements.push_back(22);

-  coboundaryElements.push_back(22);

-  coboundaryElements.push_back(24);

-  coboundaryElements.push_back(24);

-  coboundaryElements.push_back(26);

-  coboundaryElements.push_back(26);

-  coboundaryElements.push_back(21);

-  coboundaryElements.push_back(35);

-  coboundaryElements.push_back(29);

-  coboundaryElements.push_back(22);

-  coboundaryElements.push_back(36);

-  coboundaryElements.push_back(23);

-  coboundaryElements.push_back(37);

-  coboundaryElements.push_back(29);

-  coboundaryElements.push_back(31);

-  coboundaryElements.push_back(24);

-  coboundaryElements.push_back(38);

-  coboundaryElements.push_back(25);

-  coboundaryElements.push_back(39);

-  coboundaryElements.push_back(31);

-  coboundaryElements.push_back(33);

-  coboundaryElements.push_back(26);

-  coboundaryElements.push_back(40);

-  coboundaryElements.push_back(27);

-  coboundaryElements.push_back(41);

-  coboundaryElements.push_back(33);

-  coboundaryElements.push_back(36);

-  coboundaryElements.push_back(36);

-  coboundaryElements.push_back(38);

-  coboundaryElements.push_back(38);

-  coboundaryElements.push_back(40);

-  coboundaryElements.push_back(40);

-  coboundaryElements.push_back(35);

-  coboundaryElements.push_back(43);

-  coboundaryElements.push_back(36);

-  coboundaryElements.push_back(37);

-  coboundaryElements.push_back(43);

-  coboundaryElements.push_back(45);

-  coboundaryElements.push_back(38);

-  coboundaryElements.push_back(39);

-  coboundaryElements.push_back(45);

-  coboundaryElements.push_back(47);

-  coboundaryElements.push_back(40);

-  coboundaryElements.push_back(41);

-  coboundaryElements.push_back(47);

-  size_t number = 0;

-  for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) {

-    std::vector< size_t > bd = increasing.get_coboundary_of_a_cell(i);

-    for (size_t j = 0; j != bd.size(); ++j) {

-      BOOST_CHECK(coboundaryElements[number] == bd[j]);

-      ++number;

-    }

-

-  }

-}

-

-BOOST_AUTO_TEST_CASE(compute_boundary_test_3) {

-  std::vector< double > increasingFiltrationOfTopDimensionalCells;

-  increasingFiltrationOfTopDimensionalCells.push_back(1);

-  increasingFiltrationOfTopDimensionalCells.push_back(2);

-  increasingFiltrationOfTopDimensionalCells.push_back(3);

-  increasingFiltrationOfTopDimensionalCells.push_back(4);

-  increasingFiltrationOfTopDimensionalCells.push_back(5);

-  increasingFiltrationOfTopDimensionalCells.push_back(6);

-  increasingFiltrationOfTopDimensionalCells.push_back(7);

-  increasingFiltrationOfTopDimensionalCells.push_back(8);

-  increasingFiltrationOfTopDimensionalCells.push_back(9);

-

-  std::vector<unsigned> dimensions;

-  dimensions.push_back(3);

-  dimensions.push_back(3);

-

-  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);

-

-  std::vector<unsigned> dim;

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(0);

-

-  for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) {

-    BOOST_CHECK(increasing.get_dimension_of_a_cell(i) == dim[i]);

-  }

-}

-

-BOOST_AUTO_TEST_CASE(Filtration_simplex_iterator_test) {

-  std::vector< double > increasingFiltrationOfTopDimensionalCells;

-  increasingFiltrationOfTopDimensionalCells.push_back(1);

-  increasingFiltrationOfTopDimensionalCells.push_back(2);

-  increasingFiltrationOfTopDimensionalCells.push_back(3);

-  increasingFiltrationOfTopDimensionalCells.push_back(4);

-  increasingFiltrationOfTopDimensionalCells.push_back(5);

-  increasingFiltrationOfTopDimensionalCells.push_back(6);

-  increasingFiltrationOfTopDimensionalCells.push_back(7);

-  increasingFiltrationOfTopDimensionalCells.push_back(8);

-  increasingFiltrationOfTopDimensionalCells.push_back(9);

-

-  std::vector<unsigned> dimensions;

-  dimensions.push_back(3);

-  dimensions.push_back(3);

-

-  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);

-

-  std::vector< unsigned > dim;

-  dim.push_back(0);

-  dim.push_back(0);

-  dim.push_back(0);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-  dim.push_back(0);

-  dim.push_back(1);

-  dim.push_back(1);

-  dim.push_back(2);

-

-  std::vector<double> fil;

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(1);

-  fil.push_back(2);

-  fil.push_back(2);

-  fil.push_back(2);

-  fil.push_back(2);

-  fil.push_back(2);

-  fil.push_back(2);

-  fil.push_back(3);

-  fil.push_back(3);

-  fil.push_back(3);

-  fil.push_back(3);

-  fil.push_back(3);

-  fil.push_back(3);

-  fil.push_back(4);

-  fil.push_back(4);

-  fil.push_back(4);

-  fil.push_back(4);

-  fil.push_back(4);

-  fil.push_back(4);

-  fil.push_back(5);

-  fil.push_back(5);

-  fil.push_back(5);

-  fil.push_back(5);

-  fil.push_back(6);

-  fil.push_back(6);

-  fil.push_back(6);

-  fil.push_back(6);

-  fil.push_back(7);

-  fil.push_back(7);

-  fil.push_back(7);

-  fil.push_back(7);

-  fil.push_back(7);

-  fil.push_back(7);

-  fil.push_back(8);

-  fil.push_back(8);

-  fil.push_back(8);

-  fil.push_back(8);

-  fil.push_back(9);

-  fil.push_back(9);

-  fil.push_back(9);

-  fil.push_back(9);

-

-

-  Bitmap_cubical_complex<double>::Filtration_simplex_range range = increasing.filtration_simplex_range();

-  size_t position = 0;

-  for (Bitmap_cubical_complex<double>::Filtration_simplex_iterator it = range.begin(); it != range.end(); ++it) {

-    BOOST_CHECK(increasing.dimension(*it) == dim[position]);

-    BOOST_CHECK(increasing.filtration(*it) == fil[position]);

-    ++position;

-  }

-}

+#include <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <boost/program_options.hpp>
+
+// standard stuff
+#include <iostream>
+#include <sstream>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "cubical_complex"
+#include <boost/test/unit_test.hpp>
+
+
+using namespace std;
+using namespace Gudhi;
+using namespace Gudhi::Cubical_complex;
+using namespace Gudhi::persistent_cohomology;
+
+
+
+BOOST_AUTO_TEST_CASE(check_dimension) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells;
+  increasingFiltrationOfTopDimensionalCells.push_back(1);
+  increasingFiltrationOfTopDimensionalCells.push_back(2);
+  increasingFiltrationOfTopDimensionalCells.push_back(3);
+  increasingFiltrationOfTopDimensionalCells.push_back(4);
+  increasingFiltrationOfTopDimensionalCells.push_back(5);
+  increasingFiltrationOfTopDimensionalCells.push_back(6);
+  increasingFiltrationOfTopDimensionalCells.push_back(7);
+  increasingFiltrationOfTopDimensionalCells.push_back(8);
+  increasingFiltrationOfTopDimensionalCells.push_back(9);
+
+  std::vector<unsigned> dimensions;
+  dimensions.push_back(3);
+  dimensions.push_back(3);
+
+  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+  BOOST_CHECK(increasing.dimension() == 2);
+}
+
+BOOST_AUTO_TEST_CASE(topDimensionalCellsIterator_test) {
+  std::vector< double > expectedFiltrationValues1;
+  expectedFiltrationValues1.push_back(0);
+  expectedFiltrationValues1.push_back(0);
+  expectedFiltrationValues1.push_back(0);
+  expectedFiltrationValues1.push_back(0);
+  expectedFiltrationValues1.push_back(100);
+  expectedFiltrationValues1.push_back(0);
+  expectedFiltrationValues1.push_back(0);
+  expectedFiltrationValues1.push_back(0);
+  expectedFiltrationValues1.push_back(0);
+
+  std::vector< double > expectedFiltrationValues2;
+  expectedFiltrationValues2.push_back(1);
+  expectedFiltrationValues2.push_back(2);
+  expectedFiltrationValues2.push_back(3);
+  expectedFiltrationValues2.push_back(4);
+  expectedFiltrationValues2.push_back(5);
+  expectedFiltrationValues2.push_back(6);
+  expectedFiltrationValues2.push_back(7);
+  expectedFiltrationValues2.push_back(8);
+  expectedFiltrationValues2.push_back(9);
+
+  std::vector< double > increasingFiltrationOfTopDimensionalCells;
+  increasingFiltrationOfTopDimensionalCells.push_back(1);
+  increasingFiltrationOfTopDimensionalCells.push_back(2);
+  increasingFiltrationOfTopDimensionalCells.push_back(3);
+  increasingFiltrationOfTopDimensionalCells.push_back(4);
+  increasingFiltrationOfTopDimensionalCells.push_back(5);
+  increasingFiltrationOfTopDimensionalCells.push_back(6);
+  increasingFiltrationOfTopDimensionalCells.push_back(7);
+  increasingFiltrationOfTopDimensionalCells.push_back(8);
+  increasingFiltrationOfTopDimensionalCells.push_back(9);
+
+  std::vector< double > oneDimensionalCycle;
+  oneDimensionalCycle.push_back(0);
+  oneDimensionalCycle.push_back(0);
+  oneDimensionalCycle.push_back(0);
+  oneDimensionalCycle.push_back(0);
+  oneDimensionalCycle.push_back(100);
+  oneDimensionalCycle.push_back(0);
+  oneDimensionalCycle.push_back(0);
+  oneDimensionalCycle.push_back(0);
+  oneDimensionalCycle.push_back(0);
+
+  std::vector<unsigned> dimensions;
+  dimensions.push_back(3);
+  dimensions.push_back(3);
+
+  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+  Bitmap_cubical_complex<double> hole(dimensions, oneDimensionalCycle);
+
+
+  int i = 0;
+  for (Bitmap_cubical_complex<double>::Top_dimensional_cells_iterator 
+       it = increasing.top_dimensional_cells_begin(); it != increasing.top_dimensional_cells_end(); ++it) {
+    BOOST_CHECK(*it == expectedFiltrationValues2[i]);
+    ++i;
+  }
+  i = 0;
+  for (Bitmap_cubical_complex<double>::Top_dimensional_cells_iterator 
+       it = hole.top_dimensional_cells_begin(); it != hole.top_dimensional_cells_end(); ++it) {
+    BOOST_CHECK(*it == expectedFiltrationValues1[i]);
+    ++i;
+  }
+}
+
+BOOST_AUTO_TEST_CASE(compute_boundary_test_1) {
+
+  std::vector<double> boundary0;
+  std::vector<double> boundary1;
+  boundary1.push_back(0);
+  boundary1.push_back(2);
+  std::vector<double> boundary2;
+  std::vector<double> boundary3;
+  boundary3.push_back(2);
+  boundary3.push_back(4);
+  std::vector<double> boundary4;
+  std::vector<double> boundary5;
+  boundary5.push_back(4);
+  boundary5.push_back(6);
+  std::vector<double> boundary6;
+  std::vector<double> boundary7;
+  boundary7.push_back(0);
+  boundary7.push_back(14);
+  std::vector<double> boundary8;
+  boundary8.push_back(1);
+  boundary8.push_back(15);
+  boundary8.push_back(7);
+  boundary8.push_back(9);
+  std::vector<double> boundary9;
+  boundary9.push_back(2);
+  boundary9.push_back(16);
+  std::vector<double> boundary10;
+  boundary10.push_back(3);
+  boundary10.push_back(17);
+  boundary10.push_back(9);
+  boundary10.push_back(11);
+  std::vector<double> boundary11;
+  boundary11.push_back(4);
+  boundary11.push_back(18);
+  std::vector<double> boundary12;
+  boundary12.push_back(5);
+  boundary12.push_back(19);
+  boundary12.push_back(11);
+  boundary12.push_back(13);
+  std::vector<double> boundary13;
+  boundary13.push_back(6);
+  boundary13.push_back(20);
+  std::vector<double> boundary14;
+  std::vector<double> boundary15;
+  boundary15.push_back(14);
+  boundary15.push_back(16);
+  std::vector<double> boundary16;
+  std::vector<double> boundary17;
+  boundary17.push_back(16);
+  boundary17.push_back(18);
+  std::vector<double> boundary18;
+  std::vector<double> boundary19;
+  boundary19.push_back(18);
+  boundary19.push_back(20);
+  std::vector<double> boundary20;
+  std::vector<double> boundary21;
+  boundary21.push_back(14);
+  boundary21.push_back(28);
+  std::vector<double> boundary22;
+  boundary22.push_back(15);
+  boundary22.push_back(29);
+  boundary22.push_back(21);
+  boundary22.push_back(23);
+  std::vector<double> boundary23;
+  boundary23.push_back(16);
+  boundary23.push_back(30);
+  std::vector<double> boundary24;
+  boundary24.push_back(17);
+  boundary24.push_back(31);
+  boundary24.push_back(23);
+  boundary24.push_back(25);
+  std::vector<double> boundary25;
+  boundary25.push_back(18);
+  boundary25.push_back(32);
+  std::vector<double> boundary26;
+  boundary26.push_back(19);
+  boundary26.push_back(33);
+  boundary26.push_back(25);
+  boundary26.push_back(27);
+  std::vector<double> boundary27;
+  boundary27.push_back(20);
+  boundary27.push_back(34);
+  std::vector<double> boundary28;
+  std::vector<double> boundary29;
+  boundary29.push_back(28);
+  boundary29.push_back(30);
+  std::vector<double> boundary30;
+  std::vector<double> boundary31;
+  boundary31.push_back(30);
+  boundary31.push_back(32);
+  std::vector<double> boundary32;
+  std::vector<double> boundary33;
+  boundary33.push_back(32);
+  boundary33.push_back(34);
+  std::vector<double> boundary34;
+  std::vector<double> boundary35;
+  boundary35.push_back(28);
+  boundary35.push_back(42);
+  std::vector<double> boundary36;
+  boundary36.push_back(29);
+  boundary36.push_back(43);
+  boundary36.push_back(35);
+  boundary36.push_back(37);
+  std::vector<double> boundary37;
+  boundary37.push_back(30);
+  boundary37.push_back(44);
+  std::vector<double> boundary38;
+  boundary38.push_back(31);
+  boundary38.push_back(45);
+  boundary38.push_back(37);
+  boundary38.push_back(39);
+  std::vector<double> boundary39;
+  boundary39.push_back(32);
+  boundary39.push_back(46);
+  std::vector<double> boundary40;
+  boundary40.push_back(33);
+  boundary40.push_back(47);
+  boundary40.push_back(39);
+  boundary40.push_back(41);
+  std::vector<double> boundary41;
+  boundary41.push_back(34);
+  boundary41.push_back(48);
+  std::vector<double> boundary42;
+  std::vector<double> boundary43;
+  boundary43.push_back(42);
+  boundary43.push_back(44);
+  std::vector<double> boundary44;
+  std::vector<double> boundary45;
+  boundary45.push_back(44);
+  boundary45.push_back(46);
+  std::vector<double> boundary46;
+  std::vector<double> boundary47;
+  boundary47.push_back(46);
+  boundary47.push_back(48);
+  std::vector<double> boundary48;
+  std::vector< std::vector<double> > boundaries;
+  boundaries.push_back(boundary0);
+  boundaries.push_back(boundary1);
+  boundaries.push_back(boundary2);
+  boundaries.push_back(boundary3);
+  boundaries.push_back(boundary4);
+  boundaries.push_back(boundary5);
+  boundaries.push_back(boundary6);
+  boundaries.push_back(boundary7);
+  boundaries.push_back(boundary8);
+  boundaries.push_back(boundary9);
+  boundaries.push_back(boundary10);
+  boundaries.push_back(boundary11);
+  boundaries.push_back(boundary12);
+  boundaries.push_back(boundary13);
+  boundaries.push_back(boundary14);
+  boundaries.push_back(boundary15);
+  boundaries.push_back(boundary16);
+  boundaries.push_back(boundary17);
+  boundaries.push_back(boundary18);
+  boundaries.push_back(boundary19);
+  boundaries.push_back(boundary20);
+  boundaries.push_back(boundary21);
+  boundaries.push_back(boundary22);
+  boundaries.push_back(boundary23);
+  boundaries.push_back(boundary24);
+  boundaries.push_back(boundary25);
+  boundaries.push_back(boundary26);
+  boundaries.push_back(boundary27);
+  boundaries.push_back(boundary28);
+  boundaries.push_back(boundary29);
+  boundaries.push_back(boundary30);
+  boundaries.push_back(boundary31);
+  boundaries.push_back(boundary32);
+  boundaries.push_back(boundary33);
+  boundaries.push_back(boundary34);
+  boundaries.push_back(boundary35);
+  boundaries.push_back(boundary36);
+  boundaries.push_back(boundary37);
+  boundaries.push_back(boundary38);
+  boundaries.push_back(boundary39);
+  boundaries.push_back(boundary40);
+  boundaries.push_back(boundary41);
+  boundaries.push_back(boundary42);
+  boundaries.push_back(boundary43);
+  boundaries.push_back(boundary44);
+  boundaries.push_back(boundary45);
+  boundaries.push_back(boundary46);
+  boundaries.push_back(boundary47);
+  boundaries.push_back(boundary48);
+
+
+
+  std::vector< double > increasingFiltrationOfTopDimensionalCells;
+  increasingFiltrationOfTopDimensionalCells.push_back(1);
+  increasingFiltrationOfTopDimensionalCells.push_back(2);
+  increasingFiltrationOfTopDimensionalCells.push_back(3);
+  increasingFiltrationOfTopDimensionalCells.push_back(4);
+  increasingFiltrationOfTopDimensionalCells.push_back(5);
+  increasingFiltrationOfTopDimensionalCells.push_back(6);
+  increasingFiltrationOfTopDimensionalCells.push_back(7);
+  increasingFiltrationOfTopDimensionalCells.push_back(8);
+  increasingFiltrationOfTopDimensionalCells.push_back(9);
+
+  std::vector<unsigned> dimensions;
+  dimensions.push_back(3);
+  dimensions.push_back(3);
+
+  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+  for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) {
+    std::vector< size_t > bd = increasing.get_boundary_of_a_cell(i);
+    for (size_t j = 0; j != bd.size(); ++j) {
+      BOOST_CHECK(boundaries[i][j] == bd[j]);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(compute_boundary_test_2) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells;
+  increasingFiltrationOfTopDimensionalCells.push_back(1);
+  increasingFiltrationOfTopDimensionalCells.push_back(2);
+  increasingFiltrationOfTopDimensionalCells.push_back(3);
+  increasingFiltrationOfTopDimensionalCells.push_back(4);
+  increasingFiltrationOfTopDimensionalCells.push_back(5);
+  increasingFiltrationOfTopDimensionalCells.push_back(6);
+  increasingFiltrationOfTopDimensionalCells.push_back(7);
+  increasingFiltrationOfTopDimensionalCells.push_back(8);
+  increasingFiltrationOfTopDimensionalCells.push_back(9);
+
+  std::vector<unsigned> dimensions;
+  dimensions.push_back(3);
+  dimensions.push_back(3);
+
+  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+
+
+  std::vector<double> coboundaryElements;
+  coboundaryElements.push_back(7);
+  coboundaryElements.push_back(1);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(9);
+  coboundaryElements.push_back(1);
+  coboundaryElements.push_back(3);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(11);
+  coboundaryElements.push_back(3);
+  coboundaryElements.push_back(5);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(13);
+  coboundaryElements.push_back(5);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(7);
+  coboundaryElements.push_back(21);
+  coboundaryElements.push_back(15);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(9);
+  coboundaryElements.push_back(23);
+  coboundaryElements.push_back(15);
+  coboundaryElements.push_back(17);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(11);
+  coboundaryElements.push_back(25);
+  coboundaryElements.push_back(17);
+  coboundaryElements.push_back(19);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(13);
+  coboundaryElements.push_back(27);
+  coboundaryElements.push_back(19);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(21);
+  coboundaryElements.push_back(35);
+  coboundaryElements.push_back(29);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(23);
+  coboundaryElements.push_back(37);
+  coboundaryElements.push_back(29);
+  coboundaryElements.push_back(31);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(25);
+  coboundaryElements.push_back(39);
+  coboundaryElements.push_back(31);
+  coboundaryElements.push_back(33);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(27);
+  coboundaryElements.push_back(41);
+  coboundaryElements.push_back(33);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(35);
+  coboundaryElements.push_back(43);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(37);
+  coboundaryElements.push_back(43);
+  coboundaryElements.push_back(45);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(39);
+  coboundaryElements.push_back(45);
+  coboundaryElements.push_back(47);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(41);
+  coboundaryElements.push_back(47);
+  size_t number = 0;
+  for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) {
+    std::vector< size_t > bd = increasing.get_coboundary_of_a_cell(i);
+    for (size_t j = 0; j != bd.size(); ++j) {
+      BOOST_CHECK(coboundaryElements[number] == bd[j]);
+      ++number;
+    }
+
+  }
+}
+
+BOOST_AUTO_TEST_CASE(compute_boundary_test_3) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells;
+  increasingFiltrationOfTopDimensionalCells.push_back(1);
+  increasingFiltrationOfTopDimensionalCells.push_back(2);
+  increasingFiltrationOfTopDimensionalCells.push_back(3);
+  increasingFiltrationOfTopDimensionalCells.push_back(4);
+  increasingFiltrationOfTopDimensionalCells.push_back(5);
+  increasingFiltrationOfTopDimensionalCells.push_back(6);
+  increasingFiltrationOfTopDimensionalCells.push_back(7);
+  increasingFiltrationOfTopDimensionalCells.push_back(8);
+  increasingFiltrationOfTopDimensionalCells.push_back(9);
+
+  std::vector<unsigned> dimensions;
+  dimensions.push_back(3);
+  dimensions.push_back(3);
+
+  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+
+  std::vector<unsigned> dim;
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+
+  for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) {
+    BOOST_CHECK(increasing.get_dimension_of_a_cell(i) == dim[i]);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(Filtration_simplex_iterator_test) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells;
+  increasingFiltrationOfTopDimensionalCells.push_back(1);
+  increasingFiltrationOfTopDimensionalCells.push_back(2);
+  increasingFiltrationOfTopDimensionalCells.push_back(3);
+  increasingFiltrationOfTopDimensionalCells.push_back(4);
+  increasingFiltrationOfTopDimensionalCells.push_back(5);
+  increasingFiltrationOfTopDimensionalCells.push_back(6);
+  increasingFiltrationOfTopDimensionalCells.push_back(7);
+  increasingFiltrationOfTopDimensionalCells.push_back(8);
+  increasingFiltrationOfTopDimensionalCells.push_back(9);
+
+  std::vector<unsigned> dimensions;
+  dimensions.push_back(3);
+  dimensions.push_back(3);
+
+  Bitmap_cubical_complex<double> increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+
+  std::vector< unsigned > dim;
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+
+  std::vector<double> fil;
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(5);
+  fil.push_back(5);
+  fil.push_back(5);
+  fil.push_back(5);
+  fil.push_back(6);
+  fil.push_back(6);
+  fil.push_back(6);
+  fil.push_back(6);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(8);
+  fil.push_back(8);
+  fil.push_back(8);
+  fil.push_back(8);
+  fil.push_back(9);
+  fil.push_back(9);
+  fil.push_back(9);
+  fil.push_back(9);
+
+
+  Bitmap_cubical_complex<double>::Filtration_simplex_range range = increasing.filtration_simplex_range();
+  size_t position = 0;
+  for (Bitmap_cubical_complex<double>::Filtration_simplex_iterator it = range.begin(); it != range.end(); ++it) {
+    BOOST_CHECK(increasing.dimension(*it) == dim[position]);
+    BOOST_CHECK(increasing.filtration(*it) == fil[position]);
+    ++position;
+  }
+}