diff options
| author | pdlotko <pdlotko@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2016-02-09 13:26:47 +0000 |
|---|---|---|
| committer | pdlotko <pdlotko@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2016-02-09 13:26:47 +0000 |
| commit | 84399987baac2817e58bf9f5e18ded6aa6893b0f (patch) | |
| tree | 3a9c1f51a5ee1f4aa65e0b8061fc653ef7eb0bd9 /src | |
| parent | 3be6acc35255b52a60a254fa101aec5b11173b6d (diff) | |
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
Diffstat (limited to 'src')
24 files changed, 2532 insertions, 252 deletions
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 <em>elementary interval</em> is an interval of a form \f$ [n,n+1] \f$, or \f$[n,n]\f$, for \f$ n \in \mathcal{Z} \f$. The first one is called <em>non-degenerated</em>, while the second one is \a degenerated interval. A <em>boundary of a elementary
-*interval</em> is a chain \f$\partial [n,n+1] = [n+1,n+1]-[n,n] \f$ in case of non-degenerated elementary interval and \f$\partial [n,n] = 0 \f$ in case of degenerated elementary interval. An <em>elementary cube</em> \f$ C \f$ is a
+* An <em>elementary interval</em> is an interval of a form \f$ [n,n+1] \f$, or \f$[n,n]\f$, for \f$ n \in \mathcal{Z} \f$. The first one is called <em>non-degenerate</em>, while the second one is \a degenerate interval. A <em>boundary of a elementary
+*interval</em> is a chain \f$\partial [n,n+1] = [n+1,n+1]-[n,n] \f$ in case of non-degenerate elementary interval and \f$\partial [n,n] = 0 \f$ in case of degenerate elementary interval. An <em>elementary cube</em> \f$ C \f$ is a
-*product of elementary intervals, \f$C=I_1 \times \ldots \times I_n\f$. <em>Embedding dimension</em> of a cube is n, the number of elementary intervals (degenerated or not) in the product. A <em>dimension of a cube</em> \f$C=I_1 \times ... \times I_n\f$ is the
-*number of non degenerated elementary intervals in the product. A <em>boundary of a cube</em> \f$C=I_1 \times \ldots \times I_n\f$ is a chain obtained in the following way:
+*product of elementary intervals, \f$C=I_1 \times \ldots \times I_n\f$. <em>Embedding dimension</em> of a cube is n, the number of elementary intervals (degenerate or not) in the product. A <em>dimension of a cube</em> \f$C=I_1 \times ... \times I_n\f$ is the
+*number of non degenerate elementary intervals in the product. A <em>boundary of a cube</em> \f$C=I_1 \times \ldots \times I_n\f$ is a chain obtained in the following way:
*\f[\partial C = (\partial I_1 \times \ldots \times I_n) + (I_1 \times \partial I_2 \times \ldots \times I_n) + \ldots + (I_1 \times I_2 \times \ldots \times \partial I_n).\f]
*A <em>cubical complex</em> \f$\mathcal{K}\f$ is a collection of cubes closed under operation of taking boundary (i.e. boundary of every cube from the collection is in the collection). A cube \f$C\f$ in cubical complex \f$\mathcal{K}\f$ is <em>maximal</em> if it is not in
*a boundary of any other cube in \f$\mathcal{K}\f$. A \a support of a cube \f$C\f$ is the set in \f$\mathbb{R}^n\f$ occupied by \f$C\f$ (\f$n\f$ is the embedding dimension of \f$C\f$).
*
*Cubes may be equipped with a filtration values in which case we have filtered cubical complex. All the cubical complexes considered in this implementation are filtered cubical complexes (although, the range of a filtration may be a set of two elements).
*
-*For further details and theory of cubical complexes, please consult \cite kaczynski2004computational .
-*
-*as well as the following paper \cite peikert2012topological .
+*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 Binary files differindex f8207473..069c6eb2 100644 --- a/src/Bitmap_cubical_complex/doc/exampleBitmap.png +++ b/src/Bitmap_cubical_complex/doc/exampleBitmap.png 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<T>::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<T>(this) ); + is_before_in_filtration<T>(this) );
+
+ //we still need to deal here with a key_associated_to_simplex:
+ for ( size_t i = 0 ; i != simplex_associated_to_key.size() ; ++i )
+ {
+ this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i;
+ } } 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 <fstream> #include <algorithm> #include <iterator> -#include <limits> +#include <limits>
+#include <ctime> #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<unsigned>, 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 <typename K> - friend ostream& operator << ( ostream & os , const Bitmap_cubical_complex_base<K>& b ); + friend ostream& operator << ( ostream & os , const Bitmap_cubical_complex_base<K>& 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 <typename T>
+void Bitmap_cubical_complex_base<T>::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 <typename T> +void Bitmap_cubical_complex_base<T>::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 <typename T>
+std::pair< T ,T > Bitmap_cubical_complex_base<T>::min_max_filtration()
+{
+ std::pair< T ,T > min_max( std::numeric_limits<T>::max() , std::numeric_limits<T>::min() );
+ for ( size_t i = 0 ; i != this->data.size() ; ++i )
+ {
+ if ( this->data[i] < min_max.first )min_max.first = this->data[i];
+ if ( this->data[i] > min_max.second )min_max.second = this->data[i];
+ }
+ return min_max;
+}
template <typename K> @@ -422,7 +508,7 @@ void Bitmap_cubical_complex_base<T>::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<unsigned> 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 <cmath>
-#include "Bitmap_cubical_complex_base.h"
-
+#pragma once +#include <cmath> +#include "Bitmap_cubical_complex_base.h" + using namespace std; namespace Gudhi { namespace Cubical_complex -{
-
-//in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic boundary conditions). But, we set up the iterators and the procedures
-//to compute boundary and coboundary in the way that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all the cells on the
-//right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base
-
-template <typename T>
-class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T>
-{
-public:
- //constructors that take an extra parameter:
- Bitmap_cubical_complex_periodic_boundary_conditions_base(){};
- Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> 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<unsigned> dimensions , std::vector<T> 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 <typename T> +class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T> +{ +public: + //constructors that take an extra parameter: + Bitmap_cubical_complex_periodic_boundary_conditions_base(){}; + Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> 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<unsigned> dimensions , std::vector<T> 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<unsigned>& sizes ) - {
+ { unsigned multiplier = 1; for ( size_t i = 0 ; i != sizes.size() ; ++i ) { this->sizes.push_back(sizes[i]); this->multipliers.push_back(multiplier); -
- if ( directions_in_which_periodic_b_cond_are_to_be_imposed[i] )
- {
- multiplier *= 2*sizes[i];
- }
- else
- {
- multiplier *= 2*sizes[i]+1;
+ + if ( directions_in_which_periodic_b_cond_are_to_be_imposed[i] ) + { + multiplier *= 2*sizes[i]; + } + else + { + multiplier *= 2*sizes[i]+1; } } //std::reverse( this->sizes.begin() , this->sizes.end() ); this->data = std::vector<T>(multiplier,std::numeric_limits<T>::max()); this->total_number_of_cells = multiplier; - }
- Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> sizes );
- Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells );
- void construct_complex_based_on_top_dimensional_cells( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed );
-};
-
-template <typename T>
-void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells( std::vector<unsigned> dimensions , std::vector<T> 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<T>::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 <typename T>
-Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> 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 <typename T>
-Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseus_style_file )
-{
- //for Perseus style files:
- bool dbg = false;
-
+ } + Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> sizes ); + Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells ); + void construct_complex_based_on_top_dimensional_cells( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); +}; + +template <typename T> +void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells( std::vector<unsigned> dimensions , std::vector<T> 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<T>::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 <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> 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 <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseus_style_file ) +{ + + //for Perseus style files: + bool dbg = false; ifstream inFiltration; inFiltration.open( perseus_style_file ); unsigned dimensionOfData; - inFiltration >> dimensionOfData;
-
- this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>( dimensionOfData , false );
-
- std::vector<unsigned> sizes;
+ inFiltration >> dimensionOfData; + + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>( dimensionOfData , false ); + + std::vector<unsigned> sizes; sizes.reserve( dimensionOfData ); for ( size_t i = 0 ; i != dimensionOfData ; ++i ) { int size_in_this_dimension; - inFiltration >> size_in_this_dimension;
- if ( size_in_this_dimension < 0 )
- {
- this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] = true;
+ 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<T>::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<T>::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<T> data(slice,slice+w*h*d);
- delete[] slice;
- std::vector< unsigned > sizes;
- sizes.push_back(w);
- sizes.push_back(w);
- sizes.push_back(w);
-
- std::vector< bool > directions;
- directions.push_back( true );
- directions.push_back( true );
- directions.push_back( true );
- Bitmap_cubical_complex_periodic_boundary_conditions_base<T> b( sizes, data, directions );
- *this = b;
- */
-}
-
-template <typename T>
-Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> sizes )
-{
- this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>( sizes.size() , false );
- this->set_up_containers( sizes );
-}
-
-template <typename T>
-Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells )
-{
- std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>( dimensions.size() , false );
- this->construct_complex_based_on_top_dimensional_cells( dimensions , topDimensionalCells , directions_in_which_periodic_b_cond_are_to_be_imposed );
-}
-
-
-
-
-
-template <typename T>
-Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> 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 <typename T>
-std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell( size_t cell )const
-{
- bool dbg = false;
- if ( dbg ){cerr << "Computations of boundary of a cell : " << cell << endl;}
-
+ this->impose_lower_star_filtration(); + +/* + char* filename = (char*)perseus_style_file; + //char* filename = "combustionWithPeriodicBoundaryConditions/v0/tV0_000000.float"; + ifstream file( filename , ios::binary | ios::ate ); + unsigned realSizeOfFile = file.tellg(); + file.close(); + realSizeOfFile = realSizeOfFile/sizeof(T); + + unsigned w, h, d; + + w = h = d = ceil(pow( realSizeOfFile , (double)(1/(double)3) )); + + T* slice = new T[w*h*d]; + if (slice == NULL) + { + cerr << "Allocation error, cannot allocate " << w*h*d*sizeof(T) << " bytes to store the data from the file. The program will now terminate \n"; + exit(EXIT_FAILURE); + } + + FILE* fp; + if ((fp=fopen( filename, "rb" )) == NULL ) + { + cerr << "Cannot open the file: " << filename << ". The program will now terminate \n"; + exit(1); + } + + clock_t read_begin = clock(); + fread( slice,4,w*h*d,fp ); + fclose(fp); + cerr << "Time of reading the file : " << double(clock() - read_begin) / CLOCKS_PER_SEC << endl; + + + clock_t begin_creation_bitap = clock(); + std::vector<T> data(slice,slice+w*h*d); + delete[] slice; + std::vector< unsigned > sizes; + sizes.push_back(w); + sizes.push_back(w); + sizes.push_back(w); + + this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true ); + this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true ); + this->directions_in_which_periodic_b_cond_are_to_be_imposed.push_back( true ); + this->set_up_containers( sizes ); + + size_t i = 0; + for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator it = this->top_dimensional_cells_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 <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> sizes ) +{ + this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>( sizes.size() , false ); + this->set_up_containers( sizes ); +} + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells ) +{ + std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>( dimensions.size() , false ); + this->construct_complex_based_on_top_dimensional_cells( dimensions , topDimensionalCells , directions_in_which_periodic_b_cond_are_to_be_imposed ); +} + + + + + +template <typename T> +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> 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 <typename T> +std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell( size_t cell )const +{ + bool dbg = false; + if ( dbg ){cerr << "Computations of boundary of a cell : " << cell << endl;} + std::vector< size_t > boundary_elements; size_t cell1 = cell; for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) { - unsigned position = cell1/this->multipliers[i-1];
- //this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction.
+ 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 <typename T> std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell( size_t cell )const { @@ -254,12 +267,12 @@ std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T size_t cell1 = cell; for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) { - unsigned position = cell1/this->multipliers[i-1];
- //if the cell has zero length in this direction, then it will have cbd in this direction.
+ 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_base<T if ( (counter[i-1] != 2*this->sizes[i-1]) && (cell + this->multipliers[i-1] < this->data.size()) ) { coboundary_elements.push_back( cell + this->multipliers[i-1] ); - }
- }
- else
- {
- //we want to have periodic boundary conditions in this direction
- if ( counter[i-1] != 0 )
- {
- coboundary_elements.push_back( cell - this->multipliers[i-1] );
- coboundary_elements.push_back( cell + this->multipliers[i-1] );
- }
- else
- {
- //in this case counter[i-1] == 0.
- coboundary_elements.push_back( cell + this->multipliers[i-1] );
- coboundary_elements.push_back( cell + (2*this->sizes[ i-1 ]-1)*this->multipliers[i-1] );
- }
+ } } - }
+ 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 <http://www.gnu.org/licenses/>. + */ +
+#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 <typename K>
+void writeBettiNumbersAndPersistenceIntervalsToFile( char* prefix , std::pair< std::vector<std::vector< K > > , std::vector< std::vector< std::pair<K,K> > > > 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 <typename T , typename K>
+class Compute_persistence_with_phat
+{
+public:
+ Compute_persistence_with_phat( T* data_structure_ );
+ std::pair< std::vector< std::vector<K> > , std::vector< std::vector< std::pair<K,K> > > > 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 <typename T , typename K>
+void Compute_persistence_with_phat<T,K>::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 <typename T , typename K>
+phat::persistence_pairs Compute_persistence_with_phat<T,K>::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 <typename T , typename K>
+phat::persistence_pairs Compute_persistence_with_phat<T,K>::compute_persistence_pairs_twist_reduction()
+{
+ phat::persistence_pairs pairs;
+ phat::compute_persistence_pairs< phat::twist_reduction >( pairs, this->boundary_matrix );
+ return pairs;
+}
+
+template <typename T , typename K>
+phat::persistence_pairs Compute_persistence_with_phat<T,K>::compute_persistence_pairs_standard_reduction()
+{
+ phat::persistence_pairs pairs;
+ phat::compute_persistence_pairs< phat::standard_reduction >( pairs, this->boundary_matrix );
+ return pairs;
+}
+
+//template <typename T , typename K>
+//phat::persistence_pairs Compute_persistence_with_phat<T,K>::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 <typename T , typename K>
+Compute_persistence_with_phat<T,K>::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 <typename T , typename K>
+std::pair< std::vector< std::vector<K> > , std::vector< std::vector< std::pair<K,K> > > > Compute_persistence_with_phat<T,K>::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<bool> 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<K,K> > > 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<K> > birthTimesOfInfinitePersistnceClasses(this->data_structure->dimension()+1 );
+ for ( size_t i = 0 ; i != this->data_structure->dimension()+1 ; ++i )
+ {
+ std::vector<K> 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 <http://www.gnu.org/licenses/>. */
+
+#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<index> 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<index> 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<index>& 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 <http://www.gnu.org/licenses/>. */
+
+#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 <http://www.gnu.org/licenses/>. */
+
+#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 <http://www.gnu.org/licenses/>. */
+
+#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 <http://www.gnu.org/licenses/>. */
+
+#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>, 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 <http://www.gnu.org/licenses/>. */ + +#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 <http://www.gnu.org/licenses/>. */
+
+#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 <http://www.gnu.org/licenses/>. */
+
+#pragma once
+
+// STL includes
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <set>
+#include <list>
+#include <map>
+#include <algorithm>
+#include <queue>
+#include <cassert>
+#include <sstream>
+#include <algorithm>
+#include <iomanip>
+#include <cmath>
+#include <cstdlib>
+
+// 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 <stdint.h>
+#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 <omp.h>
+#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 <time.h>
+ #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 <http://www.gnu.org/licenses/>. */
+
+#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 <http://www.gnu.org/licenses/>. */
+
+#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 <http://www.gnu.org/licenses/>. */
+
+#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<vector<index>> 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 <http://www.gnu.org/licenses/>. */
+
+#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_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 <http://www.gnu.org/licenses/>. */
+
+#pragma once
+
+#include <phat/helpers/misc.h>
+#include <phat/representations/abstract_pivot_column.h>
+
+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 <http://www.gnu.org/licenses/>. */
+
+#pragma once
+
+#include <phat/helpers/misc.h>
+#include <phat/representations/abstract_pivot_column.h>
+
+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 <http://www.gnu.org/licenses/>. */ + +#pragma once + +#include <phat/helpers/misc.h> + +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 <http://www.gnu.org/licenses/>. */ + +#pragma once + +#include <phat/helpers/misc.h> + +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 <http://www.gnu.org/licenses/>. */
+
+#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<double> > 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) ); } |
