diff options
| -rw-r--r-- | .github/for_maintainers/new_gudhi_version_creation.md | 1 | ||||
| -rw-r--r-- | biblio/bibliography.bib | 26 | ||||
| -rw-r--r-- | src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h | 25 | ||||
| -rw-r--r-- | src/python/CMakeLists.txt | 7 | ||||
| -rw-r--r-- | src/python/doc/rips_complex_ref.rst | 13 | ||||
| -rw-r--r-- | src/python/doc/rips_complex_sum.inc | 3 | ||||
| -rw-r--r-- | src/python/doc/rips_complex_user.rst | 51 | ||||
| -rw-r--r-- | src/python/doc/wasserstein_distance_sum.inc | 20 | ||||
| -rw-r--r-- | src/python/gudhi/cubical_complex.pyx | 53 | ||||
| -rw-r--r-- | src/python/gudhi/periodic_cubical_complex.pyx | 53 | ||||
| -rw-r--r-- | src/python/gudhi/weighted_rips_complex.py | 59 | ||||
| -rw-r--r-- | src/python/include/Persistent_cohomology_interface.h | 44 | ||||
| -rwxr-xr-x | src/python/test/test_cubical_complex.py | 10 | ||||
| -rw-r--r-- | src/python/test/test_weighted_rips.py | 63 |
14 files changed, 417 insertions, 11 deletions
diff --git a/.github/for_maintainers/new_gudhi_version_creation.md b/.github/for_maintainers/new_gudhi_version_creation.md index f176d392..8674222b 100644 --- a/.github/for_maintainers/new_gudhi_version_creation.md +++ b/.github/for_maintainers/new_gudhi_version_creation.md @@ -16,6 +16,7 @@ rm -rf data/points/COIL_database/lucky_cat.off_dist data/points/COIL_database/lu Checkin the modifications, build and test the version: ```bash +git submodule update --init mkdir build cd build cmake -DCGAL_DIR=/your/path/to/CGAL -DWITH_GUDHI_EXAMPLE=ON -DWITH_GUDHI_BENCHMARK=ON -DUSER_VERSION_DIR=gudhi.@GUDHI_VERSION@ -DPython_ADDITIONAL_VERSIONS=3 .. diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index 3ea2f59f..677ed4df 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -1253,3 +1253,29 @@ year = "2011" year={2014}, publisher={Springer} } + +@inproceedings{dtmfiltrations, + author = {Hirokazu Anai and + Fr{\'{e}}d{\'{e}}ric Chazal and + Marc Glisse and + Yuichi Ike and + Hiroya Inakoshi and + Rapha{\"{e}}l Tinarrage and + Yuhei Umeda}, + editor = {Gill Barequet and + Yusu Wang}, + title = {DTM-Based Filtrations}, + booktitle = {35th International Symposium on Computational Geometry, SoCG 2019, + June 18-21, 2019, Portland, Oregon, {USA}}, + series = {LIPIcs}, + volume = {129}, + pages = {58:1--58:15}, + publisher = {Schloss Dagstuhl - Leibniz-Zentrum f{\"{u}}r Informatik}, + year = {2019}, + url = {https://doi.org/10.4230/LIPIcs.SoCG.2019.58}, + doi = {10.4230/LIPIcs.SoCG.2019.58}, + timestamp = {Tue, 11 Feb 2020 15:52:14 +0100}, + biburl = {https://dblp.org/rec/conf/compgeom/AnaiCGIITU19.bib}, + bibsource = {dblp computer science bibliography, https://dblp.org} +} + diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index e6a78a6d..f8f80ded 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -13,6 +13,8 @@ #include <gudhi/Bitmap_cubical_complex/counter.h> +#include <boost/config.hpp> + #include <iostream> #include <vector> #include <string> @@ -110,6 +112,16 @@ class Bitmap_cubical_complex_base { virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const; /** + * This function finds a top-dimensional cell that is incident to the input cell and has + * the same filtration value. In case several cells are suitable, an arbitrary one is + * returned. Note that the input parameter can be a cell of any dimension (vertex, edge, etc). + * On the other hand, the output is always indicating the position of + * a top-dimensional cube in the data structure. + * \pre The filtration values are assigned as per `impose_lower_star_filtration()`. + **/ + inline size_t get_top_dimensional_coface_of_a_cell(size_t splx); + + /** * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of * dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence * between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$. @@ -603,6 +615,19 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells } template <typename T> +size_t Bitmap_cubical_complex_base<T>::get_top_dimensional_coface_of_a_cell(size_t splx) { + if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;} + else{ + for (auto v : this->get_coboundary_of_a_cell(splx)){ + if(this->get_cell_data(v) == this->get_cell_data(splx)){ + return this->get_top_dimensional_coface_of_a_cell(v); + } + } + } + BOOST_UNREACHABLE_RETURN(-2); +} + +template <typename T> Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) { this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells); diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 976a8b52..ab08cd6d 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -57,6 +57,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'representations', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'wasserstein', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'point_cloud', ") + set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'weighted_rips_complex', ") add_gudhi_debug_info("Python version ${PYTHON_VERSION_STRING}") add_gudhi_debug_info("Cython version ${CYTHON_VERSION}") @@ -232,6 +233,7 @@ if(PYTHONINTERP_FOUND) file(COPY "gudhi/representations" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/") file(COPY "gudhi/wasserstein" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") + file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") add_custom_command( OUTPUT gudhi.so @@ -488,6 +490,11 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_dtm) endif() + # Weighted Rips + if(SCIPY_FOUND) + add_gudhi_py_test(test_weighted_rips) + endif() + # Set missing or not modules set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") else(CYTHON_FOUND) diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index 22b5616c..5f3e46c1 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -12,3 +12,16 @@ Rips complex reference manual :show-inheritance: .. automethod:: gudhi.RipsComplex.__init__ + +.. _weighted-rips-complex-reference-manual: + +====================================== +Weighted Rips complex reference manual +====================================== + +.. autoclass:: gudhi.weighted_rips_complex.WeightedRipsComplex + :members: + :undoc-members: + :show-inheritance: + + .. automethod:: gudhi.weighted_rips_complex.WeightedRipsComplex.__init__ diff --git a/src/python/doc/rips_complex_sum.inc b/src/python/doc/rips_complex_sum.inc index 6feb74cd..f7580714 100644 --- a/src/python/doc/rips_complex_sum.inc +++ b/src/python/doc/rips_complex_sum.inc @@ -11,6 +11,9 @@ | | | | | | This complex can be built from a point cloud and a distance function, | | | | or from a distance matrix. | | + | | | | + | | Weighted Rips complex constructs a simplicial complex from a distance | | + | | matrix and weights on vertices. | | +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+ | * :doc:`rips_complex_user` | * :doc:`rips_complex_ref` | +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/rips_complex_user.rst b/src/python/doc/rips_complex_user.rst index 8efb12e6..819568be 100644 --- a/src/python/doc/rips_complex_user.rst +++ b/src/python/doc/rips_complex_user.rst @@ -347,3 +347,54 @@ until dimension 1 - one skeleton graph in other words), the output is: points in the persistence diagram will be under the diagonal, and bottleneck distance and persistence graphical tool will not work properly, this is a known issue. + +Weighted Rips Complex +--------------------- + +`WeightedRipsComplex <rips_complex_ref.html#weighted-rips-complex-reference-manual>`_ builds a simplicial complex from a distance matrix and weights on vertices. + + +Example from a distance matrix and weights +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The following example computes the weighted Rips filtration associated with a distance matrix and weights on vertices. + +.. testcode:: + + from gudhi.weighted_rips_complex import WeightedRipsComplex + dist = [[], [1]] + weights = [1, 100] + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights) + st = w_rips.create_simplex_tree(max_dimension=2) + print(list(st.get_filtration())) + +The output is: + +.. testoutput:: + + [([0], 2.0), ([1], 200.0), ([0, 1], 200.0)] + +Example from a point cloud combined with DistanceToMeasure +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Combining with DistanceToMeasure, one can compute the DTM-filtration of a point set, as in `this notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-DTM-filtrations.ipynb>`_. + +.. testcode:: + + import numpy as np + from scipy.spatial.distance import cdist + from gudhi.point_cloud.dtm import DistanceToMeasure + from gudhi.weighted_rips_complex import WeightedRipsComplex + pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]]) + dist = cdist(pts,pts) + dtm = DistanceToMeasure(2, q=2, metric="precomputed") + r = dtm.fit_transform(dist) + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r) + st = w_rips.create_simplex_tree(max_dimension=2) + print(st.persistence()) + +The output is: + +.. testoutput:: + + [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))] diff --git a/src/python/doc/wasserstein_distance_sum.inc b/src/python/doc/wasserstein_distance_sum.inc index f9308e5e..c41de017 100644 --- a/src/python/doc/wasserstein_distance_sum.inc +++ b/src/python/doc/wasserstein_distance_sum.inc @@ -1,14 +1,12 @@ .. table:: :widths: 30 40 30 - +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+ - | .. figure:: | The q-Wasserstein distance measures the similarity between two | :Author: Theo Lacombe | - | ../../doc/Bottleneck_distance/perturb_pd.png | persistence diagrams using the sum of all edges lengths (instead of | | - | :figclass: align-center | the maximum). It allows to define sophisticated objects such as | :Since: GUDHI 3.1.0 | - | | barycenters of a family of persistence diagrams. | | - | | | :License: MIT | - | | | | - | | | :Requires: Python Optimal Transport (POT) :math:`\geq` 0.5.1 | - +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+ - | * :doc:`wasserstein_distance_user` | | - +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------+ + +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------+ + | .. figure:: | The q-Wasserstein distance measures the similarity between two | :Author: Theo Lacombe, Marc Glisse | + | ../../doc/Bottleneck_distance/perturb_pd.png | persistence diagrams using the sum of all edges lengths (instead of | | + | :figclass: align-center | the maximum). It allows to define sophisticated objects such as | :Since: GUDHI 3.1.0 | + | | barycenters of a family of persistence diagrams. | | + | | | :License: MIT, BSD-3-Clause | + +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------+ + | * :doc:`wasserstein_distance_user` | | + +-----------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 007abcb6..ca979eda 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -37,6 +37,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) void compute_persistence(int homology_coeff_field, double min_persistence) vector[pair[int, pair[double, double]]] get_persistence() + vector[vector[int]] cofaces_of_cubical_persistence_pairs() vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) @@ -170,6 +171,58 @@ cdef class CubicalComplex: self.compute_persistence(homology_coeff_field, min_persistence) return self.pcohptr.get_persistence() + def cofaces_of_persistence_pairs(self): + """A persistence interval is described by a pair of cells, one that creates the + feature and one that kills it. The filtration values of those 2 cells give coordinates + for a point in a persistence diagram, or a bar in a barcode. Structurally, in the + cubical complexes provided here, the filtration value of any cell is the minimum of the + filtration values of the maximal cells that contain it. Connecting persistence diagram + coordinates to the corresponding value in the input (i.e. the filtration values of + the top-dimensional cells) is useful for differentiation purposes. + + This function returns a list of pairs of top-dimensional cells corresponding to + the persistence birth and death cells of the filtration. The cells are represented by + their indices in the input list of top-dimensional cells (and not their indices in the + internal datastructure that includes non-maximal cells). Note that when two adjacent + top-dimensional cells have the same filtration value, we arbitrarily return one of the two + when calling the function on one of their common faces. + + :returns: The top-dimensional cells/cofaces of the positive and negative cells, + together with the corresponding homological dimension, in two lists of numpy arrays of integers. + The first list contains the regular persistence pairs, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 2]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell, + index of negative top-dimensional cell). + The second list contains the essential features, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 1]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell). + """ + + assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()" + + cdef vector[vector[int]] persistence_result + output = [[],[]] + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + pr = np.array(persistence_result) + + ess_ind = np.argwhere(pr[:,2] == -1)[:,0] + ess = pr[ess_ind] + max_h = max(ess[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(ess[:,0] == h)[:,0] + output[1].append(ess[hidxs][:,1]) + + reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) + reg = pr[reg_ind] + max_h = max(reg[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(reg[:,0] == h)[:,0] + output[0].append(reg[hidxs][:,1:]) + + return output + def betti_numbers(self): """This function returns the Betti numbers of the complex. diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 246a3a02..06309772 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -34,6 +34,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max) void compute_persistence(int homology_coeff_field, double min_persistence) vector[pair[int, pair[double, double]]] get_persistence() + vector[vector[int]] cofaces_of_cubical_persistence_pairs() vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) @@ -175,6 +176,58 @@ cdef class PeriodicCubicalComplex: self.compute_persistence(homology_coeff_field, min_persistence) return self.pcohptr.get_persistence() + def cofaces_of_persistence_pairs(self): + """A persistence interval is described by a pair of cells, one that creates the + feature and one that kills it. The filtration values of those 2 cells give coordinates + for a point in a persistence diagram, or a bar in a barcode. Structurally, in the + cubical complexes provided here, the filtration value of any cell is the minimum of the + filtration values of the maximal cells that contain it. Connecting persistence diagram + coordinates to the corresponding value in the input (i.e. the filtration values of + the top-dimensional cells) is useful for differentiation purposes. + + This function returns a list of pairs of top-dimensional cells corresponding to + the persistence birth and death cells of the filtration. The cells are represented by + their indices in the input list of top-dimensional cells (and not their indices in the + internal datastructure that includes non-maximal cells). Note that when two adjacent + top-dimensional cells have the same filtration value, we arbitrarily return one of the two + when calling the function on one of their common faces. + + :returns: The top-dimensional cells/cofaces of the positive and negative cells, + together with the corresponding homological dimension, in two lists of numpy arrays of integers. + The first list contains the regular persistence pairs, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 2]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell, + index of negative top-dimensional cell). + The second list contains the essential features, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 1]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell). + """ + cdef vector[vector[int]] persistence_result + if self.pcohptr != NULL: + output = [[],[]] + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + pr = np.array(persistence_result) + + ess_ind = np.argwhere(pr[:,2] == -1)[:,0] + ess = pr[ess_ind] + max_h = max(ess[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(ess[:,0] == h)[:,0] + output[1].append(ess[hidxs][:,1]) + + reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) + reg = pr[reg_ind] + max_h = max(reg[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(reg[:,0] == h)[:,0] + output[0].append(reg[hidxs][:,1:]) + else: + print("cofaces_of_persistence_pairs function requires persistence function" + " to be launched first.") + return output + def betti_numbers(self): """This function returns the Betti numbers of the complex. diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py new file mode 100644 index 00000000..0541572b --- /dev/null +++ b/src/python/gudhi/weighted_rips_complex.py @@ -0,0 +1,59 @@ +# 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): Raphaël Tinarrage, Yuichi Ike, Masatoshi Takenouchi +# +# Copyright (C) 2020 Inria, Copyright (C) 2020 FUjitsu Laboratories Ltd. +# +# Modification(s): +# - YYYY/MM Author: Description of the modification + +from gudhi import SimplexTree + +class WeightedRipsComplex: + """ + Class to generate a weighted Rips complex from a distance matrix and weights on vertices, + in the way described in :cite:`dtmfiltrations`. + Remark that all the filtration values are doubled compared to the definition in the paper + for the consistency with RipsComplex. + """ + def __init__(self, + distance_matrix, + weights=None, + max_filtration=float('inf')): + """ + Args: + distance_matrix (Sequence[Sequence[float]]): distance matrix (full square or lower triangular). + weights (Sequence[float]): (one half of) weight for each vertex. + max_filtration (float): specifies the maximal filtration value to be considered. + """ + self.distance_matrix = distance_matrix + if weights is not None: + self.weights = weights + else: + self.weights = [0] * len(distance_matrix) + self.max_filtration = max_filtration + + def create_simplex_tree(self, max_dimension): + """ + Args: + max_dimension (int): graph expansion until this given dimension. + """ + dist = self.distance_matrix + F = self.weights + num_pts = len(dist) + + st = SimplexTree() + + for i in range(num_pts): + if 2*F[i] <= self.max_filtration: + st.insert([i], 2*F[i]) + for i in range(num_pts): + for j in range(i): + value = max(2*F[i], 2*F[j], dist[i][j] + F[i] + F[j]) + # max is needed when F is not 1-Lipschitz + if value <= self.max_filtration: + st.insert([i,j], filtration=value) + + st.expansion(max_dimension) + return st + diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 0de9bd5c..e5a3dfba 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -13,9 +13,11 @@ #include <gudhi/Persistent_cohomology.h> +#include <cstdlib> #include <vector> #include <utility> // for std::pair #include <algorithm> // for sort +#include <unordered_map> namespace Gudhi { @@ -67,6 +69,48 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol return persistence; } + // This function computes the top-dimensional cofaces associated to the positive and negative + // simplices of a cubical complex. The output format is a vector of vectors of three integers, + // which are [homological dimension, index of top-dimensional coface of positive simplex, + // index of top-dimensional coface of negative simplex]. If the topological feature is essential, + // then the index of top-dimensional coface of negative simplex is arbitrarily set to -1. + std::vector<std::vector<int>> cofaces_of_cubical_persistence_pairs() { + + // Warning: this function is meant to be used with CubicalComplex only!! + + auto&& pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex, + persistent_cohomology::Field_Zp>::get_persistent_pairs(); + + // Gather all top-dimensional cells and store their simplex handles + std::vector<std::size_t> max_splx; + for (auto splx : stptr_->top_dimensional_cells_range()) + max_splx.push_back(splx); + // Sort these simplex handles and compute the ordering function + // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data + std::unordered_map<std::size_t, int> order; + //std::sort(max_splx.begin(), max_splx.end()); + for (unsigned int i = 0; i < max_splx.size(); i++) order.emplace(max_splx[i], i); + + std::vector<std::vector<int>> persistence_pairs; + for (auto pair : pairs) { + int h = stptr_->dimension(get<0>(pair)); + // Recursively get the top-dimensional cell / coface associated to the persistence generator + std::size_t face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair)); + // Retrieve the index of the corresponding top-dimensional cell in the input data + int splx0 = order[face0]; + + int splx1 = -1; + if (get<1>(pair) != stptr_->null_simplex()){ + // Recursively get the top-dimensional cell / coface associated to the persistence generator + std::size_t face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair)); + // Retrieve the index of the corresponding top-dimensional cell in the input data + splx1 = order[face1]; + } + persistence_pairs.push_back({ h, splx0, splx1 }); + } + return persistence_pairs; + } + std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() { std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs; auto const& pairs = Base::get_persistent_pairs(); diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index fce4875c..5c59db8f 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -147,3 +147,13 @@ def test_connected_sublevel_sets(): periodic_dimensions = periodic_dimensions) assert cub.persistence() == [(0, (2.0, float("inf")))] assert cub.betti_numbers() == [1, 0, 0] + +def test_cubical_generators(): + cub = CubicalComplex(top_dimensional_cells = [[0, 0, 0], [0, 1, 0], [0, 0, 0]]) + cub.persistence() + g = cub.cofaces_of_persistence_pairs() + assert len(g[0]) == 2 + assert len(g[1]) == 1 + assert np.array_equal(g[0][0], np.empty(shape=[0,2])) + assert np.array_equal(g[0][1], np.array([[7, 4]])) + assert np.array_equal(g[1][0], np.array([8])) diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips.py new file mode 100644 index 00000000..7ef48333 --- /dev/null +++ b/src/python/test/test_weighted_rips.py @@ -0,0 +1,63 @@ +""" 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): Yuichi Ike and Masatoshi Takenouchi + + Copyright (C) 2020 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +from gudhi.weighted_rips_complex import WeightedRipsComplex +from gudhi.point_cloud.dtm import DistanceToMeasure +import numpy as np +from math import sqrt +from scipy.spatial.distance import cdist +import pytest + +def test_non_dtm_rips_complex(): + dist = [[], [1]] + weights = [1, 100] + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights) + st = w_rips.create_simplex_tree(max_dimension=2) + assert st.filtration([0,1]) == pytest.approx(200.0) + +def test_compatibility_with_rips(): + distance_matrix = [[0], [1, 0], [1, sqrt(2), 0], [sqrt(2), 1, 1, 0]] + w_rips = WeightedRipsComplex(distance_matrix=distance_matrix,max_filtration=42) + st = w_rips.create_simplex_tree(max_dimension=1) + assert list(st.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([2], 0.0), + ([3], 0.0), + ([0, 1], 1.0), + ([0, 2], 1.0), + ([1, 3], 1.0), + ([2, 3], 1.0), + ([1, 2], sqrt(2)), + ([0, 3], sqrt(2)), + ] + +def test_compatibility_with_filtered_rips(): + distance_matrix = [[0], [1, 0], [1, sqrt(2), 0], [sqrt(2), 1, 1, 0]] + w_rips = WeightedRipsComplex(distance_matrix=distance_matrix,max_filtration=1.0) + st = w_rips.create_simplex_tree(max_dimension=1) + + assert st.__is_defined() == True + assert st.__is_persistence_defined() == False + + assert st.num_simplices() == 8 + assert st.num_vertices() == 4 + +def test_dtm_rips_complex(): + pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]]) + dist = cdist(pts,pts) + dtm = DistanceToMeasure(2, q=2, metric="precomputed") + r = dtm.fit_transform(dist) + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r) + st = w_rips.create_simplex_tree(max_dimension=2) + st.persistence() + persistence_intervals0 = st.persistence_intervals_in_dimension(0) + assert persistence_intervals0 == pytest.approx(np.array([[3.16227766, 5.39834564],[3.16227766, 5.39834564], [3.16227766, float("inf")]])) + |