diff options
38 files changed, 4939 insertions, 14 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index dc6de64d..ecd1f7fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ project(GUDHIdev) include(CMakeGUDHIVersion.txt) # Generate GUDHI official version file -configure_file(GUDHIVersion.cmake.in "${PROJECT_BINARY_DIR}/GUDHIVersion.cmake" @ONLY) +configure_file(GUDHIVersion.cmake.in "${CMAKE_SOURCE_DIR}/GUDHIVersion.cmake" @ONLY) find_package(Boost REQUIRED COMPONENTS system filesystem unit_test_framework chrono timer program_options thread REQUIRED) @@ -93,12 +93,14 @@ endif() include_directories(src/common/include/) include_directories(src/Alpha_complex/include/) + include_directories(src/Bitmap_cubical_complex/include/) include_directories(src/Bottleneck/include/) include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) include_directories(src/Persistent_cohomology/include/) include_directories(src/Simplex_tree/include/) include_directories(src/Skeleton_blocker/include/) + include_directories(src/Witness_complex/include/) add_subdirectory(src/common/example) add_subdirectory(src/common/test) @@ -109,6 +111,10 @@ endif() add_subdirectory(src/Skeleton_blocker/test) add_subdirectory(src/Skeleton_blocker/example) add_subdirectory(src/Contraction/example) + add_subdirectory(src/Witness_complex/test) + add_subdirectory(src/Witness_complex/example) + add_subdirectory(src/Bitmap_cubical_complex/test) + add_subdirectory(src/Bitmap_cubical_complex/example) add_subdirectory(src/Alpha_complex/example) add_subdirectory(src/Alpha_complex/test) diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index 3fd1c10a..9fc01a5d 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -927,3 +927,30 @@ misc{jplex_cite, year = "2009", note = "http://comptop.stanford.edu/programs/jplex/" } + +@book{kaczynski2004computational, + title={Computational Homology}, + author={Kaczynski, T. and Mischaikow, K. and Mrozek, M.}, + isbn={9780387408538}, + lccn={03061109}, + series={Applied Mathematical Sciences}, + url={https://books.google.fr/books?id=AShKtpi3GecC}, + year={2004}, + publisher={Springer New York} +} + +@ARTICLE{peikert2012topological, +year={2012}, +isbn={978-3-642-23174-2}, +booktitle={Topological Methods in Data Analysis and Visualization II}, +series={Mathematics and Visualization}, +editor={Peikert, Ronald and Hauser, Helwig and Carr, Hamish and Fuchs, Raphael}, +doi={10.1007/978-3-642-23175-9_7}, +title={Efficient Computation of Persistent Homology for Cubical Data}, +url={http://dx.doi.org/10.1007/978-3-642-23175-9_7}, +publisher={Springer Berlin Heidelberg}, +author={Wagner, Hubert and Chen, Chao and Vucini, Erald}, +pages={91-106}, +language={English} +} + diff --git a/data/bitmap/2d_torus.txt b/data/bitmap/2d_torus.txt new file mode 100644 index 00000000..8ccd45c6 --- /dev/null +++ b/data/bitmap/2d_torus.txt @@ -0,0 +1,12 @@ +2
+-3
+-3
+0
+0
+0
+0
+10
+0
+0
+0
+0
diff --git a/data/bitmap/3d_torus.txt b/data/bitmap/3d_torus.txt new file mode 100644 index 00000000..896e6c61 --- /dev/null +++ b/data/bitmap/3d_torus.txt @@ -0,0 +1,31 @@ +3
+-3
+-3
+-3
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+10
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
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/scripts/generate_version.sh b/scripts/generate_version.sh index 323396dc..6f66a638 100755 --- a/scripts/generate_version.sh +++ b/scripts/generate_version.sh @@ -44,8 +44,6 @@ echo $VERSION_DIR mkdir "$VERSION_DIR" # TOP LEVEL FILE COPY -cp $VERSION_FILE $VERSION_DIR -cp $ROOT_DIR/CMakeGUDHIVersion.txt $VERSION_DIR cp $ROOT_DIR/README $VERSION_DIR cp $ROOT_DIR/Conventions.txt $VERSION_DIR cp $ROOT_DIR/COPYING $VERSION_DIR @@ -55,11 +53,11 @@ cp $ROOT_DIR/src/Doxyfile $VERSION_DIR cp -R $ROOT_DIR/biblio $VERSION_DIR cp $ROOT_DIR/src/GUDHIConfigVersion.cmake.in $VERSION_DIR cp $ROOT_DIR/src/GUDHIConfig.cmake.in $VERSION_DIR +cp $ROOT_DIR/CMakeGUDHIVersion.txt $VERSION_DIR cp $ROOT_DIR/GUDHIVersion.cmake.in $VERSION_DIR # PACKAGE LEVEL COPY PACKAGE_INC_DIR="/include" -#PACKAGE_SRC_DIR="/source" PACKAGE_EX_DIR="/example" PACKAGE_CONCEPT_DIR="/concept" PACKAGE_DOC_DIR="/doc" diff --git a/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h new file mode 100644 index 00000000..cde0b2fc --- /dev/null +++ b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h @@ -0,0 +1,156 @@ +/* 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 DOC_GUDHI_CUBICAL_COMPLEX_COMPLEX_H_
+#define DOC_GUDHI_CUBICAL_COMPLEX_COMPLEX_H_
+
+namespace Gudhi {
+
+namespace Cubical_complex {
+
+/** \defgroup cubical_complex Cubical complex
+ *
+ * \author Pawel Dlotko
+ *
+ * @{
+ *
+
+ * Bitmap_cubical_complex is an example of a structured complex useful in computational mathematics (specially rigorous
+ * numerics) and image analysis. The presented implementation of cubical complexes is based on the following
+ * definition.
+ *
+ * An <em>elementary interval</em> is an interval of a form \f$ [n,n+1] \f$, or \f$[n,n]\f$, for \f$ n \in \mathcal{Z}
+ * \f$. The first one is called <em>non-degenerate</em>, while the second one is \a degenerate interval. A
+ * <em>boundary of a elementary interval</em> is a chain \f$\partial [n,n+1] = [n+1,n+1]-[n,n] \f$ in case of
+ * non-degenerated elementary interval and \f$\partial [n,n] = 0 \f$ in case of degenerate elementary interval. An
+ * <em>elementary cube</em> \f$ C \f$ is a product of elementary intervals, \f$C=I_1 \times \ldots \times I_n\f$.
+ * <em>Embedding dimension</em> of a cube is n, the number of elementary intervals (degenerate or not) in the product.
+ * A <em>dimension of a cube</em> \f$C=I_1 \times ... \times I_n\f$ is the number of non degenerate elementary
+ * intervals in the product. A <em>boundary of a cube</em> \f$C=I_1 \times \ldots \times I_n\f$ is a chain obtained
+ * in the following way:
+ * \f[\partial C = (\partial I_1 \times \ldots \times I_n) + (I_1 \times \partial I_2 \times \ldots \times I_n) +
+ * \ldots + (I_1 \times I_2 \times \ldots \times \partial I_n).\f]
+ * A <em>cubical complex</em> \f$\mathcal{K}\f$ is a collection of cubes closed under operation of taking boundary
+ * (i.e. boundary of every cube from the collection is in the collection). A cube \f$C\f$ in cubical complex
+ * \f$\mathcal{K}\f$ is <em>maximal</em> if it is not in a boundary of any other cube in \f$\mathcal{K}\f$. A \a
+ * support of a cube \f$C\f$ is the set in \f$\mathbb{R}^n\f$ occupied by \f$C\f$ (\f$n\f$ is the embedding dimension
+ * of \f$C\f$).
+ *
+ * Cubes may be equipped with a filtration values in which case we have filtered cubical complex. All the cubical
+ * complexes considered in this implementation are filtered cubical complexes (although, the range of a filtration may
+ * be a set of two elements).
+ *
+ * For further details and theory of cubical complexes, please consult \cite kaczynski2004computational as well as the
+ * following paper \cite peikert2012topological .
+ *
+ * \section datastructure Data structure.
+ *
+ * The implementation of Cubical complex provides a representation of complexes that occupy a rectangular region in
+ * \f$\mathbb{R}^n\f$. This extra assumption allows for a memory efficient way of storing cubical complexes in a form
+ * of so called bitmaps. Let \f$R = [b_1,e_1] \times \ldots \times [b_n,e_n]\f$, for \f$b_1,...b_n,e_1,...,e_n \in
+ * \mathbb{Z}\f$, \f$b_i \leq d_i\f$ be the considered rectangular region and let \f$\mathcal{K}\f$ be a filtered
+ * cubical complex having the rectangle \f$R\f$ as its support. Note that the structure of the coordinate system gives
+ * a way a lexicographical ordering of cells of \f$\mathcal{K}\f$. This ordering is a base of the presented
+ * bitmap-based implementation. In this implementation, the whole cubical complex is stored as a vector of the values
+ * of filtration. This, together with dimension of \f$\mathcal{K}\f$ and the sizes of \f$\mathcal{K}\f$ in all
+ * directions, allows to determine, dimension, neighborhood, boundary and coboundary of every cube \f$C \in
+ * \mathcal{K}\f$.
+ *
+ * \image html "bitmapAllCubes.png" "Cubical complex.
+ *
+ * Note that the cubical complex in the figure above is, in a natural way, a product of one dimensional cubical
+ * complexes in \f$\mathbb{R}\f$. The number of all cubes in each direction is equal \f$2n+1\f$, where \f$n\f$ is the
+ * number of maximal cubes in the considered direction. Let us consider a cube at the position \f$k\f$ in the bitmap.
+ * Knowing the sizes of the bitmap, by a series of modulo operation, we can determine which elementary intervals are
+ * present in the product that gives the cube \f$C\f$. In a similar way, we can compute boundary and the coboundary of
+ * each cube. Further details can be found in the literature.
+ *
+ * \section inputformat Input Format.
+ *
+ * In the current implantation, filtration is given at the maximal cubes, and it is then extended by the lower star
+ * filtration to all cubes. There are a number of constructors that can be used to construct cubical complex by users
+ * who want to use the code directly. They can be found in the \a Bitmap_cubical_complex class.
+ * Currently one input from a text file is used. It uses a format used already in Perseus software
+ * (http://www.sas.upenn.edu/~vnanda/perseus/) by Vidit Nanda.
+ * Below we are providing a description of the format. The first line contains a number d begin the dimension of the
+ * bitmap (2 in the example below). Next d lines are the numbers of top dimensional cubes in each dimensions (3 and 3
+ * in the example below). Next, in lexicographical order, the filtration of top dimensional cubes is given (1 4 6 8
+ * 20 4 7 6 5 in the example below).
+ *
+ *
+ * \image html "exampleBitmap.png" "Example of a input data."
+ *
+ * The input file for the following complex is:
+ * \verbatim
+2
+3
+3
+1
+4
+6
+8
+20
+4
+7
+6
+5
+\endverbatim
+
+ * \section PeriodicBoundaryConditions Periodic boundary conditions
+ * Often one would like to impose periodic boundary conditions to the cubical complex. Let \f$ I_1\times ... \times
+ * I_n \f$ be a box that is decomposed with a cubical complex \f$ \mathcal{K} \f$. Imposing periodic boundary
+ * conditions in the direction i, means that the left and the right side of a complex \f$ \mathcal{K} \f$ are
+ * considered the same. In particular, if for a bitmap \f$ \mathcal{K} \f$ periodic boundary conditions are imposed
+ * in all directions, then complex \f$ \mathcal{K} \f$ became n-dimensional torus. One can use various constructors
+ * from the file Bitmap_cubical_complex_periodic_boundary_conditions_base.h to construct cubical complex with periodic
+ * boundary conditions. One can also use Perseus style input files. To indicate periodic boundary conditions in a
+ * given direction, then number of top dimensional cells in this direction have to be multiplied by -1. For instance:
+
+ *\verbatim
+2
+-3
+3
+1
+4
+6
+8
+20
+4
+7
+6
+5
+\endverbatim
+
+ * Indicate that we have imposed periodic boundary conditions in the direction x, but not in the direction y.
+
+ * \section BitmapExamples Examples
+ * End user programs are available in example/Bitmap_cubical_complex folder.
+
+ */
+/** @} */ // end defgroup cubical_complex
+
+} // namespace Cubical_complex
+
+} // namespace Gudhi
+
+#endif // DOC_GUDHI_CUBICAL_COMPLEX_COMPLEX_H_
diff --git a/src/Bitmap_cubical_complex/doc/bitmapAllCubes.png b/src/Bitmap_cubical_complex/doc/bitmapAllCubes.png Binary files differnew file mode 100644 index 00000000..77167b13 --- /dev/null +++ b/src/Bitmap_cubical_complex/doc/bitmapAllCubes.png diff --git a/src/Bitmap_cubical_complex/doc/exampleBitmap.png b/src/Bitmap_cubical_complex/doc/exampleBitmap.png Binary files differnew file mode 100644 index 00000000..069c6eb2 --- /dev/null +++ b/src/Bitmap_cubical_complex/doc/exampleBitmap.png 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..7c7e8ac5 --- /dev/null +++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp @@ -0,0 +1,72 @@ +/* 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/>. + */ + + +#include <gudhi/reader_utils.h> +#include <gudhi/Bitmap_cubical_complex.h> +#include <gudhi/Persistent_cohomology.h> + +// standard stuff +#include <iostream> +#include <sstream> +#include <vector> + +int main(int argc, char** argv) { + std::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.\n" << std::endl; + + int p = 2; + double min_persistence = 0; + + if (argc != 2) { + std::cerr << "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; + } + + typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base; + typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex; + typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp; + typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology; + + Bitmap_cubical_complex b(argv[1]); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(b); + pcoh.init_coefficients(p); // initializes the coefficient field for homology + pcoh.compute_persistent_cohomology(min_persistence); + + std::stringstream ss; + ss << argv[1] << "_persistence"; + std::ofstream out(ss.str().c_str()); + pcoh.output_diagram(out); + out.close(); + + std::cout << "Result in file: " << ss.str().c_str() << "\n"; + + return 0; +} + diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp new file mode 100644 index 00000000..2c5d7fd3 --- /dev/null +++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp @@ -0,0 +1,74 @@ +/* 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/>. + */ + + +#include <gudhi/reader_utils.h> +#include <gudhi/Bitmap_cubical_complex.h> +#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h> +#include <gudhi/Persistent_cohomology.h> + +// standard stuff +#include <iostream> +#include <sstream> +#include <vector> + +int main(int argc, char** argv) { + std::cout << "This program computes persistent homology, by using " << + "Bitmap_cubical_complex_periodic_boundary_conditions 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.\n" << std::endl; + + int p = 2; + double min_persistence = 0; + + if (argc != 2) { + std::cerr << "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; + } + + typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double> Bitmap_base; + typedef Gudhi::Cubical_complex::Bitmap_cubical_complex< Bitmap_base > Bitmap_cubical_complex; + + Bitmap_cubical_complex b(argv[1]); + + typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp; + typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology; + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(b, true); + pcoh.init_coefficients(p); // initializes the coefficient field for homology + pcoh.compute_persistent_cohomology(min_persistence); + + std::stringstream ss; + ss << argv[1] << "_persistence"; + std::ofstream out(ss.str().c_str()); + pcoh.output_diagram(out); + out.close(); + + std::cout << "Result in file: " << ss.str().c_str() << "\n"; + + 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..8f9cfa80 --- /dev/null +++ b/src/Bitmap_cubical_complex/example/CMakeLists.txt @@ -0,0 +1,17 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIBitmap) + +add_executable ( Bitmap_cubical_complex Bitmap_cubical_complex.cpp ) +target_link_libraries(Bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY}) +add_test(Bitmap_cubical_complex_one_sphere ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt) +add_test(Bitmap_cubical_complex_two_sphere ${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) + +add_executable ( Bitmap_cubical_complex_periodic_boundary_conditions Bitmap_cubical_complex_periodic_boundary_conditions.cpp ) +target_link_libraries(Bitmap_cubical_complex_periodic_boundary_conditions ${Boost_SYSTEM_LIBRARY}) +add_test(Bitmap_cubical_complex_periodic_2d_torus ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex_periodic_boundary_conditions ${CMAKE_SOURCE_DIR}/data/bitmap/2d_torus.txt) +add_test(Bitmap_cubical_complex_periodic_3d_torus ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex_periodic_boundary_conditions ${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt) + 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..416ad8f2 --- /dev/null +++ b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp @@ -0,0 +1,83 @@ +/* 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> + +// standard stuff +#include <iostream> +#include <sstream> +#include <vector> + +int main(int argc, char** argv) { + srand(time(0)); + + std::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." << std::endl; + + int p = 2; + double min_persistence = 0; + + if (argc < 3) { + std::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() / static_cast<double>(RAND_MAX)); + } + + typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base; + typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex; + Bitmap_cubical_complex b(sizes, data); + + // Compute the persistence diagram of the complex + typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp; + typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology; + Persistent_cohomology pcoh(b); + pcoh.init_coefficients(p); // initializes the coefficient field for homology + pcoh.compute_persistent_cohomology(min_persistence); + + std::stringstream ss; + ss << "randomComplex_persistence"; + std::ofstream out(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..67e1fed3 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -0,0 +1,596 @@ +/* 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 <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h> + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_sort.h> +#endif + +#include <limits> +#include <utility> // for pair<> +#include <algorithm> // for sort +#include <vector> +#include <numeric> // for iota + +namespace Gudhi { + +namespace Cubical_complex { + +// global variable, was used just for debugging. +const bool globalDbg = false; + +template <typename T> class is_before_in_filtration; + +/** + * This is a Bitmap_cubical_complex class. It joints a functionalities of Bitmap_cubical_complex_base and + * Bitmap_cubical_complex_periodic_boundary_conditions_base classes into + * Gudhi persistent homology engine. It is a template class that inherit from its template parameter. The template + * parameter is supposed to be either Bitmap_cubical_complex_base or + * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + **/ + +/** + *@class Bitmap_cubical_complex + *@brief Cubical complex represented as a bitmap. + *@ingroup cubical_complex + */ +template <typename T> +class Bitmap_cubical_complex : public T { + public: + //*********************************************// + // Typedefs and typenames + //*********************************************// + typedef size_t Simplex_key; + typedef typename T::filtration_type Filtration_value; + typedef Simplex_key Simplex_handle; + + + //*********************************************// + // Constructors + //*********************************************// + // Over here we need to define various input types. I am proposing the following ones: + // Perseus style + // TODO(PD) H5 files? + // TODO(PD) binary files with little endiangs / big endians ? + // TODO(PD) constructor from a vector of elements of a type T. ? + + /** + * Constructor form a Perseus-style file. + **/ + Bitmap_cubical_complex(const char* perseus_style_file) : + T(perseus_style_file), key_associated_to_simplex(this->total_number_of_cells + 1) { + if (globalDbg) { + std::cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n"; + } + for (size_t i = 0; i != this->total_number_of_cells; ++i) { + this->key_associated_to_simplex[i] = i; + } + // we initialize this only once, in each constructor, when the bitmap is constructed. + // If the user decide to change some elements of the bitmap, then this procedure need + // to be called again. + this->initialize_simplex_associated_to_key(); + } + + /** + * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells + * in the following directions and vector of element of a type T + * with filtration on top dimensional cells. + **/ + Bitmap_cubical_complex(const std::vector<unsigned>& dimensions, + const std::vector<typename T::filtration_type>& top_dimensional_cells) : + T(dimensions, top_dimensional_cells), + key_associated_to_simplex(this->total_number_of_cells + 1) { + for (size_t i = 0; i != this->total_number_of_cells; ++i) { + this->key_associated_to_simplex[i] = i; + } + // we initialize this only once, in each constructor, when the bitmap is constructed. + // If the user decide to change some elements of the bitmap, then this procedure need + // to be called again. + this->initialize_simplex_associated_to_key(); + } + + /** + * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells + * in the following directions and vector of element of a type T::filtration_type + * with filtration on top dimensional cells. The last parameter of the constructor is a vector of boolean of a length + * equal to the dimension of cubical complex. + * If the position i on this vector is true, then we impose periodic boundary conditions in this direction. + **/ + Bitmap_cubical_complex(const std::vector<unsigned>& dimensions, + const std::vector<typename T::filtration_type>& top_dimensional_cells, + std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed) : + T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed), + key_associated_to_simplex(this->total_number_of_cells + 1) { + for (size_t i = 0; i != this->total_number_of_cells; ++i) { + this->key_associated_to_simplex[i] = i; + } + // we initialize this only once, in each constructor, when the bitmap is constructed. + // If the user decide to change some elements of the bitmap, then this procedure need + // to be called again. + this->initialize_simplex_associated_to_key(); + } + + /** + * Destructor of the Bitmap_cubical_complex class. + **/ + virtual ~Bitmap_cubical_complex() {} + + //*********************************************// + // Other 'easy' functions + //*********************************************// + + /** + * Returns number of all cubes in the complex. + **/ + size_t num_simplices()const { + return this->total_number_of_cells; + } + + /** + * Returns a Simplex_handle to a cube that do not exist in this complex. + **/ + static Simplex_handle null_simplex() { + if (globalDbg) { + std::cerr << "Simplex_handle null_simplex()\n"; + } + return std::numeric_limits<Simplex_handle>::max(); + } + + /** + * Returns dimension of the complex. + **/ + inline size_t dimension()const { + return this->sizes.size(); + } + + /** + * Return dimension of a cell pointed by the Simplex_handle. + **/ + inline unsigned dimension(Simplex_handle sh)const { + if (globalDbg) { + std::cerr << "unsigned dimension(const Simplex_handle& sh)\n"; + } + if (sh != std::numeric_limits<Simplex_handle>::max()) return this->get_dimension_of_a_cell(sh); + return -1; + } + + /** + * Return the filtration of a cell pointed by the Simplex_handle. + **/ + typename T::filtration_type filtration(Simplex_handle sh) { + if (globalDbg) { + std::cerr << "T::filtration_type filtration(const Simplex_handle& sh)\n"; + } + // Returns the filtration value of a simplex. + if (sh != std::numeric_limits<Simplex_handle>::max()) return this->data[sh]; + return std::numeric_limits<Simplex_handle>::max(); + } + + /** + * Return a key which is not a key of any cube in the considered data structure. + **/ + static Simplex_key null_key() { + if (globalDbg) { + std::cerr << "Simplex_key null_key()\n"; + } + return std::numeric_limits<Simplex_handle>::max(); + } + + /** + * Return the key of a cube pointed by the Simplex_handle. + **/ + Simplex_key key(Simplex_handle sh)const { + if (globalDbg) { + std::cerr << "Simplex_key key(const Simplex_handle& sh)\n"; + } + if (sh != std::numeric_limits<Simplex_handle>::max()) { + return this->key_associated_to_simplex[sh]; + } + return this->null_key(); + } + + /** + * Return the Simplex_handle given the key of the cube. + **/ + Simplex_handle simplex(Simplex_key key) { + if (globalDbg) { + std::cerr << "Simplex_handle simplex(Simplex_key key)\n"; + } + if (key != std::numeric_limits<Simplex_handle>::max()) { + return this->simplex_associated_to_key[ key ]; + } + return null_simplex(); + } + + /** + * Assign key to a cube pointed by the Simplex_handle + **/ + void assign_key(Simplex_handle sh, Simplex_key key) { + if (globalDbg) { + std::cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n"; + } + if (key == std::numeric_limits<Simplex_handle>::max()) return; + this->key_associated_to_simplex[sh] = key; + this->simplex_associated_to_key[key] = sh; + } + + /** + * Function called from a constructor. It is needed for Filtration_simplex_iterator to work. + **/ + void initialize_simplex_associated_to_key(); + + //*********************************************// + // Iterators + //*********************************************// + + /** + * Boundary_simplex_range class provides ranges for boundary iterators. + **/ + typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator; + typedef typename std::vector< Simplex_handle > Boundary_simplex_range; + + /** + * Filtration_simplex_iterator class provides an iterator though the whole structure in the order of filtration. + * Secondary criteria for filtration are: + * (1) Dimension of a cube (lower dimensional comes first). + * (2) Position in the data structure (the ones that are earlies in the data structure comes first). + **/ + class Filtration_simplex_range; + + class Filtration_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > { + // Iterator over all simplices of the complex in the order of the indexing scheme. + // 'value_type' must be 'Simplex_handle'. + public: + Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) { } + + Filtration_simplex_iterator() : b(NULL), position(0) { } + + 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; + return (*this); + } + + bool operator==(const Filtration_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n"; + } + return ( this->position == rhs.position); + } + + bool operator!=(const Filtration_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator != ( const Filtration_simplex_iterator& rhs )\n"; + } + return !(*this == rhs); + } + + Simplex_handle operator*() { + if (globalDbg) { + std::cerr << "Simplex_handle operator*()\n"; + } + return this->b->simplex_associated_to_key[ this->position ]; + } + + friend class Filtration_simplex_range; + + private: + Bitmap_cubical_complex<T>* b; + size_t position; + }; + + /** + * Filtration_simplex_range provides the ranges for Filtration_simplex_iterator. + **/ + class Filtration_simplex_range { + // Range over the simplices of the complex in the order of the filtration. + // .begin() and .end() return type Filtration_simplex_iterator. + public: + typedef Filtration_simplex_iterator const_iterator; + typedef Filtration_simplex_iterator iterator; + + Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) { } + + Filtration_simplex_iterator begin() { + if (globalDbg) { + std::cerr << "Filtration_simplex_iterator begin() \n"; + } + 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->simplex_associated_to_key.size(); + return it; + } + + private: + Bitmap_cubical_complex<T>* b; + }; + + + + //*********************************************// + // Methods to access iterators from the container: + + /** + * boundary_simplex_range creates an object of a Boundary_simplex_range class + * that provides ranges for the Boundary_simplex_iterator. + **/ + Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { + return this->get_boundary_of_a_cell(sh); + } + + /** + * filtration_simplex_range creates an object of a Filtration_simplex_range class + * that provides ranges for the Filtration_simplex_iterator. + **/ + Filtration_simplex_range filtration_simplex_range() { + if (globalDbg) { + std::cerr << "Filtration_simplex_range filtration_simplex_range()\n"; + } + // Returns a range over the simplices of the complex in the order of the filtration + return Filtration_simplex_range(this); + } + //*********************************************// + + + + //*********************************************// + // Elements which are in Gudhi now, but I (and in all the cases I asked also Marc) do not understand why they are + // there. + // TODO(PD) the file IndexingTag.h in the Gudhi library contains an empty structure, so + // I understand that this is something that was planned (for simplicial maps?) + // but was never finished. The only idea I have here is to use the same empty structure from + // IndexingTag.h file, but only if the compiler needs it. If the compiler + // do not need it, then I would rather not add here elements which I do not understand. + // typedef Indexing_tag + + /** + * Function needed for compatibility with Gudhi. Not useful for other purposes. + **/ + std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) { + std::vector< size_t > bdry = this->get_boundary_of_a_cell(sh); + if (globalDbg) { + std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n"; + std::cerr << "bdry.size() : " << bdry.size() << std::endl; + } + // this method returns two first elements from the boundary of sh. + if (bdry.size() < 2) + throw("Error in endpoints in Bitmap_cubical_complex class. The cell have less than two elements in the " + "boundary."); + return std::make_pair(bdry[0], bdry[1]); + } + + + /** + * Class needed for compatibility with Gudhi. Not useful for other purposes. + **/ + class Skeleton_simplex_range; + + class Skeleton_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > { + // Iterator over all simplices of the complex in the order of the indexing scheme. + // 'value_type' must be 'Simplex_handle'. + public: + Skeleton_simplex_iterator(Bitmap_cubical_complex* b, size_t d) : b(b), dimension(d) { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , size_t d )\n"; + } + // find the position of the first simplex of a dimension d + this->position = 0; + while ( + (this->position != b->data.size()) && + (this->b->get_dimension_of_a_cell(this->position) != this->dimension) + ) { + ++this->position; + } + } + + Skeleton_simplex_iterator() : b(NULL), position(0), dimension(0) { } + + Skeleton_simplex_iterator operator++() { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator operator++()\n"; + } + // increment the position as long as you did not get to the next element of the dimension dimension. + ++this->position; + while ( + (this->position != this->b->data.size()) && + (this->b->get_dimension_of_a_cell(this->position) != this->dimension) + ) { + ++this->position; + } + return (*this); + } + + Skeleton_simplex_iterator operator++(int) { + Skeleton_simplex_iterator result = *this; + ++(*this); + return result; + } + + Skeleton_simplex_iterator& operator=(const Skeleton_simplex_iterator& rhs) { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator operator =\n"; + } + this->b = rhs.b; + this->position = rhs.position; + this->dimension = rhs.dimension; + return (*this); + } + + bool operator==(const Skeleton_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator ==\n"; + } + return ( this->position == rhs.position); + } + + bool operator!=(const Skeleton_simplex_iterator& rhs)const { + if (globalDbg) { + std::cerr << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n"; + } + return !(*this == rhs); + } + + Simplex_handle operator*() { + if (globalDbg) { + std::cerr << "Simplex_handle operator*() \n"; + } + return this->position; + } + + friend class Skeleton_simplex_range; + private: + Bitmap_cubical_complex<T>* b; + size_t position; + unsigned dimension; + }; + + /** + * Class needed for compatibility with Gudhi. Not useful for other purposes. + **/ + class Skeleton_simplex_range { + // Range over the simplices of the complex in the order of the filtration. + // .begin() and .end() return type Filtration_simplex_iterator. + public: + typedef Skeleton_simplex_iterator const_iterator; + typedef Skeleton_simplex_iterator iterator; + + Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, unsigned dimension) : b(b), dimension(dimension) { } + + Skeleton_simplex_iterator begin() { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator begin()\n"; + } + return Skeleton_simplex_iterator(this->b, this->dimension); + } + + Skeleton_simplex_iterator end() { + if (globalDbg) { + std::cerr << "Skeleton_simplex_iterator end()\n"; + } + Skeleton_simplex_iterator it(this->b, this->dimension); + it.position = this->b->data.size(); + return it; + } + + private: + Bitmap_cubical_complex<T>* b; + unsigned dimension; + }; + + /** + * Function needed for compatibility with Gudhi. Not useful for other purposes. + **/ + Skeleton_simplex_range skeleton_simplex_range(unsigned dimension) { + if (globalDbg) { + std::cerr << "Skeleton_simplex_range skeleton_simplex_range( unsigned dimension )\n"; + } + return Skeleton_simplex_range(this, dimension); + } + + friend class is_before_in_filtration<T>; + + protected: + std::vector< size_t > key_associated_to_simplex; + std::vector< size_t > simplex_associated_to_key; +}; // Bitmap_cubical_complex + +template <typename T> +void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key() { + if (globalDbg) { + std::cerr << "void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration() \n"; + } + this->simplex_associated_to_key = std::vector<size_t>(this->data.size()); + std::iota(std::begin(simplex_associated_to_key), std::end(simplex_associated_to_key), 0); +#ifdef GUDHI_USE_TBB + tbb::parallel_sort(simplex_associated_to_key, is_before_in_filtration<T>(this)); +#else + std::sort(simplex_associated_to_key.begin(), simplex_associated_to_key.end(), is_before_in_filtration<T>(this)); +#endif + + // we still need to deal here with a key_associated_to_simplex: + for ( size_t i = 0 ; i != simplex_associated_to_key.size() ; ++i ) { + this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i; + } +} + +template <typename T> +class is_before_in_filtration { + public: + explicit is_before_in_filtration(Bitmap_cubical_complex<T> * CC) + : CC_(CC) { } + + bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1, + const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const { + // Not using st_->filtration(sh1) because it uselessly tests for null_simplex. + typename T::filtration_type fil1 = CC_->data[sh1]; + typename T::filtration_type fil2 = CC_->data[sh2]; + if (fil1 != fil2) { + return fil1 < fil2; + } + // in this case they are on the same filtration level, so the dimension decide. + size_t dim1 = CC_->get_dimension_of_a_cell(sh1); + size_t dim2 = CC_->get_dimension_of_a_cell(sh2); + if (dim1 != dim2) { + return dim1 < dim2; + } + // in this case both filtration and dimensions of the considered cubes are the same. To have stable sort, we simply + // compare their positions in the bitmap: + return sh1 < sh2; + } + + protected: + Bitmap_cubical_complex<T>* CC_; +}; + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_H_ diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h new file mode 100644 index 00000000..266ce051 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h @@ -0,0 +1,143 @@ +/* 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_COUNTER_H_ +#define BITMAP_CUBICAL_COMPLEX_COUNTER_H_ + +#include <iostream> +#include <vector> + +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(const std::vector<unsigned>& endd) : begin(endd.size(), 0), end(endd), current(endd.size(), 0) { } + + /** + * 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(const std::vector< unsigned >& beginn, const std::vector< unsigned >& endd) : begin(beginn), end(endd), current(endd.size(), 0) { + if (beginn.size() != endd.size()) + throw "In constructor of a counter, begin and end vectors do not have the same size. Program terminate"; + } + + /** + * 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(const 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) { + // std::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; +}; + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_COUNTER_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..7294da98 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -0,0 +1,819 @@ +/* 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/Bitmap_cubical_complex/counter.h> + +#include <iostream> +#include <vector> +#include <string> +#include <fstream> +#include <algorithm> +#include <iterator> +#include <limits> +#include <utility> // for pair<> + +namespace Gudhi { + +namespace Cubical_complex { + +/** + * @class Bitmap_cubical_complex_base + * @brief Cubical complex represented as a bitmap, class with basic implementation. + * @ingroup cubical_complex + */ + +/** + * This is a class implementing a basic bitmap data structure to store cubical complexes. + * It implements only the most basic subroutines. + * 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: + typedef T filtration_type; + + /** + *Default constructor + **/ + Bitmap_cubical_complex_base() : + total_number_of_cells(0) { } + /** + * There are a few constructors of a Bitmap_cubical_complex_base class. + * First one, that takes vector<unsigned>, creates an empty bitmap of a dimension equal + * the number of elements in the + * input vector and size in the i-th dimension equal the number in the position i-of the input vector. + */ + Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes); + /** + * The second constructor takes as a input a Perseus style file. For more details, + * please consult the documentations of + * Perseus software as well as examples attached to this + * implementation. + **/ + Bitmap_cubical_complex_base(const char* perseus_style_file); + /** + * The last constructor of a Bitmap_cubical_complex_base class accepts vector of dimensions (as the first one) + * together with vector of filtration values of top dimensional cells. + **/ + Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells); + + /** + * Destructor of the Bitmap_cubical_complex_base class. + **/ + virtual ~Bitmap_cubical_complex_base() { } + + /** + * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell + * and get_cell_data are the basic + * functions that compute boundary / coboundary / dimension and the filtration + * value form a position of a cell in the structure of a bitmap. The input parameter of all of those function is a + * non-negative integer, indicating a position of a cube in the data structure. + * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to + * the positions of (co)boundary element of the input cell. + */ + virtual inline std::vector< size_t > get_boundary_of_a_cell(size_t cell)const; + /** + * The functions get_coboundary_of_a_cell, get_coboundary_of_a_cell, + * get_dimension_of_a_cell and get_cell_data are the basic + * functions that compute boundary / coboundary / dimension and the filtration + * value form a position of a cell in the structure of a bitmap. + * The input parameter of all of those function is a non-negative integer, + * indicating a position of a cube in the data structure. + * In the case of functions that compute (co)boundary, the output is a vector if + * non-negative integers pointing to the + * positions of (co)boundary element of the input cell. + **/ + virtual inline std::vector< size_t > get_coboundary_of_a_cell(size_t cell)const; + /** + * In the case of get_dimension_of_a_cell function, the output is a non-negative integer + * indicating the dimension of a cell. + **/ + inline unsigned get_dimension_of_a_cell(size_t cell)const; + /** + * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position. + * This allows reading and changing the value of filtration. Note that if the value of a filtration is changed, the + * code do not check if we have a filtration or not. i.e. it do not check if the value of a filtration of a cell is + * not smaller than the value of a filtration of its boundary and not greater than the value of its coboundary. + **/ + inline T& get_cell_data(size_t cell); + + + /** + * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells. + * Then, there are a few ways one can pick the filtration of lower dimensional + * cells. The most typical one is by so called lower star filtration. This function is always called by any + * constructor which takes the top dimensional cells. If you use such a constructor, + * then there is no need to call this function. Call it only if you are putting the filtration + * of the cells by your own (for instance by using Top_dimensional_cells_iterator). + **/ + void impose_lower_star_filtration(); // assume that top dimensional cells are already set. + + /** + * Returns dimension of a complex. + **/ + inline unsigned dimension()const { + return sizes.size(); + } + + /** + * Returns number of all cubes in the data structure. + **/ + inline unsigned size()const { + return this->data.size(); + } + + /** + * Writing to stream operator. By using it we get the values T of cells in order in which they are stored in the + * structure. This procedure is used for debugging purposes. + **/ + template <typename K> + friend std::ostream& operator<<(std::ostream & os, const Bitmap_cubical_complex_base<K>& b); + + /** + * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values + * equally distributed in the range of data. + * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute + * persistence gets worst. When dealing with this type of data, one may want to put different values on cells to + * some number of bins. The function put_data_to_bins( size_t number_of_bins ) is designed for that purpose. + * The parameter of the function is the number of bins (distinct values) we want to have in the cubical complex. + **/ + void put_data_to_bins(size_t number_of_bins); + + /** + * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values + * equally distributed in the range of data. + * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute + * persistence gets worst. When dealing with this type of data, one may want to put different values on cells to + * some number of bins. The function put_data_to_bins( T diameter_of_bin ) is designed for that purpose. + * The parameter of it is the diameter of each bin. Note that the bottleneck distance between the persistence + * diagram of the cubical complex before and after using such a function will be bounded by the parameter + * diameter_of_bin. + **/ + void put_data_to_bins(T diameter_of_bin); + + /** + * Functions to find min and max values of filtration. + **/ + std::pair< T, T > min_max_filtration(); + + // ITERATORS + + /** + * Iterator through all cells in the complex (in order they appear in the structure -- i.e. + * in lexicographical order). + **/ + class All_cells_iterator : std::iterator< std::input_iterator_tag, T > { + public: + All_cells_iterator() { + this->counter = 0; + } + + All_cells_iterator operator++() { + // first find first element of the counter that can be increased: + ++this->counter; + return *this; + } + + All_cells_iterator operator++(int) { + All_cells_iterator result = *this; + ++(*this); + return result; + } + + All_cells_iterator& operator=(const All_cells_iterator& rhs) { + this->counter = rhs.counter; + return *this; + } + + bool operator==(const All_cells_iterator& rhs)const { + if (this->counter != rhs.counter)return false; + return true; + } + + bool operator!=(const All_cells_iterator& rhs)const { + return !(*this == rhs); + } + + /* + * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as + * an argument of the following functions: + * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell + * boundary and coboundary and dimension + * and in function get_cell_data to get a filtration of a cell. + */ + size_t operator*() { + return this->counter; + } + friend class Bitmap_cubical_complex_base; + protected: + size_t counter; + }; + + /** + * Function returning a All_cells_iterator to the first cell of the bitmap. + **/ + All_cells_iterator all_cells_iterator_begin() { + All_cells_iterator a; + return a; + } + + /** + * Function returning a All_cells_iterator to the last cell of the bitmap. + **/ + All_cells_iterator all_cells_iterator_end() { + All_cells_iterator a; + a.counter = this->data.size(); + return a; + } + + /** + * All_cells_range class provides ranges for All_cells_iterator + **/ + class All_cells_range { + public: + All_cells_range(Bitmap_cubical_complex_base* b) : b(b) { } + + All_cells_iterator begin() { + return b->all_cells_iterator_begin(); + } + + All_cells_iterator end() { + return b->all_cells_iterator_end(); + } + private: + Bitmap_cubical_complex_base<T>* b; + }; + + All_cells_range all_cells_range() { + return All_cells_range(this); + } + + + /** + * Boundary_range class provides ranges for boundary iterators. + **/ + typedef typename std::vector< size_t >::const_iterator Boundary_iterator; + typedef typename std::vector< size_t > Boundary_range; + + /** + * boundary_simplex_range creates an object of a Boundary_simplex_range class + * that provides ranges for the Boundary_simplex_iterator. + **/ + Boundary_range boundary_range(size_t sh) { + return this->get_boundary_of_a_cell(sh); + } + + /** + * Coboundary_range class provides ranges for boundary iterators. + **/ + typedef typename std::vector< size_t >::const_iterator Coboundary_iterator; + typedef typename std::vector< size_t > Coboundary_range; + + /** + * boundary_simplex_range creates an object of a Boundary_simplex_range class + * that provides ranges for the Boundary_simplex_iterator. + **/ + Coboundary_range coboundary_range(size_t sh) { + return this->get_coboundary_of_a_cell(sh); + } + + /** + * Iterator through top dimensional cells of the complex. The cells appear in order they are stored + * in the structure (i.e. in lexicographical order) + **/ + class Top_dimensional_cells_iterator : std::iterator< std::input_iterator_tag, T > { + public: + Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) { + this->counter = std::vector<size_t>(b.dimension()); + // std::fill( this->counter.begin() , this->counter.end() , 0 ); + } + + Top_dimensional_cells_iterator operator++() { + // first find first element of the counter that can be increased: + size_t dim = 0; + while ((dim != this->b.dimension()) && (this->counter[dim] == this->b.sizes[dim] - 1))++dim; + + if (dim != this->b.dimension()) { + ++this->counter[dim]; + for (size_t i = 0; i != dim; ++i) { + this->counter[i] = 0; + } + } else { + ++this->counter[0]; + } + return *this; + } + + Top_dimensional_cells_iterator operator++(int) { + Top_dimensional_cells_iterator result = *this; + ++(*this); + return result; + } + + Top_dimensional_cells_iterator& operator=(const Top_dimensional_cells_iterator& rhs) { + this->counter = rhs.counter; + this->b = rhs.b; + return *this; + } + + bool operator==(const Top_dimensional_cells_iterator& rhs)const { + if (&this->b != &rhs.b)return false; + if (this->counter.size() != rhs.counter.size())return false; + for (size_t i = 0; i != this->counter.size(); ++i) { + if (this->counter[i] != rhs.counter[i])return false; + } + return true; + } + + bool operator!=(const Top_dimensional_cells_iterator& rhs)const { + return !(*this == rhs); + } + + /* + * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as + * an argument of the following functions: + * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell + * boundary and coboundary and dimension + * and in function get_cell_data to get a filtration of a cell. + */ + size_t operator*() { + return this->compute_index_in_bitmap(); + } + + size_t compute_index_in_bitmap()const { + size_t index = 0; + for (size_t i = 0; i != this->counter.size(); ++i) { + index += (2 * this->counter[i] + 1) * this->b.multipliers[i]; + } + return index; + } + + void print_counter()const { + for (size_t i = 0; i != this->counter.size(); ++i) { + std::cout << this->counter[i] << " "; + } + } + friend class Bitmap_cubical_complex_base; + protected: + std::vector< size_t > counter; + Bitmap_cubical_complex_base& b; + }; + + /** + * Function returning a Top_dimensional_cells_iterator to the first top dimensional cell of the bitmap. + **/ + Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin() { + Top_dimensional_cells_iterator a(*this); + return a; + } + + /** + * Function returning a Top_dimensional_cells_iterator to the last top dimensional cell of the bitmap. + **/ + Top_dimensional_cells_iterator top_dimensional_cells_iterator_end() { + Top_dimensional_cells_iterator a(*this); + for (size_t i = 0; i != this->dimension(); ++i) { + a.counter[i] = this->sizes[i] - 1; + } + a.counter[0]++; + return a; + } + + /** + * Top_dimensional_cells_iterator_range class provides ranges for Top_dimensional_cells_iterator_range + **/ + class Top_dimensional_cells_range { + public: + Top_dimensional_cells_range(Bitmap_cubical_complex_base* b) : b(b) { } + + Top_dimensional_cells_iterator begin() { + return b->top_dimensional_cells_iterator_begin(); + } + + Top_dimensional_cells_iterator end() { + return b->top_dimensional_cells_iterator_end(); + } + private: + Bitmap_cubical_complex_base<T>* b; + }; + + Top_dimensional_cells_range top_dimensional_cells_range() { + return Top_dimensional_cells_range(this); + } + + + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + + inline size_t number_cells()const { + return this->total_number_of_cells; + } + + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + //****************************************************************************************************************// + + protected: + std::vector<unsigned> sizes; + std::vector<unsigned> multipliers; + std::vector<T> data; + size_t total_number_of_cells; + + void set_up_containers(const std::vector<unsigned>& sizes) { + unsigned multiplier = 1; + for (size_t i = 0; i != sizes.size(); ++i) { + this->sizes.push_back(sizes[i]); + this->multipliers.push_back(multiplier); + multiplier *= 2 * sizes[i] + 1; + } + this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max()); + this->total_number_of_cells = multiplier; + } + + size_t compute_position_in_bitmap(const std::vector< unsigned >& counter) { + size_t position = 0; + for (size_t i = 0; i != this->multipliers.size(); ++i) { + position += this->multipliers[i] * counter[i]; + } + return position; + } + + std::vector<unsigned> compute_counter_for_given_cell(size_t cell)const { + std::vector<unsigned> counter; + counter.reserve(this->sizes.size()); + for (size_t dim = this->sizes.size(); dim != 0; --dim) { + counter.push_back(cell / this->multipliers[dim - 1]); + cell = cell % this->multipliers[dim - 1]; + } + std::reverse(counter.begin(), counter.end()); + return counter; + } + void read_perseus_style_file(const char* perseus_style_file); + void setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions, + const std::vector<T>& top_dimensional_cells); + Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector<bool> directions); + Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, std::vector<bool> directions); + Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& top_dimensional_cells, + std::vector<bool> directions); +}; + +template <typename T> +void Bitmap_cubical_complex_base<T>::put_data_to_bins(size_t number_of_bins) { + bool bdg = false; + + std::pair< T, T > min_max = this->min_max_filtration(); + T dx = (min_max.second - min_max.first) / (T) number_of_bins; + + // now put the data into the appropriate bins: + for (size_t i = 0; i != this->data.size(); ++i) { + if (bdg) { + std::cerr << "Before binning : " << this->data[i] << std::endl; + } + this->data[i] = min_max.first + dx * (this->data[i] - min_max.first) / number_of_bins; + if (bdg) { + std::cerr << "After binning : " << this->data[i] << std::endl; + getchar(); + } + } +} + +template <typename T> +void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) { + bool bdg = false; + std::pair< T, T > min_max = this->min_max_filtration(); + + size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin; + // now put the data into the appropriate bins: + for (size_t i = 0; i != this->data.size(); ++i) { + if (bdg) { + std::cerr << "Before binning : " << this->data[i] << std::endl; + } + this->data[i] = min_max.first + diameter_of_bin * (this->data[i] - min_max.first) / number_of_bins; + if (bdg) { + std::cerr << "After binning : " << this->data[i] << std::endl; + getchar(); + } + } +} + +template <typename T> +std::pair< T, T > Bitmap_cubical_complex_base<T>::min_max_filtration() { + std::pair< T, T > min_max(std::numeric_limits<T>::max(), std::numeric_limits<T>::min()); + for (size_t i = 0; i != this->data.size(); ++i) { + if (this->data[i] < min_max.first)min_max.first = this->data[i]; + if (this->data[i] > min_max.second)min_max.second = this->data[i]; + } + return min_max; +} + +template <typename K> +std::ostream& operator<<(std::ostream & out, const Bitmap_cubical_complex_base<K>& b) { + for (typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator + it = b.all_cells_const_begin(); it != b.all_cells_const_end(); ++it) { + out << *it << " "; + } + return out; +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base +(const std::vector<unsigned>& sizes) { + this->set_up_containers(sizes); +} + +template <typename T> +void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions, + const std::vector<T>& top_dimensional_cells) { + this->set_up_containers(sizes_in_following_directions); + + size_t number_of_top_dimensional_elements = 1; + for (size_t i = 0; i != sizes_in_following_directions.size(); ++i) { + number_of_top_dimensional_elements *= sizes_in_following_directions[i]; + } + if (number_of_top_dimensional_elements != top_dimensional_cells.size()) { + std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<size_t> sizes_in_following_directions" + << ", std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " + << "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector." + << std::endl; + throw("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions," + "std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " + "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector."); + } + + Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + size_t index = 0; + for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) { + this->get_cell_data(*it) = top_dimensional_cells[index]; + ++index; + } + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base +(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) { + this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells); +} + +template <typename T> +void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus_style_file) { + bool dbg = false; + std::ifstream inFiltration; + inFiltration.open(perseus_style_file); + unsigned dimensionOfData; + inFiltration >> dimensionOfData; + + if (dbg) { + std::cerr << "dimensionOfData : " << dimensionOfData << std::endl; + getchar(); + } + + std::vector<unsigned> sizes; + sizes.reserve(dimensionOfData); + for (size_t i = 0; i != dimensionOfData; ++i) { + unsigned size_in_this_dimension; + inFiltration >> size_in_this_dimension; + sizes.push_back(size_in_this_dimension); + if (dbg) { + std::cerr << "size_in_this_dimension : " << size_in_this_dimension << std::endl; + } + } + this->set_up_containers(sizes); + + Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + it = this->top_dimensional_cells_iterator_begin(); + + while (!inFiltration.eof()) { + T filtrationLevel; + inFiltration >> filtrationLevel; + if (dbg) { + std::cerr << "Cell of an index : " + << it.compute_index_in_bitmap() + << " and dimension: " + << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) + << " get the value : " << filtrationLevel << std::endl; + } + this->get_cell_data(*it) = filtrationLevel; + ++it; + } + inFiltration.close(); + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file, + std::vector<bool> directions) { + // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary + // conditions. + // It ignores the last parameter of the function. + this->read_perseus_style_file(perseus_style_file); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, + std::vector<bool> directions) { + // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary + // conditions. + // It ignores the last parameter of the function. + this->set_up_containers(sizes); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& top_dimensional_cells, + std::vector<bool> directions) { + // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary + // conditions. + // It ignores the last parameter of the function. + this->setup_bitmap_based_on_top_dimensional_cells_list(dimensions, top_dimensional_cells); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file) { + this->read_perseus_style_file(perseus_style_file); +} + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell)const { + std::vector< size_t > boundary_elements; + + // Speed traded of for memory. Check if it is better in practice. + boundary_elements.reserve(this->dimension()*2); + + size_t cell1 = cell; + for (size_t i = this->multipliers.size(); i != 0; --i) { + unsigned position = cell1 / this->multipliers[i - 1]; + if (position % 2 == 1) { + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + } + cell1 = cell1 % this->multipliers[i - 1]; + } + return boundary_elements; +} + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell)const { + std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell); + std::vector< size_t > coboundary_elements; + size_t cell1 = cell; + for (size_t i = this->multipliers.size(); i != 0; --i) { + unsigned position = cell1 / this->multipliers[i - 1]; + if (position % 2 == 0) { + if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) { + coboundary_elements.push_back(cell - this->multipliers[i - 1]); + } + if ( + (cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) { + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + } + } + cell1 = cell1 % this->multipliers[i - 1]; + } + return coboundary_elements; +} + +template <typename T> +unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell)const { + bool dbg = false; + if (dbg) std::cerr << "\n\n\n Computing position o a cell of an index : " << cell << std::endl; + unsigned dimension = 0; + for (size_t i = this->multipliers.size(); i != 0; --i) { + unsigned position = cell / this->multipliers[i - 1]; + + if (dbg) { + std::cerr << "i-1 :" << i - 1 << std::endl; + std::cerr << "cell : " << cell << std::endl; + std::cerr << "position : " << position << std::endl; + std::cerr << "multipliers[" << i - 1 << "] = " << this->multipliers[i - 1] << std::endl; + getchar(); + } + + if (position % 2 == 1) { + if (dbg) std::cerr << "Nonzero length in this direction \n"; + dimension++; + } + cell = cell % this->multipliers[i - 1]; + } + return dimension; +} + +template <typename T> +inline T& Bitmap_cubical_complex_base<T>::get_cell_data(size_t cell) { + return this->data[cell]; +} + +template <typename T> +void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { + bool dbg = false; + + // this vector will be used to check which elements have already been taken care of in imposing lower star filtration + std::vector<bool> is_this_cell_considered(this->data.size(), false); + + size_t size_to_reserve = 1; + for (size_t i = 0; i != this->multipliers.size(); ++i) { + size_to_reserve *= (size_t) ((this->multipliers[i] - 1) / 2); + } + + std::vector<size_t> indices_to_consider; + indices_to_consider.reserve(size_to_reserve); + // we assume here that we already have a filtration on the top dimensional cells and + // we have to extend it to lower ones. + typename Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) { + indices_to_consider.push_back(it.compute_index_in_bitmap()); + } + + while (indices_to_consider.size()) { + if (dbg) { + std::cerr << "indices_to_consider in this iteration \n"; + for (size_t i = 0; i != indices_to_consider.size(); ++i) { + std::cout << indices_to_consider[i] << " "; + } + getchar(); + } + std::vector<size_t> new_indices_to_consider; + for (size_t i = 0; i != indices_to_consider.size(); ++i) { + std::vector<size_t> bd = this->get_boundary_of_a_cell(indices_to_consider[i]); + for (size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) { + if (dbg) { + std::cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[ bd[boundaryIt] ] + << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[ indices_to_consider[i] ] + << std::endl; + getchar(); + } + if (this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ]) { + this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ]; + if (dbg) { + std::cerr << "Setting the value of a cell : " << bd[boundaryIt] << " to : " + << this->data[ indices_to_consider[i] ] << std::endl; + getchar(); + } + } + 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; + } +} + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_BASE_H_ diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h new file mode 100644 index 00000000..a446c0e8 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h @@ -0,0 +1,309 @@ +/* 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_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_ +#define BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_ + +#include <gudhi/Bitmap_cubical_complex_base.h> + +#include <cmath> +#include <limits> // for numeric_limits<> +#include <vector> + +namespace Gudhi { + +namespace Cubical_complex { + +// in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic +// boundary conditions). But, we set up the iterators and the procedures to compute boundary and coboundary in the way +// that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all +// the cells on the right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base + +/** + * @class Bitmap_cubical_complex_periodic_boundary_conditions_base + * @brief Cubical complex with periodic boundary conditions represented as a bitmap. + * @ingroup cubical_complex + */ +/** + * This is a class implementing a bitmap data structure with periodic boundary conditions. Most of the functions are + * identical to the functions from Bitmap_cubical_complex_base. + * The ones that needed to be updated are the constructors and get_boundary_of_a_cell and get_coboundary_of_a_cell. + */ +template <typename T> +class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T> { + public: + // constructors that take an extra parameter: + + /** + * Default constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base() { } + /** + * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following + * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of booleans. If + * at i-th position of this vector there is true value, that means that periodic boundary conditions are to be + * imposed in this direction. In case of false, the periodic boundary conditions will not be imposed in the direction + * i. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); + /** + * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the name of Perseus + * style file as an input. Please consult the documentation about the specification of the file. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseusStyleFile); + /** + * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following + * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of top dimensional + * cells (ordered lexicographically) and (3) vector of booleans. If at i-th position of this vector there is true + * value, that means that periodic boundary conditions are to be imposed in this direction. In case of false, the + * periodic boundary conditions will not be imposed in the direction i. + */ + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector< bool >& directions_in_which_periodic_b_cond_are_to_be_imposed); + + /** + * Destructor of the Bitmap_cubical_complex_periodic_boundary_conditions_base class. + **/ + virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base() {} + + // overwritten methods co compute boundary and coboundary + /** + * A version of a function that return boundary of a given cell for an object of + * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + */ + virtual std::vector< size_t > get_boundary_of_a_cell(size_t cell) const; + + /** + * A version of a function that return coboundary of a given cell for an object of + * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + */ + virtual std::vector< size_t > get_coboundary_of_a_cell(size_t cell) const; + + protected: + std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed; + + void set_up_containers(const std::vector<unsigned>& sizes) { + unsigned multiplier = 1; + for (size_t i = 0; i != sizes.size(); ++i) { + this->sizes.push_back(sizes[i]); + this->multipliers.push_back(multiplier); + + if (directions_in_which_periodic_b_cond_are_to_be_imposed[i]) { + multiplier *= 2 * sizes[i]; + } else { + multiplier *= 2 * sizes[i] + 1; + } + } + // std::reverse( this->sizes.begin() , this->sizes.end() ); + this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max()); + this->total_number_of_cells = multiplier; + } + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes); + Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells); + void construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); +}; + +template <typename T> +void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { + this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed; + this->set_up_containers(dimensions); + + size_t i = 0; + for (auto it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) { + this->get_cell_data(*it) = topDimensionalCells[i]; + ++i; + } + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { + this->directions_in_which_periodic_b_cond_are_to_be_imposed(directions_in_which_periodic_b_cond_are_to_be_imposed); + this->set_up_containers(sizes); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseus_style_file) { + // for Perseus style files: + bool dbg = false; + + std::ifstream inFiltration; + inFiltration.open(perseus_style_file); + unsigned dimensionOfData; + inFiltration >> dimensionOfData; + + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData, false); + + std::vector<unsigned> sizes; + sizes.reserve(dimensionOfData); + for (size_t i = 0; i != dimensionOfData; ++i) { + int size_in_this_dimension; + inFiltration >> size_in_this_dimension; + if (size_in_this_dimension < 0) { + this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] = true; + } + sizes.push_back(abs(size_in_this_dimension)); + } + this->set_up_containers(sizes); + + typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator it(*this); + it = this->top_dimensional_cells_iterator_begin(); + + while (!inFiltration.eof()) { + double filtrationLevel; + inFiltration >> filtrationLevel; + if (inFiltration.eof())break; + + if (dbg) { + std::cerr << "Cell of an index : " + << it.compute_index_in_bitmap() + << " and dimension: " + << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) + << " get the value : " << filtrationLevel << std::endl; + } + this->get_cell_data(*it) = filtrationLevel; + ++it; + } + inFiltration.close(); + this->impose_lower_star_filtration(); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes) { + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(), false); + this->set_up_containers(sizes); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells) { + std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(), false); + this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells, + directions_in_which_periodic_b_cond_are_to_be_imposed); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>:: +Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, + const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { + this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells, + directions_in_which_periodic_b_cond_are_to_be_imposed); +} + +// ***********************Methods************************ // + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell(size_t cell) const { + bool dbg = false; + if (dbg) { + std::cerr << "Computations of boundary of a cell : " << cell << std::endl; + } + + std::vector< size_t > boundary_elements; + size_t cell1 = cell; + for (size_t i = this->multipliers.size(); i != 0; --i) { + unsigned position = cell1 / this->multipliers[i - 1]; + // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction. + + if (position % 2 == 1) { + // if there are no periodic boundary conditions in this direction, we do not have to do anything. + if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) { + // std::cerr << "A\n"; + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + if (dbg) { + std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " "; + } + } else { + // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end. + if (position != 2 * this->sizes[ i - 1 ] - 1) { + // std::cerr << "B\n"; + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + if (dbg) { + std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " "; + } + } else { + // std::cerr << "C\n"; + boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); + boundary_elements.push_back(cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ]); + if (dbg) { + std::cerr << cell - this->multipliers[ i - 1 ] << " " << + cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ] << " "; + } + } + } + } + cell1 = cell1 % this->multipliers[i - 1]; + } + return boundary_elements; +} + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell(size_t cell) const { + std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell); + std::vector< size_t > coboundary_elements; + size_t cell1 = cell; + for (size_t i = this->multipliers.size(); i != 0; --i) { + unsigned position = cell1 / this->multipliers[i - 1]; + // if the cell has zero length in this direction, then it will have cbd in this direction. + if (position % 2 == 0) { + if (!this->directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) { + // no periodic boundary conditions in this direction + if ((counter[i - 1] != 0) && (cell > this->multipliers[i - 1])) { + coboundary_elements.push_back(cell - this->multipliers[i - 1]); + } + if ((counter[i - 1] != 2 * this->sizes[i - 1]) && (cell + this->multipliers[i - 1] < this->data.size())) { + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + } + } else { + // we want to have periodic boundary conditions in this direction + if (counter[i - 1] != 0) { + coboundary_elements.push_back(cell - this->multipliers[i - 1]); + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + } else { + // in this case counter[i-1] == 0. + coboundary_elements.push_back(cell + this->multipliers[i - 1]); + coboundary_elements.push_back(cell + (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[i - 1]); + } + } + } + + cell1 = cell1 % this->multipliers[i - 1]; + } + return coboundary_elements; +} + +} // namespace Cubical_complex + +} // namespace Gudhi + +#endif // BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_ 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..a9162cee --- /dev/null +++ b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp @@ -0,0 +1,1378 @@ +/* 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/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "cubical_complex" +#include <boost/test/unit_test.hpp> + +#include <gudhi/reader_utils.h> +#include <gudhi/Bitmap_cubical_complex.h> +#include <gudhi/Persistent_cohomology.h> + +// standard stuff +#include <iostream> +#include <sstream> +#include <vector> + + +typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base; +typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex; + +typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double> +Bitmap_cubical_complex_periodic_boundary_conditions_base; +typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_periodic_boundary_conditions_base> +Bitmap_cubical_complex_periodic_boundary_conditions; + +BOOST_AUTO_TEST_CASE(check_dimension) { + std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + + std::vector<unsigned> dimensions({3, 3}); + + Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells); + BOOST_CHECK(increasing.dimension() == 2); +} + +BOOST_AUTO_TEST_CASE(topDimensionalCellsIterator_test) { + std::vector< double > expectedFiltrationValues1({0, 0, 0, 0, 100, 0, 0, 0, 0}); + + std::vector< double > expectedFiltrationValues2({1, 2, 3, 4, 5, 6, 7, 8, 9}); + + std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + + std::vector< double > oneDimensionalCycle({0, 0, 0, 0, 100, 0, 0, 0, 0}); + + std::vector<unsigned> dimensions({3, 3}); + + Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells); + Bitmap_cubical_complex hole(dimensions, oneDimensionalCycle); + + + int i = 0; + for (Bitmap_cubical_complex::Top_dimensional_cells_iterator + it = increasing.top_dimensional_cells_iterator_begin(); it != increasing.top_dimensional_cells_iterator_end(); ++it) { + BOOST_CHECK(increasing.get_cell_data(*it) == expectedFiltrationValues2[i]); + ++i; + } + i = 0; + for (Bitmap_cubical_complex::Top_dimensional_cells_iterator + it = hole.top_dimensional_cells_iterator_begin(); it != hole.top_dimensional_cells_iterator_end(); ++it) { + BOOST_CHECK(hole.get_cell_data(*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({1, 2, 3, 4, 5, 6, 7, 8, 9}); + + std::vector<unsigned> dimensions({3, 3}); + + Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells); + for (size_t i = 0; i != increasing.size(); ++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({1, 2, 3, 4, 5, 6, 7, 8, 9}); + + std::vector<unsigned> dimensions({3, 3}); + + Bitmap_cubical_complex 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(); ++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({1, 2, 3, 4, 5, 6, 7, 8, 9}); + + std::vector<unsigned> dimensions({3, 3}); + + Bitmap_cubical_complex 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(); ++i) { + BOOST_CHECK(increasing.get_dimension_of_a_cell(i) == dim[i]); + } +} + +BOOST_AUTO_TEST_CASE(Filtration_simplex_iterator_test) { + std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + + std::vector<unsigned> dimensions({3, 3}); + + Bitmap_cubical_complex 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::Filtration_simplex_range range = increasing.filtration_simplex_range(); + size_t position = 0; + for (Bitmap_cubical_complex::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; + } +} + +BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) { + std::vector< double > filtration({0, 0, 0, 0}); + + std::vector<unsigned> dimensions({2, 2}); + + std::vector<bool> periodic_directions({true, true}); + + Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtration, periodic_directions); + BOOST_CHECK(cmplx.dimension() == 2); + + + 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(0); + std::vector<double> boundary4; + boundary4.push_back(0); + boundary4.push_back(8); + std::vector<double> boundary5; + boundary5.push_back(1); + boundary5.push_back(9); + boundary5.push_back(4); + boundary5.push_back(6); + std::vector<double> boundary6; + boundary6.push_back(2); + boundary6.push_back(10); + std::vector<double> boundary7; + boundary7.push_back(3); + boundary7.push_back(11); + boundary7.push_back(6); + boundary7.push_back(4); + std::vector<double> boundary8; + std::vector<double> boundary9; + boundary9.push_back(8); + boundary9.push_back(10); + std::vector<double> boundary10; + std::vector<double> boundary11; + boundary11.push_back(10); + boundary11.push_back(8); + std::vector<double> boundary12; + boundary12.push_back(8); + boundary12.push_back(0); + std::vector<double> boundary13; + boundary13.push_back(9); + boundary13.push_back(1); + boundary13.push_back(12); + boundary13.push_back(14); + std::vector<double> boundary14; + boundary14.push_back(10); + boundary14.push_back(2); + std::vector<double> boundary15; + boundary15.push_back(11); + boundary15.push_back(3); + boundary15.push_back(14); + boundary15.push_back(12); + + 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); + + for (size_t i = 0; i != cmplx.size(); ++i) { + std::vector< size_t > bd = cmplx.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(coboudary_operator_2d_bitmap_with_periodic_bcond) { + std::vector< double > filtration({0, 0, 0, 0}); + + std::vector<unsigned> dimensions({2, 2}); + + std::vector<bool> periodic_directions({true, true}); + + Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtration, periodic_directions); + BOOST_CHECK(cmplx.dimension() == 2); + + + std::vector<double> coboundary0; + coboundary0.push_back(4); + coboundary0.push_back(12); + coboundary0.push_back(1); + coboundary0.push_back(3); + std::vector<double> coboundary1; + coboundary1.push_back(5); + coboundary1.push_back(13); + std::vector<double> coboundary2; + coboundary2.push_back(6); + coboundary2.push_back(14); + coboundary2.push_back(1); + coboundary2.push_back(3); + std::vector<double> coboundary3; + coboundary3.push_back(7); + coboundary3.push_back(15); + std::vector<double> coboundary4; + coboundary4.push_back(5); + coboundary4.push_back(7); + std::vector<double> coboundary5; + std::vector<double> coboundary6; + coboundary6.push_back(5); + coboundary6.push_back(7); + std::vector<double> coboundary7; + std::vector<double> coboundary8; + coboundary8.push_back(4); + coboundary8.push_back(12); + coboundary8.push_back(9); + coboundary8.push_back(11); + std::vector<double> coboundary9; + coboundary9.push_back(5); + coboundary9.push_back(13); + std::vector<double> coboundary10; + coboundary10.push_back(6); + coboundary10.push_back(14); + coboundary10.push_back(9); + coboundary10.push_back(11); + std::vector<double> coboundary11; + coboundary11.push_back(7); + coboundary11.push_back(15); + std::vector<double> coboundary12; + coboundary12.push_back(13); + coboundary12.push_back(15); + std::vector<double> coboundary13; + std::vector<double> coboundary14; + coboundary14.push_back(13); + coboundary14.push_back(15); + std::vector<double> coboundary15; + + std::vector< std::vector<double> > coboundaries; + coboundaries.push_back(coboundary0); + coboundaries.push_back(coboundary1); + coboundaries.push_back(coboundary2); + coboundaries.push_back(coboundary3); + coboundaries.push_back(coboundary4); + coboundaries.push_back(coboundary5); + coboundaries.push_back(coboundary6); + coboundaries.push_back(coboundary7); + coboundaries.push_back(coboundary8); + coboundaries.push_back(coboundary9); + coboundaries.push_back(coboundary10); + coboundaries.push_back(coboundary11); + coboundaries.push_back(coboundary12); + coboundaries.push_back(coboundary13); + coboundaries.push_back(coboundary14); + coboundaries.push_back(coboundary15); + + for (size_t i = 0; i != cmplx.size(); ++i) { + std::vector< size_t > cbd = cmplx.get_coboundary_of_a_cell(i); + for (size_t j = 0; j != cbd.size(); ++j) { + BOOST_CHECK(coboundaries[i][j] == cbd[j]); + } + } +} + +BOOST_AUTO_TEST_CASE(bitmap_2d_with_periodic_bcond_filtration) { + std::vector< double > filtrationOrg({0, 1, 2, 3}); + + std::vector<unsigned> dimensions({2, 2}); + + std::vector<bool> periodic_directions({true, true}); + + Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtrationOrg, periodic_directions); + BOOST_CHECK(cmplx.dimension() == 2); + + + std::vector<double> filtration; + filtration.push_back(0); // 0 + filtration.push_back(0); // 1 + filtration.push_back(0); // 2 + filtration.push_back(1); // 3 + filtration.push_back(0); // 4 + filtration.push_back(0); // 5 + filtration.push_back(0); // 6 + filtration.push_back(1); // 7 + filtration.push_back(0); // 8 + filtration.push_back(0); // 9 + filtration.push_back(0); // 10 + filtration.push_back(1); // 11 + filtration.push_back(2); // 12 + filtration.push_back(2); // 13 + filtration.push_back(2); // 14 + filtration.push_back(3); // 15 + + + for (size_t i = 0; i != cmplx.size(); ++i) { + BOOST_CHECK(filtration[i] == cmplx.get_cell_data(i)); + } +} +BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check) +{ + std::vector< double > expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + std::vector<unsigned> expected_dimension; + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + + std::vector<size_t> expected_boundary; + expected_boundary.push_back(0); + expected_boundary.push_back(2); + expected_boundary.push_back(2); + expected_boundary.push_back(4); + expected_boundary.push_back(0); + expected_boundary.push_back(10); + expected_boundary.push_back(1); + expected_boundary.push_back(11); + expected_boundary.push_back(5); + expected_boundary.push_back(7); + expected_boundary.push_back(2); + expected_boundary.push_back(12); + expected_boundary.push_back(3); + expected_boundary.push_back(13); + expected_boundary.push_back(7); + expected_boundary.push_back(9); + expected_boundary.push_back(4); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(12); + expected_boundary.push_back(12); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(20); + expected_boundary.push_back(11); + expected_boundary.push_back(21); + expected_boundary.push_back(15); + expected_boundary.push_back(17); + expected_boundary.push_back(12); + expected_boundary.push_back(22); + expected_boundary.push_back(13); + expected_boundary.push_back(23); + expected_boundary.push_back(17); + expected_boundary.push_back(19); + expected_boundary.push_back(14); + expected_boundary.push_back(24); + expected_boundary.push_back(20); + expected_boundary.push_back(22); + expected_boundary.push_back(22); + expected_boundary.push_back(24); + + + std::vector<size_t> expected_coboundary; + expected_coboundary.push_back(5); + expected_coboundary.push_back(1); + expected_coboundary.push_back(6); + expected_coboundary.push_back(7); + expected_coboundary.push_back(1); + expected_coboundary.push_back(3); + expected_coboundary.push_back(8); + expected_coboundary.push_back(9); + expected_coboundary.push_back(3); + expected_coboundary.push_back(6); + expected_coboundary.push_back(6); + expected_coboundary.push_back(8); + expected_coboundary.push_back(8); + expected_coboundary.push_back(5); + expected_coboundary.push_back(15); + expected_coboundary.push_back(11); + expected_coboundary.push_back(6); + expected_coboundary.push_back(16); + expected_coboundary.push_back(7); + expected_coboundary.push_back(17); + expected_coboundary.push_back(11); + expected_coboundary.push_back(13); + expected_coboundary.push_back(8); + expected_coboundary.push_back(18); + expected_coboundary.push_back(9); + expected_coboundary.push_back(19); + expected_coboundary.push_back(13); + expected_coboundary.push_back(16); + expected_coboundary.push_back(16); + expected_coboundary.push_back(18); + expected_coboundary.push_back(18); + expected_coboundary.push_back(15); + expected_coboundary.push_back(21); + expected_coboundary.push_back(16); + expected_coboundary.push_back(17); + expected_coboundary.push_back(21); + expected_coboundary.push_back(23); + expected_coboundary.push_back(18); + expected_coboundary.push_back(19); + expected_coboundary.push_back(23); + + + + std::vector< unsigned > sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector< double > data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba( sizes , data ); + int i = 0; + int bd_it = 0; + int cbd_it = 0; + for ( Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it ) + { + BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) ); + BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) ); + + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd ) + { + BOOST_CHECK( expected_boundary[bd_it] == *bd ); + ++bd_it; + } + + Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); + for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd ) + { + BOOST_CHECK( expected_coboundary[cbd_it] == *cbd ); + ++cbd_it; + } + ++i; + } +} + + + + + + + + +BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check_2) +{ + std::vector< double > expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + std::vector<unsigned> expected_dimension; + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + + std::vector<size_t> expected_boundary; + expected_boundary.push_back(0); + expected_boundary.push_back(2); + expected_boundary.push_back(2); + expected_boundary.push_back(4); + expected_boundary.push_back(0); + expected_boundary.push_back(10); + expected_boundary.push_back(1); + expected_boundary.push_back(11); + expected_boundary.push_back(5); + expected_boundary.push_back(7); + expected_boundary.push_back(2); + expected_boundary.push_back(12); + expected_boundary.push_back(3); + expected_boundary.push_back(13); + expected_boundary.push_back(7); + expected_boundary.push_back(9); + expected_boundary.push_back(4); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(12); + expected_boundary.push_back(12); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(20); + expected_boundary.push_back(11); + expected_boundary.push_back(21); + expected_boundary.push_back(15); + expected_boundary.push_back(17); + expected_boundary.push_back(12); + expected_boundary.push_back(22); + expected_boundary.push_back(13); + expected_boundary.push_back(23); + expected_boundary.push_back(17); + expected_boundary.push_back(19); + expected_boundary.push_back(14); + expected_boundary.push_back(24); + expected_boundary.push_back(20); + expected_boundary.push_back(22); + expected_boundary.push_back(22); + expected_boundary.push_back(24); + + + std::vector<size_t> expected_coboundary; + expected_coboundary.push_back(5); + expected_coboundary.push_back(1); + expected_coboundary.push_back(6); + expected_coboundary.push_back(7); + expected_coboundary.push_back(1); + expected_coboundary.push_back(3); + expected_coboundary.push_back(8); + expected_coboundary.push_back(9); + expected_coboundary.push_back(3); + expected_coboundary.push_back(6); + expected_coboundary.push_back(6); + expected_coboundary.push_back(8); + expected_coboundary.push_back(8); + expected_coboundary.push_back(5); + expected_coboundary.push_back(15); + expected_coboundary.push_back(11); + expected_coboundary.push_back(6); + expected_coboundary.push_back(16); + expected_coboundary.push_back(7); + expected_coboundary.push_back(17); + expected_coboundary.push_back(11); + expected_coboundary.push_back(13); + expected_coboundary.push_back(8); + expected_coboundary.push_back(18); + expected_coboundary.push_back(9); + expected_coboundary.push_back(19); + expected_coboundary.push_back(13); + expected_coboundary.push_back(16); + expected_coboundary.push_back(16); + expected_coboundary.push_back(18); + expected_coboundary.push_back(18); + expected_coboundary.push_back(15); + expected_coboundary.push_back(21); + expected_coboundary.push_back(16); + expected_coboundary.push_back(17); + expected_coboundary.push_back(21); + expected_coboundary.push_back(23); + expected_coboundary.push_back(18); + expected_coboundary.push_back(19); + expected_coboundary.push_back(23); + + + + std::vector< unsigned > sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector< double > data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba( sizes , data ); + int i = 0; + int bd_it = 0; + int cbd_it = 0; + + Bitmap_cubical_complex_base::All_cells_range range(&ba); + for ( Bitmap_cubical_complex_base::All_cells_iterator it = range.begin() ; it != range.end() ; ++it ) + { + BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) ); + BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) ); + + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd ) + { + BOOST_CHECK( expected_boundary[bd_it] == *bd ); + ++bd_it; + } + + Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); + for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd ) + { + BOOST_CHECK( expected_coboundary[cbd_it] == *cbd ); + ++cbd_it; + } + ++i; + } +} + + + + + + +BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check) +{ + std::vector< double > expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + std::vector<unsigned> expected_dimension; + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + + std::vector<size_t> expected_boundary; + expected_boundary.push_back(0); + expected_boundary.push_back(2); + expected_boundary.push_back(2); + expected_boundary.push_back(4); + expected_boundary.push_back(0); + expected_boundary.push_back(10); + expected_boundary.push_back(1); + expected_boundary.push_back(11); + expected_boundary.push_back(5); + expected_boundary.push_back(7); + expected_boundary.push_back(2); + expected_boundary.push_back(12); + expected_boundary.push_back(3); + expected_boundary.push_back(13); + expected_boundary.push_back(7); + expected_boundary.push_back(9); + expected_boundary.push_back(4); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(12); + expected_boundary.push_back(12); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(20); + expected_boundary.push_back(11); + expected_boundary.push_back(21); + expected_boundary.push_back(15); + expected_boundary.push_back(17); + expected_boundary.push_back(12); + expected_boundary.push_back(22); + expected_boundary.push_back(13); + expected_boundary.push_back(23); + expected_boundary.push_back(17); + expected_boundary.push_back(19); + expected_boundary.push_back(14); + expected_boundary.push_back(24); + expected_boundary.push_back(20); + expected_boundary.push_back(22); + expected_boundary.push_back(22); + expected_boundary.push_back(24); + + + std::vector<size_t> expected_coboundary; + expected_coboundary.push_back(5); + expected_coboundary.push_back(1); + expected_coboundary.push_back(6); + expected_coboundary.push_back(7); + expected_coboundary.push_back(1); + expected_coboundary.push_back(3); + expected_coboundary.push_back(8); + expected_coboundary.push_back(9); + expected_coboundary.push_back(3); + expected_coboundary.push_back(6); + expected_coboundary.push_back(6); + expected_coboundary.push_back(8); + expected_coboundary.push_back(8); + expected_coboundary.push_back(5); + expected_coboundary.push_back(15); + expected_coboundary.push_back(11); + expected_coboundary.push_back(6); + expected_coboundary.push_back(16); + expected_coboundary.push_back(7); + expected_coboundary.push_back(17); + expected_coboundary.push_back(11); + expected_coboundary.push_back(13); + expected_coboundary.push_back(8); + expected_coboundary.push_back(18); + expected_coboundary.push_back(9); + expected_coboundary.push_back(19); + expected_coboundary.push_back(13); + expected_coboundary.push_back(16); + expected_coboundary.push_back(16); + expected_coboundary.push_back(18); + expected_coboundary.push_back(18); + expected_coboundary.push_back(15); + expected_coboundary.push_back(21); + expected_coboundary.push_back(16); + expected_coboundary.push_back(17); + expected_coboundary.push_back(21); + expected_coboundary.push_back(23); + expected_coboundary.push_back(18); + expected_coboundary.push_back(19); + expected_coboundary.push_back(23); + + + + std::vector< unsigned > sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector< double > data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba( sizes , data ); + int i = 0; + int bd_it = 0; + int cbd_it = 0; + + Bitmap_cubical_complex_base::All_cells_range range = ba.all_cells_range(); + for ( Bitmap_cubical_complex_base::All_cells_iterator it = range.begin() ; it != range.end() ; ++it ) + { + BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) ); + BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) ); + + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd ) + { + BOOST_CHECK( expected_boundary[bd_it] == *bd ); + ++bd_it; + } + + Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); + for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd ) + { + BOOST_CHECK( expected_coboundary[cbd_it] == *cbd ); + ++cbd_it; + } + ++i; + } +} + +BOOST_AUTO_TEST_CASE(Top_dimensional_cells_iterator_range_check) +{ + std::vector< double > expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + + std::vector< unsigned > sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector< double > data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba( sizes , data ); + int i = 0; + + Bitmap_cubical_complex_base::Top_dimensional_cells_range range = ba.top_dimensional_cells_range(); + for ( Bitmap_cubical_complex_base::Top_dimensional_cells_iterator it = range.begin() ; it != range.end() ; ++it ) + { + BOOST_CHECK( data[i] == ba.get_cell_data( *it ) ); + BOOST_CHECK( ba.get_dimension_of_a_cell( *it ) == 2 ); + ++i; + } +} + 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 21554729..f5dccbc1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,8 @@ cmake_minimum_required(VERSION 2.6) project(GUDHI) include("CMakeGUDHIVersion.txt") +# Generate GUDHI official version file +configure_file(GUDHIVersion.cmake.in "${CMAKE_SOURCE_DIR}/GUDHIVersion.cmake" @ONLY) enable_testing() @@ -77,6 +79,8 @@ else() add_subdirectory(example/Persistent_cohomology) add_subdirectory(example/Skeleton_blocker) add_subdirectory(example/Contraction) + add_subdirectory(example/Bitmap_cubical_complex) + add_subdirectory(example/Witness_complex) add_subdirectory(example/Alpha_complex) # data points generator @@ -93,8 +97,9 @@ else() # (this registers the build-tree with a global CMake-registry) export(PACKAGE GUDHI) + message("++ make install will install ${PROJECT_NAME} in the following directory : ${CMAKE_INSTALL_PREFIX}") # Create the GUDHIConfig.cmake and GUDHIConfigVersion files - set(CONF_INCLUDE_DIRS "${PROJECT_SOURCE_DIR}/include") + set(CONF_INCLUDE_DIRS "${CMAKE_INSTALL_PREFIX}/include") configure_file(GUDHIConfig.cmake.in "${PROJECT_BINARY_DIR}/GUDHIConfig.cmake" @ONLY) configure_file(GUDHIConfigVersion.cmake.in "${PROJECT_BINARY_DIR}/GUDHIConfigVersion.cmake" @ONLY) diff --git a/src/Doxyfile b/src/Doxyfile index 9cd73444..b0afb876 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -838,8 +838,9 @@ EXAMPLE_RECURSIVE = NO IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Alpha_complex/ \ doc/common/ \ - doc/Contraction/ - + doc/Contraction/ \ + doc/Witness_complex/ \ + doc/Bitmap_cubical_complex/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/src/Persistent_cohomology/concept/FilteredComplex.h b/src/Persistent_cohomology/concept/FilteredComplex.h index 949aafc2..c19698df 100644 --- a/src/Persistent_cohomology/concept/FilteredComplex.h +++ b/src/Persistent_cohomology/concept/FilteredComplex.h @@ -43,7 +43,7 @@ struct FilteredComplex * is model of IndexingTag. */ typedef unspecified Indexing_tag; -/** Returns a Simplex_hanlde that is different from all simplex handles +/** Returns a Simplex_handle that is different from all simplex handles * of the simplices. */ Simplex_handle null_simplex(); /** \brief Returns the number of simplices in the complex. @@ -61,12 +61,13 @@ struct FilteredComplex /** \brief Returns a key that is different from the keys associated * to the simplices. */ Simplex_key null_key (); -/** \brief Returns the key associated to a simplex. */ +/** \brief Returns the key associated to a simplex. + * + * This is never called on null_simplex(). */ Simplex_key key ( Simplex_handle sh ); -/** \brief Returns the simplex associated to a key. - * - * If key is different from null_key(), returns the simplex that - * has index idx in the filtration. */ +/** \brief Returns the simplex that has index idx in the filtration. + * + * This is never called on null_key(). */ Simplex_handle simplex ( Simplex_key idx ); /** \brief Assign a key to a simplex. */ void assign_key(Simplex_handle sh, Simplex_key key); diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index f60e5c72..9388fedb 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -591,7 +591,7 @@ class Simplex_tree { bool contiguous_vertices() const { if (root_.members_.empty()) return true; if (root_.members_.begin()->first != 0) return false; - if (std::prev(root_.members_.end())->first != root_.members_.size()-1) return false; + if (std::prev(root_.members_.end())->first != static_cast<Vertex_handle>(root_.members_.size() - 1)) return false; return true; } diff --git a/src/Witness_complex/concept/Simplicial_complex_for_witness.h b/src/Witness_complex/concept/Simplicial_complex_for_witness.h new file mode 100644 index 00000000..caaf0db6 --- /dev/null +++ b/src/Witness_complex/concept/Simplicial_complex_for_witness.h @@ -0,0 +1,87 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (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 CONCEPT_WITNESS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_WITNESS_H_ +#define CONCEPT_WITNESS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_WITNESS_H_ + +namespace Gudhi { + +namespace witness_complex { + +/** \brief The concept Simplicial_Complex describes the requirements + * for a type to implement a simplicial complex, + * used for example to build a 'Witness_complex'. + */ +struct SimplicialComplexForWitness { + /** Handle to specify a simplex. */ + typedef unspecified Simplex_handle; + /** Handle to specify a vertex. Must be a non-negative integer. */ + typedef unspecified Vertex_handle; + + /** Returns a Simplex_hanlde that is different from all simplex handles + * of the simplices. */ + Simplex_handle null_simplex(); + + /** \brief Iterator over the simplices of the complex, + * in an arbitrary order. + * + * 'value_type' must be 'Simplex_handle'.*/ + typedef unspecified Complex_simplex_range; + + /** + * \brief Returns a range over all the simplices of a + * complex. + */ + Complex_simplex_range complex_simplex_range(); + + /** \brief Iterator over vertices of a simplex. + * + * 'value type' must be 'Vertex_handle'.*/ + typedef unspecified Simplex_vertex_range; + + /** \brief Returns a range over vertices of a given + * simplex. */ + Simplex_vertex_range simplex_vertex_range(Simplex_handle const & simplex); + + /** \brief Return type of an insertion of a simplex + */ + typedef unspecified Insertion_result_type; + + /** \brief Inserts a simplex with vertices from a given range + * 'vertex_range' in the simplicial complex. + * */ + template< typedef Input_vertex_range > + Insertion_result_type insert_simplex(Input_vertex_range const & vertex_range); + + /** \brief Finds a simplex with vertices given by a range + * + * If a simplex exists, its Simplex_handle is returned. + * Otherwise null_simplex() is returned. */ + template< typedef Input_vertex_range > + Simplex_handle find(Input_vertex_range const & vertex_range); +}; + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // CONCEPT_WITNESS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_WITNESS_H_ diff --git a/src/Witness_complex/doc/Witness_complex_doc.h b/src/Witness_complex/doc/Witness_complex_doc.h new file mode 100644 index 00000000..e9f78170 --- /dev/null +++ b/src/Witness_complex/doc/Witness_complex_doc.h @@ -0,0 +1,38 @@ +#ifndef WITNESS_COMPLEX_DOC_H_ +#define WITNESS_COMPLEX_DOC_H_ + +/** + \defgroup witness_complex Witness complex + + \author Siargey Kachanovich + + \section Definitions + + Witness complex \f$ Wit(W,L) \f$ is a simplicial complex defined on two sets of points in \f$\mathbb{R}^D\f$: + + \li \f$W\f$ set of **witnesses** and + \li \f$L \subseteq W\f$ set of **landmarks**. + + The simplices are based on landmarks + and a simplex belongs to the witness complex if and only if it is witnessed, that is: + + \f$ \sigma \subset L \f$ is witnessed if there exists a point \f$w \in W\f$ such that + w is closer to the vertices of \f$ \sigma \f$ than other points in \f$ L \f$ and all of its faces are witnessed as well. + + \section Implementation + + The principal class of this module is Gudhi::Witness_complex. + + In both cases, the constructor for this class takes a {witness}x{closest_landmarks} table, where each row represents a witness and consists of landmarks sorted by distance to this witness. + This table can be constructed by two additional classes Landmark_choice_by_furthest_point and Landmark_choice_by_random_point also included in the module. + + *\image html "bench_Cy8.png" "Running time as function on number of landmarks" width=10cm + *\image html "bench_sphere.png" "Running time as function on number of witnesses for |L|=300" width=10cm + + + \copyright GNU General Public License v3. + + + */ + +#endif // WITNESS_COMPLEX_DOC_H_ diff --git a/src/Witness_complex/doc/bench_Cy8.png b/src/Witness_complex/doc/bench_Cy8.png Binary files differnew file mode 100644 index 00000000..d9045294 --- /dev/null +++ b/src/Witness_complex/doc/bench_Cy8.png diff --git a/src/Witness_complex/doc/bench_sphere.png b/src/Witness_complex/doc/bench_sphere.png Binary files differnew file mode 100644 index 00000000..ba6bb381 --- /dev/null +++ b/src/Witness_complex/doc/bench_sphere.png diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt new file mode 100644 index 00000000..b304479e --- /dev/null +++ b/src/Witness_complex/example/CMakeLists.txt @@ -0,0 +1,44 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIWitnessComplex) + +# A simple example + add_executable( witness_complex_from_file witness_complex_from_file.cpp ) + add_test( witness_complex_from_bunny ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CXX_FLAGS are overwritten. + # cf. http://doc.cgal.org/latest/Manual/installation.html#title40 + # A workaround is to add "-std=c++11" again. + # A fix would be to use https://cmake.org/cmake/help/v3.1/prop_gbl/CMAKE_CXX_KNOWN_FEATURES.html + # or even better https://cmake.org/cmake/help/v3.1/variable/CMAKE_CXX_STANDARD.html + # but it implies to use cmake version 3.1 at least. + if(NOT MSVC) + include(CheckCXXCompilerFlag) + CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11) + if(COMPILER_SUPPORTS_CXX11) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") + endif() + endif() + # - End of workaround + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + message(STATUS "Eigen3 use file: ${EIGEN3_USE_FILE}.") + include_directories (BEFORE "../../include") + + add_executable ( witness_complex_sphere witness_complex_sphere.cpp ) + target_link_libraries(witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test( witness_complex_sphere_10 ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere 10) + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for witness_complex_sphere example.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile witness_complex_sphere example. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Witness_complex/example/generators.h b/src/Witness_complex/example/generators.h new file mode 100644 index 00000000..ac445261 --- /dev/null +++ b/src/Witness_complex/example/generators.h @@ -0,0 +1,147 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ +#define EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> + +#include <fstream> +#include <string> +#include <vector> + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::FT FT; +typedef K::Point_d Point_d; +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; + +/** + * \brief Rock age method of reading off file + * + */ +inline void +off_reader_cust(std::string file_name, std::vector<Point_d> & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + // Line OFF. No need in it + if (!getline(in_file, line)) { + std::cerr << "No line OFF\n"; + return; + } + // Line with 3 numbers. No need + if (!getline(in_file, line)) { + std::cerr << "No line with 3 numbers\n"; + return; + } + // Reading points + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + points.push_back(Point_d(point)); + } + in_file.close(); +} + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust(std::string file_name, Point_Vector & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +/** \brief Generate points on a grid in a cube of side 2 + * having {+-1}^D as vertices and insert them in W. + * The grid has "width" points on each side. + * If torus is true then it is supposed that the cube represents + * a flat torus, hence the opposite borders are associated. + * The points on border in this case are not placed twice. + */ +void generate_points_grid(Point_Vector& W, int width, int D, bool torus) { + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) { + std::vector<double> point; + int cell_i = i; + for (int l = 0; l < D; ++l) { + if (torus) + point.push_back(-1 + (2.0 / (width - 1))*(cell_i % width)); + else + point.push_back(-1 + (2.0 / width)*(cell_i % width)); + // attention: the bottom and the right are covered too! + cell_i /= width; + } + W.push_back(point); + } +} + +/** \brief Generate nbP points uniformly in a cube of side 2 + * having {+-1}^dim as its vertices and insert them in W. + */ +void generate_points_random_box(Point_Vector& W, int nbP, int dim) { + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) { + W.push_back(*rp++); + } +} + +/** \brief Generate nbP points uniformly on a (dim-1)-sphere + * and insert them in W. + */ +void generate_points_sphere(Point_Vector& W, int nbP, int dim) { + CGAL::Random_points_on_sphere_d<Point_d> rp(dim, 1); + for (int i = 0; i < nbP; i++) + W.push_back(*rp++); +} + +#endif // EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ diff --git a/src/Witness_complex/example/witness_complex_from_file.cpp b/src/Witness_complex/example/witness_complex_from_file.cpp new file mode 100644 index 00000000..53207ad2 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -0,0 +1,100 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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/>. + */ + +#include <sys/types.h> +#include <sys/stat.h> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/reader_utils.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <string> +#include <vector> + +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< std::vector <double> > Point_Vector; + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + if (point.size() != 1) + points.push_back(point); + } + in_file.close(); +} + +int main(int argc, char * const argv[]) { + if (argc != 3) { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file nbL \n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]); + clock_t start, end; + + // Construct the Simplex Tree + Gudhi::Simplex_tree<> simplex_tree; + + // Read the point file + Point_Vector point_vector; + read_points_cust(file_name, point_vector); + std::cout << "Successfully read " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + // Choose landmarks + start = clock(); + std::vector<std::vector< int > > knn; + Gudhi::witness_complex::landmark_choice_by_random_point(point_vector, nbL, knn); + end = clock(); + std::cout << "Landmark choice for " << nbL << " landmarks took " + << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; + + // Compute witness complex + start = clock(); + Gudhi::witness_complex::witness_complex(knn, nbL, point_vector[0].size(), simplex_tree); + end = clock(); + std::cout << "Witness complex took " + << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; +} diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp new file mode 100644 index 00000000..b26c9f36 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -0,0 +1,90 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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/>. + */ +#define BOOST_PARAMETER_MAX_ARITY 12 + + +#include <sys/types.h> +#include <sys/stat.h> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/reader_utils.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <utility> +#include <string> +#include <vector> + +#include "generators.h" + +/** Write a gnuplot readable file. + * Data range is a random access range of pairs (arg, value) + */ +template < typename Data_range > +void write_data(Data_range & data, std::string filename) { + std::ofstream ofs(filename, std::ofstream::out); + for (auto entry : data) + ofs << entry.first << ", " << entry.second << "\n"; + ofs.close(); +} + +int main(int argc, char * const argv[]) { + if (argc != 2) { + std::cerr << "Usage: " << argv[0] + << " number_of_landmarks \n"; + return 0; + } + + int number_of_landmarks = atoi(argv[1]); + clock_t start, end; + + // Construct the Simplex Tree + Gudhi::Simplex_tree<> simplex_tree; + + std::vector< std::pair<int, double> > l_time; + + // Read the point file + for (int nbP = 500; nbP < 10000; nbP += 500) { + Point_Vector point_vector; + generate_points_sphere(point_vector, nbP, 4); + std::cout << "Successfully generated " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + // Choose landmarks + start = clock(); + std::vector<std::vector< int > > knn; + Gudhi::witness_complex::landmark_choice_by_random_point(point_vector, number_of_landmarks, knn); + + // Compute witness complex + Gudhi::witness_complex::witness_complex(knn, number_of_landmarks, point_vector[0].size(), simplex_tree); + end = clock(); + double time = static_cast<double>(end - start) / CLOCKS_PER_SEC; + std::cout << "Witness complex for " << number_of_landmarks << " landmarks took " + << time << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; + l_time.push_back(std::make_pair(nbP, time)); + } + write_data(l_time, "w_time.dat"); +} diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h new file mode 100644 index 00000000..df93155b --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h @@ -0,0 +1,105 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ +#define LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ + +#include <boost/range/size.hpp> + +#include <limits> // for numeric_limits<> +#include <iterator> +#include <algorithm> // for sort +#include <vector> + +namespace Gudhi { + +namespace witness_complex { + + typedef std::vector<int> typeVectorVertex; + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by iteratively adding the furthest witness from the + * current landmark set as the new landmark. + * \details It chooses nbL landmarks from a random access range `points` and + * writes {witness}*{closest landmarks} matrix in `knn`. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges + * + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_furthest_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nb_points = boost::size(points); + assert(nb_points >= nbL); + // distance matrix witness x landmarks + std::vector<std::vector<double>> wit_land_dist(nb_points, std::vector<double>()); + // landmark list + typeVectorVertex chosen_landmarks; + + knn = KNearestNeighbours(nb_points, std::vector<int>()); + int current_number_of_landmarks = 0; // counter for landmarks + double curr_max_dist = 0; // used for defining the furhest point from L + const double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nb_points, infty); // vector of current distances to L from points + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + // or better yet std::uniform_int_distribution + int rand_int = rand() % nb_points; + int curr_max_w = rand_int; // For testing purposes a pseudo-random number is used here + + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + // curr_max_w at this point is the next landmark + chosen_landmarks.push_back(curr_max_w); + unsigned i = 0; + for (auto& p : points) { + double curr_dist = euclidean_distance(p, *(std::begin(points) + chosen_landmarks[current_number_of_landmarks])); + wit_land_dist[i].push_back(curr_dist); + knn[i].push_back(current_number_of_landmarks); + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + ++i; + } + curr_max_dist = 0; + for (i = 0; i < dist_to_L.size(); i++) + if (dist_to_L[i] > curr_max_dist) { + curr_max_dist = dist_to_L[i]; + curr_max_w = i; + } + } + for (int i = 0; i < nb_points; ++i) + std::sort(std::begin(knn[i]), + std::end(knn[i]), + [&wit_land_dist, i](int a, int b) { + return wit_land_dist[i][a] < wit_land_dist[i][b]; }); + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h new file mode 100644 index 00000000..ebf6aad1 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h @@ -0,0 +1,96 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 LANDMARK_CHOICE_BY_RANDOM_POINT_H_ +#define LANDMARK_CHOICE_BY_RANDOM_POINT_H_ + +#include <boost/range/size.hpp> + +#include <queue> // for priority_queue<> +#include <utility> // for pair<> +#include <iterator> +#include <vector> +#include <set> + +namespace Gudhi { + +namespace witness_complex { + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by taking random vertices for landmarks. + * \details It chooses nbL distinct landmarks from a random access range `points` + * and outputs a matrix {witness}*{closest landmarks} in knn. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges and + * Vertex_handle is the label type of a vertex in a simplicial complex. + * Closest_landmark_range needs to have push_back operation. + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_random_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nbP = boost::size(points); + assert(nbP >= nbL); + std::set<int> landmarks; + int current_number_of_landmarks = 0; // counter for landmarks + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + int chosen_landmark = rand() % nbP; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + while (landmarks.find(chosen_landmark) != landmarks.end()) + chosen_landmark = rand() % nbP; + landmarks.insert(chosen_landmark); + } + + int dim = boost::size(*std::begin(points)); + typedef std::pair<double, int> dist_i; + typedef bool (*comp)(dist_i, dist_i); + knn = KNearestNeighbours(nbP); + for (int points_i = 0; points_i < nbP; points_i++) { + std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([](dist_i j1, dist_i j2) { + return j1.first > j2.first; + }); + std::set<int>::iterator landmarks_it; + int landmarks_i = 0; + for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end(); + ++landmarks_it, landmarks_i++) { + dist_i dist = std::make_pair(euclidean_distance(points[points_i], points[*landmarks_it]), landmarks_i); + l_heap.push(dist); + } + for (int i = 0; i < dim + 1; i++) { + dist_i dist = l_heap.top(); + knn[points_i].push_back(dist.second); + l_heap.pop(); + } + } + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_RANDOM_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h new file mode 100644 index 00000000..489cdf11 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -0,0 +1,265 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 WITNESS_COMPLEX_H_ +#define WITNESS_COMPLEX_H_ + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +#include <boost/range/size.hpp> + +#include <gudhi/distance_functions.h> + +#include <algorithm> +#include <utility> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <ctime> +#include <iostream> + +namespace Gudhi { + +namespace witness_complex { + +// /* +// * \private +// \class Witness_complex +// \brief Constructs the witness complex for the given set of witnesses and landmarks. +// \ingroup witness_complex +// */ +template< class SimplicialComplex> +class Witness_complex { + private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } + }; + + private: + typedef typename SimplicialComplex::Simplex_handle Simplex_handle; + typedef typename SimplicialComplex::Vertex_handle Vertex_handle; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL_; // Number of landmarks + SimplicialComplex& sc_; // Simplicial complex + + public: + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /* @name Constructor + */ + + //@{ + + // Witness_range<Closest_landmark_range<Vertex_handle>> + + /* + * \brief Iterative construction of the witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn of + * nearest neighbours of the form {witnesses}x{landmarks}. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Constructor takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + template< typename KNearestNeighbors > + Witness_complex(KNearestNeighbors const & knn, + int nbL, + int dim, + SimplicialComplex & sc) : nbL_(nbL), sc_(sc) { + // Construction of the active witness list + int nbW = boost::size(knn); + typeVectorVertex vv; + int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + ActiveWitnessList active_w; // = new ActiveWitnessList(); + for (Vertex_handle i = 0; i != nbL_; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + counter++; + vv = {i}; + sc_.insert_simplex(vv); + // TODO(SK) Error if not inserted : normally no need here though + } + int k = 1; /* current dimension in iterative construction */ + for (int i = 0; i != nbW; ++i) + active_w.push_back(i); + while (!active_w.empty() && k < dim) { + typename ActiveWitnessList::iterator it = active_w.begin(); + while (it != active_w.end()) { + typeVectorVertex simplex_vector; + /* THE INSERTION: Checking if all the subfaces are in the simplex tree*/ + bool ok = all_faces_in(knn, *it, k); + if (ok) { + for (int i = 0; i != k + 1; ++i) + simplex_vector.push_back(knn[*it][i]); + sc_.insert_simplex(simplex_vector); + // TODO(SK) Error if not inserted : normally no need here though + ++it; + } else { + active_w.erase(it++); // First increase the iterator and then erase the previous element + } + } + k++; + } + } + + //@} + + private: + /* \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + template <typename KNearestNeighbors> + bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) { + std::vector< Vertex_handle > facet; + // CHECK ALL THE FACETS + for (int i = 0; i != k + 1; ++i) { + facet = {}; + for (int j = 0; j != k + 1; ++j) { + if (j != i) { + facet.push_back(knn[witness_id][j]); + } + } // endfor + if (sc_.find(facet) == sc_.null_simplex()) + return false; + } // endfor + return true; + } + + template <typename T> + static void print_vector(const std::vector<T>& v) { + std::cout << "["; + if (!v.empty()) { + std::cout << *(v.begin()); + for (auto it = v.begin() + 1; it != v.end(); ++it) { + std::cout << ","; + std::cout << *it; + } + } + std::cout << "]"; + } + + public: + // /* + // * \brief Verification if every simplex in the complex is witnessed by witnesses in knn. + // * \param print_output =true will print the witnesses for each simplex + // * \remark Added for debugging purposes. + // */ + template< class KNearestNeighbors > + bool is_witness_complex(KNearestNeighbors const & knn, bool print_output) { + for (Simplex_handle sh : sc_.complex_simplex_range()) { + bool is_witnessed = false; + typeVectorVertex simplex; + int nbV = 0; // number of verticed in the simplex + for (Vertex_handle v : sc_.simplex_vertex_range(sh)) + simplex.push_back(v); + nbV = simplex.size(); + for (typeVectorVertex w : knn) { + bool has_vertices = true; + for (Vertex_handle v : simplex) + if (std::find(w.begin(), w.begin() + nbV, v) == w.begin() + nbV) { + has_vertices = false; + } + if (has_vertices) { + is_witnessed = true; + if (print_output) { + std::cout << "The simplex "; + print_vector(simplex); + std::cout << " is witnessed by the witness "; + print_vector(w); + std::cout << std::endl; + } + break; + } + } + if (!is_witnessed) { + if (print_output) { + std::cout << "The following simplex is not witnessed "; + print_vector(simplex); + std::cout << std::endl; + } + assert(is_witnessed); + return false; + } + } + return true; + } +}; + + /** + * \ingroup witness_complex + * \brief Iterative construction of the witness complex. + * \details The witness complex is written in simplicial complex sc_ + * basing on a matrix knn of + * nearest neighbours of the form {witnesses}x{landmarks}. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Procedure takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + template <class KNearestNeighbors, class SimplicialComplexForWitness> + void witness_complex(KNearestNeighbors const & knn, + int nbL, + int dim, + SimplicialComplexForWitness & sc) { + Witness_complex<SimplicialComplexForWitness>(knn, nbL, dim, sc); + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // WITNESS_COMPLEX_H_ diff --git a/src/Witness_complex/test/CMakeLists.txt b/src/Witness_complex/test/CMakeLists.txt new file mode 100644 index 00000000..9422c152 --- /dev/null +++ b/src/Witness_complex/test/CMakeLists.txt @@ -0,0 +1,34 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIWitnessComplexUT) + +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 ( simple_witness_complexUT simple_witness_complex.cpp ) +target_link_libraries(simple_witness_complexUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +# Unitary tests definition and xml result file generation +add_test(NAME simple_witness_complexUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/simple_witness_complexUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/simple_witness_complexUT.xml --log_level=test_suite --report_level=no) + +add_executable ( witness_complex_pointsUT witness_complex_points.cpp ) +target_link_libraries(witness_complex_pointsUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +# Unitary tests definition and xml result file generation +add_test(NAME witness_complex_pointsUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_pointsUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/witness_complex_pointsUT.xml --log_level=test_suite --report_level=no) + diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp new file mode 100644 index 00000000..03df78ee --- /dev/null +++ b/src/Witness_complex/test/simple_witness_complex.cpp @@ -0,0 +1,59 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "simple_witness_complex" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> + +#include <iostream> +#include <ctime> +#include <vector> + +typedef Gudhi::Simplex_tree<> Simplex_tree; +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef Gudhi::witness_complex::Witness_complex<Simplex_tree> WitnessComplex; + +BOOST_AUTO_TEST_CASE(simple_witness_complex) { + Simplex_tree complex; + std::vector< typeVectorVertex > knn; + + knn.push_back({1, 0, 5, 2, 6, 3, 4}); + knn.push_back({2, 6, 4, 5, 0, 1, 3}); + knn.push_back({3, 4, 2, 1, 5, 6, 0}); + knn.push_back({4, 2, 1, 3, 5, 6, 0}); + knn.push_back({5, 1, 6, 0, 2, 3, 4}); + knn.push_back({6, 0, 5, 2, 1, 3, 4}); + knn.push_back({0, 5, 6, 1, 2, 3, 4}); + knn.push_back({2, 6, 4, 5, 3, 1, 0}); + knn.push_back({1, 2, 5, 4, 3, 6, 0}); + knn.push_back({3, 4, 0, 6, 5, 1, 2}); + knn.push_back({5, 0, 1, 3, 6, 2, 4}); + knn.push_back({5, 6, 1, 0, 2, 3, 4}); + knn.push_back({1, 6, 0, 5, 2, 3, 4}); + WitnessComplex witnessComplex(knn, 7, 7, complex); + + BOOST_CHECK(witnessComplex.is_witness_complex(knn, false)); +} diff --git a/src/Witness_complex/test/witness_complex_points.cpp b/src/Witness_complex/test/witness_complex_points.cpp new file mode 100644 index 00000000..bd3df604 --- /dev/null +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -0,0 +1,64 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "witness_complex_points" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/Landmark_choice_by_furthest_point.h> + +#include <iostream> +#include <vector> + +typedef std::vector<double> Point; +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef Gudhi::Simplex_tree<> Simplex_tree; +typedef Gudhi::witness_complex::Witness_complex<Simplex_tree> WitnessComplex; + +BOOST_AUTO_TEST_CASE(witness_complex_points) { + std::vector< typeVectorVertex > knn; + std::vector< Point > points; + // Add grid points as witnesses + for (double i = 0; i < 10; i += 1.0) + for (double j = 0; j < 10; j += 1.0) + for (double k = 0; k < 10; k += 1.0) + points.push_back(Point({i, j, k})); + + bool b_print_output = false; + // First test: random choice + Simplex_tree complex1; + Gudhi::witness_complex::landmark_choice_by_random_point(points, 100, knn); + assert(!knn.empty()); + WitnessComplex witnessComplex1(knn, 100, 3, complex1); + BOOST_CHECK(witnessComplex1.is_witness_complex(knn, b_print_output)); + + // Second test: furthest choice + knn.clear(); + Simplex_tree complex2; + Gudhi::witness_complex::landmark_choice_by_furthest_point(points, 100, knn); + WitnessComplex witnessComplex2(knn, 100, 3, complex2); + BOOST_CHECK(witnessComplex2.is_witness_complex(knn, b_print_output)); +} |
