diff options
-rw-r--r-- | .github/next_release.md | 3 | ||||
-rw-r--r-- | biblio/how_to_cite_gudhi.bib.in | 178 | ||||
-rw-r--r-- | src/Cech_complex/doc/Intro_cech_complex.h | 2 | ||||
-rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex.h | 3 | ||||
-rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 23 | ||||
-rw-r--r-- | src/common/doc/main_page.md | 2 | ||||
-rw-r--r-- | src/python/doc/representations_sum.inc | 22 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 1 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 92 | ||||
-rw-r--r-- | src/python/include/Simplex_tree_interface.h | 26 | ||||
-rwxr-xr-x | src/python/test/test_simplex_tree.py | 365 |
11 files changed, 510 insertions, 207 deletions
diff --git a/.github/next_release.md b/.github/next_release.md index d5fcef1c..929a7ce6 100644 --- a/.github/next_release.md +++ b/.github/next_release.md @@ -9,6 +9,9 @@ Below is a list of changes made since GUDHI 3.6.0: - [Module](link) - ... +- [Simplex tree](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html) + - New functions to initialize from a matrix or insert batches of simplices of the same dimension. + - [Rips complex](https://gudhi.inria.fr/python/latest/rips_complex_user.html) - Construction now rejects positional arguments, you need to specify `points=X`. diff --git a/biblio/how_to_cite_gudhi.bib.in b/biblio/how_to_cite_gudhi.bib.in index 579dbf41..a1a971fc 100644 --- a/biblio/how_to_cite_gudhi.bib.in +++ b/biblio/how_to_cite_gudhi.bib.in @@ -7,6 +7,16 @@ , url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/" } +@incollection{gudhi:CubicalComplex +, author = "Pawel Dlotko" +, title = "Cubical complex" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__cubical__complex.html" +, year = @GUDHI_VERSION_YEAR@ +} + @incollection{gudhi:FilteredComplexes , author = "Cl\'ement Maria" , title = "Filtered Complexes" @@ -17,33 +27,33 @@ , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:PersistentCohomology -, author = "Cl\'ement Maria" -, title = "Persistent Cohomology" +@incollection{gudhi:ToplexMap +, author = "Fran{{\c{c}}ois Godi" +, title = "Toplex map" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__persistent__cohomology.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__toplex__map.html" , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:Contraction +@incollection{gudhi:SkeletonBlocker , author = "David Salinas" -, title = "Contraction" +, title = "Skeleton-Blocker" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__contr.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__skbl.html" , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:SkeletonBlocker +@incollection{gudhi:Contraction , author = "David Salinas" -, title = "Skeleton-Blocker" +, title = "Contraction" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__skbl.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__contr.html" , year = @GUDHI_VERSION_YEAR@ } @@ -57,13 +67,33 @@ , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:CubicalComplex -, author = "Pawel Dlotko" -, title = "Cubical complex" +@incollection{gudhi:CechComplex +, author = "Vincent Rouvreau, Hind Montassif" +, title = "\v{C}ech complex" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__cubical__complex.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__cech__complex.html" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:RipsComplex +, author = "Cl\'ement Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse" +, title = "Rips complex" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__rips__complex.html" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:Collapse +, author = "Siddharth Pritam, Marc Glisse" +, title = "Edge collapse" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__edge__collapse.html" , year = @GUDHI_VERSION_YEAR@ } @@ -77,23 +107,23 @@ , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:SubSampling -, author = "Cl\'ement Jamin and Siargey Kachanovich" -, title = "Subsampling" +@incollection{gudhi:CoverComplex +, author = "Mathieu Carri\`ere" +, title = "Cover complex" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__subsampling.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__cover__complex.html" , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:SpatialSearching -, author = "Cl\'ement Jamin" -, title = "Spatial searching" +@incollection{gudhi:CoxeterTriangulation +, author = "Siargey Kachanovich" +, title = "Coxeter triangulation" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__spatial__searching.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__cover__complex.html" , year = @GUDHI_VERSION_YEAR@ } @@ -107,13 +137,13 @@ , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:RipsComplex -, author = "Cl\'ement Maria and Pawel Dlotko and Vincent Rouvreau and Marc Glisse" -, title = "Rips complex" +@incollection{gudhi:PersistentCohomology +, author = "Cl\'ement Maria" +, title = "Persistent Cohomology" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__rips__complex.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__persistent__cohomology.html" , year = @GUDHI_VERSION_YEAR@ } @@ -127,42 +157,106 @@ , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:cython -, author = "Vincent Rouvreau" -, title = "Cython interface" +@incollection{gudhi:PersistenceRepresentations +, author = "Pawel Dlotko" +, title = "Persistence representations" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group___persistence__representations.html" , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:CoverComplex -, author = "Mathieu Carri\`ere" -, title = "Cover complex" +@incollection{gudhi:SubSampling +, author = "Cl\'ement Jamin, Siargey Kachanovich" +, title = "Subsampling" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__cover__complex.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__subsampling.html" , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:PersistenceRepresentations -, author = "Pawel Dlotko" -, title = "Persistence representations" +@incollection{gudhi:SpatialSearching +, author = "Cl\'ement Jamin" +, title = "Spatial searching" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group___persistence__representations.html" +, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__spatial__searching.html" , year = @GUDHI_VERSION_YEAR@ } -@incollection{gudhi:Collapse -, author = "Siddharth Pritam and Marc Glisse" -, title = "Edge collapse" +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Python specific gudhi modules +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +@incollection{gudhi:WeightedRipsComplex +, author = "Rapha\"el Tinarrage, Yuichi Ike, Masatoshi Takenouchi" +, title = "Weighted Rips Complex" , publisher = "{GUDHI Editorial Board}" , edition = "{@GUDHI_VERSION@}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "https://gudhi.inria.fr/doc/@GUDHI_VERSION@/group__edge__collapse.html" +, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/rips_complex_user.html#weighted-rips-complex" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:DTMRipsComplex +, author = "Yuichi Ike" +, title = "DTM Rips Complex" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/rips_complex_user.html#dtm-rips-complex" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:WassersteinDistance +, author = "Theo Lacombe, Marc Glisse" +, title = "Wasserstein distance" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/wasserstein_distance_user.html" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:PersistenceRepresentationsScikitlearnInterface +, author = "Mathieu Carri\`ere, Gard Spreemann, Wojciech Reise" +, title = "Persistence representations scikit-learn like interface" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/representations.html" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:Atol +, author = "Martin Royer" +, title = "Measure Vectorization for Automatic Topologically-Oriented Learning" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/representations.html#gudhi.representations.vector_methods.Atol" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:DistanceToMeasure +, author = "Marc Glisse" +, title = "Distance to measure" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/point_cloud.html#module-gudhi.point_cloud.knn" +, year = @GUDHI_VERSION_YEAR@ +} + +@incollection{gudhi:PersistenceBasedClustering +, author = "Marc Glisse" +, title = "persistence-based clustering" +, publisher = "{GUDHI Editorial Board}" +, edition = "{@GUDHI_VERSION@}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "https://gudhi.inria.fr/python/@GUDHI_VERSION@/clustering.html" , year = @GUDHI_VERSION_YEAR@ } diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h index 595fb64b..73093c07 100644 --- a/src/Cech_complex/doc/Intro_cech_complex.h +++ b/src/Cech_complex/doc/Intro_cech_complex.h @@ -17,7 +17,7 @@ namespace cech_complex { /** \defgroup cech_complex Čech complex * - * \author Vincent Rouvreau + * \author Vincent Rouvreau, Hind montassif * * @{ * diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 625f7c9c..dbdf5e93 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -1,11 +1,12 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Vincent Rouvreau + * Author(s): Vincent Rouvreau, Hind Montassif * * Copyright (C) 2018 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification + * - 2022/02 Hind Montassif : Replace MiniBall with Sphere_circumradius */ #ifndef CECH_COMPLEX_H_ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 9059219c..ef9f8428 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -24,6 +24,7 @@ #include <boost/iterator/transform_iterator.hpp> #include <boost/graph/adjacency_list.hpp> #include <boost/range/adaptor/reversed.hpp> +#include <boost/range/size.hpp> #include <boost/container/static_vector.hpp> #ifdef GUDHI_USE_TBB @@ -702,10 +703,10 @@ class Simplex_tree { return true; } - private: - /** \brief Inserts a simplex represented by a vector of vertex. - * @param[in] simplex vector of Vertex_handles, representing the vertices of the new simplex. The vector must be - * sorted by increasing vertex handle order. + protected: + /** \brief Inserts a simplex represented by a range of vertex. + * @param[in] simplex range of Vertex_handles, representing the vertices of the new simplex. The range must be + * sorted by increasing vertex handle order, and not empty. * @param[in] filtration the filtration value assigned to the new simplex. * @return If the new simplex is inserted successfully (i.e. it was not in the * simplicial complex yet) the bool is set to true and the Simplex_handle is the handle assigned @@ -717,12 +718,13 @@ class Simplex_tree { * null_simplex. * */ - std::pair<Simplex_handle, bool> insert_vertex_vector(const std::vector<Vertex_handle>& simplex, + template <class RandomVertexHandleRange = std::initializer_list<Vertex_handle>> + std::pair<Simplex_handle, bool> insert_simplex_raw(const RandomVertexHandleRange& simplex, Filtration_value filtration) { Siblings * curr_sib = &root_; std::pair<Simplex_handle, bool> res_insert; auto vi = simplex.begin(); - for (; vi != simplex.end() - 1; ++vi) { + for (; vi != std::prev(simplex.end()); ++vi) { GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex"); res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration)); if (!(has_children(res_insert.first))) { @@ -743,9 +745,10 @@ class Simplex_tree { return std::pair<Simplex_handle, bool>(null_simplex(), false); } // otherwise the insertion has succeeded - size is a size_type - if (static_cast<int>(simplex.size()) - 1 > dimension_) { + int dim = static_cast<int>(boost::size(simplex)) - 1; + if (dim > dimension_) { // Update dimension if needed - dimension_ = static_cast<int>(simplex.size()) - 1; + dimension_ = dim; } return res_insert; } @@ -786,7 +789,7 @@ class Simplex_tree { // Copy before sorting std::vector<Vertex_handle> copy(first, last); std::sort(std::begin(copy), std::end(copy)); - return insert_vertex_vector(copy, filtration); + return insert_simplex_raw(copy, filtration); } /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of @@ -1598,7 +1601,7 @@ class Simplex_tree { Simplex_tree st_copy = *this; // Add point for coning the simplicial complex - this->insert_simplex({maxvert}, -3); + this->insert_simplex_raw({maxvert}, -3); // For each simplex std::vector<Vertex_handle> vr; diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md index ce903405..9b7c2853 100644 --- a/src/common/doc/main_page.md +++ b/src/common/doc/main_page.md @@ -178,7 +178,7 @@ The set of all simplices is filtered by the radius of their minimal enclosing ball. </td> <td width="15%"> - <b>Author:</b> Vincent Rouvreau<br> + <b>Author:</b> Vincent Rouvreau, Hind Montassif<br> <b>Introduced in:</b> GUDHI 2.2.0<br> <b>Copyright:</b> MIT [(LGPL v3)](../../licensing/)<br> <b>Requires:</b> \ref cgal diff --git a/src/python/doc/representations_sum.inc b/src/python/doc/representations_sum.inc index 4298aea9..9515f044 100644 --- a/src/python/doc/representations_sum.inc +++ b/src/python/doc/representations_sum.inc @@ -1,14 +1,14 @@ .. table:: :widths: 30 40 30 - +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------+ - | .. figure:: | Vectorizations, distances and kernels that work on persistence | :Author: Mathieu Carrière, Martin Royer | - | img/sklearn-tda.png | diagrams, compatible with scikit-learn. | | - | | | :Since: GUDHI 3.1.0 | - | | | | - | | | :License: MIT | - | | | | - | | | :Requires: `Scikit-learn <installation.html#scikit-learn>`_ | - +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------+ - | * :doc:`representations` | - +------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------------+ + +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------------------+ + | .. figure:: | Vectorizations, distances and kernels that work on persistence | :Author: Mathieu Carrière, Martin Royer, Gard Spreemann, Wojciech Reise | + | img/sklearn-tda.png | diagrams, compatible with scikit-learn. | | + | | | :Since: GUDHI 3.1.0 | + | | | | + | | | :License: MIT | + | | | | + | | | :Requires: `Scikit-learn <installation.html#scikit-learn>`_ | + +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------------------+ + | * :doc:`representations` | + +------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 5642f82d..f86f1232 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -56,6 +56,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": int upper_bound_dimension() nogil bool find_simplex(vector[int] simplex) nogil bool insert(vector[int] simplex, double filtration) nogil + void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) nogil vector[pair[vector[int], double]] get_star(vector[int] simplex) nogil vector[pair[vector[int], double]] get_cofaces(vector[int] simplex, int dimension) nogil void expansion(int max_dim) nogil except + diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index ec18b708..18215d2f 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -8,14 +8,23 @@ # - YYYY/MM Author: Description of the modification from cython.operator import dereference, preincrement -from libc.stdint cimport intptr_t +from libc.stdint cimport intptr_t, int32_t, int64_t import numpy as np cimport gudhi.simplex_tree +cimport cython __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" +ctypedef fused some_int: + int32_t + int64_t + +ctypedef fused some_float: + float + double + cdef bool callback(vector[int] simplex, void *blocker_func): return (<object>blocker_func)(simplex) @@ -228,6 +237,87 @@ cdef class SimplexTree: """ return self.get_ptr().insert(simplex, <double>filtration) + @staticmethod + @cython.boundscheck(False) + def create_from_array(filtrations, double max_filtration=np.inf): + """Creates a new, empty complex and inserts vertices and edges. The vertices are numbered from 0 to n-1, and + the filtration values are encoded in the array, with the diagonal representing the vertices. It is the + caller's responsibility to ensure that this defines a filtration, which can be achieved with either:: + + filtrations[np.diag_indices_from(filtrations)] = filtrations.min(axis=1) + + or:: + + diag = filtrations.diagonal() + filtrations = np.fmax(np.fmax(filtrations, diag[:, None]), diag[None, :]) + + :param filtrations: the filtration values of the vertices and edges to insert. The matrix is assumed to be symmetric. + :type filtrations: numpy.ndarray of shape (n,n) + :param max_filtration: only insert vertices and edges with filtration values no larger than max_filtration + :type max_filtration: float + :returns: the new complex + :rtype: SimplexTree + """ + # TODO: document which half of the matrix is actually read? + filtrations = np.asanyarray(filtrations, dtype=float) + cdef double[:,:] F = filtrations + ret = SimplexTree() + cdef int n = F.shape[0] + assert n == F.shape[1], 'create_from_array() expects a square array' + with nogil: + ret.get_ptr().insert_matrix(&F[0,0], n, F.strides[0], F.strides[1], max_filtration) + return ret + + def insert_edges_from_coo_matrix(self, edges): + """Inserts edges given by a sparse matrix in `COOrdinate format + <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html>`_. + If an edge is repeated, the smallest filtration value is used. Missing entries are not inserted. + Diagonal entries are currently interpreted as vertices, although we do not guarantee this behavior + in the future, and this is only useful if you want to insert vertices with a smaller filtration value + than the smallest edge containing it, since vertices are implicitly inserted together with the edges. + + :param edges: the edges to insert and their filtration values. + :type edges: scipy.sparse.coo_matrix of shape (n,n) + + .. seealso:: :func:`insert_batch` + """ + # TODO: optimize this? + for edge in zip(edges.row, edges.col, edges.data): + self.get_ptr().insert((edge[0], edge[1]), edge[2]) + + @cython.boundscheck(False) + @cython.wraparound(False) + def insert_batch(self, some_int[:,:] vertex_array, some_float[:] filtrations): + """Inserts k-simplices given by a sparse array in a format similar + to `torch.sparse <https://pytorch.org/docs/stable/sparse.html>`_. + The n-th simplex has vertices `vertex_array[0,n]`, ..., + `vertex_array[k,n]` and filtration value `filtrations[n]`. + If a simplex is repeated, the smallest filtration value is used. + Simplices with a repeated vertex are currently interpreted as lower + dimensional simplices, but we do not guarantee this behavior in the + future. Any time a simplex is inserted, its faces are inserted as well + if needed to preserve a simplicial complex. + + :param vertex_array: the k-simplices to insert. + :type vertex_array: numpy.array of shape (k+1,n) + :param filtrations: the filtration values. + :type filtrations: numpy.array of shape (n,) + """ + # This may be slow if we end up inserting vertices in a bad order (flat_map). + # We could first insert the vertices from np.unique(vertex_array), or leave it to the caller. + cdef Py_ssize_t k = vertex_array.shape[0] + cdef Py_ssize_t n = vertex_array.shape[1] + assert filtrations.shape[0] == n, 'inconsistent sizes for vertex_array and filtrations' + cdef Py_ssize_t i + cdef Py_ssize_t j + cdef vector[int] v + with nogil: + for i in range(n): + for j in range(k): + v.push_back(vertex_array[j, i]) + self.get_ptr().insert(v, filtrations[i]) + v.clear() + def get_simplices(self): """This function returns a generator with simplices and their given filtration values. diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 3848c5ad..0317ea39 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -40,6 +40,8 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; using Extended_filtration_data = typename Base::Extended_filtration_data; using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator; + using Siblings = typename Base::Siblings; + using Node = typename Base::Node; typedef bool (*blocker_func_t)(Simplex simplex, void *user_data); public: @@ -62,6 +64,30 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return (result.second); } + void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) { + // We could delegate to insert_graph, but wrapping the matrix in a graph interface is too much work, + // and this is a bit more efficient. + auto& rm = this->root()->members_; + for(int i=0; i<n; ++i) { + char* p = reinterpret_cast<char*>(filtrations) + i * stride0; + double fv = *reinterpret_cast<double*>(p + i * stride1); + if(fv > max_filtration) continue; + auto sh = rm.emplace_hint(rm.end(), i, Node(this->root(), fv)); + Siblings* children = nullptr; + // Should we make a first pass to count the number of edges so we can reserve the right space? + for(int j=i+1; j<n; ++j) { + double fe = *reinterpret_cast<double*>(p + j * stride1); + if(fe > max_filtration) continue; + if(!children) { + children = new Siblings(this->root(), i); + sh->second.assign_children(children); + } + children->members().emplace_hint(children->members().end(), j, Node(children, fe)); + } + } + + } + // Do not interface this function, only used in alpha complex interface for complex creation bool insert_simplex(const Simplex& simplex, Filtration_value filtration = 0) { Insertion_result result = Base::insert_simplex(simplex, filtration); diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 54bafed5..2ccbfbf5 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -249,6 +249,7 @@ def test_make_filtration_non_decreasing(): assert st.filtration([3, 4]) == 2.0 assert st.filtration([4, 5]) == 2.0 + def test_extend_filtration(): # Inserted simplex: @@ -257,86 +258,87 @@ def test_extend_filtration(): # / \ / # o o # /2\ /3 - # o o - # 1 0 - - st = SimplexTree() - st.insert([0,2]) - st.insert([1,2]) - st.insert([0,3]) - st.insert([2,5]) - st.insert([3,4]) - st.insert([3,5]) - st.assign_filtration([0], 1.) - st.assign_filtration([1], 2.) - st.assign_filtration([2], 3.) - st.assign_filtration([3], 4.) - st.assign_filtration([4], 5.) - st.assign_filtration([5], 6.) - - assert list(st.get_filtration()) == [ - ([0, 2], 0.0), - ([1, 2], 0.0), - ([0, 3], 0.0), - ([3, 4], 0.0), - ([2, 5], 0.0), - ([3, 5], 0.0), - ([0], 1.0), - ([1], 2.0), - ([2], 3.0), - ([3], 4.0), - ([4], 5.0), - ([5], 6.0) + # o o + # 1 0 + + st = SimplexTree() + st.insert([0, 2]) + st.insert([1, 2]) + st.insert([0, 3]) + st.insert([2, 5]) + st.insert([3, 4]) + st.insert([3, 5]) + st.assign_filtration([0], 1.0) + st.assign_filtration([1], 2.0) + st.assign_filtration([2], 3.0) + st.assign_filtration([3], 4.0) + st.assign_filtration([4], 5.0) + st.assign_filtration([5], 6.0) + + assert list(st.get_filtration()) == [ + ([0, 2], 0.0), + ([1, 2], 0.0), + ([0, 3], 0.0), + ([3, 4], 0.0), + ([2, 5], 0.0), + ([3, 5], 0.0), + ([0], 1.0), + ([1], 2.0), + ([2], 3.0), + ([3], 4.0), + ([4], 5.0), + ([5], 6.0), ] - + st.extend_filtration() - - assert list(st.get_filtration()) == [ - ([6], -3.0), - ([0], -2.0), - ([1], -1.8), - ([2], -1.6), - ([0, 2], -1.6), - ([1, 2], -1.6), - ([3], -1.4), - ([0, 3], -1.4), - ([4], -1.2), - ([3, 4], -1.2), - ([5], -1.0), - ([2, 5], -1.0), - ([3, 5], -1.0), - ([5, 6], 1.0), - ([4, 6], 1.2), - ([3, 6], 1.4), + + assert list(st.get_filtration()) == [ + ([6], -3.0), + ([0], -2.0), + ([1], -1.8), + ([2], -1.6), + ([0, 2], -1.6), + ([1, 2], -1.6), + ([3], -1.4), + ([0, 3], -1.4), + ([4], -1.2), + ([3, 4], -1.2), + ([5], -1.0), + ([2, 5], -1.0), + ([3, 5], -1.0), + ([5, 6], 1.0), + ([4, 6], 1.2), + ([3, 6], 1.4), ([3, 4, 6], 1.4), - ([3, 5, 6], 1.4), - ([2, 6], 1.6), - ([2, 5, 6], 1.6), - ([1, 6], 1.8), - ([1, 2, 6], 1.8), - ([0, 6], 2.0), - ([0, 2, 6], 2.0), - ([0, 3, 6], 2.0) + ([3, 5, 6], 1.4), + ([2, 6], 1.6), + ([2, 5, 6], 1.6), + ([1, 6], 1.8), + ([1, 2, 6], 1.8), + ([0, 6], 2.0), + ([0, 2, 6], 2.0), + ([0, 3, 6], 2.0), ] - dgms = st.extended_persistence(min_persistence=-1.) + dgms = st.extended_persistence(min_persistence=-1.0) assert len(dgms) == 4 # Sort by (death-birth) descending - we are only interested in those with the longest life span for idx in range(4): - dgms[idx] = sorted(dgms[idx], key=lambda x:(-abs(x[1][0]-x[1][1]))) + dgms[idx] = sorted(dgms[idx], key=lambda x: (-abs(x[1][0] - x[1][1]))) + + assert dgms[0][0][1][0] == pytest.approx(2.0) + assert dgms[0][0][1][1] == pytest.approx(3.0) + assert dgms[1][0][1][0] == pytest.approx(5.0) + assert dgms[1][0][1][1] == pytest.approx(4.0) + assert dgms[2][0][1][0] == pytest.approx(1.0) + assert dgms[2][0][1][1] == pytest.approx(6.0) + assert dgms[3][0][1][0] == pytest.approx(6.0) + assert dgms[3][0][1][1] == pytest.approx(1.0) - assert dgms[0][0][1][0] == pytest.approx(2.) - assert dgms[0][0][1][1] == pytest.approx(3.) - assert dgms[1][0][1][0] == pytest.approx(5.) - assert dgms[1][0][1][1] == pytest.approx(4.) - assert dgms[2][0][1][0] == pytest.approx(1.) - assert dgms[2][0][1][1] == pytest.approx(6.) - assert dgms[3][0][1][0] == pytest.approx(6.) - assert dgms[3][0][1][1] == pytest.approx(1.) def test_simplices_iterator(): st = SimplexTree() - + assert st.insert([0, 1, 2], filtration=4.0) == True assert st.insert([2, 3, 4], filtration=2.0) == True @@ -346,9 +348,10 @@ def test_simplices_iterator(): print("filtration is: ", simplex[1]) assert st.filtration(simplex[0]) == simplex[1] + def test_collapse_edges(): st = SimplexTree() - + assert st.insert([0, 1], filtration=1.0) == True assert st.insert([1, 2], filtration=1.0) == True assert st.insert([2, 3], filtration=1.0) == True @@ -360,31 +363,33 @@ def test_collapse_edges(): st.collapse_edges() assert st.num_simplices() == 9 - assert st.find([0, 2]) == False # [1, 3] would be fine as well + assert st.find([0, 2]) == False # [1, 3] would be fine as well for simplex in st.get_skeleton(0): - assert simplex[1] == 1. + assert simplex[1] == 1.0 + def test_reset_filtration(): st = SimplexTree() - - assert st.insert([0, 1, 2], 3.) == True - assert st.insert([0, 3], 2.) == True - assert st.insert([3, 4, 5], 3.) == True - assert st.insert([0, 1, 6, 7], 4.) == True + + assert st.insert([0, 1, 2], 3.0) == True + assert st.insert([0, 3], 2.0) == True + assert st.insert([3, 4, 5], 3.0) == True + assert st.insert([0, 1, 6, 7], 4.0) == True # Guaranteed by construction for simplex in st.get_simplices(): - assert st.filtration(simplex[0]) >= 2. - + assert st.filtration(simplex[0]) >= 2.0 + # dimension until 5 even if simplex tree is of dimension 3 to test the limits for dimension in range(5, -1, -1): - st.reset_filtration(0., dimension) + st.reset_filtration(0.0, dimension) for simplex in st.get_skeleton(3): print(simplex) if len(simplex[0]) < (dimension) + 1: - assert st.filtration(simplex[0]) >= 2. + assert st.filtration(simplex[0]) >= 2.0 else: - assert st.filtration(simplex[0]) == 0. + assert st.filtration(simplex[0]) == 0.0 + def test_boundaries_iterator(): st = SimplexTree() @@ -400,16 +405,17 @@ def test_boundaries_iterator(): list(st.get_boundaries([])) with pytest.raises(RuntimeError): - list(st.get_boundaries([0, 4])) # (0, 4) does not exist + list(st.get_boundaries([0, 4])) # (0, 4) does not exist with pytest.raises(RuntimeError): - list(st.get_boundaries([6])) # (6) does not exist + list(st.get_boundaries([6])) # (6) does not exist + def test_persistence_intervals_in_dimension(): # Here is our triangulation of a 2-torus - taken from https://dioscuri-tda.org/Paris_TDA_Tutorial_2021.html # 0-----3-----4-----0 # | \ | \ | \ | \ | - # | \ | \ | \| \ | + # | \ | \ | \| \ | # 1-----8-----7-----1 # | \ | \ | \ | \ | # | \ | \ | \ | \ | @@ -418,50 +424,52 @@ def test_persistence_intervals_in_dimension(): # | \ | \ | \ | \ | # 0-----3-----4-----0 st = SimplexTree() - st.insert([0,1,8]) - st.insert([0,3,8]) - st.insert([3,7,8]) - st.insert([3,4,7]) - st.insert([1,4,7]) - st.insert([0,1,4]) - st.insert([1,2,5]) - st.insert([1,5,8]) - st.insert([5,6,8]) - st.insert([6,7,8]) - st.insert([2,6,7]) - st.insert([1,2,7]) - st.insert([0,2,3]) - st.insert([2,3,5]) - st.insert([3,4,5]) - st.insert([4,5,6]) - st.insert([0,4,6]) - st.insert([0,2,6]) + st.insert([0, 1, 8]) + st.insert([0, 3, 8]) + st.insert([3, 7, 8]) + st.insert([3, 4, 7]) + st.insert([1, 4, 7]) + st.insert([0, 1, 4]) + st.insert([1, 2, 5]) + st.insert([1, 5, 8]) + st.insert([5, 6, 8]) + st.insert([6, 7, 8]) + st.insert([2, 6, 7]) + st.insert([1, 2, 7]) + st.insert([0, 2, 3]) + st.insert([2, 3, 5]) + st.insert([3, 4, 5]) + st.insert([4, 5, 6]) + st.insert([0, 4, 6]) + st.insert([0, 2, 6]) st.compute_persistence(persistence_dim_max=True) - + H0 = st.persistence_intervals_in_dimension(0) - assert np.array_equal(H0, np.array([[ 0., float("inf")]])) + assert np.array_equal(H0, np.array([[0.0, float("inf")]])) H1 = st.persistence_intervals_in_dimension(1) - assert np.array_equal(H1, np.array([[ 0., float("inf")], [ 0., float("inf")]])) + assert np.array_equal(H1, np.array([[0.0, float("inf")], [0.0, float("inf")]])) H2 = st.persistence_intervals_in_dimension(2) - assert np.array_equal(H2, np.array([[ 0., float("inf")]])) + assert np.array_equal(H2, np.array([[0.0, float("inf")]])) # Test empty case assert st.persistence_intervals_in_dimension(3).shape == (0, 2) + def test_equality_operator(): st1 = SimplexTree() st2 = SimplexTree() assert st1 == st2 - st1.insert([1,2,3], 4.) + st1.insert([1, 2, 3], 4.0) assert st1 != st2 - st2.insert([1,2,3], 4.) + st2.insert([1, 2, 3], 4.0) assert st1 == st2 + def test_simplex_tree_deep_copy(): st = SimplexTree() - st.insert([1, 2, 3], 0.) + st.insert([1, 2, 3], 0.0) # compute persistence only on the original st.compute_persistence() @@ -480,14 +488,15 @@ def test_simplex_tree_deep_copy(): for a_splx in a_filt_list: assert a_splx in st_filt_list - + # test double free del st del st_copy + def test_simplex_tree_deep_copy_constructor(): st = SimplexTree() - st.insert([1, 2, 3], 0.) + st.insert([1, 2, 3], 0.0) # compute persistence only on the original st.compute_persistence() @@ -506,56 +515,132 @@ def test_simplex_tree_deep_copy_constructor(): for a_splx in a_filt_list: assert a_splx in st_filt_list - + # test double free del st del st_copy + def test_simplex_tree_constructor_exception(): with pytest.raises(TypeError): - st = SimplexTree(other = "Construction from a string shall raise an exception") + st = SimplexTree(other="Construction from a string shall raise an exception") + + +def test_create_from_array(): + a = np.array([[1, 4, 13, 6], [4, 3, 11, 5], [13, 11, 10, 12], [6, 5, 12, 2]]) + st = SimplexTree.create_from_array(a, max_filtration=5.0) + assert list(st.get_filtration()) == [([0], 1.0), ([3], 2.0), ([1], 3.0), ([0, 1], 4.0), ([1, 3], 5.0)] + + +def test_insert_edges_from_coo_matrix(): + try: + from scipy.sparse import coo_matrix + from scipy.spatial import cKDTree + except ImportError: + print("Skipping, no SciPy") + return + + st = SimplexTree() + st.insert([1, 2, 7], 7) + row = np.array([2, 5, 3]) + col = np.array([1, 4, 6]) + dat = np.array([1, 2, 3]) + edges = coo_matrix((dat, (row, col))) + st.insert_edges_from_coo_matrix(edges) + assert list(st.get_filtration()) == [ + ([1], 1.0), + ([2], 1.0), + ([1, 2], 1.0), + ([4], 2.0), + ([5], 2.0), + ([4, 5], 2.0), + ([3], 3.0), + ([6], 3.0), + ([3, 6], 3.0), + ([7], 7.0), + ([1, 7], 7.0), + ([2, 7], 7.0), + ([1, 2, 7], 7.0), + ] + + pts = np.random.rand(100, 2) + tree = cKDTree(pts) + edges = tree.sparse_distance_matrix(tree, max_distance=0.15, output_type="coo_matrix") + st = SimplexTree() + st.insert_edges_from_coo_matrix(edges) + assert 100 < st.num_simplices() < 1000 + + +def test_insert_batch(): + st = SimplexTree() + # vertices + st.insert_batch(np.array([[6, 1, 5]]), np.array([-5.0, 2.0, -3.0])) + # triangles + st.insert_batch(np.array([[2, 10], [5, 0], [6, 11]]), np.array([4.0, 0.0])) + # edges + st.insert_batch(np.array([[1, 5], [2, 5]]), np.array([1.0, 3.0])) + + assert list(st.get_filtration()) == [ + ([6], -5.0), + ([5], -3.0), + ([0], 0.0), + ([10], 0.0), + ([0, 10], 0.0), + ([11], 0.0), + ([0, 11], 0.0), + ([10, 11], 0.0), + ([0, 10, 11], 0.0), + ([1], 1.0), + ([2], 1.0), + ([1, 2], 1.0), + ([2, 5], 4.0), + ([2, 6], 4.0), + ([5, 6], 4.0), + ([2, 5, 6], 4.0), + ] + def test_expansion_with_blocker(): - st=SimplexTree() - st.insert([0,1],0) - st.insert([0,2],1) - st.insert([0,3],2) - st.insert([1,2],3) - st.insert([1,3],4) - st.insert([2,3],5) - st.insert([2,4],6) - st.insert([3,6],7) - st.insert([4,5],8) - st.insert([4,6],9) - st.insert([5,6],10) - st.insert([6],10) + st = SimplexTree() + st.insert([0, 1], 0) + st.insert([0, 2], 1) + st.insert([0, 3], 2) + st.insert([1, 2], 3) + st.insert([1, 3], 4) + st.insert([2, 3], 5) + st.insert([2, 4], 6) + st.insert([3, 6], 7) + st.insert([4, 5], 8) + st.insert([4, 6], 9) + st.insert([5, 6], 10) + st.insert([6], 10) def blocker(simplex): try: # Block all simplices that contain vertex 6 simplex.index(6) - print(simplex, ' is blocked') + print(simplex, " is blocked") return True except ValueError: - print(simplex, ' is accepted') - st.assign_filtration(simplex, st.filtration(simplex) + 1.) + print(simplex, " is accepted") + st.assign_filtration(simplex, st.filtration(simplex) + 1.0) return False st.expansion_with_blocker(2, blocker) assert st.num_simplices() == 22 assert st.dimension() == 2 - assert st.find([4,5,6]) == False - assert st.filtration([0,1,2]) == 4. - assert st.filtration([0,1,3]) == 5. - assert st.filtration([0,2,3]) == 6. - assert st.filtration([1,2,3]) == 6. + assert st.find([4, 5, 6]) == False + assert st.filtration([0, 1, 2]) == 4.0 + assert st.filtration([0, 1, 3]) == 5.0 + assert st.filtration([0, 2, 3]) == 6.0 + assert st.filtration([1, 2, 3]) == 6.0 st.expansion_with_blocker(3, blocker) assert st.num_simplices() == 23 assert st.dimension() == 3 - assert st.find([4,5,6]) == False - assert st.filtration([0,1,2]) == 4. - assert st.filtration([0,1,3]) == 5. - assert st.filtration([0,2,3]) == 6. - assert st.filtration([1,2,3]) == 6. - assert st.filtration([0,1,2,3]) == 7. + assert st.find([4, 5, 6]) == False + assert st.filtration([0, 1, 2]) == 4.0 + assert st.filtration([0, 1, 3]) == 5.0 + assert st.filtration([0, 2, 3]) == 6.0 + assert st.filtration([1, 2, 3]) == 6.0 + assert st.filtration([0, 1, 2, 3]) == 7.0 |