diff options
Diffstat (limited to 'ripser.cpp')
| -rw-r--r-- | ripser.cpp | 69 |
1 files changed, 39 insertions, 30 deletions
@@ -10,6 +10,7 @@ //#define USE_GOOGLE_HASHMAP +#include <algorithm> #include <iostream> #include <fstream> #include <cassert> @@ -362,11 +363,10 @@ template <typename Heap> diameter_entry_t pop_pivot(Heap& column, coefficient_t column.pop(); if (coefficient == 0) { - if (column.empty()) { + if (column.empty()) return diameter_entry_t(-1); - } else { + else pivot = column.top(); - } } } while (!column.empty() && get_index(column.top()) == get_index(pivot)); if (get_index(pivot) != -1) { set_coefficient(pivot, coefficient); } @@ -401,6 +401,11 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce columns_to_reduce.clear(); +#ifdef INDICATE_PROGRESS + std::cout << "\033[K" + << "assembling " << num_simplices << " columns" << std::flush << "\r"; +#endif + for (index_t index = 0; index < num_simplices; ++index) { if (pivot_column_index.find(index) == pivot_column_index.end()) { value_t diameter = comp.diameter(index); @@ -408,8 +413,16 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce } } +#ifdef INDICATE_PROGRESS + std::cout << "\033[K" + << "sorting " << num_simplices << " columns" << std::flush << "\r"; +#endif + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), greater_diameter_or_smaller_index<diameter_index_t>()); +#ifdef INDICATE_PROGRESS + std::cout << "\033[K"; +#endif } template <typename ValueType> class compressed_sparse_matrix { @@ -496,9 +509,10 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, value_t diameter = get_diameter(column_to_reduce); #ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " - << diameter << ")" << std::flush << "\r"; + if ((i + 1) % 1000 == 0) + std::cout << "\033[K" + << "reducing column " << i + 1 << "/" << columns_to_reduce.size() + << " (diameter " << diameter << ")" << std::flush << "\r"; #endif index_t j = i; @@ -666,6 +680,12 @@ bool is_prime(const coefficient_t n) { return true; } +template <typename T> T read(std::ifstream& s) { + T result; + s.read(reinterpret_cast<char*>(&result), sizeof(T)); + return result; // on little endian: boost::endian::little_to_native(result); +} + void print_usage_and_exit(int exit_code) { std::cerr << "Usage: " << "ripser " @@ -674,7 +694,7 @@ void print_usage_and_exit(int exit_code) { std::cerr << "Options:" << std::endl; std::cerr << std::endl; std::cerr << " --help print this screen" << std::endl; - std::cerr << " --top_dim <k> compute persistent homology up to dimension <k>" << std::endl; + std::cerr << " --dim <k> compute persistent homology up to dimension <k>" << std::endl; std::cerr << " --threshold <t> compute Rips complexes up to diameter <t>" << std::endl; #ifdef USE_COEFFICIENTS std::cerr << " --modulus <p> compute homology with coefficients in the prime field Z/<p>Z" @@ -703,7 +723,7 @@ int main(int argc, char** argv) { const std::string arg(argv[i]); if (arg == "--help") { print_usage_and_exit(0); - } else if (arg == "--top_dim") { + } else if (arg == "--dim") { std::string parameter = std::string(argv[++i]); size_t next_pos; dim_max = std::stol(parameter, &next_pos); @@ -735,33 +755,24 @@ int main(int argc, char** argv) { std::vector<value_t> diameters; #ifdef FILE_FORMAT_DIPHA - int64_t magic_number; - input_stream.read(reinterpret_cast<char*>(&magic_number), sizeof(int64_t)); - if (magic_number != 8067171840) { + + if (read<int64_t>(input_stream) != 8067171840) { std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; exit(-1); } - int64_t file_type; - input_stream.read(reinterpret_cast<char*>(&file_type), sizeof(int64_t)); - if (file_type != 7) { + if (read<int64_t>(input_stream) != 7) { std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl; exit(-1); } - int64_t n; - input_stream.read(reinterpret_cast<char*>(&n), sizeof(int64_t)); + index_t n = read<int64_t>(input_stream); distance_matrix dist_full; dist_full.distances = std::vector<std::vector<value_t>>(n, std::vector<value_t>(n)); - for (int i = 0; i < n; ++i) { - for (int j = 0; j < n; ++j) { - double val; - input_stream.read(reinterpret_cast<char*>(&val), sizeof(double)); - dist_full.distances[i][j] = val; - } - } + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) dist_full.distances[i][j] = read<double>(input_stream); std::cout << "distance matrix with " << n << " points" << std::endl; compressed_lower_distance_matrix_adapter dist(diameters, dist_full); @@ -770,9 +781,8 @@ int main(int argc, char** argv) { #ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV std::vector<value_t> distances; - std::string value_string; - while (std::getline(input_stream, value_string, ',')) - distances.push_back(std::stod(value_string)); + value_t value; + while ((input_stream >> value).ignore()) distances.push_back(value); compressed_upper_distance_matrix_adapter dist_upper(distances); @@ -785,9 +795,8 @@ int main(int argc, char** argv) { #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV std::vector<value_t>& distances = diameters; - std::string value_string; - while (std::getline(input_stream, value_string, ',')) - distances.push_back(std::stod(value_string)); + value_t value; + while ((input_stream >> value).ignore()) distances.push_back(value); compressed_lower_distance_matrix_adapter dist(distances); @@ -799,7 +808,7 @@ int main(int argc, char** argv) { auto result = std::minmax_element(dist.distances.begin(), dist.distances.end()); std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl; - assert(dim_max <= n - 2); + dim_max = std::min(dim_max, n - 2); binomial_coeff_table binomial_coeff(n, dim_max + 2); std::vector<coefficient_t> multiplicative_inverse(multiplicative_inverse_vector(modulus)); |