From bf9bf990219174c93e86132fe088713c5760f93b Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 28 Jul 2024 19:41:12 +0200 Subject: [PATCH] Full_distance_matrix (unused) --- src/Ripser/include/gudhi/ripser.h | 31 ++++++++++++++++++++++++++++++- src/Ripser/utility/ripser.cc | 11 ++++++----- src/Ripser/utility/test_rips | 5 ++++- 3 files changed, 40 insertions(+), 7 deletions(-) diff --git a/src/Ripser/include/gudhi/ripser.h b/src/Ripser/include/gudhi/ripser.h index 1831f639eb..c2b1f4bdff 100644 --- a/src/Ripser/include/gudhi/ripser.h +++ b/src/Ripser/include/gudhi/ripser.h @@ -53,7 +53,7 @@ // - once cohomology is computed for some dim, assemble the relevant simplices (don't forget the essential ones) and reduce them homology-style. I think we have 2 choices: reduce the birth and check which columns we add (only possibility for essential classes, but do we really need to do cohomology first if we are going to do that?), or reduce the death and look at the column after reduction (Ripser(er) chose this). // * check out the persistence image branch // -// * allow non-0 filtration value on vertices, so we can handle all flag-type filtrations, not just plain Rips +// * allow non-0 filtration value on vertices, so we can handle all flag-type filtrations, not just plain Rips. cf scikit-tda // // * assert -> GUDHI_* // INDICATE_PROGRESS -> GUDHI_INDICATE_PROGRESS @@ -198,6 +198,34 @@ class Tag_dense {}; // Use directly, iterate on vertices class Tag_sparse {}; // Use directly, iterate on edges class Tag_other {}; // Could use directly as dense, but prefer converting it. +template +struct Full_distance_matrix { + public: + typedef Tag_dense Category; + typedef typename Params::vertex_t vertex_t; + typedef typename Params::value_t value_t; + std::vector distances; // TODO: private + private: + vertex_t n; + public: + //Full_distance_matrix(std::vector&& _distances) + // : distances(std::move(_distances)), n(std::sqrt(_distances.size())) {} + + template + Full_distance_matrix(const DistanceMatrix& mat) + : distances(static_cast(mat.size()) * mat.size()), n(mat.size()) { // vertex_t is meant for vertices. Using it for edges could be unsafe, so we cast to size_t. + for (vertex_t i = 0; i < size(); ++i) + for (vertex_t j = 0; j < size(); ++j) + distances[i * n + j] = mat(i, j); // If mat::operator() involves computations, should we try to take advantage of the symmetry? + } + + // Confusing: why is the order j,i significantly faster than i,j? + value_t operator()(const vertex_t j, const vertex_t i) const { + return distances[i * n + j]; + } + vertex_t size() const { return n; } +}; + enum Compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR }; template @@ -734,6 +762,7 @@ template st assert(k != -1); } value_t cofacet_diameter = get_diameter(simplex); + // The order of j and i matters for performance for (vertex_t i : vertices) cofacet_diameter = std::max(cofacet_diameter, dist(j, i)); simplex_t cofacet_index = idx_above + simplex_encoding(j--, k + 1) + idx_below; coefficient_t cofacet_coefficient = parent.get_coefficient(simplex); diff --git a/src/Ripser/utility/ripser.cc b/src/Ripser/utility/ripser.cc index d21c93bebb..6eb39a9864 100644 --- a/src/Ripser/utility/ripser.cc +++ b/src/Ripser/utility/ripser.cc @@ -51,7 +51,7 @@ struct Params1 { typedef float value_t; typedef int8_t dimension_t; // Does it need to be signed? Experimentally no. typedef int vertex_t; // Currently needs to be signed for Simplex_coboundary_enumerator::has_next. Reducing to int16_t helps perf a bit. - typedef unsigned __int128 simplex_t; // FIXME: don't try to use that on windows... + // typedef Gudhi::numbers::uint128_t simplex_t; // We could introduce a smaller edge_t, but it is not trivial to separate and probably not worth it typedef uint16_t coefficient_storage_t; // For the table of multiplicative inverses typedef uint_least32_t coefficient_t; // We need x * y % z to work, but promotion from uint16_t to int is not enough @@ -69,10 +69,11 @@ struct Ripser_all { typedef typename Params::value_t value_t; typedef typename Params::dimension_t dimension_t; typedef typename Params::vertex_t vertex_t; - typedef typename Params::simplex_t simplex_t; + //typedef typename Params::simplex_t simplex_t; typedef typename Params::coefficient_t coefficient_t; - static constexpr bool use_coefficients = Params::use_coefficients; + //static constexpr bool use_coefficients = Params::use_coefficients; + typedef Gudhi::ripser::Full_distance_matrix Full_distance_matrix; typedef Compressed_distance_matrix Compressed_lower_distance_matrix; typedef Compressed_distance_matrix Compressed_upper_distance_matrix; typedef Gudhi::ripser::Sparse_distance_matrix Sparse_distance_matrix; @@ -370,7 +371,7 @@ struct Ripser_all { size_t num_edges = 0; value_t enclosing_radius = std::numeric_limits::infinity(); - if (threshold == std::numeric_limits::max()) { + if (threshold >= std::numeric_limits::max()) { for (vertex_t i = 0; i < dist.size(); ++i) { value_t r_i = -std::numeric_limits::infinity(); for (vertex_t j = 0; j < dist.size(); ++j) r_i = std::max(r_i, dist(i, j)); @@ -386,7 +387,7 @@ struct Ripser_all { } std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl; - if (threshold == std::numeric_limits::max()) { + if (threshold >= std::numeric_limits::max()) { std::cout << "distance matrix with " << dist.size() << " points, using threshold at enclosing radius " << enclosing_radius << std::endl; diff --git a/src/Ripser/utility/test_rips b/src/Ripser/utility/test_rips index f43846127d..0dda95a369 100755 --- a/src/Ripser/utility/test_rips +++ b/src/Ripser/utility/test_rips @@ -1,7 +1,10 @@ #!/usr/bin/zsh #set -e #with -O3 -march=native and g++ 12 or 13, the result slightly differs (order, rounding). That's suspicious, but ok if the difference is just that. -g++-snapshot ripser.cc -O2 -DNDEBUG -Wall -g `gudhi-includes` `boost-includes` || exit -1 +#g++-snapshot ripser.cc -O2 -DNDEBUG -Wall -g `gudhi-includes` `boost-includes` || exit -1 +export PATH=/usr/lib/ccache:"$PATH" +g++-14 ripser.cc -c -O2 -DNDEBUG -Wall -g `gudhi-includes` `boost-includes` || exit -1 +g++-14 ripser.o -O2 -g || exit -1 #g++-snapshot ripser.cc -O2 -Wall -g `gudhi-includes` `boost-includes` || exit -1 # push time ./a.out --format point-cloud ~/repos/ripser/examples/o3_1024.txt >! /tmp/new1