Full_distance_matrix (unused)

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 <class Params>
+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<value_t> distances; // TODO: private
+  private:
+    vertex_t n;
+  public:
+    //Full_distance_matrix(std::vector<value_t>&& _distances)
+    //  : distances(std::move(_distances)), n(std::sqrt(_distances.size())) {}
+
+    template <typename DistanceMatrix>
+      Full_distance_matrix(const DistanceMatrix& mat)
+      : distances(static_cast<std::size_t>(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 <class Params, Compressed_matrix_layout Layout>
@@ -734,6 +762,7 @@ template <typename DistanceMatrix, typename SimplexEncoding, typename Params> 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);

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<Compressed_lower_distance_matrix>::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<Params> Full_distance_matrix;
   typedef Compressed_distance_matrix<Params, LOWER_TRIANGULAR> Compressed_lower_distance_matrix;
   typedef Compressed_distance_matrix<Params, UPPER_TRIANGULAR> Compressed_upper_distance_matrix;
   typedef Gudhi::ripser::Sparse_distance_matrix<Params> Sparse_distance_matrix;
@@ -370,7 +371,7 @@ struct Ripser_all {
       size_t num_edges = 0;
 
       value_t enclosing_radius = std::numeric_limits<value_t>::infinity();
-      if (threshold == std::numeric_limits<value_t>::max()) {
+      if (threshold >= std::numeric_limits<value_t>::max()) {
         for (vertex_t i = 0; i < dist.size(); ++i) {
           value_t r_i = -std::numeric_limits<value_t>::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<value_t>::max()) {
+      if (threshold >= std::numeric_limits<value_t>::max()) {
         std::cout << "distance matrix with " << dist.size()
           << " points, using threshold at enclosing radius " << enclosing_radius
           << std::endl;

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