From 8e4f3d151818b78a29d11cdc6ca171947bfd6dd9 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 3 Mar 2020 15:33:17 +0100 Subject: update wasserstein distance with pot so that it can return optimal matching now! --- src/python/doc/wasserstein_distance_user.rst | 24 ++++++++++ src/python/gudhi/wasserstein.py | 69 ++++++++++++++++++++++------ src/python/test/test_wasserstein_distance.py | 31 +++++++++---- 3 files changed, 102 insertions(+), 22 deletions(-) (limited to 'src') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 94b454e2..d3daa318 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -47,3 +47,27 @@ The output is: .. testoutput:: Wasserstein distance value = 1.45 + +We can also have access to the optimal matching by letting `matching=True`. +It is encoded as a list of indices (i,j), meaning that the i-th point in X +is mapped to the j-th point in Y. +An index of -1 represents the diagonal. + +.. testcode:: + + import gudhi.wasserstein + import numpy as np + + diag1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) + diag2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1]]) + cost, matching = gudhi.wasserstein.wasserstein_distance(diag1, diag2, matching=True, order=1., internal_p=2.) + + message = "Wasserstein distance value = %.2f, optimal matching: %s" %(cost, matching) + print(message) + +The output is: + +.. testoutput:: + + Wasserstein distance value = 2.15, optimal matching: [(0, 0), (1, 2), (2, -1), (-1, 1)] + diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 13102094..ba0f7343 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -62,14 +62,39 @@ def _perstot(X, order, internal_p): return (np.sum(np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**order))**(1./order) -def wasserstein_distance(X, Y, order=2., internal_p=2.): +def _clean_match(match, n, m): ''' - :param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate). + :param match: a list of the form [(i,j) ...] + :param n: int, size of the first dgm + :param m: int, size of the second dgm + :return: a modified version of match where indices greater than n, m are replaced by -1, encoding the diagonal. + and (-1, -1) are removed + ''' + new_match = [] + for i,j in match: + if i >= n: + if j < m: + new_match.append((-1, j)) + elif j >= m: + if i < n: + new_match.append((i,-1)) + else: + new_match.append((i,j)) + return new_match + + +def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): + ''' + :param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points + (i.e. with infinite coordinate). :param Y: (m x 2) numpy.array encoding the second diagram. + :param matching: if True, computes and returns the optimal matching between X and Y, encoded as... :param order: exponent for Wasserstein; Default value is 2. - :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); Default value is 2 (Euclidean norm). - :returns: the Wasserstein distance of order q (1 <= q < infinity) between persistence diagrams with respect to the internal_p-norm as ground metric. - :rtype: float + :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); + Default value is 2 (Euclidean norm). + :returns: the Wasserstein distance of order q (1 <= q < infinity) between persistence diagrams with + respect to the internal_p-norm as ground metric. + If matching is set to True, also returns the optimal matching between X and Y. ''' n = len(X) m = len(Y) @@ -77,21 +102,39 @@ def wasserstein_distance(X, Y, order=2., internal_p=2.): # handle empty diagrams if X.size == 0: if Y.size == 0: - return 0. + if not matching: + return 0. + else: + return 0., [] else: - return _perstot(Y, order, internal_p) + if not matching: + return _perstot(Y, order, internal_p) + else: + return _perstot(Y, order, internal_p), [(-1, j) for j in range(m)] elif Y.size == 0: - return _perstot(X, order, internal_p) + if not matching: + return _perstot(X, order, internal_p) + else: + return _perstot(X, order, internal_p), [(i, -1) for i in range(n)] M = _build_dist_matrix(X, Y, order=order, internal_p=internal_p) - a = np.full(n+1, 1. / (n + m) ) # weight vector of the input diagram. Uniform here. - a[-1] = a[-1] * m # normalized so that we have a probability measure, required by POT - b = np.full(m+1, 1. / (n + m) ) # weight vector of the input diagram. Uniform here. - b[-1] = b[-1] * n # so that we have a probability measure, required by POT + a = np.ones(n+1) # weight vector of the input diagram. Uniform here. + a[-1] = m + b = np.ones(m+1) # weight vector of the input diagram. Uniform here. + b[-1] = n + + if matching: + P = ot.emd(a=a,b=b,M=M, numItermax=2000000) + ot_cost = np.sum(np.multiply(P,M)) + P[P < 0.5] = 0 # trick to avoid numerical issue, could it be improved? + match = np.argwhere(P) + # Now we turn to -1 points encoding the diagonal + match = _clean_match(match, n, m) + return ot_cost ** (1./order) , match # Comptuation of the otcost using the ot.emd2 library. # Note: it is the Wasserstein distance to the power q. # The default numItermax=100000 is not sufficient for some examples with 5000 points, what is a good value? - ot_cost = (n+m) * ot.emd2(a, b, M, numItermax=2000000) + ot_cost = ot.emd2(a, b, M, numItermax=2000000) return ot_cost ** (1./order) diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index 6a6b217b..02a1d2c9 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -51,14 +51,27 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True): assert wasserstein_distance(diag3, diag4, internal_p=1., order=2.) == approx(np.sqrt(5)) assert wasserstein_distance(diag3, diag4, internal_p=4.5, order=2.) == approx(np.sqrt(5)) - if(not test_infinity): - return + if test_infinity: + diag5 = np.array([[0, 3], [4, np.inf]]) + diag6 = np.array([[7, 8], [4, 6], [3, np.inf]]) - diag5 = np.array([[0, 3], [4, np.inf]]) - diag6 = np.array([[7, 8], [4, 6], [3, np.inf]]) + assert wasserstein_distance(diag4, diag5) == np.inf + assert wasserstein_distance(diag5, diag6, order=1, internal_p=np.inf) == approx(4.) + + + if test_matching: + match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1] + assert match == [] + match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] + assert match == [] + match = wasserstein_distance(emptydiag, diag2, matching=True, internal_p=np.inf, order=2.)[1] + assert match == [(-1, 0), (-1, 1)] + match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] + assert match == [(0, -1), (1, -1)] + match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] + assert match == [(0, 0), (1, 1), (2, -1)] + - assert wasserstein_distance(diag4, diag5) == np.inf - assert wasserstein_distance(diag5, diag6, order=1, internal_p=np.inf) == approx(4.) def hera_wrap(delta): def fun(*kargs,**kwargs): @@ -66,8 +79,8 @@ def hera_wrap(delta): return fun def test_wasserstein_distance_pot(): - _basic_wasserstein(pot, 1e-15, test_infinity=False) + _basic_wasserstein(pot, 1e-15, test_infinity=False, test_matching=True) def test_wasserstein_distance_hera(): - _basic_wasserstein(hera_wrap(1e-12), 1e-12) - _basic_wasserstein(hera_wrap(.1), .1) + _basic_wasserstein(hera_wrap(1e-12), 1e-12, test_matching=False) + _basic_wasserstein(hera_wrap(.1), .1, test_matching=False) -- cgit v1.2.3 From 2141ef8adfee531f3eaf822cf4076b9b010e6f94 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 3 Mar 2020 16:22:48 +0100 Subject: correction missing arg in test_wasserstein_distance --- src/python/test/test_wasserstein_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index 02a1d2c9..d0f0323c 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -17,7 +17,7 @@ __author__ = "Theo Lacombe" __copyright__ = "Copyright (C) 2019 Inria" __license__ = "MIT" -def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True): +def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_matching=True): diag1 = np.array([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]]) diag2 = np.array([[2.8, 4.45], [9.5, 14.1]]) diag3 = np.array([[0, 2], [4, 6]]) -- cgit v1.2.3 From 570a9b83eb3f714bc52735dae289a5195874bf41 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Thu, 5 Mar 2020 15:40:45 +0100 Subject: completed as... --- src/python/gudhi/wasserstein.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index ba0f7343..aab0cb3c 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -30,7 +30,9 @@ def _build_dist_matrix(X, Y, order=2., internal_p=2.): :param order: exponent for the Wasserstein metric. :param internal_p: Ground metric (i.e. norm L^p). :returns: (n+1) x (m+1) np.array encoding the cost matrix C. - For 1 <= i <= n, 1 <= j <= m, C[i,j] encodes the distance between X[i] and Y[j], while C[i, m+1] (resp. C[n+1, j]) encodes the distance (to the p) between X[i] (resp Y[j]) and its orthogonal proj onto the diagonal. + For 1 <= i <= n, 1 <= j <= m, C[i,j] encodes the distance between X[i] and Y[j], + while C[i, m+1] (resp. C[n+1, j]) encodes the distance (to the p) between X[i] (resp Y[j]) + and its orthogonal proj onto the diagonal. note also that C[n+1, m+1] = 0 (it costs nothing to move from the diagonal to the diagonal). ''' Xdiag = _proj_on_diag(X) @@ -88,7 +90,9 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): :param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate). :param Y: (m x 2) numpy.array encoding the second diagram. - :param matching: if True, computes and returns the optimal matching between X and Y, encoded as... + :param matching: if True, computes and returns the optimal matching between X and Y, encoded as + a list of tuple [...(i,j)...], meaning the i-th point in X is matched to + the j-th point in Y, with the convention (-1) represents the diagonal. :param order: exponent for Wasserstein; Default value is 2. :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); Default value is 2 (Euclidean norm). -- cgit v1.2.3 From 78ccc10eb0034a4648df303f2913b6b4680b085e Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Fri, 6 Mar 2020 17:35:38 +0100 Subject: Generators for simplex tree --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 86 +++++++++++++++++++++- .../gudhi/Simplex_tree/Simplex_tree_iterators.h | 12 +-- src/Simplex_tree/test/simplex_tree_unit_test.cpp | 36 +++++++++ 3 files changed, 122 insertions(+), 12 deletions(-) (limited to 'src') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 76608008..7315bf45 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -24,6 +24,7 @@ #include #include #include +#include #ifdef GUDHI_USE_TBB #include @@ -246,8 +247,8 @@ class Simplex_tree { * which is consequenlty * equal to \f$(-1)^{\text{dim} \sigma}\f$ the canonical orientation on the simplex. */ - Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) { - assert(sh != null_simplex()); // Empty simplex + Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) const { + GUDHI_CHECK(sh != null_simplex(), "empty simplex"); return Simplex_vertex_range(Simplex_vertex_iterator(this, sh), Simplex_vertex_iterator(this)); } @@ -450,6 +451,15 @@ class Simplex_tree { return true; } + /** \brief Returns the filtration value of a simplex. + * + * Same as `filtration()`, but does not handle `null_simplex()`. + */ + static Filtration_value filtration_(Simplex_handle sh) { + GUDHI_CHECK (sh != null_simplex(), "null simplex"); + return sh->second.filtration(); + } + public: /** \brief Returns the key associated to a simplex. * @@ -827,7 +837,7 @@ class Simplex_tree { /** Returns the Siblings containing a simplex.*/ template - Siblings* self_siblings(SimplexHandle sh) { + static Siblings* self_siblings(SimplexHandle sh) { if (sh->second.children()->parent() == sh->first) return sh->second.children()->oncles(); else @@ -1465,6 +1475,76 @@ class Simplex_tree { } } + /** \brief Returns a vertex of `sh` that has the same filtration value as `sh` if it exists, and `null_vertex()` otherwise. + * + * For a lower-star filtration built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which vertex had its filtration value propagated to `sh`. + * If several vertices have the same filtration value, the one it returns is arbitrary. */ + Vertex_handle vertex_with_same_filtration(Simplex_handle sh) { + auto filt = filtration_(sh); + for(auto v : simplex_vertex_range(sh)) + if(filtration_(find_vertex(v)) == filt) + return v; + return null_vertex(); + } + + /** \brief Returns an edge of `sh` that has the same filtration value as `sh` if it exists, and `null_simplex()` otherwise. + * + * For a flag-complex built with `expansion()`, this is a way to invert the process and find out which edge had its filtration value propagated to `sh`. + * If several edges have the same filtration value, the one it returns is arbitrary. + * + * \pre `sh` must have dimension at least 1. */ + Simplex_handle edge_with_same_filtration(Simplex_handle sh) { +#if 0 + // FIXME: Only do this if dim >= 2, since we don't want to return a vertex... + // Test if we are lucky and the parent has the same filtration value. + Siblings* sib = self_siblings(sh); + Vertex_handle v_par = sib->parent(); + sib = sib->oncles(); + Simplex_handle par = sib->find(v_par); + if(filtration_(par) == filt) return edge_with_same_filtration(par); +#endif + auto&& vertices = simplex_vertex_range(sh); + auto end = std::end(vertices); + auto vi = std::begin(vertices); + GUDHI_CHECK(vi != end, "empty simplex"); + auto v0 = *vi; + ++vi; + GUDHI_CHECK(vi != end, "simplex of dimension 0"); + if(std::next(vi) == end) return sh; // shortcut for dimension 1 + boost::container::static_vector suffix; + suffix.push_back(v0); + auto filt = filtration_(sh); + do + { + Vertex_handle v = *vi; + auto&& children1 = find_vertex(v)->second.children()->members_; + for(auto w : suffix){ + // Can we take advantage of the fact that suffix is ordered? + Simplex_handle s = children1.find(w); + if(filtration_(s) == filt) + return s; + } + suffix.push_back(v); + } + while(++vi != end); + return null_simplex(); + } + + /** \brief Returns an edge of `sh` that has the same filtration value as `sh` if it exists, and `null_simplex()` otherwise. + * + * For a flag-complex built with `expansion()`, this is a way to invert the process and find out which edge had its filtration value propagated to `sh`. + * If several edges have the same filtration value, the one it returns is arbitrary. */ + Simplex_handle minimal_simplex_with_same_filtration(Simplex_handle sh) { + if(dimension(sh) == 0) // vertices are minimal + return sh; + auto filt = filtration_(sh); + // Naive implementation, it can be sped up. + for(auto b : boundary_simplex_range(sh)) + if(filtration_(b) == filt) + return minimal_simplex_with_same_filtration(b); + return sh; // None of its faces has the same filtration. + } + private: Vertex_handle null_vertex_; /** \brief Total number of simplices in the complex, without the empty simplex.*/ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h index efccf2f2..9007b6bd 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h @@ -15,9 +15,7 @@ #include #include -#if BOOST_VERSION >= 105600 -# include -#endif +#include #include @@ -42,13 +40,13 @@ class Simplex_tree_simplex_vertex_iterator : public boost::iterator_facade< typedef typename SimplexTree::Siblings Siblings; typedef typename SimplexTree::Vertex_handle Vertex_handle; - explicit Simplex_tree_simplex_vertex_iterator(SimplexTree * st) + explicit Simplex_tree_simplex_vertex_iterator(SimplexTree const* st) : // any end() iterator sib_(nullptr), v_(st->null_vertex()) { } - Simplex_tree_simplex_vertex_iterator(SimplexTree * st, Simplex_handle sh) + Simplex_tree_simplex_vertex_iterator(SimplexTree const* st, Simplex_handle sh) : sib_(st->self_siblings(sh)), v_(sh->first) { } @@ -166,15 +164,11 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< // Most of the storage should be moved to the range, iterators should be light. Vertex_handle last_; // last vertex of the simplex Vertex_handle next_; // next vertex to push in suffix_ -#if BOOST_VERSION >= 105600 // 40 seems a conservative bound on the dimension of a Simplex_tree for now, // as it would not fit on the biggest hard-drive. boost::container::static_vector suffix_; // static_vector still has some overhead compared to a trivial hand-made // version using std::aligned_storage, or compared to making suffix_ static. -#else - std::vector suffix_; -#endif Siblings * sib_; // where the next search will start from Simplex_handle sh_; // current Simplex_handle in the boundary SimplexTree * st_; // simplex containing the simplicial complex diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index 58bfa8db..2a2e2b25 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -986,5 +986,41 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(insert_duplicated_vertices, typeST, list_of_tested << " - num_simplices = " << st.num_simplices() << std::endl; BOOST_CHECK(st.dimension() == 1); BOOST_CHECK(st.num_simplices() == st.num_vertices() + 1); +} +BOOST_AUTO_TEST_CASE_TEMPLATE(generators, typeST, list_of_tested_variants) { + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST FIND GENERATORS" << std::endl; + { + typeST st; + st.insert_simplex_and_subfaces({0,1,2,3,4,5,6},0); + st.assign_filtration(st.find({0,2,4}), 10); + st.assign_filtration(st.find({1,5}), 20); + st.assign_filtration(st.find({1,2,4}), 30); + st.assign_filtration(st.find({3}), 5); + st.make_filtration_non_decreasing(); + BOOST_CHECK(st.filtration(st.find({1,2}))==0); + BOOST_CHECK(st.filtration(st.find({0,1,2,3,4}))==30); + BOOST_CHECK(st.minimal_simplex_with_same_filtration(st.find({0,1,2,3,4,5}))==st.find({1,2,4})); + BOOST_CHECK(st.minimal_simplex_with_same_filtration(st.find({0,2,3}))==st.find({3})); + auto s=st.minimal_simplex_with_same_filtration(st.find({0,2,6})); + BOOST_CHECK(s==st.find({0})||s==st.find({2})||s==st.find({6})); + BOOST_CHECK(st.vertex_with_same_filtration(st.find({2}))==2); + BOOST_CHECK(st.vertex_with_same_filtration(st.find({1,5}))==st.null_vertex()); + BOOST_CHECK(st.vertex_with_same_filtration(st.find({5,6}))>=5); + } + { + typeST st; + st.insert_simplex_and_subfaces({0,1}, 8); + st.insert_simplex_and_subfaces({0,2}, 10); + st.insert_simplex_and_subfaces({3,4}, 6); + st.insert_simplex_and_subfaces({1,2}, 5); + st.insert_simplex_and_subfaces({1,5}, 4); + st.insert_simplex_and_subfaces({0,5}, 3); + st.insert_simplex_and_subfaces({2,5}, 2); + st.insert_simplex_and_subfaces({1,3}, 9); + st.expansion(50); + BOOST_CHECK(st.edge_with_same_filtration(st.find({0,1,2,5}))==st.find({0,2})); + BOOST_CHECK(st.edge_with_same_filtration(st.find({1,5}))==st.find({1,5})); + } } -- cgit v1.2.3 From 5748647e9e89bb91e361b06c7c6cb053081618bf Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Fri, 6 Mar 2020 18:10:05 +0100 Subject: Fix doc of minimal_simplex_with_same_filtration --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 7315bf45..b6973756 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1530,13 +1530,11 @@ class Simplex_tree { return null_simplex(); } - /** \brief Returns an edge of `sh` that has the same filtration value as `sh` if it exists, and `null_simplex()` otherwise. + /** \brief Returns a minimal face of `sh` that has the same filtration value as `sh`. * - * For a flag-complex built with `expansion()`, this is a way to invert the process and find out which edge had its filtration value propagated to `sh`. - * If several edges have the same filtration value, the one it returns is arbitrary. */ + * For a complex built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which simplex had its filtration value propagated to `sh`. + * If several minimal (for inclusion) simplices have the same filtration value, the one it returns is arbitrary, and it is not guaranteed to be the one with smallest dimension. */ Simplex_handle minimal_simplex_with_same_filtration(Simplex_handle sh) { - if(dimension(sh) == 0) // vertices are minimal - return sh; auto filt = filtration_(sh); // Naive implementation, it can be sped up. for(auto b : boundary_simplex_range(sh)) -- cgit v1.2.3 From 8bd39f74f69e8fcb662873e0e045c953b814f28f Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Fri, 6 Mar 2020 18:21:40 +0100 Subject: Tweak doc. --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index b6973756..ad592a92 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1532,7 +1532,7 @@ class Simplex_tree { /** \brief Returns a minimal face of `sh` that has the same filtration value as `sh`. * - * For a complex built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which simplex had its filtration value propagated to `sh`. + * For a filtration built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which simplex had its filtration value propagated to `sh`. * If several minimal (for inclusion) simplices have the same filtration value, the one it returns is arbitrary, and it is not guaranteed to be the one with smallest dimension. */ Simplex_handle minimal_simplex_with_same_filtration(Simplex_handle sh) { auto filt = filtration_(sh); -- cgit v1.2.3 From 2eca5c75b1fbd7157e2656b875e730dc5f00f373 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 15:45:45 +0100 Subject: removed P[P < 0.5] thresholding ; as it shouldn't happen anymore. --- src/python/gudhi/wasserstein.py | 1 - 1 file changed, 1 deletion(-) (limited to 'src') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index aab0cb3c..e28c63e6 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -130,7 +130,6 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): if matching: P = ot.emd(a=a,b=b,M=M, numItermax=2000000) ot_cost = np.sum(np.multiply(P,M)) - P[P < 0.5] = 0 # trick to avoid numerical issue, could it be improved? match = np.argwhere(P) # Now we turn to -1 points encoding the diagonal match = _clean_match(match, n, m) -- cgit v1.2.3 From 967ceab26b09ad74e0cff0d84429a766af267f6b Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 16:47:09 +0100 Subject: removed _clean_match and changed matching format, it is now a (n x 2) numpy array --- src/python/gudhi/wasserstein.py | 31 ++++++------------------------- 1 file changed, 6 insertions(+), 25 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index e28c63e6..9efa946e 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -64,34 +64,13 @@ def _perstot(X, order, internal_p): return (np.sum(np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**order))**(1./order) -def _clean_match(match, n, m): - ''' - :param match: a list of the form [(i,j) ...] - :param n: int, size of the first dgm - :param m: int, size of the second dgm - :return: a modified version of match where indices greater than n, m are replaced by -1, encoding the diagonal. - and (-1, -1) are removed - ''' - new_match = [] - for i,j in match: - if i >= n: - if j < m: - new_match.append((-1, j)) - elif j >= m: - if i < n: - new_match.append((i,-1)) - else: - new_match.append((i,j)) - return new_match - - def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): ''' :param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate). :param Y: (m x 2) numpy.array encoding the second diagram. :param matching: if True, computes and returns the optimal matching between X and Y, encoded as - a list of tuple [...(i,j)...], meaning the i-th point in X is matched to + a (n x 2) np.array [...[i,j]...], meaning the i-th point in X is matched to the j-th point in Y, with the convention (-1) represents the diagonal. :param order: exponent for Wasserstein; Default value is 2. :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); @@ -114,12 +93,12 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): if not matching: return _perstot(Y, order, internal_p) else: - return _perstot(Y, order, internal_p), [(-1, j) for j in range(m)] + return _perstot(Y, order, internal_p), np.array([[-1, j] for j in range(m)]) elif Y.size == 0: if not matching: return _perstot(X, order, internal_p) else: - return _perstot(X, order, internal_p), [(i, -1) for i in range(n)] + return _perstot(X, order, internal_p), np.array([[i, -1] for i in range(n)]) M = _build_dist_matrix(X, Y, order=order, internal_p=internal_p) a = np.ones(n+1) # weight vector of the input diagram. Uniform here. @@ -130,9 +109,11 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): if matching: P = ot.emd(a=a,b=b,M=M, numItermax=2000000) ot_cost = np.sum(np.multiply(P,M)) + P[-1, -1] = 0 # Remove matching corresponding to the diagonal match = np.argwhere(P) # Now we turn to -1 points encoding the diagonal - match = _clean_match(match, n, m) + match[:,0][match[:,0] >= n] = -1 + match[:,1][match[:,1] >= m] = -1 return ot_cost ** (1./order) , match # Comptuation of the otcost using the ot.emd2 library. -- cgit v1.2.3 From 4aea5deab6ce4cbb491f4c9c2b7e9f023efbbe01 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 17:41:38 +0100 Subject: changed output of matching as a (n x 2) array, adapted tests and doc --- src/python/doc/wasserstein_distance_user.rst | 2 +- src/python/gudhi/wasserstein.py | 2 +- src/python/test/test_wasserstein_distance.py | 10 +++++----- 3 files changed, 7 insertions(+), 7 deletions(-) (limited to 'src') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index d3daa318..9519caa6 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -69,5 +69,5 @@ The output is: .. testoutput:: - Wasserstein distance value = 2.15, optimal matching: [(0, 0), (1, 2), (2, -1), (-1, 1)] + Wasserstein distance value = 2.15, optimal matching: [[0, 0], [1, 2], [2, -1], [-1, 1]] diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 9efa946e..9e4dc7d5 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -88,7 +88,7 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): if not matching: return 0. else: - return 0., [] + return 0., np.array([]) else: if not matching: return _perstot(Y, order, internal_p) diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index d0f0323c..ca9a4a61 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -61,15 +61,15 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat if test_matching: match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1] - assert match == [] + assert np.array_equal(match, np.array([])) match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] - assert match == [] + assert np.array_equal(match, np.array([])) match = wasserstein_distance(emptydiag, diag2, matching=True, internal_p=np.inf, order=2.)[1] - assert match == [(-1, 0), (-1, 1)] + assert np.array_equal(match , np.array([[-1, 0], [-1, 1]])) match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] - assert match == [(0, -1), (1, -1)] + assert np.array_equal(match , np.array([[0, -1], [1, -1]])) match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] - assert match == [(0, 0), (1, 1), (2, -1)] + assert np.array_equal(match, np.array_equal([[0, 0], [1, 1], [2, -1]])) -- cgit v1.2.3 From fc4e10863d103ee6bc22863f48548fe246a3ddd6 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 18:03:21 +0100 Subject: correction of typo in the doc --- src/python/gudhi/wasserstein.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 9e4dc7d5..12337780 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -30,10 +30,10 @@ def _build_dist_matrix(X, Y, order=2., internal_p=2.): :param order: exponent for the Wasserstein metric. :param internal_p: Ground metric (i.e. norm L^p). :returns: (n+1) x (m+1) np.array encoding the cost matrix C. - For 1 <= i <= n, 1 <= j <= m, C[i,j] encodes the distance between X[i] and Y[j], - while C[i, m+1] (resp. C[n+1, j]) encodes the distance (to the p) between X[i] (resp Y[j]) + For 0 <= i < n, 0 <= j < m, C[i,j] encodes the distance between X[i] and Y[j], + while C[i, m] (resp. C[n, j]) encodes the distance (to the p) between X[i] (resp Y[j]) and its orthogonal proj onto the diagonal. - note also that C[n+1, m+1] = 0 (it costs nothing to move from the diagonal to the diagonal). + note also that C[n, m] = 0 (it costs nothing to move from the diagonal to the diagonal). ''' Xdiag = _proj_on_diag(X) Ydiag = _proj_on_diag(Y) -- cgit v1.2.3 From 6c369a6aa566dfcb8cdb501d0c39eafb32219669 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 18:08:15 +0100 Subject: fix typo in test_wasserstein_distance --- src/python/test/test_wasserstein_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index ca9a4a61..f92208c0 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -69,7 +69,7 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] assert np.array_equal(match , np.array([[0, -1], [1, -1]])) match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] - assert np.array_equal(match, np.array_equal([[0, 0], [1, 1], [2, -1]])) + assert np.array_equal(match, np.array([[0, 0], [1, 1], [2, -1]])) -- cgit v1.2.3 From 753290475ab6e95c2de1baad97ee6f755a0ce19a Mon Sep 17 00:00:00 2001 From: Théo Lacombe Date: Tue, 10 Mar 2020 18:25:10 +0100 Subject: Update src/python/gudhi/wasserstein.py Co-Authored-By: Marc Glisse --- src/python/gudhi/wasserstein.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 12337780..83a682df 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -32,7 +32,7 @@ def _build_dist_matrix(X, Y, order=2., internal_p=2.): :returns: (n+1) x (m+1) np.array encoding the cost matrix C. For 0 <= i < n, 0 <= j < m, C[i,j] encodes the distance between X[i] and Y[j], while C[i, m] (resp. C[n, j]) encodes the distance (to the p) between X[i] (resp Y[j]) - and its orthogonal proj onto the diagonal. + and its orthogonal projection onto the diagonal. note also that C[n, m] = 0 (it costs nothing to move from the diagonal to the diagonal). ''' Xdiag = _proj_on_diag(X) -- cgit v1.2.3 From c9d6e27495c8927d736d593afb0450360b46ccc9 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 18:55:19 +0100 Subject: fix indentation in wasserstein --- src/python/gudhi/wasserstein.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 83a682df..3dd993f9 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -70,13 +70,13 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): (i.e. with infinite coordinate). :param Y: (m x 2) numpy.array encoding the second diagram. :param matching: if True, computes and returns the optimal matching between X and Y, encoded as - a (n x 2) np.array [...[i,j]...], meaning the i-th point in X is matched to - the j-th point in Y, with the convention (-1) represents the diagonal. + a (n x 2) np.array [...[i,j]...], meaning the i-th point in X is matched to + the j-th point in Y, with the convention (-1) represents the diagonal. :param order: exponent for Wasserstein; Default value is 2. :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); - Default value is 2 (Euclidean norm). + Default value is 2 (Euclidean norm). :returns: the Wasserstein distance of order q (1 <= q < infinity) between persistence diagrams with - respect to the internal_p-norm as ground metric. + respect to the internal_p-norm as ground metric. If matching is set to True, also returns the optimal matching between X and Y. ''' n = len(X) -- cgit v1.2.3 From a17a09a2c58bba79e897d0ba00aada05da556967 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Wed, 11 Mar 2020 10:41:53 +0100 Subject: clean test_wasserstein from useless np.array --- src/python/test/test_wasserstein_distance.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index f92208c0..0d70e11a 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -61,15 +61,15 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat if test_matching: match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1] - assert np.array_equal(match, np.array([])) + assert np.array_equal(match, []) match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] - assert np.array_equal(match, np.array([])) + assert np.array_equal(match, []) match = wasserstein_distance(emptydiag, diag2, matching=True, internal_p=np.inf, order=2.)[1] - assert np.array_equal(match , np.array([[-1, 0], [-1, 1]])) + assert np.array_equal(match , [[-1, 0], [-1, 1]]) match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] - assert np.array_equal(match , np.array([[0, -1], [1, -1]])) + assert np.array_equal(match , [[0, -1], [1, -1]]) match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] - assert np.array_equal(match, np.array([[0, 0], [1, 1], [2, -1]])) + assert np.array_equal(match, [[0, 0], [1, 1], [2, -1]]) -- cgit v1.2.3 From 5c55e976606b4dd020bd4e21c93ae22143ef5348 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 18:01:16 +0100 Subject: changed doc of matchings for a more explicit (and hopefully sphinx-valid) version --- src/python/doc/wasserstein_distance_user.rst | 29 ++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 9519caa6..4c3b53dd 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -58,16 +58,29 @@ An index of -1 represents the diagonal. import gudhi.wasserstein import numpy as np - diag1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) - diag2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1]]) - cost, matching = gudhi.wasserstein.wasserstein_distance(diag1, diag2, matching=True, order=1., internal_p=2.) - - message = "Wasserstein distance value = %.2f, optimal matching: %s" %(cost, matching) - print(message) + dgm1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) + dgm2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1]]) + cost, matchings = gudhi.wasserstein.wasserstein_distance(diag1, diag2, matching=True, order=1., internal_p=2.) + + message_cost = "Wasserstein distance value = %.2f" %cost + print(message_cost) + dgm1_to_diagonal = matchings[np.where(matchings[:,0] == -1)][:,1] + dgm2_to_diagonal = matchings[np.where(matchings[:,1] == -1)][:,0] + off_diagonal_match = np.delete(matchings, np.where(matchings == -1)[0], axis=0) + + for i,j in off_diagonal_match: + print("point %s in dgm1 is matched to point %s in dgm2" %(i,j)) + for i in dgm1_to_diagonal: + print("point %s in dgm1 is matched to the diagonal" %i) + for j in dgm2_to_diagonal: + print("point %s in dgm2 is matched to the diagonal" %j) The output is: .. testoutput:: - Wasserstein distance value = 2.15, optimal matching: [[0, 0], [1, 2], [2, -1], [-1, 1]] - + Wasserstein distance value = 2.15 + point 0 in dgm1 is matched to point 0 in dgm2 + point 1 in dgm1 is matched to point 2 in dgm2 + point 2 in dgm1 is matched to the diagonal + point 1 in dgm2 is matched to the diagonal -- cgit v1.2.3 From 66f0b08a8f8d5006f8d29352c169525cc53a22e6 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 19:11:30 +0100 Subject: changed typo in doc (diag --> dgm), used integer for order and internal p, simplify th ecode --- src/python/doc/wasserstein_distance_user.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 4c3b53dd..f43b2217 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -36,10 +36,10 @@ Note that persistence diagrams must be submitted as (n x 2) numpy arrays and mus import gudhi.wasserstein import numpy as np - diag1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) - diag2 = np.array([[2.8, 4.45],[9.5, 14.1]]) + dgm1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) + dgm2 = np.array([[2.8, 4.45],[9.5, 14.1]]) - message = "Wasserstein distance value = " + '%.2f' % gudhi.wasserstein.wasserstein_distance(diag1, diag2, order=1., internal_p=2.) + message = "Wasserstein distance value = " + '%.2f' % gudhi.wasserstein.wasserstein_distance(dgm1, dgm2, order=1., internal_p=2.) print(message) The output is: @@ -60,12 +60,12 @@ An index of -1 represents the diagonal. dgm1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) dgm2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1]]) - cost, matchings = gudhi.wasserstein.wasserstein_distance(diag1, diag2, matching=True, order=1., internal_p=2.) + cost, matchings = gudhi.wasserstein.wasserstein_distance(dgm1, dgm2, matching=True, order=1, internal_p=2) message_cost = "Wasserstein distance value = %.2f" %cost print(message_cost) - dgm1_to_diagonal = matchings[np.where(matchings[:,0] == -1)][:,1] - dgm2_to_diagonal = matchings[np.where(matchings[:,1] == -1)][:,0] + dgm1_to_diagonal = matching[matching[:,0] == -1, 1] + dgm2_to_diagonal = matching[matching[:,1] == -1, 0] off_diagonal_match = np.delete(matchings, np.where(matchings == -1)[0], axis=0) for i,j in off_diagonal_match: -- cgit v1.2.3 From a253c0c4f54a9a148740ed9c20457df0ea43c842 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 19:36:07 +0100 Subject: correction typo in usr.rst --- src/python/doc/wasserstein_distance_user.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index f43b2217..25e51d68 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -64,8 +64,8 @@ An index of -1 represents the diagonal. message_cost = "Wasserstein distance value = %.2f" %cost print(message_cost) - dgm1_to_diagonal = matching[matching[:,0] == -1, 1] - dgm2_to_diagonal = matching[matching[:,1] == -1, 0] + dgm1_to_diagonal = matchings[matchings[:,0] == -1, 1] + dgm2_to_diagonal = matchings[matchings[:,1] == -1, 0] off_diagonal_match = np.delete(matchings, np.where(matchings == -1)[0], axis=0) for i,j in off_diagonal_match: -- cgit v1.2.3 From 60d11e3f06e08b66e49997f389c4dc01b00b793f Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 21:17:38 +0100 Subject: correction of typo in usr.rst --- src/python/doc/wasserstein_distance_user.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 25e51d68..a9b21fa5 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -64,8 +64,8 @@ An index of -1 represents the diagonal. message_cost = "Wasserstein distance value = %.2f" %cost print(message_cost) - dgm1_to_diagonal = matchings[matchings[:,0] == -1, 1] - dgm2_to_diagonal = matchings[matchings[:,1] == -1, 0] + dgm1_to_diagonal = matchings[matchings[:,1] == -1, 0] + dgm2_to_diagonal = matchings[matchings[:,0] == -1, 1] off_diagonal_match = np.delete(matchings, np.where(matchings == -1)[0], axis=0) for i,j in off_diagonal_match: -- cgit v1.2.3 From 47e72c7250bff568735df829a40bcbea0a48f7c2 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 17 Mar 2020 13:28:35 +0100 Subject: Remove code that was commented out. --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) (limited to 'src') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index ad592a92..af711075 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1494,16 +1494,8 @@ class Simplex_tree { * * \pre `sh` must have dimension at least 1. */ Simplex_handle edge_with_same_filtration(Simplex_handle sh) { -#if 0 - // FIXME: Only do this if dim >= 2, since we don't want to return a vertex... - // Test if we are lucky and the parent has the same filtration value. - Siblings* sib = self_siblings(sh); - Vertex_handle v_par = sib->parent(); - sib = sib->oncles(); - Simplex_handle par = sib->find(v_par); - if(filtration_(par) == filt) return edge_with_same_filtration(par); -#endif - auto&& vertices = simplex_vertex_range(sh); + // See issue #251 for potential speed improvements. + auto&& vertices = simplex_vertex_range(sh); // vertices in decreasing order auto end = std::end(vertices); auto vi = std::begin(vertices); GUDHI_CHECK(vi != end, "empty simplex"); -- cgit v1.2.3