From a877a13a7abd84fef2acb902d420826d04bffaf3 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 30 Jul 2016 19:11:29 +0200 Subject: distance matrix cleanup --- README.md | 2 +- ripser.cpp | 100 ++++++++++++++++----------------------- ripser.xcodeproj/project.pbxproj | 4 +- 3 files changed, 44 insertions(+), 62 deletions(-) diff --git a/README.md b/README.md index cfe2f2f..f8b7372 100644 --- a/README.md +++ b/README.md @@ -73,7 +73,7 @@ The following options are supported at the command line: - `--dim k`: compute persistent homology up to dimension *k* - `--threshold t`: compute Rips complexes up to diameter *t* - - `--modulus p`: compute homology with coefficients in the prime field Z/*p*Z (only available when build with the option `USE_COEFFICIENTS`) + - `--modulus p`: compute homology with coefficients in the prime field Z/*p*Z (only available when built with the option `USE_COEFFICIENTS`) ### Planned features diff --git a/ripser.cpp b/ripser.cpp index bd614de..b17616f 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -283,50 +283,28 @@ public: } }; -class distance_matrix { -public: - typedef value_t value_type; - std::vector> distances; - value_t operator()(const index_t a, const index_t b) const { return distances[a][b]; } - - size_t size() const { return distances.size(); } -}; - enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR }; -template class compressed_distance_matrix_adapter { +template class compressed_distance_matrix { public: typedef value_t value_type; - std::vector& distances; + std::vector distances; std::vector rows; - void init_distances() { distances.resize(size() * (size() - 1) / 2); } - void init_rows(); - compressed_distance_matrix_adapter(std::vector& _distances) - : distances(_distances), rows((1 + std::sqrt(1 + 8 * _distances.size())) / 2) { + compressed_distance_matrix(std::vector _distances) + : distances(std::move(_distances)), rows((1 + std::sqrt(1 + 8 * distances.size())) / 2) { assert(distances.size() == size() * (size() - 1) / 2); init_rows(); } - template - compressed_distance_matrix_adapter(std::vector& _distances, const DistanceMatrix& mat) - : distances(_distances), rows(mat.size()) { - init_distances(); - init_rows(); - - for (index_t i = 1; i < size(); ++i) { - for (index_t j = 0; j < i; ++j) { rows[i][j] = mat(i, j); } - } - } - - value_t operator()(const index_t i, const index_t j) const; + value_t operator()(const index_t i, const index_t j) const; size_t size() const { return rows.size(); } }; -template <> void compressed_distance_matrix_adapter::init_rows() { +template <> void compressed_distance_matrix::init_rows() { value_t* pointer = &distances[0]; for (index_t i = 1; i < size(); ++i) { rows[i] = pointer; @@ -334,7 +312,7 @@ template <> void compressed_distance_matrix_adapter::init_rows } } -template <> void compressed_distance_matrix_adapter::init_rows() { +template <> void compressed_distance_matrix::init_rows() { value_t* pointer = &distances[0] - 1; for (index_t i = 0; i < size() - 1; ++i) { rows[i] = pointer; @@ -343,7 +321,7 @@ template <> void compressed_distance_matrix_adapter::init_rows } template <> -value_t compressed_distance_matrix_adapter::operator()(const index_t i, +value_t compressed_distance_matrix::operator()(const index_t i, const index_t j) const { if (i < j) return rows[i][j]; @@ -354,7 +332,7 @@ value_t compressed_distance_matrix_adapter::operator()(const i } template <> -value_t compressed_distance_matrix_adapter::operator()(const index_t i, +value_t compressed_distance_matrix::operator()(const index_t i, const index_t j) const { if (i > j) return rows[i][j]; @@ -364,10 +342,10 @@ value_t compressed_distance_matrix_adapter::operator()(const i return 0; } -typedef compressed_distance_matrix_adapter - compressed_lower_distance_matrix_adapter; -typedef compressed_distance_matrix_adapter - compressed_upper_distance_matrix_adapter; +typedef compressed_distance_matrix + compressed_lower_distance_matrix; +typedef compressed_distance_matrix + compressed_upper_distance_matrix; template diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) { if (column.empty()) @@ -771,10 +749,8 @@ int main(int argc, char** argv) { exit(-1); } - std::vector diameters; - #ifdef FILE_FORMAT_DIPHA - + if (read(input_stream) != 8067171840) { std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; exit(-1); @@ -787,44 +763,50 @@ int main(int argc, char** argv) { index_t n = read(input_stream); - distance_matrix dist_full; - dist_full.distances = std::vector>(n, std::vector(n)); + std::vector distances; for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) dist_full.distances[i][j] = read(input_stream); - std::cout << "distance matrix with " << n << " points" << std::endl; + for (int j = 0; j < n; ++j) + if (i > j) distances.push_back(read(input_stream)); else read(input_stream); - compressed_lower_distance_matrix_adapter dist(diameters, dist_full); - std::cout << "distance matrix transformed to lower triangular form" << std::endl; -#endif - -#ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV - std::vector distances; - value_t value; - while ((input_stream >> value).ignore()) distances.push_back(value); - - compressed_upper_distance_matrix_adapter dist_upper(distances); - - index_t n = dist_upper.size(); - std::cout << "distance matrix with " << n << " points" << std::endl; + compressed_lower_distance_matrix dist(distances); + + std::cout << "distance matrix with " << n << " points" << std::endl; - compressed_lower_distance_matrix_adapter dist(diameters, dist_upper); - std::cout << "distance matrix transformed to lower triangular form" << std::endl; #endif #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV - std::vector& distances = diameters; + + std::vector distances; value_t value; while (input_stream >> value) { distances.push_back(value); input_stream.ignore(); } - compressed_lower_distance_matrix_adapter dist(distances); + compressed_lower_distance_matrix dist(distances); index_t n = dist.size(); std::cout << "distance matrix with " << n << " points" << std::endl; + +#endif + +#ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV + + std::vector distances; + value_t value; + while (input_stream >> value) { + distances.push_back(value); + input_stream.ignore(); + } + + compressed_upper_distance_matrix dist(distances); + + index_t n = dist.size(); + + std::cout << "distance matrix with " << n << " points" << std::endl; + #endif auto result = std::minmax_element(dist.distances.begin(), dist.distances.end()); diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index 08adcbd..a14431e 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -7,7 +7,7 @@ objects = { /* Begin PBXBuildFile section */ - 55E113311BD63CF3002B6F51 /* ripser.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 55E113301BD63CF3002B6F51 /* ripser.cpp */; settings = {ASSET_TAGS = (); }; }; + 55E113311BD63CF3002B6F51 /* ripser.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 55E113301BD63CF3002B6F51 /* ripser.cpp */; }; /* End PBXBuildFile section */ /* Begin PBXCopyFilesBuildPhase section */ @@ -209,7 +209,7 @@ buildSettings = { GCC_PREPROCESSOR_DEFINITIONS = ( "$(inherited)", - FILE_FORMAT_UPPER_TRIANGULAR_CSV, + FILE_FORMAT_LOWER_TRIANGULAR_CSV, ); PRODUCT_NAME = "$(TARGET_NAME)"; }; -- cgit v1.2.3