diff options
| -rw-r--r-- | src/Rips_complex/doc/Intro_rips_complex.h | 5 | ||||
| -rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 1 | ||||
| -rw-r--r-- | src/cython/cython/rips_complex.pyx | 50 | ||||
| -rw-r--r-- | src/cython/include/Rips_complex_interface.h | 39 | ||||
| -rwxr-xr-x | src/cython/test/test_rips_complex.py | 1 |
5 files changed, 60 insertions, 36 deletions
diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 712d3b6e..3db5287e 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -53,7 +53,10 @@ namespace rips_complex { * The number of simplices in the full Rips complex is exponential in the * number of vertices, it is thus usually restricted, by excluding all the * simplices with filtration value larger than some threshold, and keeping only - * the dim_max-skeleton. + * the dim_max-skeleton. It may also be a good idea to start by making the + * point set sparser, for instance with + * `Gudhi::subsampling::sparsify_point_set`, since small clusters of points + * have a disproportionate cost without affecting the persistence diagram much. * * In order to build this complex, the algorithm first builds the graph. * The filtration value of each edge is computed from a user-given distance diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index ee96d5a2..2c4c53a3 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1035,6 +1035,7 @@ class Simplex_tree { * The Simplex_tree must contain no simplex of dimension bigger than * 1 when calling the method. */ void expansion(int max_dim) { + if (max_dim <= 1) return; dimension_ = max_dim; for (Dictionary_it root_it = root_.members_.begin(); root_it != root_.members_.end(); ++root_it) { diff --git a/src/cython/cython/rips_complex.pyx b/src/cython/cython/rips_complex.pyx index 59c16bff..620130ca 100644 --- a/src/cython/cython/rips_complex.pyx +++ b/src/cython/cython/rips_complex.pyx @@ -33,7 +33,11 @@ __license__ = "GPL v3" cdef extern from "Rips_complex_interface.h" namespace "Gudhi": cdef cppclass Rips_complex_interface "Gudhi::rips_complex::Rips_complex_interface": - Rips_complex_interface(vector[vector[double]] values, double threshold, bool euclidean) + Rips_complex_interface() + void init_points(vector[vector[double]] values, double threshold) + void init_matrix(vector[vector[double]] values, double threshold) + void init_points_sparse(vector[vector[double]] values, double threshold, double sparse) + void init_matrix_sparse(vector[vector[double]] values, double threshold, double sparse) void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, int dim_max) # RipsComplex python interface @@ -44,14 +48,14 @@ cdef class RipsComplex: function, or a distance matrix. """ - cdef Rips_complex_interface * thisptr + cdef Rips_complex_interface thisref # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=None, distance_matrix=None, max_edge_length=float('inf')): + def __init__(self, points=None, distance_matrix=None, max_edge_length=float('inf'), sparse=None): """RipsComplex constructor. :param max_edge_length: Rips value. - :type max_edge_length: int + :type max_edge_length: float :param points: A list of points in d-Dimension. :type points: list of list of double @@ -61,27 +65,31 @@ cdef class RipsComplex: :param distance_matrix: A distance matrix (full square or lower triangular). :type points: list of list of double + + And in both cases + :param sparse: If this is not None, it switches to building a sparse Rips and represents the approximation parameter epsilon. + :type sparse: float """ # The real cython constructor - def __cinit__(self, points=None, distance_matrix=None, max_edge_length=float('inf')): - if distance_matrix is not None: - self.thisptr = new Rips_complex_interface(distance_matrix, max_edge_length, False) + def __cinit__(self, points=None, distance_matrix=None, max_edge_length=float('inf'), sparse=None): + if sparse is not None: + if distance_matrix is not None: + self.thisref.init_matrix_sparse(distance_matrix, max_edge_length, sparse) + else: + if points is None: + # Empty Rips construction + points=[] + self.thisref.init_points_sparse(points, max_edge_length, sparse) else: - if points is None: - # Empty Rips construction - points=[] - self.thisptr = new Rips_complex_interface(points, max_edge_length, True) - - - def __dealloc__(self): - if self.thisptr != NULL: - del self.thisptr + if distance_matrix is not None: + self.thisref.init_matrix(distance_matrix, max_edge_length) + else: + if points is None: + # Empty Rips construction + points=[] + self.thisref.init_points(points, max_edge_length) - def __is_defined(self): - """Returns true if RipsComplex pointer is not NULL. - """ - return self.thisptr != NULL def create_simplex_tree(self, max_dimension=1): """ @@ -92,5 +100,5 @@ cdef class RipsComplex: :rtype: SimplexTree """ simplex_tree = SimplexTree() - self.thisptr.create_simplex_tree(simplex_tree.thisptr, max_dimension) + self.thisref.create_simplex_tree(simplex_tree.thisptr, max_dimension) return simplex_tree diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h index 8b6c9c35..f6fc79a1 100644 --- a/src/cython/include/Rips_complex_interface.h +++ b/src/cython/include/Rips_complex_interface.h @@ -25,8 +25,11 @@ #include <gudhi/Simplex_tree.h> #include <gudhi/Rips_complex.h> +#include <gudhi/Sparse_rips_complex.h> #include <gudhi/distance_functions.h> +#include <boost/optional.hpp> + #include "Simplex_tree_interface.h" #include <iostream> @@ -43,28 +46,38 @@ class Rips_complex_interface { using Distance_matrix = std::vector<std::vector<Simplex_tree_interface<>::Filtration_value>>; public: - Rips_complex_interface(const std::vector<std::vector<double>>& values, double threshold, bool euclidean) { - if (euclidean) { - // Rips construction where values is a vector of points - rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold, - Gudhi::Euclidean_distance()); - } else { - // Rips construction where values is a distance matrix - rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold); - } + void init_points(const std::vector<std::vector<double>>& points, double threshold) { + rips_complex_.emplace(points, threshold, Gudhi::Euclidean_distance()); + } + void init_matrix(const std::vector<std::vector<double>>& matrix, double threshold) { + rips_complex_.emplace(matrix, threshold); } - ~Rips_complex_interface() { - delete rips_complex_; + void init_points_sparse(const std::vector<std::vector<double>>& points, double threshold, double epsilon) { + sparse_rips_complex_.emplace(points, Gudhi::Euclidean_distance(), epsilon); + threshold_ = threshold; + } + void init_matrix_sparse(const std::vector<std::vector<double>>& matrix, double threshold, double epsilon) { + sparse_rips_complex_.emplace(matrix, epsilon); + threshold_ = threshold; } void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, int dim_max) { - rips_complex_->create_complex(*simplex_tree, dim_max); + if (rips_complex_) + rips_complex_->create_complex(*simplex_tree, dim_max); + else { + sparse_rips_complex_->create_complex(*simplex_tree, dim_max); + // This pruning should be done much earlier! It isn't that useful for sparse Rips, but it would be inconsistent not to do it. + simplex_tree->prune_above_filtration(threshold_); + } simplex_tree->initialize_filtration(); } private: - Rips_complex<Simplex_tree_interface<>::Filtration_value>* rips_complex_; + // std::variant would work, but we don't require C++17 yet, and boost::variant is not super convenient. Anyway, storing a graph would make more sense. Or changing the interface completely so there is no such storage. + boost::optional<Rips_complex<Simplex_tree_interface<>::Filtration_value>> rips_complex_; + boost::optional<Sparse_rips_complex<Simplex_tree_interface<>::Filtration_value>> sparse_rips_complex_; + double threshold_; }; } // namespace rips_complex diff --git a/src/cython/test/test_rips_complex.py b/src/cython/test/test_rips_complex.py index c37b5400..2eb1c61c 100755 --- a/src/cython/test/test_rips_complex.py +++ b/src/cython/test/test_rips_complex.py @@ -30,7 +30,6 @@ __license__ = "GPL v3" def test_empty_rips(): rips_complex = RipsComplex() - assert rips_complex.__is_defined() == True def test_rips_from_points(): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] |