From c817c81d9d79fe77b0729216b7db5a7eab56a23c Mon Sep 17 00:00:00 2001 From: "ulrich.bauer@gmail.com" Date: Thu, 28 Feb 2013 19:31:04 +0000 Subject: initial checkin git-svn-id: https://phat.googlecode.com/svn/trunk@2 8e3bb3c2-eed4-f18f-5264-0b6c94e6926d --- src/phat.cpp | 173 +++++++++++++++++++++++++++++++++++++++++++++++++ src/self_test.cpp | 154 +++++++++++++++++++++++++++++++++++++++++++ src/simple_example.cpp | 127 ++++++++++++++++++++++++++++++++++++ 3 files changed, 454 insertions(+) create mode 100644 src/phat.cpp create mode 100644 src/self_test.cpp create mode 100644 src/simple_example.cpp (limited to 'src') 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 . */ + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + + +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 +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 . */ + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +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 . */ + +// 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 + +// main data structure (choice affects performance) +#include + +// algorithm (choice affects performance) +#include +#include +#include +#include + +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; +} -- cgit v1.2.3