From f0ca2ba2080c8f2b01a7cd459560980dfcb964f8 Mon Sep 17 00:00:00 2001 From: pdlotko Date: Fri, 11 Dec 2015 12:44:20 +0000 Subject: Pictures in png format git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bitmap@943 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 470f82a103e17ffbda2b784fa829e5b0b11458e6 --- src/Bitmap_cubical_complex/doc/exampleBitmap.png | Bin 0 -> 9594 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 src/Bitmap_cubical_complex/doc/exampleBitmap.png (limited to 'src/Bitmap_cubical_complex/doc/exampleBitmap.png') diff --git a/src/Bitmap_cubical_complex/doc/exampleBitmap.png b/src/Bitmap_cubical_complex/doc/exampleBitmap.png new file mode 100644 index 00000000..f8207473 Binary files /dev/null and b/src/Bitmap_cubical_complex/doc/exampleBitmap.png differ -- cgit v1.2.3 From 84399987baac2817e58bf9f5e18ded6aa6893b0f Mon Sep 17 00:00:00 2001 From: pdlotko Date: Tue, 9 Feb 2016 13:26:47 +0000 Subject: adding missing partsy git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bitmap@1008 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f6de1ee4317763b50233f9a7212bdbf6587ee686 --- .../doc/Gudhi_Cubical_Complex_doc.h | 41 +- src/Bitmap_cubical_complex/doc/exampleBitmap.png | Bin 9594 -> 2549 bytes .../include/gudhi/Bitmap_cubical_complex.h | 28 +- .../include/gudhi/Bitmap_cubical_complex_base.h | 114 +++++- ...cal_complex_periodic_boundary_conditions_base.h | 451 +++++++++++---------- .../include/gudhi/Compute_persistence_with_phat.h | 242 +++++++++++ .../gudhi/phat/algorithms/chunk_reduction.h | 240 +++++++++++ .../include/gudhi/phat/algorithms/row_reduction.h | 55 +++ .../gudhi/phat/algorithms/standard_reduction.h | 45 ++ .../gudhi/phat/algorithms/twist_reduction.h | 50 +++ .../include/gudhi/phat/boundary_matrix.h | 336 +++++++++++++++ .../include/gudhi/phat/compute_persistence_pairs.h | 69 ++++ .../include/gudhi/phat/helpers/dualize.h | 63 +++ .../include/gudhi/phat/helpers/misc.h | 74 ++++ .../gudhi/phat/helpers/thread_local_storage.h | 52 +++ .../include/gudhi/phat/persistence_pairs.h | 151 +++++++ .../phat/representations/abstract_pivot_column.h | 158 ++++++++ .../phat/representations/bit_tree_pivot_column.h | 169 ++++++++ .../gudhi/phat/representations/full_pivot_column.h | 81 ++++ .../phat/representations/sparse_pivot_column.h | 62 +++ .../gudhi/phat/representations/vector_list.h | 98 +++++ .../gudhi/phat/representations/vector_set.h | 100 +++++ .../gudhi/phat/representations/vector_vector.h | 93 +++++ src/Bitmap_cubical_complex/test/Bitmap_test.cpp | 12 +- 24 files changed, 2532 insertions(+), 252 deletions(-) create mode 100644 src/Bitmap_cubical_complex/include/gudhi/Compute_persistence_with_phat.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/chunk_reduction.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/row_reduction.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/standard_reduction.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/twist_reduction.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/boundary_matrix.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/compute_persistence_pairs.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/helpers/dualize.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/helpers/misc.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/helpers/thread_local_storage.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/persistence_pairs.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/representations/abstract_pivot_column.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/representations/bit_tree_pivot_column.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/representations/full_pivot_column.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/representations/sparse_pivot_column.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_list.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_set.h create mode 100644 src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_vector.h (limited to 'src/Bitmap_cubical_complex/doc/exampleBitmap.png') diff --git a/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h index 1a6310fb..c06678a1 100644 --- a/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h +++ b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h @@ -38,20 +38,18 @@ namespace Cubical_complex *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 elementary interval 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 non-degenerated, while the second one is \a degenerated interval. A boundary of a elementary -*interval 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 degenerated elementary interval. An elementary cube \f$ C \f$ is a +* An elementary interval 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 non-degenerate, while the second one is \a degenerate interval. A boundary of a elementary +*interval is a chain \f$\partial [n,n+1] = [n+1,n+1]-[n,n] \f$ in case of non-degenerate elementary interval and \f$\partial [n,n] = 0 \f$ in case of degenerate elementary interval. An elementary cube \f$ C \f$ is a -*product of elementary intervals, \f$C=I_1 \times \ldots \times I_n\f$. Embedding dimension of a cube is n, the number of elementary intervals (degenerated or not) in the product. A dimension of a cube \f$C=I_1 \times ... \times I_n\f$ is the -*number of non degenerated elementary intervals in the product. A boundary of a cube \f$C=I_1 \times \ldots \times I_n\f$ is a chain obtained in the following way: +*product of elementary intervals, \f$C=I_1 \times \ldots \times I_n\f$. Embedding dimension of a cube is n, the number of elementary intervals (degenerate or not) in the product. A dimension of a cube \f$C=I_1 \times ... \times I_n\f$ is the +*number of non degenerate elementary intervals in the product. A boundary of a cube \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 cubical complex \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 maximal 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 . +*For further details and theory of cubical complexes, please consult \cite kaczynski2004computational as well as the following paper \cite peikert2012topological . * *\section datastructure Data structure. * @@ -73,7 +71,8 @@ namespace Cubical_complex *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. +*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." @@ -84,6 +83,29 @@ namespace Cubical_complex 3 3 1 +4 +6 +8 +20 +4 +7 +6 +5 +\endverbatim + +\section 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 2 3 8 @@ -94,6 +116,9 @@ namespace Cubical_complex 5 \endverbatim +Indicate that we have imposed periodic boundary conditions in the direction x, but not in the direction y. + + */ /** @} */ // end defgroup cubical_complex diff --git a/src/Bitmap_cubical_complex/doc/exampleBitmap.png b/src/Bitmap_cubical_complex/doc/exampleBitmap.png index f8207473..069c6eb2 100644 Binary files a/src/Bitmap_cubical_complex/doc/exampleBitmap.png and b/src/Bitmap_cubical_complex/doc/exampleBitmap.png differ diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index 63edcadd..82ea8672 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -1,4 +1,4 @@ - /* This file is part of the Gudhi Library. The Gudhi library +/* This file is part of the Gudhi Library. The Gudhi library * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * @@ -65,7 +65,8 @@ public: **/ Bitmap_cubical_complex( const char* perseus_style_file ): T(perseus_style_file),key_associated_to_simplex(this->total_number_of_cells+1) - { + { + //clock_t begin = clock(); if ( globalDbg ){cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n";} for ( size_t i = 0 ; i != this->total_number_of_cells ; ++i ) { @@ -74,7 +75,8 @@ public: //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(); + this->initialize_simplex_associated_to_key(); + //cerr << "Time of running Bitmap_cubical_complex( const char* perseus_style_file ) constructor : " << double(clock() - begin) / CLOCKS_PER_SEC << endl; } @@ -115,7 +117,8 @@ public: //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(); - } + } + //*********************************************// //Other 'easy' functions @@ -328,6 +331,15 @@ public: **/ Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { + /* + std::vector< size_t > bdry = this->get_boundary_of_a_cell(sh); + Boundary_simplex_range result( bdry.size() ); + for ( size_t i = 0 ; i != bdry.size() ; ++i ) + { + result[i] = this->simplex_associated_to_key[ bdry[i] ]; + } + return result; + */ return this->get_boundary_of_a_cell(sh); } @@ -500,7 +512,13 @@ void Bitmap_cubical_complex::initialize_simplex_associated_to_key() std::iota (std::begin(simplex_associated_to_key), std::end(simplex_associated_to_key), 0); std::sort( simplex_associated_to_key.begin() , simplex_associated_to_key.end() , - is_before_in_filtration(this) ); + is_before_in_filtration(this) ); + + //we still need to deal here with a key_associated_to_simplex: + for ( size_t i = 0 ; i != simplex_associated_to_key.size() ; ++i ) + { + this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i; + } } diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 600f250d..22b703a9 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -29,7 +29,8 @@ #include #include #include -#include +#include +#include #include "Bitmap_cubical_complex/counter.h" @@ -68,7 +69,9 @@ public: /** *Default constructor **/ - Bitmap_cubical_complex_base(){} + Bitmap_cubical_complex_base() + { + } /** * There are a few constructors of a Bitmap_cubical_complex_base class. * First one, that takes vector, creates an empty bitmap of a dimension equal @@ -110,7 +113,7 @@ public: * non-negative integers pointing to the * positions of (co)boundary element of the input cell. **/ - inline std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const; + 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. @@ -140,7 +143,7 @@ public: /** * Returns number of all cubes in the data structure. **/ - inline unsigned size_of_bitmap()const + inline unsigned size()const { return this->data.size(); } @@ -149,7 +152,19 @@ public: * Writing to stream operator. **/ template - friend ostream& operator << ( ostream & os , const Bitmap_cubical_complex_base& b ); + friend ostream& operator << ( ostream & os , const Bitmap_cubical_complex_base& b ); + + + /** + * Functions that put the input data to bins. + **/ + void put_data_toBins( size_t number_of_bins ); + void put_data_toBins( T diameter_of_bin ); + + /** + * Functions to find min and max values of filtration. + **/ + std::pair< T ,T > min_max_filtration(); //ITERATORS @@ -157,22 +172,41 @@ public: * Iterator through all cells in the complex (in order they appear in the structure -- i.e. * in lexicographical order). **/ - typedef typename std::vector< T >::iterator all_cells_iterator; - all_cells_iterator all_cells_begin()const + typedef typename std::vector< T >::iterator all_cells_iterator; + + /** + * Function returning an iterator to the first cell of the bitmap. + **/ + all_cells_iterator all_cells_begin() { return this->data.begin(); - } + } + + /** + * Function returning an iterator to the last cell of the bitmap. + **/ all_cells_iterator all_cells_end()const { return this->data.end(); } - - typedef typename std::vector< T >::const_iterator all_cells_const_iterator; + /** + * Constant iterator through all cells in the complex (in order they appear in the structure -- i.e. + * in lexicographical order). + **/ + typedef typename std::vector< T >::const_iterator all_cells_const_iterator; + + /** + * Function returning a constant iterator to the first cell of the bitmap. + **/ all_cells_const_iterator all_cells_const_begin()const { return this->data.begin(); - } + } + + /** + * Function returning a constant iterator to the last cell of the bitmap. + **/ all_cells_const_iterator all_cells_const_end()const { return this->data.end(); @@ -269,12 +303,20 @@ public: protected: std::vector< size_t > counter; Bitmap_cubical_complex_base& b; - }; + }; + + /** + * Function returning a Top_dimensional_cells_iterator to the first top dimensional cell cell of the bitmap. + **/ Top_dimensional_cells_iterator top_dimensional_cells_begin() { Top_dimensional_cells_iterator a(*this); return a; - } + } + + /** + * Function returning a Top_dimensional_cells_iterator to the last top dimensional cell cell of the bitmap. + **/ Top_dimensional_cells_iterator top_dimensional_cells_end() { Top_dimensional_cells_iterator a(*this); @@ -351,6 +393,50 @@ protected: }; +template +void Bitmap_cubical_complex_base::put_data_toBins( size_t number_of_bins ) +{ + bool bdg = false; + + std::pair< T ,T > min_max = this->min_max_filtration(); + T dx = (min_max.second-min_max.first)/(T)number_of_bins; + + //now put the data into the appropriate bins: + for ( size_t i = 0 ; i != this->data.size() ; ++i ) + { + if ( bdg ){cerr << "Before binning : " << this->data[i] << endl;} + this->data[i] = min_max.first + dx*(this->data[i]-min_max.first)/number_of_bins; + if ( bdg ){cerr << "After binning : " << this->data[i] << endl;getchar();} + } +} + +template +void Bitmap_cubical_complex_base::put_data_toBins( T diameter_of_bin ) +{ + bool bdg = false; + std::pair< T ,T > min_max = this->min_max_filtration(); + + size_t number_of_bins = (min_max.second - min_max.first)/diameter_of_bin; + //now put the data into the appropriate bins: + for ( size_t i = 0 ; i != this->data.size() ; ++i ) + { + if ( bdg ){cerr << "Before binning : " << this->data[i] << endl;} + this->data[i] = min_max.first + diameter_of_bin*(this->data[i]-min_max.first)/number_of_bins; + if ( bdg ){cerr << "After binning : " << this->data[i] << endl;getchar();} + } +} + +template +std::pair< T ,T > Bitmap_cubical_complex_base::min_max_filtration() +{ + std::pair< T ,T > min_max( std::numeric_limits::max() , std::numeric_limits::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 @@ -422,7 +508,7 @@ void Bitmap_cubical_complex_base::read_perseus_style_file( const char* perseu unsigned dimensionOfData; inFiltration >> dimensionOfData; - if (dbg){cerr << "dimensionOfData : " << dimensionOfData << endl;} + if (dbg){cerr << "dimensionOfData : " << dimensionOfData << endl;getchar();} std::vector sizes; sizes.reserve( dimensionOfData ); diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h index 28911ea8..38c218dc 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h @@ -1,120 +1,120 @@ -#pragma once -#include -#include "Bitmap_cubical_complex_base.h" - +#pragma once +#include +#include "Bitmap_cubical_complex_base.h" + using namespace std; namespace Gudhi { namespace Cubical_complex -{ - -//in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic boundary conditions). But, we set up the iterators and the procedures -//to compute boundary and coboundary in the way that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all the cells on the -//right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base - -template -class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base -{ -public: - //constructors that take an extra parameter: - Bitmap_cubical_complex_periodic_boundary_conditions_base(){}; - Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); - Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseusStyleFile ); - Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); - - //overwritten methods co compute boundary and coboundary - virtual std::vector< size_t > get_boundary_of_a_cell( size_t cell )const; - std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const; - //inline unsigned get_dimension_of_a_cell( size_t cell )const; - -protected: - std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed; +{ + +//in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic boundary conditions). But, we set up the iterators and the procedures +//to compute boundary and coboundary in the way that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all the cells on the +//right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base + +template +class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base +{ +public: + //constructors that take an extra parameter: + Bitmap_cubical_complex_periodic_boundary_conditions_base(){}; + Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); + Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseusStyleFile ); + Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); + + //overwritten methods co compute boundary and coboundary + virtual std::vector< size_t > get_boundary_of_a_cell( size_t cell )const; + std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const; + //inline unsigned get_dimension_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& 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; + + 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(multiplier,std::numeric_limits::max()); this->total_number_of_cells = multiplier; - } - Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes ); - Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells ); - void construct_complex_based_on_top_dimensional_cells( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); -}; - -template -void Bitmap_cubical_complex_periodic_boundary_conditions_base::construct_complex_based_on_top_dimensional_cells( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ) -{ - this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed; - this->set_up_containers( dimensions ); - - size_t i = 0; - for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base::Top_dimensional_cells_iterator it = this->top_dimensional_cells_begin() ; it != this->top_dimensional_cells_end() ; ++it ) - { - *it = topDimensionalCells[i]; - ++i; - } - this->impose_lower_star_filtration(); -} - -template -Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes , 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 -Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseus_style_file ) -{ - //for Perseus style files: - bool dbg = false; - + } + Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes ); + Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells ); + void construct_complex_based_on_top_dimensional_cells( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); +}; + +template +void Bitmap_cubical_complex_periodic_boundary_conditions_base::construct_complex_based_on_top_dimensional_cells( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ) +{ + this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed; + this->set_up_containers( dimensions ); + + size_t i = 0; + for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base::Top_dimensional_cells_iterator it = this->top_dimensional_cells_begin() ; it != this->top_dimensional_cells_end() ; ++it ) + { + *it = topDimensionalCells[i]; + ++i; + } + this->impose_lower_star_filtration(); +} + +template +Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes , 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 +Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseus_style_file ) +{ + + //for Perseus style files: + bool dbg = false; ifstream inFiltration; inFiltration.open( perseus_style_file ); unsigned dimensionOfData; - inFiltration >> dimensionOfData; - - this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( dimensionOfData , false ); - - std::vector sizes; + inFiltration >> dimensionOfData; + + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( dimensionOfData , false ); + + std::vector 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; + 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 ); - + this->set_up_containers( sizes ); + typename Bitmap_cubical_complex_periodic_boundary_conditions_base::Top_dimensional_cells_iterator it(*this); it = this->top_dimensional_cells_begin(); while ( !inFiltration.eof() ) { double filtrationLevel; - inFiltration >> filtrationLevel; - if ( inFiltration.eof() )break; + inFiltration >> filtrationLevel; + if ( inFiltration.eof() )break; if ( dbg ) { @@ -128,124 +128,137 @@ Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_comp ++it; } inFiltration.close(); - this->impose_lower_star_filtration(); - - /* - char* filename = (char*)perseus_style_file; - //char* filename = "combustionWithPeriodicBoundaryConditions/v0/tV0_000000.float"; - ifstream file( filename , ios::binary | ios::ate ); - unsigned realSizeOfFile = file.tellg(); - file.close(); - realSizeOfFile = realSizeOfFile/sizeof(T); - - unsigned w, h, d; - - w = h = d = ceil(pow( realSizeOfFile , (double)(1/(double)3) )); - - T* slice = new T[w*h*d]; - if (slice == NULL) - { - cerr << "Allocation error, cannot allocate " << w*h*d*sizeof(T) << " bytes to store the data from the file. The program will now terminate \n"; - exit(EXIT_FAILURE); - } - - FILE* fp; - if ((fp=fopen( filename, "rb" )) == NULL ) - { - cerr << "Cannot open the file: " << filename << ". The program will now terminate \n"; - exit(1); - } - fread( slice,4,w*h*d,fp ); - fclose(fp); - std::vector data(slice,slice+w*h*d); - delete[] slice; - std::vector< unsigned > sizes; - sizes.push_back(w); - sizes.push_back(w); - sizes.push_back(w); - - std::vector< bool > directions; - directions.push_back( true ); - directions.push_back( true ); - directions.push_back( true ); - Bitmap_cubical_complex_periodic_boundary_conditions_base b( sizes, data, directions ); - *this = b; - */ -} - -template -Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes ) -{ - this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( sizes.size() , false ); - this->set_up_containers( sizes ); -} - -template -Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells ) -{ - std::vector directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( dimensions.size() , false ); - this->construct_complex_based_on_top_dimensional_cells( dimensions , topDimensionalCells , directions_in_which_periodic_b_cond_are_to_be_imposed ); -} - - - - - -template -Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells , 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 -std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base::get_boundary_of_a_cell( size_t cell )const -{ - bool dbg = false; - if ( dbg ){cerr << "Computations of boundary of a cell : " << cell << endl;} - + this->impose_lower_star_filtration(); + +/* + char* filename = (char*)perseus_style_file; + //char* filename = "combustionWithPeriodicBoundaryConditions/v0/tV0_000000.float"; + ifstream file( filename , ios::binary | ios::ate ); + unsigned realSizeOfFile = file.tellg(); + file.close(); + realSizeOfFile = realSizeOfFile/sizeof(T); + + unsigned w, h, d; + + w = h = d = ceil(pow( realSizeOfFile , (double)(1/(double)3) )); + + T* slice = new T[w*h*d]; + if (slice == NULL) + { + cerr << "Allocation error, cannot allocate " << w*h*d*sizeof(T) << " bytes to store the data from the file. The program will now terminate \n"; + exit(EXIT_FAILURE); + } + + FILE* fp; + if ((fp=fopen( filename, "rb" )) == NULL ) + { + cerr << "Cannot open the file: " << filename << ". The program will now terminate \n"; + exit(1); + } + + clock_t read_begin = clock(); + fread( slice,4,w*h*d,fp ); + fclose(fp); + cerr << "Time of reading the file : " << double(clock() - read_begin) / CLOCKS_PER_SEC << endl; + + + clock_t begin_creation_bitap = clock(); + std::vector data(slice,slice+w*h*d); + delete[] slice; + std::vector< unsigned > sizes; + sizes.push_back(w); + sizes.push_back(w); + sizes.push_back(w); + + this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true ); + this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true ); + this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true ); + this->set_up_containers( sizes ); + + size_t i = 0; + for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base::Top_dimensional_cells_iterator it = this->top_dimensional_cells_begin() ; it != this->top_dimensional_cells_end() ; ++it ) + { + *it = data[i]; + ++i; + } + this->impose_lower_star_filtration(); + cerr << "Time of creation of a bitmap : " << double(clock() - begin_creation_bitap ) / CLOCKS_PER_SEC << endl; +*/ +} + +template +Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes ) +{ + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( sizes.size() , false ); + this->set_up_containers( sizes ); +} + +template +Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells ) +{ + std::vector directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( dimensions.size() , false ); + this->construct_complex_based_on_top_dimensional_cells( dimensions , topDimensionalCells , directions_in_which_periodic_b_cond_are_to_be_imposed ); +} + + + + + +template +Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells , 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 +std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base::get_boundary_of_a_cell( size_t cell )const +{ + bool dbg = false; + if ( dbg ){cerr << "Computations of boundary of a cell : " << cell << endl;} + std::vector< size_t > boundary_elements; size_t cell1 = cell; for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) { - unsigned position = cell1/this->multipliers[i-1]; - //this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction. + unsigned position = cell1/this->multipliers[i-1]; + //this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction. if ( position%2 == 1 ) - { - //if there are no periodic boundary conditions in this direction, we do not have to do anything. - if ( !directions_in_which_periodic_b_cond_are_to_be_imposed[i-1] ) - { - //cerr << "A\n"; - boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); - boundary_elements.push_back( cell + this->multipliers[ i-1 ] ); - if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";} - } - else - { - //in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end. - if ( position != 2*this->sizes[ i-1 ]-1 ) - { - //cerr << "B\n"; - boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); - boundary_elements.push_back( cell + this->multipliers[ i-1 ] ); - if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";} - } - else - { - //cerr << "C\n"; - boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); - boundary_elements.push_back( cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] ); - if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] << " ";} - } + { + //if there are no periodic boundary conditions in this direction, we do not have to do anything. + if ( !directions_in_which_periodic_b_cond_are_to_be_imposed[i-1] ) + { + //cerr << "A\n"; + boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); + boundary_elements.push_back( cell + this->multipliers[ i-1 ] ); + if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";} + } + else + { + //in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end. + if ( position != 2*this->sizes[ i-1 ]-1 ) + { + //cerr << "B\n"; + boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); + boundary_elements.push_back( cell + this->multipliers[ i-1 ] ); + if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";} + } + else + { + //cerr << "C\n"; + boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); + boundary_elements.push_back( cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] ); + if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] << " ";} + } } } cell1 = cell1%this->multipliers[i-1]; } - return boundary_elements; -} - + return boundary_elements; +} + template std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base::get_coboundary_of_a_cell( size_t cell )const { @@ -254,12 +267,12 @@ std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_basemultipliers.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. + 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] ) - { + { + 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]) ) { @@ -268,31 +281,31 @@ std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_basesizes[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] ); - } + } } - } + else + { + //we want to have periodic boundary conditions in this direction + if ( counter[i-1] != 0 ) + { + coboundary_elements.push_back( cell - this->multipliers[i-1] ); + coboundary_elements.push_back( cell + this->multipliers[i-1] ); + } + else + { + //in this case counter[i-1] == 0. + coboundary_elements.push_back( cell + this->multipliers[i-1] ); + coboundary_elements.push_back( cell + (2*this->sizes[ i-1 ]-1)*this->multipliers[i-1] ); + } + } + } cell1 = cell1%this->multipliers[i-1]; } return coboundary_elements; -} - - - -}//Cubical_complex -}//namespace Gudhi +} + + + +}//Cubical_complex +}//namespace Gudhi diff --git a/src/Bitmap_cubical_complex/include/gudhi/Compute_persistence_with_phat.h b/src/Bitmap_cubical_complex/include/gudhi/Compute_persistence_with_phat.h new file mode 100644 index 00000000..9f4ada45 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/Compute_persistence_with_phat.h @@ -0,0 +1,242 @@ +/* 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 . + */ + +#pragma once + + +#include "phat/compute_persistence_pairs.h" +#include "phat/representations/vector_vector.h" +#include "phat/algorithms/standard_reduction.h" +#include "phat/algorithms/chunk_reduction.h" +#include "phat/algorithms/row_reduction.h" +#include "phat/algorithms/twist_reduction.h" + + +namespace Gudhi +{ + + +//the only aim of this class is to have a ability to compute persistence with phat. +template +void writeBettiNumbersAndPersistenceIntervalsToFile( char* prefix , std::pair< std::vector > , std::vector< std::vector< std::pair > > > resutsFromPhat ) +{ + std::ostringstream filenameStr; + filenameStr << prefix << "_bettiNumbers"; + std::string str = filenameStr.str(); + const char* filename = str.c_str(); + ofstream out; + out.open( filename ); + for ( size_t dim = 0 ; dim != resutsFromPhat.first.size() ; ++dim ) + { + out << "Dimension : " << dim << endl; + for ( size_t i = 0 ; i != resutsFromPhat.first[dim].size() ; ++i ) + { + out << resutsFromPhat.first[dim][i] << endl; + } + out << endl; + } + out.close(); + + + cerr << "Write persistence to file \n"; + for ( size_t dim = 0 ; dim != resutsFromPhat.second.size() ; ++dim ) + { + cerr << "resutsFromPhat.second[dim].size() : " << resutsFromPhat.second[dim].size() << endl; + if ( resutsFromPhat.second[dim].size() == 0 )continue; + std::ostringstream filenameStr; + filenameStr << prefix << "_persistence_" << dim; + std::string str = filenameStr.str(); + const char* filename = str.c_str(); + ofstream out1; + out1.open( filename ); + for ( size_t i = 0 ; i != resutsFromPhat.second[dim].size() ; ++i ) + { + out1 << resutsFromPhat.second[dim][i].first << " " << resutsFromPhat.second[dim][i].second << endl; + } + out1.close(); + } +}//writeBettiNumbersAndPersistenceIntervalsToFile + + +template +class Compute_persistence_with_phat +{ +public: + Compute_persistence_with_phat( T* data_structure_ ); + std::pair< std::vector< std::vector > , std::vector< std::vector< std::pair > > > get_the_intervals( phat::persistence_pairs pairs ); + + phat::persistence_pairs compute_persistence_pairs_dualized_chunk_reduction(); + phat::persistence_pairs compute_persistence_pairs_twist_reduction(); + phat::persistence_pairs compute_persistence_pairs_standard_reduction(); + //phat::persistence_pairs compute_persistence_pairs_spectral_sequence_reduction(); +private: + void print_bd_matrix(); + phat::boundary_matrix< phat::vector_vector > boundary_matrix; + T* data_structure; +}; + +template +void Compute_persistence_with_phat::print_bd_matrix() +{ + std::cout << "The boundary matrix has " << this->boundary_matrix.get_num_cols() << " columns: " << std::endl; + for( phat::index col_idx = 0; col_idx < this->boundary_matrix.get_num_cols(); col_idx++ ) { + std::cout << "Colum " << col_idx << " represents a cell of dimension " << (int)this->boundary_matrix.get_dim( col_idx ) << ". "; + if( !this->boundary_matrix.is_empty( col_idx ) ) { + std::vector< phat::index > temp_col; + this->boundary_matrix.get_col( col_idx, temp_col ); + std::cout << "Its boundary consists of the cells"; + for( phat::index idx = 0; idx < (phat::index)temp_col.size(); idx++ ) + std::cout << " " << temp_col[ idx ]; + } + std::cout << std::endl; + } +} + +template +phat::persistence_pairs Compute_persistence_with_phat::compute_persistence_pairs_dualized_chunk_reduction() +{ + phat::persistence_pairs pairs; + phat::compute_persistence_pairs_dualized< phat::chunk_reduction >( pairs, this->boundary_matrix ); + return pairs; +} + +template +phat::persistence_pairs Compute_persistence_with_phat::compute_persistence_pairs_twist_reduction() +{ + phat::persistence_pairs pairs; + phat::compute_persistence_pairs< phat::twist_reduction >( pairs, this->boundary_matrix ); + return pairs; +} + +template +phat::persistence_pairs Compute_persistence_with_phat::compute_persistence_pairs_standard_reduction() +{ + phat::persistence_pairs pairs; + phat::compute_persistence_pairs< phat::standard_reduction >( pairs, this->boundary_matrix ); + return pairs; +} + +//template +//phat::persistence_pairs Compute_persistence_with_phat::compute_persistence_pairs_spectral_sequence_reduction() +//{ +// phat::persistence_pairs pairs; +// phat::compute_persistence_pairs< phat::spectral_sequence_reduction >( pairs, this->boundary_matrix ); +// return pairs; +//} + +template +Compute_persistence_with_phat::Compute_persistence_with_phat( T* data_structure_ ):data_structure( data_structure_ ) +{ + bool dbg = false; + this->boundary_matrix.set_num_cols( this->data_structure->num_simplices() ); + + //setting up the dimensions of cells: + for ( size_t i = 0 ; i != this->data_structure->num_simplices() ; ++i ) + { + this->boundary_matrix.set_dim( i, this->data_structure->dimension( this->data_structure->simplex(i) ) ); + } + + + //now it is time to set up the boundary matrix: + typename T::Filtration_simplex_range range = this->data_structure->filtration_simplex_range(); + std::vector< phat::index > temp_col; + for ( typename T::Filtration_simplex_iterator it = range.begin() ; it != range.end() ; ++it ) + { + typename T::Boundary_simplex_range boundary_range = this->data_structure->boundary_simplex_range( *it ); + for ( typename T::Boundary_simplex_iterator bd = boundary_range.begin() ; bd != boundary_range.end() ; ++bd ) + { + temp_col.push_back( this->data_structure->key( *bd ) ); + } + //we do not know if the boundary elements are sorted according to filtration, that is why I am enforcing it here: + this->boundary_matrix.set_col( this->data_structure->key( *it ) , temp_col ); + temp_col.clear(); + } +} + +template +std::pair< std::vector< std::vector > , std::vector< std::vector< std::pair > > > Compute_persistence_with_phat::get_the_intervals( phat::persistence_pairs pairs ) +{ + bool dbg = false; + //in order to find the birth times of the infinite homology classes, we need to know which elements are not paired. To search for them, we will use this vector: + std::vector isTheElementPaired( this->data_structure->num_simplices() , false ); + + //now it is time to recover the finite persistence pairs and the Betti numbers: + std::vector< std::vector< std::pair > > finitePersistencePairs( this->data_structure->dimension() ); + for( phat::index idx = 0; idx < pairs.get_num_pairs(); idx++ ) + { + typename T::Simplex_key positionOfBeginOfInterval = pairs.get_pair( idx ).first; + typename T::Simplex_key positionOfEndOfInterval = pairs.get_pair( idx ).second; + + typename T::Simplex_handle first_simplex = this->data_structure->simplex(positionOfBeginOfInterval); + typename T::Simplex_handle second_simplex = this->data_structure->simplex(positionOfEndOfInterval); + + typename T::Filtration_value valueFirst = this->data_structure->filtration( first_simplex ); + typename T::Filtration_value valueSecond = this->data_structure->filtration( second_simplex ); + + if ( valueFirst > valueSecond ){std::swap( valueFirst , valueSecond );} + + unsigned dimFirst = this->data_structure->dimension(first_simplex); + unsigned dimSecond = this->data_structure->dimension(second_simplex); + unsigned dim = std::min( dimFirst , dimSecond ); + + + //we are ignoring trivial barcodes + if ( valueFirst != valueSecond ) + { + finitePersistencePairs[ dim ].push_back( std::make_pair(valueFirst , valueSecond) ); + if ( dbg ){cerr << "Adding barcode : " << valueFirst << "," << valueSecond << endl;} + } + + //isTheElementPaired[ positionOfBeginOfIntervalInBitmap ] = true; + //isTheElementPaired[ positionOfEndOfIntervalInBitmap ] = true; + isTheElementPaired[ pairs.get_pair( idx ).first ] = true; + isTheElementPaired[ pairs.get_pair( idx ).second ] = true; + } + + + std::vector< std::vector > birthTimesOfInfinitePersistnceClasses(this->data_structure->dimension()+1 ); + for ( size_t i = 0 ; i != this->data_structure->dimension()+1 ; ++i ) + { + std::vector v; + birthTimesOfInfinitePersistnceClasses[i] = v; + } + for ( size_t i = 0 ; i != isTheElementPaired.size() ; ++i ) + { + if ( isTheElementPaired[i] == false ) + { + //i-th element is not paired, therefore it gives an infinite class + typename T::Simplex_handle simplex = this->data_structure->simplex(i); + birthTimesOfInfinitePersistnceClasses[ this->data_structure->dimension( simplex ) ].push_back( this->data_structure->filtration(simplex) ); + } + } + + //sorting finite persistence pairs: + for ( size_t dim = 0 ; dim != finitePersistencePairs.size() ; ++dim ) + { + std::sort( finitePersistencePairs[dim].begin() , finitePersistencePairs[dim].end() ); + } + return std::make_pair( birthTimesOfInfinitePersistnceClasses , finitePersistencePairs ); +}//Compute_persistence_with_phat + + + +}//namespace Gudhi diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/chunk_reduction.h b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/chunk_reduction.h new file mode 100644 index 00000000..352392a8 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/chunk_reduction.h @@ -0,0 +1,240 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "../helpers/misc.h" +#include "../boundary_matrix.h" + +namespace phat { + class chunk_reduction { + public: + enum column_type { GLOBAL + , LOCAL_POSITIVE + , LOCAL_NEGATIVE }; + + public: + template< typename Representation > + void operator() ( boundary_matrix< Representation >& boundary_matrix ) { + + const index nr_columns = boundary_matrix.get_num_cols(); + const dimension max_dim = boundary_matrix.get_max_dim(); + + std::vector< index > lowest_one_lookup( nr_columns, -1 ); + std::vector < column_type > column_type( nr_columns, GLOBAL ); + std::vector< char > is_active( nr_columns, false ); + + std::vector chunk_boundaries; + _get_chunks( boundary_matrix, chunk_boundaries ); + + // Phase 1: Reduce chunks locally -- 1st pass + #pragma omp parallel for schedule( guided, 1 ) + for( index chunk_id = 0; chunk_id < (index) chunk_boundaries.size() - 2; chunk_id += 2 ) + _local_chunk_reduction( boundary_matrix, lowest_one_lookup, column_type, max_dim, + chunk_boundaries[chunk_id], chunk_boundaries[chunk_id+2] - 1 ); + boundary_matrix.sync(); + + // Phase 1: Reduce chunks locally -- 2nd pass + #pragma omp parallel for schedule( guided, 1 ) + for( index chunk_id = 1; chunk_id < (index) chunk_boundaries.size() - 2; chunk_id += 2 ) + _local_chunk_reduction( boundary_matrix, lowest_one_lookup, column_type, max_dim, + chunk_boundaries[chunk_id], chunk_boundaries[chunk_id+2] - 1 ); + boundary_matrix.sync(); + + // get global columns + std::vector< index > global_columns; + for( index cur_col_idx = 0; cur_col_idx < nr_columns; cur_col_idx++ ) + if( column_type[ cur_col_idx ] == GLOBAL ) + global_columns.push_back( cur_col_idx ); + + // get active columns + #pragma omp parallel for + for( index idx = 0; idx < (index)global_columns.size(); idx++ ) + is_active[ global_columns[ idx ] ] = true; + _get_active_columns( boundary_matrix, lowest_one_lookup, column_type, global_columns, is_active ); + + // Phase 2+3: Simplify columns and reduce them + for( dimension cur_dim = max_dim; cur_dim >= 1; cur_dim-- ) { + // Phase 2: Simplify columns + std::vector< index > temp_col; + #pragma omp parallel for schedule( guided, 1 ), private( temp_col ) + for( index idx = 0; idx < (index)global_columns.size(); idx++ ) + if( boundary_matrix.get_dim( global_columns[ idx ] ) == cur_dim ) + _global_column_simplification( global_columns[ idx ], boundary_matrix, lowest_one_lookup, column_type, is_active, temp_col ); + boundary_matrix.sync(); + + // Phase 3: Reduce columns + for( index idx = 0; idx < (index)global_columns.size(); idx++ ) { + index cur_col = global_columns[ idx ]; + if( boundary_matrix.get_dim( cur_col ) == cur_dim && column_type[ cur_col ] == GLOBAL ) { + index lowest_one = boundary_matrix.get_max_index( cur_col ); + while( lowest_one != -1 && lowest_one_lookup[ lowest_one ] != -1 ) { + boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); + lowest_one = boundary_matrix.get_max_index( cur_col ); + } + if( lowest_one != -1 ) { + lowest_one_lookup[ lowest_one ] = cur_col; + boundary_matrix.clear( lowest_one ); + } + } + } + } + + boundary_matrix.sync(); + } + + protected: + template< typename Representation > + void _get_chunks( const boundary_matrix< Representation >& boundary_matrix + , std::vector< index >& chunk_boundaries) + { + chunk_boundaries.clear(); + std::vector temp_chunk_boundaries; + const index nr_columns = boundary_matrix.get_num_cols(); + + // size of chuks = sqrt(N) + const index chunk_size = (index) sqrt( (float)nr_columns ); + + // size of chunks = N / num_threads + //const index chunk_size = nr_columns / omp_get_max_threads(); + + for ( index cur_col = 0; cur_col < nr_columns; cur_col++ ) + if( cur_col % chunk_size == 0 ) + temp_chunk_boundaries.push_back( cur_col ); + temp_chunk_boundaries.push_back( nr_columns ); + + // subdivide chunks for interleaved 2 pass appraoch + for( index chunk_id = 0; chunk_id < (index) temp_chunk_boundaries.size(); chunk_id ++ ) { + chunk_boundaries.push_back( temp_chunk_boundaries[ chunk_id ] ); + if( chunk_id < (index) temp_chunk_boundaries.size() - 1 ) { + index midPoint = ( temp_chunk_boundaries[ chunk_id ] + temp_chunk_boundaries[ chunk_id + 1 ] ) / 2; + chunk_boundaries.push_back( midPoint ); + } + } + } + + template< typename Representation > + void _local_chunk_reduction( boundary_matrix< Representation >& boundary_matrix + , std::vector& lowest_one_lookup + , std::vector< column_type >& column_type + , const dimension max_dim + , const index chunk_begin + , const index chunk_end ) { + for( dimension cur_dim = max_dim; cur_dim >= 1; cur_dim-- ) { + for( index cur_col = chunk_begin; cur_col <= chunk_end; cur_col++ ) { + if( column_type[ cur_col ] == GLOBAL && boundary_matrix.get_dim( cur_col ) == cur_dim ) { + index lowest_one = boundary_matrix.get_max_index( cur_col ); + while( lowest_one != -1 && lowest_one >= chunk_begin && lowest_one_lookup[ lowest_one ] != -1 ) { + boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); + lowest_one = boundary_matrix.get_max_index( cur_col ); + } + if( lowest_one >= chunk_begin ) { + lowest_one_lookup[ lowest_one ] = cur_col; + column_type[ cur_col ] = LOCAL_NEGATIVE; + column_type[ lowest_one ] = LOCAL_POSITIVE; + boundary_matrix.clear( lowest_one ); + } + } + } + } + } + + template< typename Representation > + void _get_active_columns( const boundary_matrix< Representation >& boundary_matrix + , const std::vector< index >& lowest_one_lookup + , const std::vector< column_type >& column_type + , const std::vector< index >& global_columns + , std::vector< char >& is_active ) { + + const index nr_columns = boundary_matrix.get_num_cols(); + std::vector< char > finished( nr_columns, false ); + + std::vector< std::pair < index, index > > stack; + std::vector< index > cur_col_values; + #pragma omp parallel for schedule( guided, 1 ), private( stack, cur_col_values ) + for( index idx = 0; idx < (index)global_columns.size(); idx++ ) { + bool pop_next = false; + index start_col = global_columns[ idx ]; + stack.push_back( std::pair< index, index >( start_col, -1 ) ); + while( !stack.empty() ) { + index cur_col = stack.back().first; + index prev_col = stack.back().second; + if( pop_next ) { + stack.pop_back(); + pop_next = false; + if( prev_col != -1 ) { + if( is_active[ cur_col ] ) { + is_active[ prev_col ] = true; + } + if( prev_col == stack.back().first ) { + finished[ prev_col ] = true; + pop_next = true; + } + } + } else { + pop_next = true; + boundary_matrix.get_col( cur_col, cur_col_values ); + for( index idx = 0; idx < (index) cur_col_values.size(); idx++ ) { + index cur_row = cur_col_values[ idx ]; + if( ( column_type[ cur_row ] == GLOBAL ) ) { + is_active[ cur_col ] = true; + } else if( column_type[ cur_row ] == LOCAL_POSITIVE ) { + index next_col = lowest_one_lookup[ cur_row ]; + if( next_col != cur_col && !finished[ cur_col ] ) { + stack.push_back( std::make_pair( next_col, cur_col ) ); + pop_next = false; + } + } + } + } + } + } + } + + template< typename Representation > + void _global_column_simplification( const index col_idx + , boundary_matrix< Representation >& boundary_matrix + , const std::vector< index >& lowest_one_lookup + , const std::vector< column_type >& column_type + , const std::vector< char >& is_active + , std::vector< index >& temp_col ) + { + temp_col.clear(); + while( !boundary_matrix.is_empty( col_idx ) ) { + index cur_row = boundary_matrix.get_max_index( col_idx ); + switch( column_type[ cur_row ] ) { + case GLOBAL: + temp_col.push_back( cur_row ); + boundary_matrix.remove_max( col_idx ); + break; + case LOCAL_NEGATIVE: + boundary_matrix.remove_max( col_idx ); + break; + case LOCAL_POSITIVE: + if( is_active[ lowest_one_lookup[ cur_row ] ] ) + boundary_matrix.add_to( lowest_one_lookup[ cur_row ], col_idx ); + else + boundary_matrix.remove_max( col_idx ); + break; + } + } + std::reverse( temp_col.begin(), temp_col.end() ); + boundary_matrix.set_col( col_idx, temp_col ); + } + }; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/row_reduction.h b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/row_reduction.h new file mode 100644 index 00000000..2047cafd --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/row_reduction.h @@ -0,0 +1,55 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "../helpers/misc.h" +#include "../boundary_matrix.h" + +namespace phat { + class row_reduction { + public: + template< typename Representation > + void operator() ( boundary_matrix< Representation >& boundary_matrix ) { + + const index nr_columns = boundary_matrix.get_num_cols(); + std::vector< std::vector< index > > lowest_one_lookup( nr_columns ); + + for( index cur_col = nr_columns - 1; cur_col >= 0; cur_col-- ) { + if( !boundary_matrix.is_empty( cur_col ) ) + lowest_one_lookup[ boundary_matrix.get_max_index( cur_col ) ].push_back( cur_col ); + + if( !lowest_one_lookup[ cur_col ].empty() ) { + boundary_matrix.clear( cur_col ); + std::vector< index >& cols_with_cur_lowest = lowest_one_lookup[ cur_col ]; + index source = *min_element( cols_with_cur_lowest.begin(), cols_with_cur_lowest.end() ); + for( index idx = 0; idx < (index)cols_with_cur_lowest.size(); idx++ ) { + index target = cols_with_cur_lowest[ idx ]; + if( target != source && !boundary_matrix.is_empty( target ) ) { + boundary_matrix.add_to( source, target ); + if( !boundary_matrix.is_empty( target ) ) { + index lowest_one_of_target = boundary_matrix.get_max_index( target ); + lowest_one_lookup[ lowest_one_of_target ].push_back( target ); + } + } + } + } + } + } + }; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/standard_reduction.h b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/standard_reduction.h new file mode 100644 index 00000000..b2c91a85 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/standard_reduction.h @@ -0,0 +1,45 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "../helpers/misc.h" +#include "../boundary_matrix.h" + +namespace phat { + class standard_reduction { + public: + template< typename Representation > + void operator() ( boundary_matrix< Representation >& boundary_matrix ) { + + const index nr_columns = boundary_matrix.get_num_cols(); + std::vector< index > lowest_one_lookup( nr_columns, -1 ); + + for( index cur_col = 0; cur_col < nr_columns; cur_col++ ) { + index lowest_one = boundary_matrix.get_max_index( cur_col ); + while( lowest_one != -1 && lowest_one_lookup[ lowest_one ] != -1 ) { + boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); + lowest_one = boundary_matrix.get_max_index( cur_col ); + } + if( lowest_one != -1 ) { + lowest_one_lookup[ lowest_one ] = cur_col; + } + } + } + }; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/twist_reduction.h b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/twist_reduction.h new file mode 100644 index 00000000..1bdd8de2 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/algorithms/twist_reduction.h @@ -0,0 +1,50 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "../helpers/misc.h" +#include "../boundary_matrix.h" + +namespace phat { + class twist_reduction { + public: + template< typename Representation > + void operator () ( boundary_matrix< Representation >& boundary_matrix ) { + + const index nr_columns = boundary_matrix.get_num_cols(); + std::vector< index > lowest_one_lookup( nr_columns, -1 ); + + for( index cur_dim = boundary_matrix.get_max_dim(); cur_dim >= 1 ; cur_dim-- ) { + for( index cur_col = 0; cur_col < nr_columns; cur_col++ ) { + if( boundary_matrix.get_dim( cur_col ) == cur_dim ) { + index lowest_one = boundary_matrix.get_max_index( cur_col ); + while( lowest_one != -1 && lowest_one_lookup[ lowest_one ] != -1 ) { + boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); + lowest_one = boundary_matrix.get_max_index( cur_col ); + } + if( lowest_one != -1 ) { + lowest_one_lookup[ lowest_one ] = cur_col; + boundary_matrix.clear( lowest_one ); + } + } + } + } + } + }; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/boundary_matrix.h b/src/Bitmap_cubical_complex/include/gudhi/phat/boundary_matrix.h new file mode 100644 index 00000000..941537da --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/boundary_matrix.h @@ -0,0 +1,336 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "helpers/misc.h" +#include "representations/bit_tree_pivot_column.h" + +// interface class for the main data structure -- implementations of the interface can be found in ./representations +namespace phat { + template< class Representation = bit_tree_pivot_column > + class boundary_matrix + { + + protected: + Representation rep; + + // interface functions -- actual implementation and complexity depends on chosen @Representation template + public: + // get overall number of columns in boundary_matrix + index get_num_cols() const { return rep._get_num_cols(); } + + // set overall number of columns in boundary_matrix + void set_num_cols( index nr_of_columns ) { rep._set_num_cols( nr_of_columns ); } + + // get dimension of given index + dimension get_dim( index idx ) const { return rep._get_dim( idx ); } + + // set dimension of given index + void set_dim( index idx, dimension dim ) { rep._set_dim( idx, dim ); } + + // replaces content of @col with boundary of given index + void get_col( index idx, column& col ) const { rep._get_col( idx, col ); } + + // set column @idx to the values contained in @col + void set_col( index idx, const column& col ) { rep._set_col( idx, col ); } + + // true iff boundary of given column is empty + bool is_empty( index idx ) const { return rep._is_empty( idx ); } + + // largest index of given column (new name for lowestOne()) + index get_max_index( index idx ) const { return rep._get_max_index( idx ); } + + // removes maximal index from given column + void remove_max( index idx ) { rep._remove_max( idx ); } + + // adds column @source to column @target' + void add_to( index source, index target ) { rep._add_to( source, target ); } + + // clears given column + void clear( index idx ) { rep._clear( idx ); } + + // syncronizes all internal data structures -- has to be called before and after any multithreaded access! + void sync() { rep._sync(); } + + // info functions -- independent of chosen 'Representation' + public: + // maximal dimension + dimension get_max_dim() const { + dimension cur_max_dim = 0; + for( index idx = 0; idx < get_num_cols(); idx++ ) + cur_max_dim = get_dim( idx ) > cur_max_dim ? get_dim( idx ) : cur_max_dim; + return cur_max_dim; + } + + // number of nonzero rows for given column @idx + index get_num_rows( index idx ) const { + column cur_col; + get_col( idx, cur_col ); + return cur_col.size(); + } + + // maximal number of nonzero rows of all columns + index get_max_col_entries() const { + index max_col_entries = -1; + const index nr_of_columns = get_num_cols(); + for( index idx = 0; idx < nr_of_columns; idx++ ) + max_col_entries = get_num_rows( idx ) > max_col_entries ? get_num_rows( idx ) : max_col_entries; + return max_col_entries; + } + + // maximal number of nonzero cols of all rows + index get_max_row_entries() const { + size_t max_row_entries = 0; + const index nr_of_columns = get_num_cols(); + std::vector< std::vector< index > > transposed_matrix( nr_of_columns ); + column temp_col; + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) { + get_col( cur_col, temp_col ); + for( index idx = 0; idx < (index)temp_col.size(); idx++) + transposed_matrix[ temp_col[ idx ] ].push_back( cur_col ); + } + for( index idx = 0; idx < nr_of_columns; idx++ ) + max_row_entries = transposed_matrix[ idx ].size() > max_row_entries ? transposed_matrix[ idx ].size() : max_row_entries; + return max_row_entries; + } + + // overall number of entries in the matrix + index get_num_entries() const { + index number_of_nonzero_entries = 0; + const index nr_of_columns = get_num_cols(); + for( index idx = 0; idx < nr_of_columns; idx++ ) + number_of_nonzero_entries += get_num_rows( idx ); + return number_of_nonzero_entries; + } + + // operators / constructors + public: + boundary_matrix() {}; + + template< class OtherRepresentation > + boundary_matrix( const boundary_matrix< OtherRepresentation >& other ) { + *this = other; + } + + template< typename OtherRepresentation > + bool operator==( const boundary_matrix< OtherRepresentation >& other_boundary_matrix ) const { + const index number_of_columns = this->get_num_cols(); + + if( number_of_columns != other_boundary_matrix.get_num_cols() ) + return false; + + column temp_col; + column other_temp_col; + for( index idx = 0; idx < number_of_columns; idx++ ) { + this->get_col( idx, temp_col ); + other_boundary_matrix.get_col( idx, other_temp_col ); + if( temp_col != other_temp_col || this->get_dim( idx ) != other_boundary_matrix.get_dim( idx ) ) + return false; + } + return true; + } + + template< typename OtherRepresentation > + bool operator!=( const boundary_matrix< OtherRepresentation >& other_boundary_matrix ) const { + return !( *this == other_boundary_matrix ); + } + + template< typename OtherRepresentation > + boundary_matrix< Representation >& operator=( const boundary_matrix< OtherRepresentation >& other ) + { + const index nr_of_columns = other.get_num_cols(); + this->set_num_cols( nr_of_columns ); + column temp_col; + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) { + this->set_dim( cur_col, other.get_dim( cur_col ) ); + other.get_col( cur_col, temp_col ); + this->set_col( cur_col, temp_col ); + } + + // by convention, always return *this + return *this; + } + + // I/O -- independent of chosen 'Representation' + public: + + // initializes boundary_matrix from (vector, vector) pair -- untested + template< typename index_type, typename dimemsion_type > + void load_vector_vector( const std::vector< std::vector< index_type > >& input_matrix, const std::vector< dimemsion_type >& input_dims ) { + const index nr_of_columns = (index)input_matrix.size(); + this->set_num_cols( nr_of_columns ); + column temp_col; + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) { + this->set_dim( cur_col, (dimension)input_dims[ cur_col ] ); + + index num_rows = input_matrix[ cur_col ].size(); + temp_col.resize( num_rows ); + for( index cur_row = 0; cur_row < num_rows; cur_row++ ) + temp_col[ cur_row ] = (index)input_matrix[ cur_col ][ cur_row ]; + this->set_col( cur_col, temp_col ); + } + } + + template< typename index_type, typename dimemsion_type > + void save_vector_vector( std::vector< std::vector< index_type > >& output_matrix, std::vector< dimemsion_type >& output_dims ) { + const index nr_of_columns = get_num_cols(); + output_matrix.resize( nr_of_columns ); + output_dims.resize( nr_of_columns ); + column temp_col; + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) { + output_dims[ cur_col ] = (dimemsion_type)get_dim( cur_col ); + get_col( cur_col, temp_col ); + index num_rows = temp_col.size(); + output_matrix[ cur_col ].clear(); + output_matrix[ cur_col ].resize( num_rows ); + for( index cur_row = 0; cur_row < num_rows; cur_row++ ) + output_matrix[ cur_col ][ cur_row ] = (index_type)temp_col[ cur_row ]; + } + } + + // Loads the boundary_matrix from given file in ascii format + // Format: each line represents a column, first number is dimension, other numbers are the content of the column. + // Ignores empty lines and lines starting with a '#'. + bool load_ascii( std::string filename ) { + // first count number of columns: + std::string cur_line; + std::ifstream dummy( filename .c_str() ); + if( dummy.fail() ) + return false; + + index number_of_columns = 0; + while( getline( dummy, cur_line ) ) { + cur_line.erase(cur_line.find_last_not_of(" \t\n\r\f\v") + 1); + if( cur_line != "" && cur_line[ 0 ] != '#' ) + number_of_columns++; + + } + this->set_num_cols( number_of_columns ); + dummy.close(); + + std::ifstream input_stream( filename.c_str() ); + if( input_stream.fail() ) + return false; + + column temp_col; + index cur_col = -1; + while( getline( input_stream, cur_line ) ) { + cur_line.erase(cur_line.find_last_not_of(" \t\n\r\f\v") + 1); + if( cur_line != "" && cur_line[ 0 ] != '#' ) { + cur_col++; + std::stringstream ss( cur_line ); + + int64_t temp_dim; + ss >> temp_dim; + this->set_dim( cur_col, (dimension) temp_dim ); + + int64_t temp_index; + temp_col.clear(); + while( ss.good() ) { + ss >> temp_index; + temp_col.push_back( (index)temp_index ); + } + std::sort( temp_col.begin(), temp_col.end() ); + this->set_col( cur_col, temp_col ); + } + } + + input_stream.close(); + return true; + } + + // Saves the boundary_matrix to given file in ascii format + // Format: each line represents a column, first number is dimension, other numbers are the content of the column + bool save_ascii( std::string filename ) { + std::ofstream output_stream( filename.c_str() ); + if( output_stream.fail() ) + return false; + + const index nr_columns = this->get_num_cols(); + column tempCol; + for( index cur_col = 0; cur_col < nr_columns; cur_col++ ) { + output_stream << (int64_t)this->get_dim( cur_col ); + this->get_col( cur_col, tempCol ); + for( index cur_row_idx = 0; cur_row_idx < (index)tempCol.size(); cur_row_idx++ ) + output_stream << " " << tempCol[ cur_row_idx ]; + output_stream << std::endl; + } + + output_stream.close(); + return true; + } + + // Loads boundary_matrix from given file + // Format: nr_columns % dim1 % N1 % row1 row2 % ...% rowN1 % dim2 % N2 % ... + bool load_binary( std::string filename ) { + std::ifstream input_stream( filename.c_str(), std::ios_base::binary | std::ios_base::in ); + if( input_stream.fail() ) + return false; + + int64_t nr_columns; + input_stream.read( (char*)&nr_columns, sizeof( int64_t ) ); + this->set_num_cols( (index)nr_columns ); + + column temp_col; + for( index cur_col = 0; cur_col < nr_columns; cur_col++ ) { + int64_t cur_dim; + input_stream.read( (char*)&cur_dim, sizeof( int64_t ) ); + this->set_dim( cur_col, (dimension) cur_dim ); + int64_t nr_rows; + input_stream.read( (char*)&nr_rows, sizeof( int64_t ) ); + temp_col.resize( (std::size_t)nr_rows ); + for( index idx = 0; idx < nr_rows; idx++ ) { + int64_t cur_row; + input_stream.read( (char*)&cur_row, sizeof( int64_t ) ); + temp_col[ idx ] = (index)cur_row; + } + this->set_col( cur_col, temp_col ); + } + + input_stream.close(); + return true; + } + + // Saves the boundary_matrix to given file in binary format + // Format: nr_columns % dim1 % N1 % row1 row2 % ...% rowN1 % dim2 % N2 % ... + bool save_binary( std::string filename ) { + std::ofstream output_stream( filename.c_str(), std::ios_base::binary | std::ios_base::out ); + if( output_stream.fail() ) + return false; + + const int64_t nr_columns = this->get_num_cols(); + output_stream.write( (char*)&nr_columns, sizeof( int64_t ) ); + column tempCol; + for( index cur_col = 0; cur_col < nr_columns; cur_col++ ) { + int64_t cur_dim = this->get_dim( cur_col ); + output_stream.write( (char*)&cur_dim, sizeof( int64_t ) ); + this->get_col( cur_col, tempCol ); + int64_t cur_nr_rows = tempCol.size(); + output_stream.write( (char*)&cur_nr_rows, sizeof( int64_t ) ); + for( index cur_row_idx = 0; cur_row_idx < (index)tempCol.size(); cur_row_idx++ ) { + int64_t cur_row = tempCol[ cur_row_idx ]; + output_stream.write( (char*)&cur_row, sizeof( int64_t ) ); + } + } + + output_stream.close(); + return true; + } + }; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/compute_persistence_pairs.h b/src/Bitmap_cubical_complex/include/gudhi/phat/compute_persistence_pairs.h new file mode 100644 index 00000000..f5b04d5a --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/compute_persistence_pairs.h @@ -0,0 +1,69 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "persistence_pairs.h" +#include "boundary_matrix.h" +#include "helpers/dualize.h" +#include "algorithms/twist_reduction.h" + +namespace phat { + + template< typename ReductionAlgorithm, typename Representation > + void compute_persistence_pairs( persistence_pairs& pairs, boundary_matrix< Representation >& boundary_matrix ) { + ReductionAlgorithm reduce; + reduce( boundary_matrix ); + pairs.clear(); + for( index idx = 0; idx < boundary_matrix.get_num_cols(); idx++ ) { + if( !boundary_matrix.is_empty( idx ) ) { + index birth = boundary_matrix.get_max_index( idx ); + index death = idx; + pairs.append_pair( birth, death ); + } + } + } + + template< typename ReductionAlgorithm, typename Representation > + void compute_persistence_pairs_dualized( persistence_pairs& pairs, boundary_matrix< Representation >& boundary_matrix ) { + ReductionAlgorithm reduce; + const index nr_columns = boundary_matrix.get_num_cols(); + dualize( boundary_matrix ); + reduce( boundary_matrix ); + pairs.clear(); + for( index idx = 0; idx < nr_columns; idx++ ) { + if( !boundary_matrix.is_empty( idx ) ) { + index death = nr_columns - 1 - boundary_matrix.get_max_index( idx ); + index birth = nr_columns - 1 - idx; + pairs.append_pair( birth, death ); + } + } + } + + template< typename Representation > + void compute_persistence_pairs( persistence_pairs& pairs, boundary_matrix< Representation >& boundary_matrix ) { + phat::compute_persistence_pairs< twist_reduction >( pairs, boundary_matrix ); + } + + + template< typename Representation > + void compute_persistence_pairs_dualized( persistence_pairs& pairs, boundary_matrix< Representation >& boundary_matrix ) { + compute_persistence_pairs_dualized< twist_reduction >( pairs, boundary_matrix ); + } + +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/dualize.h b/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/dualize.h new file mode 100644 index 00000000..60caa05c --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/dualize.h @@ -0,0 +1,63 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "misc.h" +#include "../boundary_matrix.h" + +namespace phat { + template< typename Representation > + void dualize( boundary_matrix< Representation >& boundary_matrix ) { + + std::vector< dimension > dual_dims; + std::vector< std::vector< index > > dual_matrix; + + index nr_of_columns = boundary_matrix.get_num_cols(); + dual_matrix.resize( nr_of_columns ); + dual_dims.resize( nr_of_columns ); + + std::vector< index > dual_sizes( nr_of_columns, 0 ); + + column temp_col; + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) { + boundary_matrix.get_col( cur_col, temp_col ); + for( index idx = 0; idx < (index)temp_col.size(); idx++) + dual_sizes[ nr_of_columns - 1 - temp_col[ idx ] ]++; + } + + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) { + dual_matrix[cur_col].reserve(dual_sizes[cur_col]); + } + + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) { + boundary_matrix.get_col( cur_col, temp_col ); + for( index idx = 0; idx < (index)temp_col.size(); idx++) + dual_matrix[ nr_of_columns - 1 - temp_col[ idx ] ].push_back( nr_of_columns - 1 - cur_col ); + } + + const dimension max_dim = boundary_matrix.get_max_dim(); + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) + dual_dims[ nr_of_columns - 1 - cur_col ] = max_dim - boundary_matrix.get_dim( cur_col ); + + for( index cur_col = 0; cur_col < nr_of_columns; cur_col++ ) + std::reverse( dual_matrix[ cur_col ].begin(), dual_matrix[ cur_col ].end() ); + + boundary_matrix.load_vector_vector( dual_matrix, dual_dims ); + } +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/misc.h b/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/misc.h new file mode 100644 index 00000000..abbf8d53 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/misc.h @@ -0,0 +1,74 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +// STL includes +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// VS2008 and below unfortunately do not support stdint.h +#if defined(_MSC_VER)&& _MSC_VER < 1600 + typedef __int8 int8_t; + typedef unsigned __int8 uint8_t; + typedef __int16 int16_t; + typedef unsigned __int16 uint16_t; + typedef __int32 int32_t; + typedef unsigned __int32 uint32_t; + typedef __int64 int64_t; + typedef unsigned __int64 uint64_t; +#else + #include +#endif + +// basic types. index can be changed to int32_t to save memory on small instances +namespace phat { + typedef int64_t index; + typedef int8_t dimension; + typedef std::vector< index > column; +} + +// OpenMP (proxy) functions +#if defined _OPENMP + #include +#else + #define omp_get_thread_num() 0 + #define omp_get_max_threads() 1 + #define omp_get_num_threads() 1 + void omp_set_num_threads( int ) {}; + #include + #define omp_get_wtime() (float)clock() / (float)CLOCKS_PER_SEC +#endif + +#include "thread_local_storage.h" + + + diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/thread_local_storage.h b/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/thread_local_storage.h new file mode 100644 index 00000000..06e95c20 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/helpers/thread_local_storage.h @@ -0,0 +1,52 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "misc.h" + +// should ideally be equal to the cache line size of the CPU +#define PHAT_TLS_SPACING_FACTOR 64 + +// ThreadLocalStorage with some spacing to avoid "false sharing" (see wikipedia) +template< typename T > +class thread_local_storage +{ +public: + + thread_local_storage() : per_thread_storage( omp_get_max_threads() * PHAT_TLS_SPACING_FACTOR ) {}; + + T& operator()() { + return per_thread_storage[ omp_get_thread_num() * PHAT_TLS_SPACING_FACTOR ]; + } + + const T& operator()() const { + return per_thread_storage[ omp_get_thread_num() * PHAT_TLS_SPACING_FACTOR ]; + } + + T& operator[]( int tid ) { + return per_thread_storage[ tid * PHAT_TLS_SPACING_FACTOR ]; + } + + const T& operator[]( int tid ) const { + return per_thread_storage[ tid * PHAT_TLS_SPACING_FACTOR ]; + } + +protected: + std::vector< T > per_thread_storage; +}; diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/persistence_pairs.h b/src/Bitmap_cubical_complex/include/gudhi/phat/persistence_pairs.h new file mode 100644 index 00000000..f8006353 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/persistence_pairs.h @@ -0,0 +1,151 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "helpers/misc.h" + +namespace phat { + class persistence_pairs { + + protected: + std::vector< std::pair< index, index > > pairs; + + public: + index get_num_pairs() const { + return (index)pairs.size(); + } + + void append_pair( index birth, index death ) { + pairs.push_back( std::make_pair( birth, death ) ); + } + + std::pair< index, index > get_pair( index idx ) const { + return pairs[ idx ]; + } + + void clear() { + pairs.clear(); + } + + void sort() { + std::sort( pairs.begin(), pairs.end() ); + } + + // Loads the persistence pairs from given file in asci format + // Format: nr_pairs % newline % birth1 % death1 % newline % birth2 % death2 % newline ... + bool load_ascii( std::string filename ) { + std::ifstream input_stream( filename.c_str() ); + if( input_stream.fail() ) + return false; + + int64_t nr_pairs; + input_stream >> nr_pairs; + pairs.clear(); + for( index idx = 0; idx < nr_pairs; idx++ ) { + int64_t birth; + input_stream >> birth; + int64_t death; + input_stream >> death; + append_pair( (index)birth, (index)death ); + } + + input_stream.close(); + return true; + } + + // Saves the persistence pairs to given file in binary format + // Format: nr_pairs % newline % birth1 % death1 % newline % birth2 % death2 % newline ... + bool save_ascii( std::string filename ) { + std::ofstream output_stream( filename.c_str() ); + if( output_stream.fail() ) + return false; + + this->sort(); + output_stream << get_num_pairs() << std::endl; + for( std::size_t idx = 0; idx < pairs.size(); idx++ ) { + output_stream << pairs[idx].first << " " << pairs[idx].second << std::endl; + } + + output_stream.close(); + return true; + } + + // Loads the persistence pairs from given file in binary format + // Format: nr_pairs % birth1 % death1 % birth2 % death2 ... + bool load_binary( std::string filename ) { + std::ifstream input_stream( filename.c_str(), std::ios_base::binary | std::ios_base::in ); + if( input_stream.fail() ) + return false; + + int64_t nr_pairs; + input_stream.read( (char*)&nr_pairs, sizeof( int64_t ) ); + for( index idx = 0; idx < nr_pairs; idx++ ) { + int64_t birth; + input_stream.read( (char*)&birth, sizeof( int64_t ) ); + int64_t death; + input_stream.read( (char*)&death, sizeof( int64_t ) ); + append_pair( (index)birth, (index)death ); + } + + input_stream.close(); + return true; + } + + // Saves the persistence pairs to given file in binary format + // Format: nr_pairs % birth1 % death1 % birth2 % death2 ... + bool save_binary( std::string filename ) { + std::ofstream output_stream( filename.c_str(), std::ios_base::binary | std::ios_base::out ); + if( output_stream.fail() ) + return false; + + this->sort(); + int64_t nr_pairs = get_num_pairs(); + output_stream.write( (char*)&nr_pairs, sizeof( int64_t ) ); + for( std::size_t idx = 0; idx < pairs.size(); idx++ ) { + int64_t birth = pairs[ idx ].first; + output_stream.write( (char*)&birth, sizeof( int64_t ) ); + int64_t death = pairs[ idx ].second; + output_stream.write( (char*)&death, sizeof( int64_t ) ); + } + + output_stream.close(); + return true; + } + + bool operator==( persistence_pairs& other_pairs ) { + this->sort(); + other_pairs.sort(); + if( pairs.size() != (std::size_t)other_pairs.get_num_pairs() ) + return false; + + for( index idx = 0; idx < (index)pairs.size(); idx++ ) + if( get_pair( idx ) != other_pairs.get_pair( idx ) ) + return false; + + return true; + } + + bool operator!=( persistence_pairs& other_pairs ) { + return !( *this == other_pairs ); + } + }; + + + +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/representations/abstract_pivot_column.h b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/abstract_pivot_column.h new file mode 100644 index 00000000..41104108 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/abstract_pivot_column.h @@ -0,0 +1,158 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "../helpers/misc.h" +#include "../representations/vector_vector.h" + +namespace phat { + + // Note: We could even make the rep generic in the underlying Const representation + // But I cannot imagine that anything else than vector> would + // make sense + template< typename PivotColumn > + class abstract_pivot_column : public vector_vector { + public: + + protected: + typedef vector_vector Base; + typedef PivotColumn pivot_column; + + // For parallization purposes, it could be more than one full column + mutable thread_local_storage< pivot_column > _pivot_columns; + mutable thread_local_storage< index > pos_of_pivot_columns; + + pivot_column& get_pivot_column() const { + return _pivot_columns(); + } + + bool is_represented_by_pivot_column( index idx ) const { + return pos_of_pivot_columns() == idx; + } + + void unset_pos_of_pivot_column() { + index idx = pos_of_pivot_columns(); + if( idx != -1 ) { + _pivot_columns().get_column_and_clear( this->matrix[ idx ] ); + } + pos_of_pivot_columns() = -1; + } + + void represent_by_pivot_column( index idx ) { + pos_of_pivot_columns() = idx; + get_pivot_column().add_column( matrix[ idx ] ); + } + + public: + + void _set_num_cols( index nr_of_columns ) { + #pragma omp parallel for + for( int tid = 0; tid < omp_get_num_threads(); tid++ ) { + _pivot_columns[ tid ].init( nr_of_columns ); + pos_of_pivot_columns[ tid ] = -1; + } + Base::_set_num_cols( nr_of_columns ); + } + // replaces(!) content of 'col' with boundary of given index + void _get_col( index idx, column& col ) const { + col.clear(); + if( is_represented_by_pivot_column( idx ) ) { + pivot_column& pivot_column = get_pivot_column(); + pivot_column.get_column_and_clear( col ); + pivot_column.add_column( col ); + } else { + Base::_get_col( idx, col ); + } + } + + // true iff boundary of given idx is empty + bool _is_empty( index idx ) const { + return is_represented_by_pivot_column( idx ) ? get_pivot_column().empty() : Base::_is_empty( idx ); + } + + // largest row index of given column idx (new name for lowestOne()) + index _get_max_index( index idx ) const { + if( is_represented_by_pivot_column( idx ) ) { + pivot_column& pivot_column = get_pivot_column(); + if( pivot_column.empty() ) { + return -1; + } else { + return pivot_column.max_index(); + } + } else { + return Base::_get_max_index( idx ); + } + } + + // adds column 'source' to column 'target' + void _add_to( index source, index target ) { + if( !is_represented_by_pivot_column( target ) ) { + unset_pos_of_pivot_column(); + represent_by_pivot_column( target ); + } + get_pivot_column().add_column( matrix[source] ); + } + + // clears given column + void _clear( index idx ) { + if( is_represented_by_pivot_column( idx ) ) { + column dummy; + get_pivot_column().get_column_and_clear(dummy); + } else { + Base::_clear( idx ); + } + } + + void _set_col( index idx, const column& col ) { + if( is_represented_by_pivot_column( idx ) ) { + column dummy; + pivot_column& pivot_column = get_pivot_column(); + pivot_column.get_column_and_clear( dummy ); + pivot_column.add_column( col ); + } else { + Base::_set_col( idx, col ); + } + } + + // removes the maximal index of a column + void _remove_max( index idx ) { + _toggle( idx, _get_max_index( idx ) ); + } + + //// toggles given index pair + void _toggle( index col_idx, index row_idx ) { + if( !is_represented_by_pivot_column( col_idx ) ) { + unset_pos_of_pivot_column(); + represent_by_pivot_column( col_idx ); + } + get_pivot_column().add_index( row_idx ); + } + + // syncronizes all data structures (essential for openmp stuff) + // has to be called before and after any multithreaded access! + void _sync() { + #pragma omp parallel for + for( int tid = 0; tid < omp_get_num_threads(); tid++ ) + unset_pos_of_pivot_column(); + } + + }; +} + + diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/representations/bit_tree_pivot_column.h b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/bit_tree_pivot_column.h new file mode 100644 index 00000000..34366d6a --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/bit_tree_pivot_column.h @@ -0,0 +1,169 @@ +/* Copyright 2013 IST Austria + Contributed by: Hubert Wagner + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "../helpers/misc.h" +#include "../representations/abstract_pivot_column.h" + +namespace phat { + + // This is a bitset indexed with a 64-ary tree. Each node in the index + // has 64 bits; i-th bit says that the i-th subtree is non-empty. + // Supports practically O(1), inplace, zero-allocation: insert, remove, max_element + // and clear in O(number of ones in the bitset). + // 'add_index' is still the real bottleneck in practice. + class bit_tree_column + { + size_t offset; // data[i + offset] = ith block of the data-bitset + typedef uint64_t block_type; + std::vector< block_type > data; + + // this static is not a problem with OMP, it's initialized just after program starts + static const size_t debrujin_magic_table[64]; + + enum { block_size_in_bits = 64 }; + enum { block_shift = 6 }; + + // Some magic: http://graphics.stanford.edu/~seander/bithacks.html + // Gets the position of the rightmost bit of 'x'. 0 means the most significant bit. + // (-x)&x isolates the rightmost bit. + // The whole method is much faster than calling log2i, and very comparable to using ScanBitForward/Reverse intrinsic, + // which should be one CPU instruction, but is not portable. + size_t rightmost_pos(const block_type value) const + { + return 64 - 1 - debrujin_magic_table[((value & (-(int64_t)value))*0x07EDD5E59A4E28C2) >> 58]; + } + + public: + + void init(index num_cols) + { + int64_t n = 1; // in case of overflow + int64_t bottom_blocks_needed = (num_cols+block_size_in_bits-1)/block_size_in_bits; + int64_t upper_blocks = 1; + + // How many blocks/nodes of index needed to index the whole bitset? + while(n * block_size_in_bits < bottom_blocks_needed) + { + n *= block_size_in_bits; + upper_blocks += n; + } + + offset = upper_blocks; + data.resize(upper_blocks + bottom_blocks_needed, 0); + } + + index max_index() const + { + if (!data[0]) + return -1; + + const size_t size = data.size(); + size_t n = 0; + + while(true) + { + const block_type val = data[n]; + const size_t index = rightmost_pos(val); + const size_t newn = (n << block_shift) + index + 1; + + if (newn >= size) + { + const size_t bottom_index = n - offset; + return (bottom_index << block_shift) + index; + } + + n = newn; + } + + return -1; + } + + bool empty() const + { + return data[0] == 0; + } + + void add_index(const size_t entry) + { + const block_type ONE = 1; + const block_type block_modulo_mask = ((ONE << block_shift) - 1); + size_t index_in_level = entry >> block_shift; + size_t address = index_in_level + offset; + size_t index_in_block = entry & block_modulo_mask; + + block_type mask = (ONE << (block_size_in_bits - index_in_block - 1)); + + while(true) + { + data[address]^=mask; + + // First we check if we reached the root. + // Also, if anyone else was in this block, we don't need to update the path up. + if (!address || (data[address] & ~mask)) + return; + + index_in_block = index_in_level & block_modulo_mask; + index_in_level >>= block_shift; + --address; + address >>= block_shift; + mask = (ONE << (block_size_in_bits - index_in_block - 1)); + } + } + + void get_column_and_clear(column &out) + { + out.clear(); + while(true) + { + index mx = this->max_index(); + if (mx == -1) + break; + out.push_back(mx); + add_index(mx); + } + + std::reverse(out.begin(), out.end()); + } + + void add_column(const column &col) + { + for (size_t i = 0; i < col.size(); ++i) + { + add_index(col[i]); + } + } + + bool empty() { + return !data[0]; + } + }; + + const size_t bit_tree_column::debrujin_magic_table[64] = { + 63, 0, 58, 1, 59, 47, 53, 2, + 60, 39, 48, 27, 54, 33, 42, 3, + 61, 51, 37, 40, 49, 18, 28, 20, + 55, 30, 34, 11, 43, 14, 22, 4, + 62, 57, 46, 52, 38, 26, 32, 41, + 50, 36, 17, 19, 29, 10, 13, 21, + 56, 45, 25, 31, 35, 16, 9, 12, + 44, 24, 15, 8, 23, 7, 6, 5}; + + typedef abstract_pivot_column bit_tree_pivot_column; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/representations/full_pivot_column.h b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/full_pivot_column.h new file mode 100644 index 00000000..9e0c348d --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/full_pivot_column.h @@ -0,0 +1,81 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include +#include + +namespace phat { + class full_column { + + protected: + std::priority_queue< index > m_history; + std::vector< char > m_isInHistory; + std::vector< char > m_data; + + public: + void init( const index total_size ) { + m_data.resize( total_size, false ); + m_isInHistory.resize( total_size, false ); + } + + void add_column( const column& col ) { + for( index idx = 0; idx < (index) col.size(); idx++ ) { + add_index( col[ idx ] ); + } + } + void add_index( const index idx ) { + if( !m_isInHistory[ idx ] ) { + m_history.push( idx ); + m_isInHistory[ idx ] = true; + } + + m_data[ idx ] = !m_data[ idx ]; + } + + index max_index() { + while( m_history.size() > 0 ) { + index topIndex = m_history.top(); + if( m_data[ topIndex ] ) { + return topIndex; + } else { + m_history.pop(); + m_isInHistory[ topIndex ] = false; + } + } + + return -1; + } + + void get_column_and_clear( column& col ) { + col.clear(); + while( !empty() ) { + col.push_back( max_index() ); + add_index( max_index() ); + } + std::reverse( col.begin(), col.end() ); + } + + bool empty() { + return (max_index() == -1); + } + }; + + typedef abstract_pivot_column< full_column > full_pivot_column; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/representations/sparse_pivot_column.h b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/sparse_pivot_column.h new file mode 100644 index 00000000..c851a2b5 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/sparse_pivot_column.h @@ -0,0 +1,62 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include +#include + +namespace phat { + class sparse_column { + + protected: + std::set< index > m_data; + + public: + void init( const index total_size ) { + m_data.clear(); + } + + void add_column( const column& col ) { + for( index idx = 0; idx < (index) col.size(); idx++ ) + add_index( col[ idx ] ); + } + + void add_index( const index idx ) { + std::pair< std::set< index >::iterator, bool > result = m_data.insert( idx ); + if( result.second == false ) + m_data.erase( result.first ); + } + + index max_index() { + return m_data.empty() ? -1 : *m_data.rbegin(); + } + + void get_column_and_clear( column& col ) { + col.clear(); + col.assign( m_data.begin(), m_data.end() ); + m_data.clear(); + } + + bool empty() { + return m_data.empty(); + } + }; + + typedef abstract_pivot_column< sparse_column > sparse_pivot_column; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_list.h b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_list.h new file mode 100644 index 00000000..38e3090c --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_list.h @@ -0,0 +1,98 @@ +/* Copyright 2013 IST Austria + Contributed by: Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include + +namespace phat { + class vector_list { + + protected: + typedef std::list< index > internal_column; + std::vector< dimension > dims; + std::vector< internal_column > matrix; + + public: + // overall number of cells in boundary_matrix + index _get_num_cols() const { + return (index)matrix.size(); + } + void _set_num_cols( index nr_of_columns ) { + dims.resize( nr_of_columns ); + matrix.resize( nr_of_columns ); + } + + // dimension of given index + dimension _get_dim( index idx ) const { + return dims[ idx ]; + } + void _set_dim( index idx, dimension dim ) { + dims[ idx ] = dim; + } + + // replaces(!) content of 'col' with boundary of given index + void _get_col( index idx, column& col ) const { + col.clear(); + col.reserve( matrix[idx].size() ); + std::copy (matrix[idx].begin(), matrix[idx].end(), std::back_inserter(col) ); + } + + void _set_col( index idx, const column& col ) { + matrix[ idx ].clear(); + matrix[ idx ].resize( col.size() ); + std::copy (col.begin(), col.end(), matrix[ idx ].begin() ); + } + + // true iff boundary of given idx is empty + bool _is_empty( index idx ) const { + return matrix[ idx ].empty(); + } + + // largest row index of given column idx (new name for lowestOne()) + index _get_max_index( index idx ) const { + return matrix[ idx ].empty() ? -1 : *matrix[ idx ].rbegin(); + } + + // removes the maximal index of a column + void _remove_max( index idx ) { + internal_column::iterator it = matrix[ idx ].end(); + it--; + matrix[ idx ].erase( it ); + } + + // clears given column + void _clear( index idx ) { + matrix[ idx ].clear(); + } + + // syncronizes all data structures (essential for openmp stuff) + void _sync() {} + + // adds column 'source' to column 'target' + void _add_to( index source, index target ) { + internal_column& source_col = matrix[ source ]; + internal_column& target_col = matrix[ target ]; + internal_column temp_col; + target_col.swap( temp_col ); + std::set_symmetric_difference( temp_col.begin(), temp_col.end(), + source_col.begin(), source_col.end(), + std::back_inserter( target_col ) ); + } + }; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_set.h b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_set.h new file mode 100644 index 00000000..dadf1b34 --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_set.h @@ -0,0 +1,100 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include + +namespace phat { + class vector_set { + + protected: + typedef std::set< index > internal_column; + std::vector< dimension > dims; + std::vector< internal_column > matrix; + + public: + // overall number of cells in boundary_matrix + index _get_num_cols() const { + return (index)matrix.size(); + } + void _set_num_cols( index nr_of_columns ) { + dims.resize( nr_of_columns ); + matrix.resize( nr_of_columns ); + } + + // dimension of given index + dimension _get_dim( index idx ) const { + return dims[ idx ]; + } + void _set_dim( index idx, dimension dim ) { + dims[ idx ] = dim; + } + + // replaces(!) content of 'col' with boundary of given index + void _get_col( index idx, column& col ) const { + col.clear(); + col.reserve( matrix[idx].size() ); + std::copy (matrix[idx].begin(), matrix[idx].end(), std::back_inserter(col) ); + } + void _set_col( index idx, const column& col ) { + matrix[ idx ].clear(); + matrix[ idx ].insert( col.begin(), col.end() ); + } + + // true iff boundary of given idx is empty + bool _is_empty( index idx ) const { + return matrix[ idx ].empty(); + } + + // largest row index of given column idx (new name for lowestOne()) + index _get_max_index( index idx ) const { + return matrix[ idx ].empty() ? -1 : *matrix[ idx ].rbegin(); + } + + // removes the maximal index of a column + void _remove_max( index idx ) { + internal_column::iterator it = matrix[ idx ].end(); + it--; + matrix[ idx ].erase( it ); + } + + // clears given column + void _clear( index idx ) { + matrix[ idx ].clear(); + } + + // syncronizes all data structures (essential for openmp stuff) + void _sync() {} + + // adds column 'source' to column 'target' + void _add_to( index source, index target ) { + for( internal_column::iterator it = matrix[ source ].begin(); it != matrix[ source ].end(); it++ ) + _toggle( target, *it ); + } + + //// toggles given index pair + void _toggle( index col_idx, index row_idx ) { + internal_column& col = matrix[ col_idx ]; + std::pair< internal_column::iterator, bool > result = col.insert( row_idx ); + if( !result.second ) { + col.erase( result.first ); + } + } + }; +} diff --git a/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_vector.h b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_vector.h new file mode 100644 index 00000000..efe8de4d --- /dev/null +++ b/src/Bitmap_cubical_complex/include/gudhi/phat/representations/vector_vector.h @@ -0,0 +1,93 @@ +/* Copyright 2013 IST Austria + Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include "../helpers/misc.h" + +namespace phat { + class vector_vector { + + protected: + std::vector< dimension > dims; + std::vector< column > matrix; + + thread_local_storage< column > temp_column_buffer; + + public: + // overall number of cells in boundary_matrix + index _get_num_cols() const { + return (index)matrix.size(); + } + void _set_num_cols( index nr_of_columns ) { + dims.resize( nr_of_columns ); + matrix.resize( nr_of_columns ); + } + + // dimension of given index + dimension _get_dim( index idx ) const { + return dims[ idx ]; + } + void _set_dim( index idx, dimension dim ) { + dims[ idx ] = dim; + } + + // replaces(!) content of 'col' with boundary of given index + void _get_col( index idx, column& col ) const { + col = matrix[ idx ]; + } + void _set_col( index idx, const column& col ) { + matrix[ idx ] = col; + } + + // true iff boundary of given idx is empty + bool _is_empty( index idx ) const { + return matrix[ idx ].empty(); + } + + // largest row index of given column idx (new name for lowestOne()) + index _get_max_index( index idx ) const { + return matrix[ idx ].empty() ? -1 : matrix[ idx ].back(); + } + + // removes the maximal index of a column + void _remove_max( index idx ) { + matrix[ idx ].pop_back(); + } + + // clears given column + void _clear( index idx ) { + matrix[ idx ].clear(); + } + + // syncronizes all data structures (essential for openmp stuff) + void _sync() {} + + // adds column 'source' to column 'target' + void _add_to( index source, index target ) { + column& source_col = matrix[ source ]; + column& target_col = matrix[ target ]; + column& temp_col = temp_column_buffer(); + temp_col.clear(); + std::set_symmetric_difference( target_col.begin(), target_col.end(), + source_col.begin(), source_col.end(), + std::back_inserter( temp_col ) ); + target_col.swap( temp_col ); + } + }; +} diff --git a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp index 4235761f..968483a3 100644 --- a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp +++ b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp @@ -310,7 +310,7 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) { dimensions.push_back(3); Bitmap_cubical_complex< Bitmap_cubical_complex_base > increasing(dimensions, increasingFiltrationOfTopDimensionalCells); - for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) { + 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]); @@ -423,7 +423,7 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_2) { coboundaryElements.push_back(41); coboundaryElements.push_back(47); size_t number = 0; - for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) { + 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]); @@ -502,7 +502,7 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_3) { dim.push_back(1); dim.push_back(0); - for (size_t i = 0; i != increasing.size_of_bitmap(); ++i) { + for (size_t i = 0; i != increasing.size(); ++i) { BOOST_CHECK(increasing.get_dimension_of_a_cell(i) == dim[i]); } } @@ -726,7 +726,7 @@ BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) { boundaries.push_back( boundary14 ); boundaries.push_back( boundary15 ); - for (size_t i = 0; i != cmplx.size_of_bitmap(); ++i) { + 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]); @@ -826,7 +826,7 @@ BOOST_AUTO_TEST_CASE(coboudary_operator_2d_bitmap_with_periodic_bcond) { coboundaries.push_back( coboundary14 ); coboundaries.push_back( coboundary15 ); - for (size_t i = 0; i != cmplx.size_of_bitmap(); ++i) { + 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]); @@ -879,7 +879,7 @@ BOOST_AUTO_TEST_CASE(bitmap_2d_with_periodic_bcond_filtration) { filtration.push_back(3);//15 - for (size_t i = 0; i != cmplx.size_of_bitmap(); ++i) + for (size_t i = 0; i != cmplx.size(); ++i) { BOOST_CHECK( filtration[i] == cmplx.get_cell_data(i) ); } -- cgit v1.2.3