initial checkin

git-svn-id: https://phat.googlecode.com/svn/trunk@2 8e3bb3c2-eed4-f18f-5264-0b6c94e6926d

3 files changed, 454 insertions, 0 deletions

diff --git a/src/phat.cpp b/src/phat.cpp
new file mode 100644
index 0000000..895d52d
--- /dev/null
+++ b/src/phat.cpp
@@ -0,0 +1,173 @@
+/*  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/>. */

+

+#include <phat/compute_persistence_pairs.h>

+

+#include <phat/representations/vector_vector.h>

+#include <phat/representations/vector_set.h>

+#include <phat/representations/sparse_pivot_column.h>

+#include <phat/representations/full_pivot_column.h> 

+

+#include <phat/algorithms/twist_reduction.h>

+#include <phat/algorithms/standard_reduction.h> 

+#include <phat/algorithms/row_reduction.h>

+#include <phat/algorithms/chunk_reduction.h>

+

+

+enum Representation_type  {VEC_VEC, VEC_SET, SPARSE_PIVOT, FULL_PIVOT};

+enum Algorithm_type  {STANDARD, TWIST, ROW, CHUNK };

+

+void print_help() {

+    std::cerr << "Usage: " << "phat " << "[options] input_filename output_filename" << std::endl;

+    std::cerr << std::endl;

+    std::cerr << "Options:" << std::endl;

+    std::cerr << std::endl;

+    std::cerr << "--ascii   --  use ascii file format" << std::endl;

+    std::cerr << "--binary  --  use binary file format (default)" << std::endl;

+    std::cerr << "--help    --  prints this screen" << std::endl;

+    std::cerr << "--verbose --  verbose output" << std::endl;

+    std::cerr << "--dualize   --  use dualization approach" << std::endl;

+    std::cerr << "--vec-vec, --vec-set, --full-pivot, --sparse-pivot  --  selects a representation data structure for boundary matrices (default is '--sparse-pivot')" << std::endl;

+    std::cerr << "--standard, --twist, --chunk, --row  --  selects a reduction algorithm (default is '--twist')" << std::endl;

+}

+

+void print_help_and_exit() {

+    print_help();

+    exit( EXIT_FAILURE );

+}   

+

+void parse_command_line( int argc, char** argv, bool& use_binary, Representation_type& rep_type, Algorithm_type& reduction,

+                         std::string& input_filename, std::string& output_filename, bool& verbose, bool& dualize ) {

+    

+    if( argc < 3 ) print_help_and_exit();

+

+    input_filename = argv[ argc - 2 ];

+    output_filename = argv[ argc - 1 ];

+    

+    for( int idx = 1; idx < argc - 2; idx++ ) {

+        const std::string option = argv[ idx ];

+

+        if( option == "--ascii" ) use_binary = false;

+        else if( option == "--binary" ) use_binary = true;

+        else if( option == "--dualize" ) dualize = true;

+        else if( option == "--vec-vec" ) rep_type = VEC_VEC;

+        else if( option == "--vec-set" ) rep_type = VEC_SET;

+        else if( option == "--full-pivot" )  rep_type = FULL_PIVOT;

+        else if( option == "--sparse-pivot" ) rep_type = SPARSE_PIVOT;

+        else if( option == "--standard" ) reduction = STANDARD;

+        else if( option == "--twist" ) reduction = TWIST;

+        else if( option == "--row" ) reduction = ROW;

+        else if( option == "--chunk" ) reduction = CHUNK;

+        else if( option == "--verbose" ) verbose = true;

+        else if( option == "--help" ) print_help_and_exit();

+        else print_help_and_exit();

+    }

+}

+

+#define LOG(msg) if( verbose ) std::cout << msg << std::endl; 

+

+template<typename Representation, typename Algorithm>

