summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorvrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb>2015-09-25 14:57:29 +0000
committervrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb>2015-09-25 14:57:29 +0000
commit62937147e40a7d2da7aa7a7a604808feeccaa75e (patch)
tree8a3bfc688d4e71c3302bbba6c8508d1f5cf313fe
parent3fd88081e98379b4e94d01987166b9884fb7db4c (diff)
Add bitmap cubical complex
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bitmap@794 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 50d9b8eb80e0fe99f871afa5bdbb853add97e25e
-rw-r--r--CMakeLists.txt3
-rw-r--r--data/bitmap/CubicalOneSphere.txt12
-rw-r--r--data/bitmap/CubicalTwoSphere.txt31
-rw-r--r--src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp69
-rw-r--r--src/Bitmap_cubical_complex/example/CMakeLists.txt12
-rw-r--r--src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp85
-rw-r--r--src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h706
-rw-r--r--src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h577
-rw-r--r--src/Bitmap_cubical_complex/include/gudhi/counter.h136
-rw-r--r--src/Bitmap_cubical_complex/test/Bitmap_test.cpp623
-rw-r--r--src/Bitmap_cubical_complex/test/CMakeLists.txt25
-rw-r--r--src/CMakeLists.txt1
12 files changed, 2280 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6bea06e2..43e0558c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -57,6 +57,7 @@ else()
message(STATUS "boost library dirs:" ${Boost_LIBRARY_DIRS})
include_directories(src/common/include/)
+ include_directories(src/Bitmap_cubical_complex/include/)
include_directories(src/Bottleneck/include/)
include_directories(src/Contraction/include/)
include_directories(src/Hasse_complex/include/)
@@ -74,6 +75,8 @@ else()
add_subdirectory(src/Hasse_complex/example)
add_subdirectory(src/Bottleneck/example)
add_subdirectory(src/Bottleneck/test)
+ add_subdirectory(src/Bitmap_cubical_complex/example)
+ add_subdirectory(src/Bitmap_cubical_complex/test)
# GudhUI
add_subdirectory(src/GudhUI)
diff --git a/data/bitmap/CubicalOneSphere.txt b/data/bitmap/CubicalOneSphere.txt
new file mode 100644
index 00000000..4dec3fc6
--- /dev/null
+++ b/data/bitmap/CubicalOneSphere.txt
@@ -0,0 +1,12 @@
+2
+3
+3
+0
+0
+0
+0
+100
+0
+0
+0
+0
diff --git a/data/bitmap/CubicalTwoSphere.txt b/data/bitmap/CubicalTwoSphere.txt
new file mode 100644
index 00000000..32f5d933
--- /dev/null
+++ b/data/bitmap/CubicalTwoSphere.txt
@@ -0,0 +1,31 @@
+3
+3
+3
+3
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+100
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
new file mode 100644
index 00000000..c0dbaf36
--- /dev/null
+++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
@@ -0,0 +1,69 @@
+/* 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>
+
+using namespace Gudhi;
+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 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); //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
new file mode 100644
index 00000000..05ef1319
--- /dev/null
+++ b/src/Bitmap_cubical_complex/example/CMakeLists.txt
@@ -0,0 +1,12 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHISimplexTreeFromFile)
+
+add_executable ( Bitmap_cubical_complex Bitmap_cubical_complex.cpp )
+target_link_libraries(Bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY})
+add_test(Bitmap_cubical_complex ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt)
+add_test(Bitmap_cubical_complex ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt)
+
+add_executable ( Random_bitmap_cubical_complex Random_bitmap_cubical_complex.cpp )
+target_link_libraries(Random_bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY})
+add_test(Random_bitmap_cubical_complex ${CMAKE_CURRENT_BINARY_DIR}/Random_bitmap_cubical_complex 2 100 100)
+
diff --git a/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp
new file mode 100644
index 00000000..de9d96e0
--- /dev/null
+++ b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp
@@ -0,0 +1,85 @@
+/* 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>
+
+using namespace Gudhi;
+using namespace Gudhi::persistent_cohomology;
+
+//standard stuff
+#include <iostream>
+#include <sstream>
+#include <vector>
+#include <cstdlib>
+#include <ctime>
+
+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); //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
new file mode 100644
index 00000000..61ae8105
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
@@ -0,0 +1,706 @@
+/* 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 <limits>
+
+#include <gudhi/Bitmap_cubical_complex_base.h>
+
+//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) {
+ 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 -- 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
+ //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(char* perseusStyleFile) : Bitmap_cubical_complex_base<T>(perseusStyleFile) {
+ if (globalDbg) {
+ 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) {
+ 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) {
+ 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) {
+ 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->keyAssociatedToSimplex[ 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->simplexAssociatedToKey[ 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->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) {
+ 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) {
+ 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->boundaryElements.clear();
+ this->boundaryElementsinsert(this->boundaryElements.end(), rhs.boundaryElements.begin(), rhs.boundaryElements.end());
+ }
+
+ bool operator==(const Boundary_simplex_iterator& rhs) {
+ if (globalDbg) {
+ 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) {
+ cerr << "bool operator != \n";
+ }
+ return !(*this == rhs);
+ }
+
+ Simplex_handle operator*() {
+ if (globalDbg) {
+ 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) {
+ 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.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) {
+ 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->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) {
+ 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->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) {
+ 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 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) {
+ 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;
+ 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) {
+ 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;
+ int dimension;
+ };
+
+ /**
+ * Function needed for compatibility with Gudhi. Not useful for other purposes.
+ **/
+ Skeleton_simplex_range skeleton_simplex_range(int dimension) {
+ if (globalDbg) {
+ 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) {
+ cerr << i << " -> " << this->simplexAssociatedToKey[i] << 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;
+ std::vector< size_t > elementsOrderedAccordingToFiltration; //needed by Filtration_simplex_iterator. If this iterator is not used, this field is not initialized.
+}; //Bitmap_cubical_complex
+
+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) {
+ 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) {
+ 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 -- 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>);
+
+ std::vector< size_t > elementsOfBitmapOrderedAccordingToFiltrationThenAccordingToDimensionThenAccordingToPositionInBitmap(this->data.size());
+ for (size_t i = 0; i != dataOfElementsFromBitmap.size(); ++i) {
+ elementsOfBitmapOrderedAccordingToFiltrationThenAccordingToDimensionThenAccordingToPositionInBitmap[i] = dataOfElementsFromBitmap[i].first;
+ }
+ this->elementsOrderedAccordingToFiltration = elementsOfBitmapOrderedAccordingToFiltrationThenAccordingToDimensionThenAccordingToPositionInBitmap;
+}
+
+
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+//****************************************************************************************************************//
+
+#endif // BITMAP_CUBICAL_COMPLEX_H_
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
new file mode 100644
index 00000000..26c97872
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -0,0 +1,577 @@
+/* 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 <iostream>
+#include <string>
+#include <vector>
+#include <string>
+#include <cstdlib>
+#include <climits>
+#include <fstream>
+#include <algorithm>
+#include <iterator>
+
+#include <gudhi/counter.h>
+
+using namespace std;
+
+/**
+ * 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 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_) {
+ 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) {
+ 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>
+ostream& operator<<(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()) {
+ 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." << 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;
+ ifstream inFiltration, inIds;
+ inFiltration.open(perseusStyleFile_);
+ unsigned dimensionOfData;
+ inFiltration >> dimensionOfData;
+
+ if (dbg) {
+ cerr << "dimensionOfData : " << dimensionOfData << 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) {
+ cerr << "sizeInThisDimension : " << sizeInThisDimension << 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.computeIndexInBitmap() << " and dimension: " << this->get_dimension_of_a_cell(it.computeIndexInBitmap()) << " 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_) {
+ 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) {
+ cerr << "dimensionsInWhichCellHasNonzeroLength : \n";
+ for (size_t i = 0; i != dimensionsInWhichCellHasNonzeroLength.size(); ++i) {
+ cerr << dimensionsInWhichCellHasNonzeroLength[i] << 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) cerr << "multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " << multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << endl;
+ if (bdg) cerr << "cell_ - multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " << cell_ - multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << endl;
+ if (bdg) cerr << "cell_ + multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " << cell_ + multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << 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) {
+ cerr << "dimensionsInWhichCellHasZeroLength : \n";
+ for (size_t i = 0; i != dimensionsInWhichCellHasZeroLength.size(); ++i) {
+ cerr << dimensionsInWhichCellHasZeroLength[i] << endl;
+ }
+ cerr << "\n counter : " << endl;
+ for (size_t i = 0; i != counter.size(); ++i) {
+ cerr << counter[i] << endl;
+ }
+ getchar();
+ }
+
+ std::vector< size_t > coboundaryElements;
+ if (dimensionsInWhichCellHasZeroLength.size() == 0)return coboundaryElements;
+ for (size_t i = 0; i != dimensionsInWhichCellHasZeroLength.size(); ++i) {
+ if (bdg) {
+ cerr << "Dimension : " << i << endl;
+ if (counter[dimensionsInWhichCellHasZeroLength[i]] == 0) {
+ 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]]) {
+ 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)cerr << "Subtracting : " << cell_ - multipliers[dimensionsInWhichCellHasZeroLength[i]] << 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)cerr << "Adding : " << cell_ + multipliers[dimensionsInWhichCellHasZeroLength[i]] << 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)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> 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) {
+ cerr << "indicesToConsider in this iteration \n";
+ for (size_t i = 0; i != indicesToConsider.size(); ++i) {
+ 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) {
+ cerr << " position : " << position << endl;
+ cerr << c << endl;
+ getchar();
+ }
+
+ c.increment();
+ }
+
+ return result;
+}
+
+#endif // BITMAP_CUBICAL_COMPLEX_BASE_H_
diff --git a/src/Bitmap_cubical_complex/include/gudhi/counter.h b/src/Bitmap_cubical_complex/include/gudhi/counter.h
new file mode 100644
index 00000000..9df819b2
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/counter.h
@@ -0,0 +1,136 @@
+/* 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>
+
+using namespace std;
+
+/**
+ * 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 endd, 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 end 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) {
+ //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< int > begin;
+ std::vector< int > end;
+ std::vector< int > current;
+};
+
+#endif // COUNTER_H_
diff --git a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
new file mode 100644
index 00000000..1c204bae
--- /dev/null
+++ b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
@@ -0,0 +1,623 @@
+#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;
+ }
+}
diff --git a/src/Bitmap_cubical_complex/test/CMakeLists.txt b/src/Bitmap_cubical_complex/test/CMakeLists.txt
new file mode 100644
index 00000000..97c374e6
--- /dev/null
+++ b/src/Bitmap_cubical_complex/test/CMakeLists.txt
@@ -0,0 +1,25 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBitmapCCUT)
+
+if (GCOVR_PATH)
+ # for gcovr to make coverage reports - Corbera Jenkins plugin
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage")
+ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage")
+ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage")
+endif()
+if (GPROF_PATH)
+ # for gprof to make coverage reports - Jenkins
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg")
+ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg")
+ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg")
+endif()
+
+add_executable ( BitmapCCUT Bitmap_test.cpp )
+target_link_libraries(BitmapCCUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
+# Unitary tests
+add_test(NAME BitmapCCUT
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BitmapCCUT
+ # XML format for Jenkins xUnit plugin
+ --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BitmapCCUT.xml --log_level=test_suite --report_level=no)
+
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 545b0b58..576ba353 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -44,6 +44,7 @@ else()
add_subdirectory(example/Contraction)
add_subdirectory(example/Hasse_complex)
add_subdirectory(example/Bottleneck)
+ add_subdirectory(example/Bitmap_cubical_complex)
# GudhUI
add_subdirectory(GudhUI)