+void generic_compute_pairing( std::string input_filename,

+                              std::string output_filename,

+	                          bool use_binary,

+	                          bool verbose,

+                              bool dualize ) {

+

+    phat::boundary_matrix< Representation > matrix;

+    bool read_successful;

+

+    double read_timer = omp_get_wtime();

+    if( use_binary ) {

+	    LOG( "Reading input file " << input_filename << " in binary mode" )

+        read_successful = matrix.load_binary( input_filename );

+    } else {

+        LOG( "Reading input file " << input_filename << " in ascii mode" )

+        read_successful = matrix.load_ascii( input_filename );

+    }

+    LOG( "Reading input file took " << omp_get_wtime() - read_timer <<"s" )

+

+    if( !read_successful ) {

+        std::cerr << "Error opening file " << input_filename << std::endl;

+        print_help_and_exit();

+    }

+

+    double pairs_timer = omp_get_wtime();

+    phat::persistence_pairs pairs;

+    LOG( "Computing persistence pairs ..." )

+    if( dualize )

+        phat::compute_persistence_pairs_dualized< Algorithm > ( pairs, matrix );

+    else 

+        phat::compute_persistence_pairs < Algorithm > ( pairs, matrix );

+    LOG( "Computing persistence pairs took " << omp_get_wtime() - pairs_timer <<"s" )

+

+    double write_timer = omp_get_wtime();

+    if( use_binary ) {

+        LOG( "Writing output file " << output_filename << " in binary mode ..." )

+        pairs.save_binary( output_filename );

+    } else {

+        LOG( "Writing output file " << output_filename << " in ascii mode ..." )

+        pairs.save_ascii( output_filename );

+    }

+    LOG( "Writing output file took " << omp_get_wtime() - write_timer <<"s" )

+    

+}

+

+#define CALL_GENERIC_CODE(rep,alg) generic_compute_pairing < rep, alg >( input_filename, output_filename, use_binary, verbose, dualize );

+

+int main( int argc, char** argv ) 

+{

+    bool use_binary = true; // interpret input as binary or ascii file

+    Representation_type rep_type = SPARSE_PIVOT; // representation class

+    Algorithm_type reduction = TWIST; // reduction algorithm

+    std::string input_filename; // name of file that contains the boundary matrix

+    std::string output_filename; // name of file that will contain the persistence pairs

+    bool verbose = false; // print timings / info

+    bool dualize = false; // toggle for dualization approach

+

+    parse_command_line( argc, argv, use_binary, rep_type, reduction, input_filename, output_filename, verbose, dualize );

+

+    switch( rep_type ) {

+    case VEC_VEC:       switch( reduction ) {

+                        case STANDARD: CALL_GENERIC_CODE(phat::vector_vector, phat::standard_reduction) break;

+                        case TWIST: CALL_GENERIC_CODE(phat::vector_vector, phat::twist_reduction) break;

+                        case ROW: CALL_GENERIC_CODE(phat::vector_vector, phat::row_reduction) break;

+                        case CHUNK: CALL_GENERIC_CODE(phat::vector_vector, phat::chunk_reduction) break;

+                        } break;

+

+    case VEC_SET:       switch( reduction ) {

+                        case STANDARD: CALL_GENERIC_CODE(phat::vector_set, phat::standard_reduction) break;

+                        case TWIST: CALL_GENERIC_CODE(phat::vector_set, phat::twist_reduction) break;

+                        case ROW: CALL_GENERIC_CODE(phat::vector_set, phat::row_reduction) break;

+                        case CHUNK: CALL_GENERIC_CODE(phat::vector_set, phat::chunk_reduction) break;

+                        } break;

+

+    case FULL_PIVOT:    switch( reduction ) {

+                        case STANDARD: CALL_GENERIC_CODE(phat::full_pivot_column, phat::standard_reduction) break;

+                        case TWIST: CALL_GENERIC_CODE(phat::full_pivot_column, phat::twist_reduction) break;

+                        case ROW: CALL_GENERIC_CODE(phat::full_pivot_column, phat::row_reduction) break;

+                        case CHUNK: CALL_GENERIC_CODE(phat::full_pivot_column, phat::chunk_reduction) break;

+                        } break;

+

+    case SPARSE_PIVOT:  switch( reduction ) {

+                        case STANDARD: CALL_GENERIC_CODE(phat::sparse_pivot_column, phat::standard_reduction) break;

+                        case TWIST: CALL_GENERIC_CODE(phat::sparse_pivot_column, phat::twist_reduction) break;

+                        case ROW: CALL_GENERIC_CODE(phat::sparse_pivot_column, phat::row_reduction) break;

+                        case CHUNK: CALL_GENERIC_CODE(phat::sparse_pivot_column, phat::chunk_reduction) break;

+                        } break;

+    }

+}

diff --git a/src/self_test.cpp b/src/self_test.cpp
new file mode 100644
index 0000000..4cda0ae
--- /dev/null
+++ b/src/self_test.cpp
@@ -0,0 +1,154 @@
+/*  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/>. */
+
+#include <phat/compute_persistence_pairs.h>
+
+#include <phat/representations/vector_vector.h>
+#include <phat/representations/vector_set.h>
+#include <phat/representations/sparse_pivot_column.h>
+#include <phat/representations/full_pivot_column.h>
+
+#include <phat/algorithms/twist_reduction.h>
+#include <phat/algorithms/standard_reduction.h>
+#include <phat/algorithms/row_reduction.h>
+#include <phat/algorithms/chunk_reduction.h>
+
+int main( int argc, char** argv )
+{
+    std::string test_data = argc > 1 ? argv[ 1 ] : "torus.bin";
+    
+    typedef phat::sparse_pivot_column Sparse;
+    typedef phat::full_pivot_column Full;
+    typedef phat::vector_vector Vec_vec;
+    typedef phat::vector_set Vec_set;
+    
+    std::cout << "Reading test data " << test_data << " in binary format ..." << std::endl;
+    phat::boundary_matrix< Full > boundary_matrix;
+    if( !boundary_matrix.load_binary( test_data ) ) {
+        std::cerr << "Error: test data " << test_data << " not found!" << std::endl;
+        return EXIT_FAILURE;
+    }
+    
+    bool error = false;
+    std::cout << "Comparing representations using Chunk algorithm ..." << std::endl;
+    {
+        std::cout << "Running Chunk - Sparse ..." << std::endl;
+        phat::persistence_pairs sparse_pairs;
+        phat::boundary_matrix< Sparse > sparse_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::chunk_reduction >( sparse_pairs, sparse_boundary_matrix );
+        
+        std::cout << "Running Chunk - Full ..." << std::endl;
+        phat::persistence_pairs full_pairs;
+        phat::boundary_matrix< Full > full_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::chunk_reduction >( full_pairs, full_boundary_matrix );
+        
+        std::cout << "Running Chunk - Vec_vec ..." << std::endl;
+        phat::persistence_pairs vec_vec_pairs;
+        phat::boundary_matrix< Vec_vec > vec_vec_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::chunk_reduction >( vec_vec_pairs, vec_vec_boundary_matrix );
+        
+        std::cout << "Running Chunk - Vec_set ..." << std::endl;
+        phat::persistence_pairs vec_set_pairs;
+        phat::boundary_matrix< Vec_set > vec_set_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::chunk_reduction >( vec_set_pairs, vec_set_boundary_matrix );
+        
+        if( sparse_pairs != full_pairs ) {
+            std::cerr << "Error: sparse and full differ!" << std::endl;
+            error = true;
+        }
+        if( full_pairs != vec_vec_pairs ) {
+            std::cerr << "Error: full and vec_vec differ!" << std::endl;
+            error = true;
+        }
+        if( vec_vec_pairs != vec_set_pairs ) {
+            std::cerr << "Error: vec_vec and vec_set differ!" << std::endl;
+            error = true;
+        }
+        if( vec_set_pairs != sparse_pairs ) {
+            std::cerr << "Error: vec_set and sparse differ!" << std::endl;
+            error = true;
+        }
+        
+        if( error ) return EXIT_FAILURE;
+        else std::cout << "All results are identical (as they should be)" << std::endl;
+    }
+    
+    std::cout << "Comparing algorithms using full_pivot_column representation ..." << std::endl;
+    {
+        std::cout << "Running Twist - Full ..." << std::endl;
+        phat::persistence_pairs twist_pairs;
+        phat::boundary_matrix< Full > twist_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::twist_reduction >( twist_pairs, twist_boundary_matrix );
+        
+        std::cout << "Running Standard - Full ..." << std::endl;
+        phat::persistence_pairs std_pairs;
+        phat::boundary_matrix< Full > std_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::standard_reduction >( std_pairs, std_boundary_matrix );
+        
+        std::cout << "Running Chunk - Full ..." << std::endl;
+        phat::persistence_pairs chunk_pairs;
+        phat::boundary_matrix< Full > chunk_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::chunk_reduction >( chunk_pairs, chunk_boundary_matrix );
+        
+        std::cout << "Running Row - Full ..." << std::endl;
+        phat::persistence_pairs row_pairs;
+        phat::boundary_matrix< Full > row_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::row_reduction >( row_pairs, row_boundary_matrix );
+        
+        if( twist_pairs != std_pairs ) {
+            std::cerr << "Error: twist and standard differ!" << std::endl;
+            error = true;
+        }
+        if( std_pairs != chunk_pairs ) {
+            std::cerr << "Error: standard and chunk differ!" << std::endl;
+            error = true;
+        }
+        if( chunk_pairs != row_pairs ) {
+            std::cerr << "Error: chunk and row differ!" << std::endl;
+            error = true;
+        }
+        if( row_pairs != twist_pairs ) {
+            std::cerr << "Error: row and twist differ!" << std::endl;
+            error = true;
+        }
+        
+        if( error ) return EXIT_FAILURE;
+        else std::cout << "All results are identical (as they should be)" << std::endl;
+    }
+    
+    std::cout << "Comparing primal and dual approach using Chunk - Full ..." << std::endl;
+    {
+        phat::persistence_pairs primal_pairs;
+        phat::boundary_matrix< Full > primal_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs< phat::chunk_reduction >( primal_pairs, primal_boundary_matrix );
+        
+        phat::persistence_pairs dual_pairs;
+        phat::boundary_matrix< Full > dual_boundary_matrix = boundary_matrix;
+        phat::compute_persistence_pairs_dualized< phat::chunk_reduction >( dual_pairs, dual_boundary_matrix );
+        
+        if( primal_pairs != dual_pairs ) {
+            std::cerr << "Error: primal and dual differ!" << std::endl;
+            error = true;
+        }
+        
+        if( error ) return EXIT_FAILURE;
+        else std::cout << "All results are identical (as they should be)" << std::endl;
+    }
+    
+    return EXIT_SUCCESS;
+}
diff --git a/src/simple_example.cpp b/src/simple_example.cpp
new file mode 100644
index 0000000..aff3ff9
--- /dev/null
+++ b/src/simple_example.cpp
@@ -0,0 +1,127 @@
+/*  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/>. */
+
+// This file contains a simple example that demonstrates the usage of the library interface
+
+// wrapper algorithm that computes the persistence pairs of a given boundary matrix using a specified algorithm
+#include <phat/compute_persistence_pairs.h>
+
+// main data structure (choice affects performance)
+#include <phat/representations/vector_vector.h>
+
+// algorithm (choice affects performance)
+#include <phat/algorithms/standard_reduction.h>
+#include <phat/algorithms/chunk_reduction.h>
+#include <phat/algorithms/row_reduction.h>
+#include <phat/algorithms/twist_reduction.h>
+
+int main( int argc, char** argv ) 
+{
+    std::cout << "We will build an ordered boundary matrix of this simplicial complex consisting of a single triangle: " << std::endl;
+    std::cout << std::endl;
+    std::cout << " 3" << std::endl;
+    std::cout << " |\\" << std::endl;
+    std::cout << " | \\" << std::endl;
+    std::cout << " |  \\" << std::endl;
+    std::cout << " |   \\ 4" << std::endl;
+    std::cout << "5|    \\" << std::endl;
+    std::cout << " |     \\" << std::endl;
+    std::cout << " |  6   \\" << std::endl;
+    std::cout << " |       \\" << std::endl;
+    std::cout << " |________\\" << std::endl;
+    std::cout << " 0    2    1" << std::endl;
+
+
+    // first define a boundary matrix with the chosen internal representation
+    phat::boundary_matrix< phat::vector_vector > boundary_matrix;
+
+    // set the number of columns (has to be 7 since we have 7 simplices)
+    boundary_matrix.set_num_cols( 7 );
+    
+    // set the dimension of the cell that a column represents:
+    boundary_matrix.set_dim( 0, 0 );
+    boundary_matrix.set_dim( 1, 0 );
+    boundary_matrix.set_dim( 2, 1 );
+    boundary_matrix.set_dim( 3, 0 );
+    boundary_matrix.set_dim( 4, 1 );
+    boundary_matrix.set_dim( 5, 1 );
+    boundary_matrix.set_dim( 6, 2 );
+
+    // set the respective columns -- the columns entries have to be sorted
+    std::vector< phat::index > temp_col;
+
+    boundary_matrix.set_col( 0, temp_col );
+
+    boundary_matrix.set_col( 1, temp_col );
+
+    temp_col.push_back( 0 );
+    temp_col.push_back( 1 );
+    boundary_matrix.set_col( 2, temp_col );
+    temp_col.clear();
+
+    boundary_matrix.set_col( 3, temp_col );
+
+    temp_col.push_back( 1 );
+    temp_col.push_back( 3 );
+    boundary_matrix.set_col( 4, temp_col );
+    temp_col.clear();
+
+    temp_col.push_back( 0 );
+    temp_col.push_back( 3 );
+    boundary_matrix.set_col( 5, temp_col );
+    temp_col.clear();
+
+    temp_col.push_back( 2 );
+    temp_col.push_back( 4 );
+    temp_col.push_back( 5 );
+    boundary_matrix.set_col( 6, temp_col );
+    temp_col.clear();
+
+    // print some information of the boundary matrix:
+    std::cout << std::endl;
+    std::cout << "The boundary matrix has " << boundary_matrix.get_num_cols() << " columns: " << std::endl;
+    for( phat::index col_idx = 0; col_idx < boundary_matrix.get_num_cols(); col_idx++ ) {
+        std::cout << "Colum " << col_idx << " represents a cell of dimension " << (int)boundary_matrix.get_dim( col_idx ) << ". ";
+        if( !boundary_matrix.is_empty( col_idx ) ) {
+            std::vector< phat::index > temp_col;
+            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;
+    }
+    std::cout << "Overall, the boundary matrix has " << boundary_matrix.get_num_entries() << " entries." << std::endl;  
+    
+
+    // define the object to hold the resulting persistence pairs
+    phat::persistence_pairs pairs;
+
+    // choose an algorithm (choice affects performance) and compute the persistence pair
+    // (modifies boundary_matrix)
+    phat::compute_persistence_pairs< phat::twist_reduction >( pairs, boundary_matrix );
+    
+    // sort the persistence pairs by birth index 
+    pairs.sort();
+
+    // print the pairs:
+    std::cout << std::endl;
+    std::cout << "There are " << pairs.get_num_pairs() << " persistence pairs: " << std::endl;
+    for( phat::index idx = 0; idx < pairs.get_num_pairs(); idx++ )
+        std::cout << "Birth: " << pairs.get_pair( idx ).first << ", Death: " << pairs.get_pair( idx ).second << std::endl;
+}