diff options
33 files changed, 2985 insertions, 736 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 005df6d7..d1a3c7bf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -123,7 +123,7 @@ else() include_directories(src/common/include/) include_directories(src/Alpha_complex/include/) include_directories(src/Bitmap_cubical_complex/include/) - include_directories(src/Bottleneck/include/) + include_directories(src/Bottleneck_distance/include/) include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) include_directories(src/Persistent_cohomology/include/) @@ -156,6 +156,8 @@ else() add_subdirectory(src/Tangential_complex/example) add_subdirectory(src/Tangential_complex/test) add_subdirectory(src/Tangential_complex/benchmark) + add_subdirectory(src/Bottleneck_distance/example) + add_subdirectory(src/Bottleneck_distance/test) # data points generator add_subdirectory(data/points/generator) diff --git a/biblio/how_to_cite_gudhi.bib b/biblio/how_to_cite_gudhi.bib index 03c05728..d88ecb75 100644 --- a/biblio/how_to_cite_gudhi.bib +++ b/biblio/how_to_cite_gudhi.bib @@ -33,7 +33,7 @@ , year = 2015 } -@incollection{gudhi:SkeletonBlocker +@incollection{gudhi:Skeleton-Blocker , author = "David Salinas" , title = "Skeleton-Blocker" , publisher = "{GUDHI Editorial Board}" @@ -78,20 +78,11 @@ , year = 2016 } -@incollection{gudhi:SpatialSearching -, author = "Cl\'ement Jamin" -, title = "Spatial searching" +@incollection{gudhi:BottleneckDistance +, author = "Fran{{\c{c}}ois Godi" +, title = "Bottleneck distance" , publisher = "{GUDHI Editorial Board}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "http://gudhi.gforge.inria.fr/doc/latest/group__spatial__searching.html" +, url = "http://gudhi.gforge.inria.fr/doc/latest/group__bottleneck__distance.html" , year = 2016 } - -@incollection{gudhi:TangentialComplex -, author = "Cl\'ement Jamin" -, title = "Tangential complex" -, publisher = "{GUDHI Editorial Board}" -, booktitle = "{GUDHI} User and Reference Manual" -, url = "http://gudhi.gforge.inria.fr/doc/latest/group__tangential__complex.html" -, year = 2016 -}
\ No newline at end of file diff --git a/src/Bottleneck/concept/Persistence_diagram.h b/src/Bottleneck/concept/Persistence_diagram.h deleted file mode 100644 index eaaf8bc5..00000000 --- a/src/Bottleneck/concept/Persistence_diagram.h +++ /dev/null @@ -1,7 +0,0 @@ -typedef typename std::pair<double,double> Diagram_point; - -struct Persistence_Diagram -{ - const_iterator<Diagram_point> cbegin() const; - const_iterator<Diagram_point> cend() const; -}; diff --git a/src/Bottleneck/example/CMakeLists.txt b/src/Bottleneck/example/CMakeLists.txt deleted file mode 100644 index 77797202..00000000 --- a/src/Bottleneck/example/CMakeLists.txt +++ /dev/null @@ -1,5 +0,0 @@ -cmake_minimum_required(VERSION 2.6) -project(Bottleneck_examples) - -add_executable ( RandomDiagrams random_diagrams.cpp ) -add_test(RandomDiagrams ${CMAKE_CURRENT_BINARY_DIR}/RandomDiagrams) diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h deleted file mode 100644 index ea47e1d5..00000000 --- a/src/Bottleneck/include/gudhi/Graph_matching.h +++ /dev/null @@ -1,197 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ - -#include <deque> -#include <list> -#include <vector> - -#include "gudhi/Layered_neighbors_finder.h" - -namespace Gudhi { - -namespace bottleneck { - -template<typename Persistence_diagram1, typename Persistence_diagram2> -double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.); - -class Graph_matching { - public: - Graph_matching(const Persistence_diagrams_graph& g); - Graph_matching& operator=(const Graph_matching& m); - bool perfect() const; - bool multi_augment(); - void set_r(double r); - - private: - const Persistence_diagrams_graph& g; - double r; - std::vector<int> v_to_u; - std::list<int> unmatched_in_u; - - Layered_neighbors_finder* layering() const; - bool augment(Layered_neighbors_finder* layered_nf, int u_start_index, int max_depth); - void update(std::deque<int>& path); -}; - -Graph_matching::Graph_matching(const Persistence_diagrams_graph& g) - : g(g), r(0), v_to_u(g.size()), unmatched_in_u() { - for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) - unmatched_in_u.emplace_back(u_point_index); -} - -Graph_matching& Graph_matching::operator=(const Graph_matching& m) { - r = m.r; - v_to_u = m.v_to_u; - unmatched_in_u = m.unmatched_in_u; - return *this; -} - -inline bool Graph_matching::perfect() const { - return unmatched_in_u.empty(); -} - -inline bool Graph_matching::multi_augment() { - if (perfect()) - return false; - Layered_neighbors_finder* layered_nf = layering(); - double rn = sqrt(g.size()); - int nblmax = layered_nf->vlayers_number()*2 + 1; - // verification of a necessary criterion - if ((unmatched_in_u.size() > rn && nblmax > rn) || nblmax == 0) - return false; - bool successful = false; - std::list<int>* tries = new std::list<int>(unmatched_in_u); - for (auto it = tries->cbegin(); it != tries->cend(); it++) - successful = successful || augment(layered_nf, *it, nblmax); - delete tries; - delete layered_nf; - return successful; -} - -inline void Graph_matching::set_r(double r) { - this->r = r; -} - -Layered_neighbors_finder* Graph_matching::layering() const { - bool end = false; - int layer = 0; - std::list<int> u_vertices(unmatched_in_u); - std::list<int> v_vertices; - Neighbors_finder nf(g, r); - Layered_neighbors_finder* layered_nf = new Layered_neighbors_finder(g, r); - for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) - nf.add(v_point_index); - while (!u_vertices.empty()) { - for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) { - std::list<int>* u_succ = nf.pull_all_near(*it); - for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) { - layered_nf->add(*it, layer); - v_vertices.emplace_back(*it); - } - delete u_succ; - } - u_vertices.clear(); - for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) { - if (v_to_u.at(*it) == null_point_index()) - end = true; - else - u_vertices.emplace_back(v_to_u.at(*it)); - } - if (end) - return layered_nf; - v_vertices.clear(); - layer++; - } - return layered_nf; -} - -bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) { - std::deque<int> path; - path.emplace_back(u_start_index); - // start is a point from U - do { - if (static_cast<int>(path.size()) > max_depth) { - path.pop_back(); - path.pop_back(); - } - if (path.empty()) - return false; - int w = path.back(); - path.emplace_back(layered_nf->pull_near(w, path.size() / 2)); - while (path.back() == null_point_index()) { - path.pop_back(); - path.pop_back(); - if (path.empty()) - return false; - path.pop_back(); - path.emplace_back(layered_nf->pull_near(path.back(), path.size() / 2)); - } - path.emplace_back(v_to_u.at(path.back())); - } while (path.back() != null_point_index()); - path.pop_back(); - update(path); - return true; -} - -void Graph_matching::update(std::deque<int>& path) { - unmatched_in_u.remove(path.front()); - for (auto it = path.cbegin(); it != path.cend(); ++it) { - int tmp = *it; - ++it; - v_to_u[*it] = tmp; - } -} - -template<typename Persistence_diagram1, typename Persistence_diagram2> -double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) { - Persistence_diagrams_graph g(diag1, diag2, e); - std::vector<double>* sd = g.sorted_distances(); - int idmin = 0; - int idmax = sd->size() - 1; - double alpha = pow(sd->size(), 0.25); - Graph_matching m(g); - Graph_matching biggest_unperfect = m; - while (idmin != idmax) { - int pas = static_cast<int>((idmax - idmin) / alpha); - m.set_r(sd->at(idmin + pas)); - while (m.multi_augment()) {} - if (m.perfect()) { - idmax = idmin + pas; - m = biggest_unperfect; - } else { - biggest_unperfect = m; - idmin = idmin + pas + 1; - } - } - double b = sd->at(idmin); - delete sd; - return b; -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h deleted file mode 100644 index de36e00b..00000000 --- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h +++ /dev/null @@ -1,74 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ - -#include <vector> - -#include "Neighbors_finder.h" - -// Layered_neighbors_finder is a data structure used to find if a query point from U has neighbors in V in a given -// vlayer of the vlayered persistence diagrams graph. V's points have to be added manually using their index. -// A neighbor returned is automatically removed. - -namespace Gudhi { - -namespace bottleneck { - -class Layered_neighbors_finder { - public: - Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r); - void add(int v_point_index, int vlayer); - int pull_near(int u_point_index, int vlayer); - int vlayers_number() const; - - private: - const Persistence_diagrams_graph& g; - const double r; - std::vector<Neighbors_finder> neighbors_finder; -}; - -Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r), neighbors_finder() { } - -inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { - for (int l = neighbors_finder.size(); l <= vlayer; l++) - neighbors_finder.emplace_back(Neighbors_finder(g, r)); - neighbors_finder.at(vlayer).add(v_point_index); -} - -inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { - if (static_cast<int> (neighbors_finder.size()) <= vlayer) - return null_point_index(); - return neighbors_finder.at(vlayer).pull_near(u_point_index); -} - -inline int Layered_neighbors_finder::vlayers_number() const { - return neighbors_finder.size(); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h deleted file mode 100644 index 98256571..00000000 --- a/src/Bottleneck/include/gudhi/Neighbors_finder.h +++ /dev/null @@ -1,96 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ - -#include <unordered_set> -#include <list> - -#include "gudhi/Planar_neighbors_finder.h" - -namespace Gudhi { - -namespace bottleneck { - -// Neighbors_finder is a data structure used to find if a query point from U has neighbors in V -// in the persistence diagrams graph. -// V's points have to be added manually using their index. A neighbor returned is automatically removed. - -class Neighbors_finder { - public: - Neighbors_finder(const Persistence_diagrams_graph& g, double r); - void add(int v_point_index); - int pull_near(int u_point_index); - std::list<int>* pull_all_near(int u_point_index); - - private: - const Persistence_diagrams_graph& g; - const double r; - Planar_neighbors_finder planar_neighbors_f; - std::unordered_set<int> projections_f; -}; - -Neighbors_finder::Neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r), planar_neighbors_f(g, r), projections_f() { } - -inline void Neighbors_finder::add(int v_point_index) { - if (g.on_the_v_diagonal(v_point_index)) - projections_f.emplace(v_point_index); - else - planar_neighbors_f.add(v_point_index); -} - -inline int Neighbors_finder::pull_near(int u_point_index) { - int v_challenger = g.corresponding_point_in_v(u_point_index); - if (planar_neighbors_f.contains(v_challenger) && g.distance(u_point_index, v_challenger) < r) { - planar_neighbors_f.remove(v_challenger); - return v_challenger; - } - if (g.on_the_u_diagonal(u_point_index)) { - auto it = projections_f.cbegin(); - if (it != projections_f.cend()) { - int tmp = *it; - projections_f.erase(it); - return tmp; - } - } else { - return planar_neighbors_f.pull_near(u_point_index); - } - return null_point_index(); -} - -inline std::list<int>* Neighbors_finder::pull_all_near(int u_point_index) { - std::list<int>* all_pull = planar_neighbors_f.pull_all_near(u_point_index); - int last_pull = pull_near(u_point_index); - while (last_pull != null_point_index()) { - all_pull->emplace_back(last_pull); - last_pull = pull_near(u_point_index); - } - return all_pull; -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h deleted file mode 100644 index 73ad940b..00000000 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ /dev/null @@ -1,147 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ - -#include <vector> -#include <set> -#include <cmath> -#include <utility> // for pair<> -#include <algorithm> // for max - -namespace Gudhi { - -namespace bottleneck { - -// Diagram_point is the type of the persistence diagram's points -typedef std::pair<double, double> Diagram_point; - -// Return the used index for encoding none of the points -int null_point_index(); - -// Persistence_diagrams_graph is the interface beetwen any external representation of the two persistence diagrams and -// the bottleneck distance computation. An interface is necessary to ensure basic functions complexity. - -class Persistence_diagrams_graph { - public: - // Persistence_diagram1 and 2 are the types of any externals representations of persistence diagrams. - // They have to have an iterator over points, which have to have fields first (for birth) and second (for death). - template<typename Persistence_diagram1, typename Persistence_diagram2> - Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.); - Persistence_diagrams_graph(); - bool on_the_u_diagonal(int u_point_index) const; - bool on_the_v_diagonal(int v_point_index) const; - int corresponding_point_in_u(int v_point_index) const; - int corresponding_point_in_v(int u_point_index) const; - double distance(int u_point_index, int v_point_index) const; - int size() const; - std::vector<double>* sorted_distances(); - - private: - std::vector<Diagram_point> u; - std::vector<Diagram_point> v; - Diagram_point get_u_point(int u_point_index) const; - Diagram_point get_v_point(int v_point_index) const; -}; - -inline int null_point_index() { - return -1; -} - -template<typename Persistence_diagram1, typename Persistence_diagram2> -Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) - : u(), v() { - for (auto it = diag1.cbegin(); it != diag1.cend(); ++it) - if (it->second - it->first > e) - u.emplace_back(*it); - for (auto it = diag2.cbegin(); it != diag2.cend(); ++it) - if (it->second - it->first > e) - v.emplace_back(*it); - if (u.size() < v.size()) - swap(u, v); -} - -Persistence_diagrams_graph::Persistence_diagrams_graph() - : u(), v() { } - -inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const { - return u_point_index >= static_cast<int> (u.size()); -} - -inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const { - return v_point_index >= static_cast<int> (v.size()); -} - -inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const { - return on_the_v_diagonal(v_point_index) ? - v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size()); -} - -inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const { - return on_the_u_diagonal(u_point_index) ? - u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size()); -} - -inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) const { - // could be optimized for the case where one point is the projection of the other - if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) - return 0; - Diagram_point p_u = get_u_point(u_point_index); - Diagram_point p_v = get_v_point(v_point_index); - return (std::max)(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second)); -} - -inline int Persistence_diagrams_graph::size() const { - return static_cast<int> (u.size() + v.size()); -} - -inline std::vector<double>* Persistence_diagrams_graph::sorted_distances() { - // could be optimized - std::set<double> sorted_distances; - for (int u_point_index = 0; u_point_index < size(); ++u_point_index) - for (int v_point_index = 0; v_point_index < size(); ++v_point_index) - sorted_distances.emplace(distance(u_point_index, v_point_index)); - return new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend()); -} - -inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const { - if (!on_the_u_diagonal(u_point_index)) - return u.at(u_point_index); - Diagram_point projector = v.at(corresponding_point_in_v(u_point_index)); - double x = (projector.first + projector.second) / 2; - return Diagram_point(x, x); -} - -inline Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const { - if (!on_the_v_diagonal(v_point_index)) - return v.at(v_point_index); - Diagram_point projector = u.at(corresponding_point_in_u(v_point_index)); - double x = (projector.first + projector.second) / 2; - return Diagram_point(x, x); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h deleted file mode 100644 index 4af672e4..00000000 --- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h +++ /dev/null @@ -1,119 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ - -#include <list> -#include <iostream> -#include <set> - -#include "Persistence_diagrams_graph.h" - -namespace Gudhi { - -namespace bottleneck { - -// Planar_neighbors_finder is a data structure used to find if a query point from U has planar neighbors in V with the -// planar distance. -// V's points have to be added manually using their index. A neighbor returned is automatically removed but we can also -// remove points manually using their index. - -class Abstract_planar_neighbors_finder { - public: - Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r); - virtual ~Abstract_planar_neighbors_finder() = 0; - virtual void add(int v_point_index) = 0; - virtual void remove(int v_point_index) = 0; - virtual bool contains(int v_point_index) const = 0; - virtual int pull_near(int u_point_index) = 0; - virtual std::list<int>* pull_all_near(int u_point_index); - - protected: - const Persistence_diagrams_graph& g; - const double r; -}; - - -// Naive_pnf is a nave implementation of Abstract_planar_neighbors_finder - -class Naive_pnf : public Abstract_planar_neighbors_finder { - public: - Naive_pnf(const Persistence_diagrams_graph& g, double r); - void add(int v_point_index); - void remove(int v_point_index); - bool contains(int v_point_index) const; - int pull_near(int u_point_index); - - private: - std::set<int> candidates; -}; - - -// Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation -typedef Naive_pnf Planar_neighbors_finder; - -Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r) { } - -inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() { } - -inline std::list<int>* Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { - std::list<int>* all_pull = new std::list<int>(); - int last_pull = pull_near(u_point_index); - while (last_pull != null_point_index()) { - all_pull->emplace_back(last_pull); - last_pull = pull_near(u_point_index); - } - return all_pull; -} - -Naive_pnf::Naive_pnf(const Persistence_diagrams_graph& g, double r) : - Abstract_planar_neighbors_finder(g, r), candidates() { } - -inline void Naive_pnf::add(int v_point_index) { - candidates.emplace(v_point_index); -} - -inline void Naive_pnf::remove(int v_point_index) { - candidates.erase(v_point_index); -} - -inline bool Naive_pnf::contains(int v_point_index) const { - return (candidates.count(v_point_index) > 0); -} - -inline int Naive_pnf::pull_near(int u_point_index) { - for (auto it = candidates.begin(); it != candidates.end(); ++it) - if (g.distance(u_point_index, *it) <= r) { - int tmp = *it; - candidates.erase(it); - return tmp; - } - return null_point_index(); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/test/CMakeLists.txt b/src/Bottleneck/test/CMakeLists.txt deleted file mode 100644 index 9d88ab25..00000000 --- a/src/Bottleneck/test/CMakeLists.txt +++ /dev/null @@ -1,21 +0,0 @@ -cmake_minimum_required(VERSION 2.6) -project(Bottleneck_tests) - -if (GCOVR_PATH) - # for gcovr to make coverage reports - Corbera Jenkins plugin - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") -endif() -if (GPROF_PATH) - # for gprof to make coverage reports - Jenkins - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") -endif() - -add_executable ( BottleneckUT bottleneck_unit_test.cpp ) -target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) - -# Unitary tests -add_test(NAME BottleneckUT - COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUT - # XML format for Jenkins xUnit plugin - --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BottleneckUT.xml --log_level=test_suite --report_level=no) - diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp deleted file mode 100644 index c60f5d8a..00000000 --- a/src/Bottleneck/test/bottleneck_unit_test.cpp +++ /dev/null @@ -1,26 +0,0 @@ -#define BOOST_TEST_DYN_LINK -#define BOOST_TEST_MODULE "bottleneck" -#include <boost/test/unit_test.hpp> - -#include "gudhi/Graph_matching.h" -#include <iostream> - -using namespace Gudhi::bottleneck; - -BOOST_AUTO_TEST_CASE(random_diagrams) { - int n = 100; - // Random construction - std::vector< std::pair<double, double> > v1, v2; - for (int i = 0; i < n; i++) { - int a = rand() % n; - v1.emplace_back(a, a + rand() % (n - a)); - int b = rand() % n; - v2.emplace_back(b, b + rand() % (n - b)); - } - // v1 and v2 are persistence diagrams containing each 100 randoms points. - double b = bottleneck_distance(v1, v2, 0); - // - std::cout << b << std::endl; - const double EXPECTED_DISTANCE = 98.5; - BOOST_CHECK(b == EXPECTED_DISTANCE); -} diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h new file mode 100644 index 00000000..e4766628 --- /dev/null +++ b/src/Bottleneck_distance/concept/Persistence_diagram.h @@ -0,0 +1,49 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi + * + * Copyright (C) 2015 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ +#define CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ + +namespace Gudhi { + +namespace bottleneck_distance { + +/** \brief Concept of Diagram_point. std::get<0>(point) must return the birth of the corresponding component and std::get<1>(point) its death. + * A valid implementation of this concept is std::pair<double,double>. + * Birth must be 0 or positive, death should be larger than birth, death can be std::numeric_limits<double>::infinity() for components which stay alive. + * + * \ingroup bottleneck_distance + */ +struct Diagram_point; + +/** \brief Concept of persistence diagram. It's a range of Diagram_point. + * std::begin(diagram) and std::end(diagram) must return corresponding iterators. + * + * \ingroup bottleneck_distance + */ +struct Persistence_Diagram; + +} // namespace bottleneck_distance + +} // namespace Gudhi + +#endif // CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h new file mode 100644 index 00000000..1f443f7f --- /dev/null +++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h @@ -0,0 +1,51 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef DOC_BOTTLENECK_DISTANCE_H_ +#define DOC_BOTTLENECK_DISTANCE_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +// needs namespace for Doxygen to link on classes +namespace bottleneck_distance { + +/** \defgroup bottleneck_distance Bottleneck distance + * + * \author François Godi + * @{ + * + * \section bottleneckdefinition Definition + * + * Bottleneck distance measures the similarity between two persistence diagrams. + * It's the shortest distance b for which there exists a perfect matching between + * the points of the two diagrams (and also all the points of the diagonal) such that + * any couple of matched points are at distance at most b. + * + * \image html perturb_pd.png Bottleneck distance is the length of the longest edge on this picture. + * + */ +/** @} */ // end defgroup bottleneck_distance + +} // namespace bottleneck_distance + +} // namespace Gudhi + diff --git a/src/Bottleneck_distance/doc/perturb_pd.png b/src/Bottleneck_distance/doc/perturb_pd.png Binary files differnew file mode 100644 index 00000000..be638de0 --- /dev/null +++ b/src/Bottleneck_distance/doc/perturb_pd.png diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt new file mode 100644 index 00000000..cd53ccfc --- /dev/null +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -0,0 +1,13 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_examples) + +# requires CGAL 4.8 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.8.0) + if (EIGEN3_FOUND) + add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) + add_executable (bottleneck_basic_example bottleneck_basic_example.cpp) + endif() + endif () +endif() diff --git a/src/Bottleneck/example/random_diagrams.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp index 71f152a6..78e00e57 100644 --- a/src/Bottleneck/example/random_diagrams.cpp +++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Authors: Francois Godi, small modifications by Pawel Dlotko * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -20,21 +20,29 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "gudhi/Graph_matching.h" +#include <gudhi/Bottleneck.h> #include <iostream> -using namespace Gudhi::bottleneck; - int main() { - int n = 100; - std::vector< std::pair<double, double> > v1, v2; - for (int i = 0; i < n; i++) { - int a = rand() % n; - v1.emplace_back(a, a + rand() % (n - a)); - int b = rand() % n; - v2.emplace_back(b, b + rand() % (n - b)); - } - // v1 and v2 are persistence diagrams containing each 100 randoms points. - double b = bottleneck_distance(v1, v2, 0); - std::cout << b << std::endl; + + std::vector< std::pair<double, double> > v1, v2; + + v1.emplace_back(2.7, 3.7); + v1.emplace_back(9.6, 14.); + v1.emplace_back(34.2, 34.974); + v1.emplace_back(3., std::numeric_limits<double>::infinity()); + + v2.emplace_back(2.8, 4.45); + v2.emplace_back(9.5, 14.1); + v2.emplace_back(3.2, std::numeric_limits<double>::infinity()); + + + double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); + + std::cout << "Bottleneck distance = " << b << std::endl; + + b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); + + std::cout << "Approx bottleneck distance = " << b << std::endl; + } diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp new file mode 100644 index 00000000..ceedccc5 --- /dev/null +++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp @@ -0,0 +1,77 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#define CGAL_HAS_THREADS + +#include <gudhi/Bottleneck.h> +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> + +std::vector< std::pair<double, double> > read_diagram_from_file( const char* filename ) +{ + std::ifstream in; + in.open( filename ); + std::vector< std::pair<double, double> > result; + if ( !in.is_open() ) + { + std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; + throw "File do not exist \n"; + } + + std::string line; + while (!in.eof()) + { + getline(in,line); + if ( line.length() != 0 ) + { + std::stringstream lineSS; + lineSS << line; + double beginn, endd; + lineSS >> beginn; + lineSS >> endd; + result.push_back( std::make_pair( beginn , endd ) ); + } + } + in.close(); + return result; +} //read_diagram_from_file + +int main( int argc , char** argv ) +{ + if ( argc < 3 ) + { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " << + "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << + " distance (set by default to zero). The program will now terminate \n"; + } + std::vector< std::pair< double , double > > diag1 = read_diagram_from_file( argv[1] ); + std::vector< std::pair< double , double > > diag2 = read_diagram_from_file( argv[2] ); + double tolerance = 0.; + if ( argc == 4 ) + { + tolerance = atof( argv[3] ); + } + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); + std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; +} diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h new file mode 100644 index 00000000..24a31ac1 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -0,0 +1,75 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef BOTTLENECK_H_ +#define BOTTLENECK_H_ + +#include <gudhi/Graph_matching.h> +#include <cmath> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams (see Concepts). + * If the last parameter e is not 0 (default value if not explicited), you get an additive e-approximation. + * + * \ingroup bottleneck_distance + */ +template<typename Persistence_diagram1, typename Persistence_diagram2> +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) { + Persistence_graph g(diag1, diag2, e); + double b = g.bottleneck_alive(); + if(b == std::numeric_limits<double>::infinity()) + return std::numeric_limits<double>::infinity(); + std::vector<double> sd; + if(e == 0.) + sd = g.sorted_distances(); + long idmin = 0; + long idmax = e==0. ? sd.size() - 1 : g.diameter_bound()/e + 1; + double alpha = std::pow(g.size(), 1./5.); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); + while (idmin != idmax) { + long step = static_cast<long>((idmax - idmin - 1)/alpha); + m.set_r(e == 0. ? sd.at(idmin + step) : e*(idmin + step)); + while (m.multi_augment()); + //The above while compute a maximum matching (according to the r setted before) + if (m.perfect()) { + idmax = idmin + step; + m = biggest_unperfect; + } else { + biggest_unperfect = m; + idmin = idmin + step + 1; + } + } + b = std::max(b, e == 0. ? sd.at(idmin) : e*(idmin)); + return b; +} + + + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // BOTTLENECK_H_ diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h new file mode 100644 index 00000000..f085b0da --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h @@ -0,0 +1,582 @@ +// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>), +// : Waqar Khan <wkhan@mpi-inf.mpg.de> + +#ifndef CGAL_KD_TREE_H +#define CGAL_KD_TREE_H + +#include "Kd_tree_node.h" + +#include <CGAL/basic.h> +#include <CGAL/assertions.h> +#include <vector> + +#include <CGAL/algorithm.h> +#include <CGAL/internal/Get_dimension_tag.h> +#include <CGAL/Search_traits.h> + + +#include <deque> +#include <boost/container/deque.hpp> +#include <boost/optional.hpp> + +#ifdef CGAL_HAS_THREADS +#include <CGAL/mutex.h> +#endif + +namespace CGAL { + +//template <class SearchTraits, class Splitter_=Median_of_rectangle<SearchTraits>, class UseExtendedNode = Tag_true > +template <class SearchTraits, class Splitter_=Sliding_midpoint<SearchTraits>, class UseExtendedNode = Tag_true > +class Kd_tree { + +public: + typedef SearchTraits Traits; + typedef Splitter_ Splitter; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Splitter::Container Point_container; + + typedef typename SearchTraits::FT FT; + typedef Kd_tree_node<SearchTraits, Splitter, UseExtendedNode > Node; + typedef Kd_tree_leaf_node<SearchTraits, Splitter, UseExtendedNode > Leaf_node; + typedef Kd_tree_internal_node<SearchTraits, Splitter, UseExtendedNode > Internal_node; + typedef Kd_tree<SearchTraits, Splitter> Tree; + typedef Kd_tree<SearchTraits, Splitter,UseExtendedNode> Self; + + typedef Node* Node_handle; + typedef const Node* Node_const_handle; + typedef Leaf_node* Leaf_node_handle; + typedef const Leaf_node* Leaf_node_const_handle; + typedef Internal_node* Internal_node_handle; + typedef const Internal_node* Internal_node_const_handle; + typedef typename std::vector<const Point_d*>::const_iterator Point_d_iterator; + typedef typename std::vector<const Point_d*>::const_iterator Point_d_const_iterator; + typedef typename Splitter::Separator Separator; + typedef typename std::vector<Point_d>::const_iterator iterator; + typedef typename std::vector<Point_d>::const_iterator const_iterator; + + typedef typename std::vector<Point_d>::size_type size_type; + + typedef typename internal::Get_dimension_tag<SearchTraits>::Dimension D; + +private: + SearchTraits traits_; + Splitter split; + + + // wokaround for https://svn.boost.org/trac/boost/ticket/9332 +#if (_MSC_VER == 1800) && (BOOST_VERSION == 105500) + std::deque<Internal_node> internal_nodes; + std::deque<Leaf_node> leaf_nodes; +#else + boost::container::deque<Internal_node> internal_nodes; + boost::container::deque<Leaf_node> leaf_nodes; +#endif + + Node_handle tree_root; + + Kd_tree_rectangle<FT,D>* bbox; + std::vector<Point_d> pts; + + // Instead of storing the points in arrays in the Kd_tree_node + // we put all the data in a vector in the Kd_tree. + // and we only store an iterator range in the Kd_tree_node. + // + std::vector<const Point_d*> data; + + + #ifdef CGAL_HAS_THREADS + mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build() + #endif + bool built_; + bool removed_; + + // protected copy constructor + Kd_tree(const Tree& tree) + : traits_(tree.traits_),built_(tree.built_) + {}; + + + // Instead of the recursive construction of the tree in the class Kd_tree_node + // we do this in the tree class. The advantage is that we then can optimize + // the allocation of the nodes. + + // The leaf node + Node_handle + create_leaf_node(Point_container& c) + { + Leaf_node node(true , static_cast<unsigned int>(c.size())); + std::ptrdiff_t tmp = c.begin() - data.begin(); + node.data = pts.begin() + tmp; + + leaf_nodes.push_back(node); + Leaf_node_handle nh = &leaf_nodes.back(); + + + return nh; + } + + + // The internal node + + Node_handle + create_internal_node(Point_container& c, const Tag_true&) + { + return create_internal_node_use_extension(c); + } + + Node_handle + create_internal_node(Point_container& c, const Tag_false&) + { + return create_internal_node(c); + } + + + + // TODO: Similiar to the leaf_init function above, a part of the code should be + // moved to a the class Kd_tree_node. + // It is not proper yet, but the goal was to see if there is + // a potential performance gain through the Compact_container + Node_handle + create_internal_node_use_extension(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + + Separator sep; + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + int cd = nh->cutting_dimension(); + if(!c_low.empty()){ + nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd); + nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd); + } + else{ + nh->lower_low_val = nh->cutting_value(); + nh->lower_high_val = nh->cutting_value(); + } + if(!c.empty()){ + nh->upper_low_val = c.tight_bounding_box().min_coord(cd); + nh->upper_high_val = c.tight_bounding_box().max_coord(cd); + } + else{ + nh->upper_low_val = nh->cutting_value(); + nh->upper_high_val = nh->cutting_value(); + } + + CGAL_assertion(nh->cutting_value() >= nh->lower_low_val); + CGAL_assertion(nh->cutting_value() <= nh->upper_high_val); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node_use_extension(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node_use_extension(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + + return nh; + } + + + // Note also that I duplicated the code to get rid if the if's for + // the boolean use_extension which was constant over the construction + Node_handle + create_internal_node(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + Separator sep; + + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + return nh; + } + + + +public: + + Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + {} + + template <class InputIterator> + Kd_tree(InputIterator first, InputIterator beyond, + Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + { + pts.insert(pts.end(), first, beyond); + } + + bool empty() const { + return pts.empty(); + } + + void + build() + { + // This function is not ready to be called when a tree already exists, one + // must call invalidate_built() first. + CGAL_assertion(!is_built()); + CGAL_assertion(!removed_); + const Point_d& p = *pts.begin(); + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object(); + int dim = static_cast<int>(std::distance(ccci(p), ccci(p,0))); + + data.reserve(pts.size()); + for(unsigned int i = 0; i < pts.size(); i++){ + data.push_back(&pts[i]); + } + Point_container c(dim, data.begin(), data.end(),traits_); + bbox = new Kd_tree_rectangle<FT,D>(c.bounding_box()); + if (c.size() <= split.bucket_size()){ + tree_root = create_leaf_node(c); + }else { + tree_root = create_internal_node(c, UseExtendedNode()); + } + + //Reorder vector for spatial locality + std::vector<Point_d> ptstmp; + ptstmp.resize(pts.size()); + for (std::size_t i = 0; i < pts.size(); ++i){ + ptstmp[i] = *data[i]; + } + for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ + std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); + leaf_nodes[i].data = ptstmp.begin() + tmp; + } + pts.swap(ptstmp); + + data.clear(); + + built_ = true; + } + +private: + //any call to this function is for the moment not threadsafe + void const_build() const { + #ifdef CGAL_HAS_THREADS + //this ensure that build() will be called once + CGAL_SCOPED_LOCK(building_mutex); + if(!is_built()) + #endif + const_cast<Self*>(this)->build(); //THIS IS NOT THREADSAFE + } +public: + + bool is_built() const + { + return built_; + } + + void invalidate_built() + { + if(removed_){ + // Walk the tree to collect the remaining points. + // Writing directly to pts would likely work, but better be safe. + std::vector<Point_d> ptstmp; + //ptstmp.resize(root()->num_items()); + root()->tree_items(std::back_inserter(ptstmp)); + pts.swap(ptstmp); + removed_=false; + CGAL_assertion(is_built()); // the rest of the cleanup must happen + } + if(is_built()){ + internal_nodes.clear(); + leaf_nodes.clear(); + data.clear(); + delete bbox; + built_ = false; + } + } + + void clear() + { + invalidate_built(); + pts.clear(); + removed_ = false; + } + + void + insert(const Point_d& p) + { + invalidate_built(); + pts.push_back(p); + } + + template <class InputIterator> + void + insert(InputIterator first, InputIterator beyond) + { + invalidate_built(); + pts.insert(pts.end(),first, beyond); + } + +private: + struct Equal_by_coordinates { + SearchTraits const* traits; + Point_d const* pp; + bool operator()(Point_d const&q) const { + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object(); + return std::equal(ccci(*pp), ccci(*pp,0), ccci(q)); + } + }; + Equal_by_coordinates equal_by_coordinates(Point_d const&p){ + Equal_by_coordinates ret = { &traits(), &p }; + return ret; + } + +public: + void + remove(const Point_d& p) + { + remove(p, equal_by_coordinates(p)); + } + + template<class Equal> + void + remove(const Point_d& p, Equal const& equal_to_p) + { +#if 0 + // This code could have quadratic runtime. + if (!is_built()) { + std::vector<Point_d>::iterator pi = std::find(pts.begin(), pts.end(), p); + // Precondition: the point must be there. + CGAL_assertion (pi != pts.end()); + pts.erase(pi); + return; + } +#endif + bool success = remove_(p, 0, false, 0, false, root(), equal_to_p); + CGAL_assertion(success); + + // Do not set the flag is the tree has been cleared. + if(is_built()) + removed_ |= success; + } +private: + template<class Equal> + bool remove_(const Point_d& p, + Internal_node_handle grandparent, bool parent_islower, + Internal_node_handle parent, bool islower, + Node_handle node, Equal const& equal_to_p) { + // Recurse to locate the point + if (!node->is_leaf()) { + Internal_node_handle newparent = static_cast<Internal_node_handle>(node); + // FIXME: This should be if(x<y) remove low; else remove up; + if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] <= newparent->cutting_value()) { + if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p)) + return true; + } + //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value()) + return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p); + + CGAL_assertion(false); // Point was not found + } + + // Actual removal + Leaf_node_handle lnode = static_cast<Leaf_node_handle>(node); + if (lnode->size() > 1) { + iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p); + // FIXME: we should ensure this never happens + if (pi == lnode->end()) return false; + iterator lasti = lnode->end() - 1; + if (pi != lasti) { + // Hack to get a non-const iterator + std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin())); + } + lnode->drop_last_point(); + } else if (!equal_to_p(*lnode->begin())) { + // FIXME: we should ensure this never happens + return false; + } else if (grandparent) { + Node_handle brother = islower ? parent->upper() : parent->lower(); + if (parent_islower) + grandparent->set_lower(brother); + else + grandparent->set_upper(brother); + } else if (parent) { + tree_root = islower ? parent->upper() : parent->lower(); + } else { + clear(); + } + return true; + } + +public: + //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known). + void reserve(size_t size) + { + pts.reserve(size); + } + + //Get the capacity of the underlying points vector. + size_t capacity() + { + return pts.capacity(); + } + + + template <class OutputIterator, class FuzzyQueryItem> + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle<FT,D> b(*bbox); + return tree_root->search(it,q,b); + } + return it; + } + + + template <class FuzzyQueryItem> + boost::optional<Point_d> + search_any_point(const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle<FT,D> b(*bbox); + return tree_root->search_any_point(q,b); + } + return boost::none; + } + + + ~Kd_tree() { + if(is_built()){ + delete bbox; + } + } + + + const SearchTraits& + traits() const + { + return traits_; + } + + Node_const_handle + root() const + { + if(! is_built()){ + const_build(); + } + return tree_root; + } + + Node_handle + root() + { + if(! is_built()){ + build(); + } + return tree_root; + } + + void + print() const + { + if(! is_built()){ + const_build(); + } + root()->print(); + } + + const Kd_tree_rectangle<FT,D>& + bounding_box() const + { + if(! is_built()){ + const_build(); + } + return *bbox; + } + + const_iterator + begin() const + { + return pts.begin(); + } + + const_iterator + end() const + { + return pts.end(); + } + + size_type + size() const + { + return pts.size(); + } + + // Print statistics of the tree. + std::ostream& + statistics(std::ostream& s) const + { + if(! is_built()){ + const_build(); + } + s << "Tree statistics:" << std::endl; + s << "Number of items stored: " + << root()->num_items() << std::endl; + s << "Number of nodes: " + << root()->num_nodes() << std::endl; + s << " Tree depth: " << root()->depth() << std::endl; + return s; + } + + +}; + +} // namespace CGAL + +#endif // CGAL_KD_TREE_H diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h new file mode 100644 index 00000000..909ee260 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h @@ -0,0 +1,586 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Authors : Hans Tangelder (<hanst@cs.uu.nl>) + +#ifndef CGAL_KD_TREE_NODE_H +#define CGAL_KD_TREE_NODE_H + +#include "CGAL/Splitters.h" + +#include <CGAL/Compact_container.h> +#include <boost/cstdint.hpp> + +namespace CGAL { + + template <class SearchTraits, class Splitter, class UseExtendedNode> + class Kd_tree; + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_node { + + friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>; + + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_handle Internal_node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_const_handle Internal_node_const_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_handle Leaf_node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_const_handle Leaf_node_const_handle; + typedef typename TreeTraits::Point_d Point_d; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Point_d_iterator Point_d_iterator; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::D D; + + bool leaf; + + public : + Kd_tree_node(bool leaf_) + :leaf(leaf_){} + + bool is_leaf() const{ + return leaf; + } + + std::size_t + num_items() const + { + if (is_leaf()){ + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + return node->size(); + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + return node->lower()->num_items() + node->upper()->num_items(); + } + } + + std::size_t + num_nodes() const + { + if (is_leaf()) return 1; + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + return node->lower()->num_nodes() + node->upper()->num_nodes(); + } + } + + int + depth(const int current_max_depth) const + { + if (is_leaf()){ + return current_max_depth; + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + return + (std::max)( node->lower()->depth(current_max_depth + 1), + node->upper()->depth(current_max_depth + 1)); + } + } + + int + depth() const + { + return depth(1); + } + + template <class OutputIterator> + OutputIterator + tree_items(OutputIterator it) const { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {*it=*i; ++it;} + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + it=node->lower()->tree_items(it); + it=node->upper()->tree_items(it); + } + return it; + } + + + boost::optional<Point_d> + any_tree_item() const { + boost::optional<Point_d> result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0){ + return boost::make_optional(*(node->begin())); + } + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + result = node->lower()->any_tree_item(); + if(! result){ + result = node->upper()->any_tree_item(); + } + } + return result; + } + + + void + indent(int d) const + { + for(int i = 0; i < d; i++){ + std::cout << " "; + } + } + + + void + print(int d = 0) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + indent(d); + std::cout << "leaf" << std::endl; + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {indent(d);std::cout << *i << std::endl;} + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + indent(d); + std::cout << "lower tree" << std::endl; + node->lower()->print(d+1); + indent(d); + std::cout << "separator: dim = " << node->cutting_dimension() << " val = " << node->cutting_value() << std::endl; + indent(d); + std::cout << "upper tree" << std::endl; + node->upper()->print(d+1); + } + } + + + template <class OutputIterator, class FuzzyQueryItem> + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q, + Kd_tree_rectangle<FT,D>& b) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + {*it++=*i;} + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle<FT,D> b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)) + it=node->lower()->tree_items(it); + else + if (q.inner_range_intersects(b)) + it=node->lower()->search(it,q,b); + if (q.outer_range_contains(b_upper)) + it=node->upper()->tree_items(it); + else + if (q.inner_range_intersects(b_upper)) + it=node->upper()->search(it,q,b_upper); + }; + return it; + } + + + template <class FuzzyQueryItem> + boost::optional<Point_d> + search_any_point(const FuzzyQueryItem& q, + Kd_tree_rectangle<FT,D>& b) const + { + boost::optional<Point_d> result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + { result = *i; break; } + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle<FT,D> b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)){ + result = node->lower()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b)){ + result = node->lower()->search_any_point(q,b); + } + } + if(result){ + return result; + } + if (q.outer_range_contains(b_upper)){ + result = node->upper()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b_upper)) + result = node->upper()->search_any_point(q,b_upper); + } + } + return result; + } + + }; + + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>; + + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator; + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename TreeTraits::Point_d Point_d; + + private: + + // private variables for leaf nodes + boost::int32_t n; // denotes number of items in a leaf node + iterator data; // iterator to data in leaf node + + + public: + + // default constructor + Kd_tree_leaf_node() + {} + + Kd_tree_leaf_node(bool leaf_ ) + : Base(leaf_) + {} + + Kd_tree_leaf_node(bool leaf_,unsigned int n_ ) + : Base(leaf_), n(n_) + {} + + // members for all nodes + + // members for leaf nodes only + inline + unsigned int + size() const + { + return n; + } + + inline + iterator + begin() const + { + return data; + } + + inline + iterator + end() const + { + return data + n; + } + + inline + void + drop_last_point() + { + --n; + } + + }; //leaf node + + + + template < class TreeTraits, class Splitter, class UseExtendedNode> + class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>; + + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator; + + private: + + // private variables for internal nodes + boost::int32_t cut_dim; + FT cut_val; + Node_handle lower_ch, upper_ch; + + + // private variables for extended internal nodes + FT upper_low_val; + FT upper_high_val; + FT lower_low_val; + FT lower_high_val; + + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + // members for extended internal node only + inline + FT + upper_low_value() const + { + return upper_low_val; + } + + inline + FT + upper_high_value() const + { + return upper_high_val; + } + + inline + FT + lower_low_value() const + { + return lower_low_val; + } + + inline + FT + lower_high_value() const + { + return lower_high_val; + } + + /*Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + template < class TreeTraits, class Splitter> + class Kd_tree_internal_node<TreeTraits,Splitter,Tag_false> : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ + + friend class Kd_tree<TreeTraits,Splitter,Tag_false>; + + typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; + typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_handle Node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Separator Separator; + + private: + + // private variables for internal nodes + boost::uint8_t cut_dim; + FT cut_val; + + Node_handle lower_ch, upper_ch; + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + /* Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + + +} // namespace CGAL +#endif // CGAL_KDTREE_NODE_H diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h b/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h new file mode 100644 index 00000000..dbe707ed --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h @@ -0,0 +1,621 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>) + +#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH +#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH + +#include "CGAL/Kd_tree.h" + +#include <cstring> +#include <list> +#include <queue> +#include <memory> +#include <CGAL/Euclidean_distance.h> +#include <CGAL/tuple.h> + +namespace CGAL { + + template <class SearchTraits, + class Distance_= typename internal::Spatial_searching_default_distance<SearchTraits>::type, + class Splitter_ = Sliding_midpoint<SearchTraits>, + class Tree_= Kd_tree<SearchTraits, Splitter_, Tag_true> > + class Orthogonal_incremental_neighbor_search { + + public: + typedef Splitter_ Splitter; + typedef Tree_ Tree; + typedef Distance_ Distance; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Distance::Query_item Query_item; + typedef typename SearchTraits::FT FT; + typedef typename Tree::Point_d_iterator Point_d_iterator; + typedef typename Tree::Node_const_handle Node_const_handle; + + typedef std::pair<Point_d,FT> Point_with_transformed_distance; + typedef CGAL::cpp11::tuple<Node_const_handle,FT,std::vector<FT> > Node_with_distance; + typedef std::vector<Node_with_distance*> Node_with_distance_vector; + typedef std::vector<Point_with_transformed_distance*> Point_with_transformed_distance_vector; + + template<class T> + struct Object_wrapper + { + T object; + Object_wrapper(const T& t):object(t){} + const T& operator* () const { return object; } + const T* operator-> () const { return &object; } + }; + + class Iterator_implementation { + SearchTraits traits; + public: + + int number_of_neighbours_computed; + int number_of_internal_nodes_visited; + int number_of_leaf_nodes_visited; + int number_of_items_visited; + + private: + + typedef std::vector<FT> Distance_vector; + + Distance_vector dists; + + Distance Orthogonal_distance_instance; + + FT multiplication_factor; + + Query_item query_point; + + FT distance_to_root; + + bool search_nearest_neighbour; + + FT rd; + + + class Priority_higher { + public: + + bool search_nearest; + + Priority_higher(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool + operator() (Node_with_distance* n1, Node_with_distance* n2) const + { + return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1)); + } + }; + + class Distance_smaller { + + public: + + bool search_nearest; + + Distance_smaller(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const + { + return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second); + } + }; + + + std::priority_queue<Node_with_distance*, Node_with_distance_vector, + Priority_higher> PriorityQueue; + + public: + std::priority_queue<Point_with_transformed_distance*, Point_with_transformed_distance_vector, + Distance_smaller> Item_PriorityQueue; + + + public: + + int reference_count; + + + + // constructor + Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr, + FT Eps=FT(0.0), bool search_nearest=true) + : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0), + number_of_leaf_nodes_visited(0), number_of_items_visited(0), + Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), + query_point(q), search_nearest_neighbour(search_nearest), + PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)), + reference_count(1) + + + { + if (tree.empty()) return; + + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object(); + int dim = static_cast<int>(std::distance(ccci(q), ccci(q,0))); + + dists.resize(dim); + for(int i=0 ; i<dim ; ++i){ + dists[i] = 0; + } + + if (search_nearest){ + distance_to_root= + Orthogonal_distance_instance.min_distance_to_rectangle(q, tree.bounding_box(),dists); + Node_with_distance *The_Root = new Node_with_distance(tree.root(), + distance_to_root, dists); + PriorityQueue.push(The_Root); + + // rd is the distance of the top of the priority queue to q + rd=CGAL::cpp11::get<1>(*The_Root); + Compute_the_next_nearest_neighbour(); + } + else{ + distance_to_root= + Orthogonal_distance_instance.max_distance_to_rectangle(q, + tree.bounding_box(), dists); + Node_with_distance *The_Root = new Node_with_distance(tree.root(), + distance_to_root, dists); + PriorityQueue.push(The_Root); + + // rd is the distance of the top of the priority queue to q + rd=CGAL::cpp11::get<1>(*The_Root); + Compute_the_next_furthest_neighbour(); + } + + + } + + // * operator + const Point_with_transformed_distance& + operator* () const + { + return *(Item_PriorityQueue.top()); + } + + // prefix operator + Iterator_implementation& + operator++() + { + Delete_the_current_item_top(); + if(search_nearest_neighbour) + Compute_the_next_nearest_neighbour(); + else + Compute_the_next_furthest_neighbour(); + return *this; + } + + // postfix operator + Object_wrapper<Point_with_transformed_distance> + operator++(int) + { + Object_wrapper<Point_with_transformed_distance> result( *(Item_PriorityQueue.top()) ); + ++*this; + return result; + } + + // Print statistics of the general priority search process. + std::ostream& + statistics (std::ostream& s) const { + s << "Orthogonal priority search statistics:" + << std::endl; + s << "Number of internal nodes visited:" + << number_of_internal_nodes_visited << std::endl; + s << "Number of leaf nodes visited:" + << number_of_leaf_nodes_visited << std::endl; + s << "Number of items visited:" + << number_of_items_visited << std::endl; + s << "Number of neighbours computed:" + << number_of_neighbours_computed << std::endl; + return s; + } + + + //destructor + ~Iterator_implementation() + { + while (!PriorityQueue.empty()) { + Node_with_distance* The_top=PriorityQueue.top(); + PriorityQueue.pop(); + delete The_top; + } + while (!Item_PriorityQueue.empty()) { + Point_with_transformed_distance* The_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_top; + } + } + + private: + + void + Delete_the_current_item_top() + { + Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_item_top; + } + + void + Compute_the_next_nearest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (multiplication_factor*rd > Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast<typename Tree::Internal_node_const_handle>(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff1; + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), new_rd, dists); + PriorityQueue.push(Upper_Child); + dists[new_cut_dim] = dst; + N=node->lower(); + + } + else { // compute new distance + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff2; + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), new_rd, dists); + PriorityQueue.push(Lower_Child); + dists[new_cut_dim] = dst; + N=node->upper(); + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast<typename Tree::Leaf_node_const_handle>(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd > + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + + + void + Compute_the_next_furthest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (rd < multiplication_factor*Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast<typename Tree::Internal_node_const_handle>(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + diff1 = val - node->upper_high_value(); + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), copy_rd, dists); + PriorityQueue.push(Lower_Child); + N=node->upper(); + dists[new_cut_dim] = diff1; + copy_rd=new_rd; + + } + else { // compute new distance + diff2 = val - node->lower_low_value(); + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), copy_rd, dists); + PriorityQueue.push(Upper_Child); + N=node->lower(); + dists[new_cut_dim] = diff2; + copy_rd=new_rd; + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast<typename Tree::Leaf_node_const_handle>(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd < + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + }; // class Iterator_implementaion + + + + + + + + + + public: + class iterator; + typedef iterator const_iterator; + + // constructor + Orthogonal_incremental_neighbor_search(const Tree& tree, + const Query_item& q, FT Eps = FT(0.0), + bool search_nearest=true, const Distance& tr=Distance()) + : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest) + {} + + iterator + begin() const + { + return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest); + } + + iterator + end() const + { + return iterator(); + } + + std::ostream& + statistics(std::ostream& s) + { + begin()->statistics(s); + return s; + } + + + + + class iterator { + + public: + + typedef std::input_iterator_tag iterator_category; + typedef Point_with_transformed_distance value_type; + typedef Point_with_transformed_distance* pointer; + typedef const Point_with_transformed_distance& reference; + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + typedef int distance_type; + + //class Iterator_implementation; + Iterator_implementation *Ptr_implementation; + + + public: + + // default constructor + iterator() + : Ptr_implementation(0) + {} + + int + the_number_of_items_visited() + { + return Ptr_implementation->number_of_items_visited; + } + + // constructor + iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0), + bool search_nearest=true) + : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest)) + {} + + // copy constructor + iterator(const iterator& Iter) + { + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + + iterator& operator=(const iterator& Iter) + { + if (Ptr_implementation != Iter.Ptr_implementation){ + if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) { + delete Ptr_implementation; + } + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + return *this; + } + + + const Point_with_transformed_distance& + operator* () const + { + return *(*Ptr_implementation); + } + + // -> operator + const Point_with_transformed_distance* + operator-> () const + { + return &*(*Ptr_implementation); + } + + // prefix operator + iterator& + operator++() + { + ++(*Ptr_implementation); + return *this; + } + + // postfix operator + Object_wrapper<Point_with_transformed_distance> + operator++(int) + { + return (*Ptr_implementation)++; + } + + + bool + operator==(const iterator& It) const + { + if ( + ((Ptr_implementation == 0) || + Ptr_implementation->Item_PriorityQueue.empty()) && + ((It.Ptr_implementation == 0) || + It.Ptr_implementation->Item_PriorityQueue.empty()) + ) + return true; + // else + return (Ptr_implementation == It.Ptr_implementation); + } + + bool + operator!=(const iterator& It) const + { + return !(*this == It); + } + + std::ostream& + statistics (std::ostream& s) + { + Ptr_implementation->statistics(s); + return s; + } + + ~iterator() + { + if (Ptr_implementation != 0) { + Ptr_implementation->reference_count--; + if (Ptr_implementation->reference_count==0) { + delete Ptr_implementation; + Ptr_implementation = 0; + } + } + } + + + }; // class iterator + + //data members + const Tree& m_tree; + Query_item m_query; + Distance m_dist; + FT m_Eps; + bool m_search_nearest; + }; // class + + template <class Traits, class Query_item, class Distance> + void swap (typename Orthogonal_incremental_neighbor_search<Traits, + Query_item, Distance>::iterator& x, + typename Orthogonal_incremental_neighbor_search<Traits, + Query_item, Distance>::iterator& y) + { + typename Orthogonal_incremental_neighbor_search<Traits, + Query_item, Distance>::iterator::Iterator_implementation + *tmp = x.Ptr_implementation; + x.Ptr_implementation = y.Ptr_implementation; + y.Ptr_implementation = tmp; + } + +} // namespace CGAL + +#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..e9f455d7 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -0,0 +1,179 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GRAPH_MATCHING_H_ +#define GRAPH_MATCHING_H_ + +#include <gudhi/Neighbors_finder.h> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph. + * + * \ingroup bottleneck_distance + */ +class Graph_matching { +public: + /** \internal \brief Constructor constructing an empty matching. */ + explicit Graph_matching(Persistence_graph &g); + /** \internal \brief Copy operator. */ + Graph_matching& operator=(const Graph_matching& m); + /** \internal \brief Is the matching perfect ? */ + bool perfect() const; + /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */ + bool multi_augment(); + /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ + void set_r(double r); + +private: + Persistence_graph& g; + double r; + /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ + std::vector<int> v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::list<int> unmatched_in_u; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + Layered_neighbors_finder layering() const; + /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ + void update(std::vector<int> & path); +}; + +inline Graph_matching::Graph_matching(Persistence_graph& g) + : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() { + for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) + unmatched_in_u.emplace_back(u_point_index); +} + +inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) { + g = m.g; + r = m.r; + v_to_u = m.v_to_u; + unmatched_in_u = m.unmatched_in_u; + return *this; +} + +inline bool Graph_matching::perfect() const { + return unmatched_in_u.empty(); +} + +inline bool Graph_matching::multi_augment() { + if (perfect()) + return false; + Layered_neighbors_finder layered_nf(layering()); + int max_depth = layered_nf.vlayers_number()*2 - 1; + double rn = sqrt(g.size()); + // verification of a necessary criterion in order to shortcut if possible + if (max_depth <0 || (unmatched_in_u.size() > rn && max_depth >= rn)) + return false; + bool successful = false; + std::list<int> tries(unmatched_in_u); + for (auto it = tries.cbegin(); it != tries.cend(); it++) + // 'augment' has side-effects which have to be always executed, don't change order + successful = augment(layered_nf, *it, max_depth) || successful; + return successful; +} + +inline void Graph_matching::set_r(double r) { + this->r = r; +} + +inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { + //V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::vector<int> path; + path.emplace_back(u_start_index); + do { + if (static_cast<int>(path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + path.emplace_back(layered_nf.pull_near(path.back(), static_cast<int>(path.size())/2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + //if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex + path.pop_back(); + update(path); + return true; +} + +inline Layered_neighbors_finder Graph_matching::layering() const { + std::list<int> u_vertices(unmatched_in_u); + std::list<int> v_vertices; + Neighbors_finder nf(g, r); + for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) + nf.add(v_point_index); + Layered_neighbors_finder layered_nf(g, r); + for(int layer = 0; !u_vertices.empty(); layer++) { + // one layer is one step in the BFS + for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { + std::vector<int> u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { + layered_nf.add(*it2, layer); + v_vertices.emplace_back(*it2); + } + } + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS + u_vertices.clear(); + bool end = false; + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) + if (v_to_u.at(*it) == null_point_index()) + // we stop when a nearest exposed V vertex (from U exposed vertices) has been found + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + // When the above for finishes, we have progress of one half-step (from V to U) in the BFS + if (end) + return layered_nf; + v_vertices.clear(); + } + return layered_nf; +} + +inline void Graph_matching::update(std::vector<int>& path) { + unmatched_in_u.remove(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + // Be careful, the iterator is incremented twice each time + int tmp = *it; + v_to_u[*(++it)] = tmp; + } +} + + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // GRAPH_MATCHING_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h new file mode 100644 index 00000000..70342d0c --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -0,0 +1,70 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef INTERNAL_POINT_H_ +#define INTERNAL_POINT_H_ + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Returns the used index for encoding none of the points */ +int null_point_index(); + +/** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ +struct Internal_point { + double vec[2]; + int point_index; + Internal_point() {} + Internal_point(double x, double y, int p_i) { vec[0]=x; vec[1]=y; point_index = p_i; } + double x() const { return vec[ 0 ]; } + double y() const { return vec[ 1 ]; } + double& x() { return vec[ 0 ]; } + double& y() { return vec[ 1 ]; } + bool operator==(const Internal_point& p) const + { + return point_index==p.point_index; + } + bool operator!=(const Internal_point& p) const { return !(*this == p); } +}; + +inline int null_point_index() { + return -1; +} + +struct Construct_coord_iterator { + typedef const double* result_type; + const double* operator()(const Internal_point& p) const + { return p.vec; } + const double* operator()(const Internal_point& p, int) const + { return p.vec+2; } +}; + +} // namespace persistence_diagram + +} // namespace Gudhi + + + + + +#endif // INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..792925b7 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -0,0 +1,169 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef NEIGHBORS_FINDER_H_ +#define NEIGHBORS_FINDER_H_ + +// Inclusion order is important for CGAL patch +#include "CGAL/Kd_tree_node.h" +#include "CGAL/Kd_tree.h" +#include "CGAL/Orthogonal_k_neighbor_search.h" + +#include <CGAL/Weighted_Minkowski_distance.h> +#include <CGAL/Search_traits.h> + +#include <gudhi/Persistence_graph.h> +#include <gudhi/Internal_point.h> + +#include <unordered_set> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief data structure used to find any point (including projections) in V near to a query point from U + * (which can be a projection). + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed. + * + * \ingroup bottleneck_distance + */ +class Neighbors_finder { + + typedef CGAL::Dimension_tag<2> D; + typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits; + typedef CGAL::Weighted_Minkowski_distance<Traits> Distance; + typedef CGAL::Orthogonal_k_neighbor_search<Traits, Distance> K_neighbor_search; + typedef K_neighbor_search::Tree Kd_tree; + +public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ + int pull_near(int u_point_index); + /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ + std::vector<int> pull_all_near(int u_point_index); + +private: + const Persistence_graph& g; + const double r; + Kd_tree kd_t; + std::unordered_set<int> projections_f; +}; + +/** \internal \brief data structure used to find any point (including projections) in V near to a query point from U + * (which can be a projection) in a layered graph layer given as parmeter. + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed. + * + * \ingroup bottleneck_distance + */ +class Layered_neighbors_finder { +public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Layered_neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index, int vlayer); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ + int pull_near(int u_point_index, int vlayer); + /** \internal \brief Returns the number of layers. */ + int vlayers_number() const; + +private: + const Persistence_graph& g; + const double r; + std::vector<std::unique_ptr<Neighbors_finder>> neighbors_finder; +}; + +inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), kd_t(), projections_f() { } + +inline void Neighbors_finder::add(int v_point_index) { + if (g.on_the_v_diagonal(v_point_index)) + projections_f.emplace(v_point_index); + else + kd_t.insert(g.get_v_point(v_point_index)); +} + +inline int Neighbors_finder::pull_near(int u_point_index) { + int tmp; + int c = g.corresponding_point_in_v(u_point_index); + if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()){ + //Any pair of projection is at distance 0 + tmp = *projections_f.cbegin(); + projections_f.erase(tmp); + } + else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)){ + //Is the query point near to its projection ? + tmp = c; + projections_f.erase(tmp); + } + else{ + //Is the query point near to a V point in the plane ? + Internal_point u_point = g.get_u_point(u_point_index); + std::array<double, 2> w = { {1., 1.} }; + K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end())); + auto it = search.begin(); + if(it==search.end() || g.distance(u_point_index, it->first.point_index) > r) + return null_point_index(); + tmp = it->first.point_index; + kd_t.remove(g.get_v_point(tmp)); + } + return tmp; +} + +inline std::vector<int> Neighbors_finder::pull_all_near(int u_point_index) { + std::vector<int> all_pull; + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull.emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; +} + +inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), neighbors_finder() { } + +inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(std::unique_ptr<Neighbors_finder>(new Neighbors_finder(g, r))); + neighbors_finder.at(vlayer)->add(v_point_index); +} + +inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { + if (static_cast<int> (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer)->pull_near(u_point_index); +} + +inline int Layered_neighbors_finder::vlayers_number() const { + return static_cast<int>(neighbors_finder.size()); +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h new file mode 100644 index 00000000..2ee55995 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -0,0 +1,179 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef PERSISTENCE_GRAPH_H_ +#define PERSISTENCE_GRAPH_H_ + +#include <vector> +#include <algorithm> +#include <gudhi/Internal_point.h> + +namespace Gudhi { + +namespace persistence_diagram { + + +/** \internal \brief Structure representing an euclidean bipartite graph containing + * the points from the two persistence diagrams (including the projections). + * + * \ingroup bottleneck_distance + */ +class Persistence_graph { +public: + /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */ + template<typename Persistence_diagram1, typename Persistence_diagram2> + Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); + /** \internal \brief Is the given point from U the projection of a point in V ? */ + bool on_the_u_diagonal(int u_point_index) const; + /** \internal \brief Is the given point from V the projection of a point in U ? */ + bool on_the_v_diagonal(int v_point_index) const; + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + int corresponding_point_in_u(int v_point_index) const; + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + int corresponding_point_in_v(int u_point_index) const; + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + double distance(int u_point_index, int v_point_index) const; + /** \internal \brief Returns size = |U| = |V|. */ + int size() const; + /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */ + double bottleneck_alive() const; + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + std::vector<double> sorted_distances() const; + /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */ + double diameter_bound() const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_u_point(int u_point_index) const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_v_point(int v_point_index) const; + +private: + std::vector<Internal_point> u; + std::vector<Internal_point> v; + double b_alive; +}; + +template<typename Persistence_diagram1, typename Persistence_diagram2> +Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, + const Persistence_diagram2 &diag2, double e) + : u(), v(), b_alive(0.) +{ + std::vector<double> u_alive; + std::vector<double> v_alive; + for (auto it = std::begin(diag1); it != std::end(diag1); ++it){ + if(std::get<1>(*it) == std::numeric_limits<double>::infinity()) + u_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size())); + } + for (auto it = std::begin(diag2); it != std::end(diag2); ++it){ + if(std::get<1>(*it) == std::numeric_limits<double>::infinity()) + v_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size())); + } + if (u.size() < v.size()) + swap(u, v); + std::sort(u_alive.begin(), u_alive.end()); + std::sort(v_alive.begin(), v_alive.end()); + if(u_alive.size() != v_alive.size()) + b_alive = std::numeric_limits<double>::infinity(); + else for(auto it_u=u_alive.cbegin(), it_v=v_alive.cbegin();it_u != u_alive.cend();++it_u, ++it_v) + b_alive = std::max(b_alive, std::fabs(*it_u - *it_v)); +} + +inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const { + return u_point_index >= static_cast<int> (u.size()); +} + +inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const { + return v_point_index >= static_cast<int> (v.size()); +} + +inline int Persistence_graph::corresponding_point_in_u(int v_point_index) const { + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size()); +} + +inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const { + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size()); +} + +inline double Persistence_graph::distance(int u_point_index, int v_point_index) const { + if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) + return 0.; + Internal_point p_u = get_u_point(u_point_index); + Internal_point p_v = get_v_point(v_point_index); + return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); +} + +inline int Persistence_graph::size() const { + return static_cast<int> (u.size() + v.size()); +} + +inline double Persistence_graph::bottleneck_alive() const{ + return b_alive; +} + +inline std::vector<double> Persistence_graph::sorted_distances() const { + std::vector<double> distances; + distances.push_back(0.); + for (int u_point_index = 0; u_point_index < size(); ++u_point_index){ + distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index))); + for (int v_point_index = 0; v_point_index < size(); ++v_point_index) + distances.push_back(distance(u_point_index, v_point_index)); + } + std::sort(distances.begin(), distances.end()); + return distances; +} + +inline Internal_point Persistence_graph::get_u_point(int u_point_index) const { + if (!on_the_u_diagonal(u_point_index)) + return u.at(u_point_index); + Internal_point projector = v.at(corresponding_point_in_v(u_point_index)); + double m = (projector.x() + projector.y()) / 2; + return Internal_point(m,m,u_point_index); +} + +inline Internal_point Persistence_graph::get_v_point(int v_point_index) const { + if (!on_the_v_diagonal(v_point_index)) + return v.at(v_point_index); + Internal_point projector = u.at(corresponding_point_in_u(v_point_index)); + double m = (projector.x() + projector.y()) / 2; + return Internal_point(m,m,v_point_index); +} + +inline double Persistence_graph::diameter_bound() const { + double max = 0.; + for(auto it = u.cbegin(); it != u.cend(); it++) + max = std::max(max, it->y()); + for(auto it = v.cbegin(); it != v.cend(); it++) + max = std::max(max, it->y()); + return max; +} + + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // PERSISTENCE_GRAPH_H_ diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt new file mode 100644 index 00000000..13213075 --- /dev/null +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -0,0 +1,29 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_tests) + + +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") +endif() + +# requires CGAL 4.8 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.8.0) + if (EIGEN3_FOUND) + add_executable ( bottleneckUT bottleneck_unit_test.cpp ) + add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) + target_link_libraries(bottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + + # Unitary tests + add_test(NAME bottleneckUT COMMAND ${CMAKE_CURRENT_BINARY_DIR}/bottleneckUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/bottleneckUT.xml --log_level=test_suite --report_level=no) + endif() + endif () +endif() diff --git a/src/Bottleneck/test/README b/src/Bottleneck_distance/test/README index 0e7b8673..0e7b8673 100644 --- a/src/Bottleneck/test/README +++ b/src/Bottleneck_distance/test/README diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp new file mode 100644 index 00000000..dcf98460 --- /dev/null +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -0,0 +1,63 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Bottleneck.h> +#include <chrono> +#include <fstream> +#include <random> + +using namespace Gudhi::persistence_diagram; + + +double upper_bound = 400.; // any real >0 + +int main(){ + std::ofstream objetfichier; + objetfichier.open("results.csv", std::ios::out); + + for(int n = 1000; n<=4000; n+=1000){ + std::uniform_real_distribution<double> unif1(0.,upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound/1000.,upper_bound/100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a,b), std::max(a,b)); + v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); + if(i%5==0) + v1.emplace_back(std::min(a,b),std::min(a,b)+x); + if(i%3==0) + v2.emplace_back(std::max(a,b),std::max(a,b)+y); + } + double epsilon = 0.0000000000000001; + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + double b = bottleneck_distance(v1,v2, epsilon); + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + typedef std::chrono::duration<int,std::milli> millisecs_t; + millisecs_t duration(std::chrono::duration_cast<millisecs_t>(end-start)); + objetfichier << n << ";" << duration.count() << ";" << b << std::endl; + } + objetfichier.close(); +} diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp new file mode 100644 index 00000000..fba1d369 --- /dev/null +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -0,0 +1,167 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "bottleneck distance" +#include <boost/test/unit_test.hpp> + +#include <random> +#include <gudhi/Bottleneck.h> + +using namespace Gudhi::persistence_diagram; + +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 406.43; // any real >0 + + +std::uniform_real_distribution<double> unif(0.,upper_bound); +std::default_random_engine re; +std::vector< std::pair<double, double> > v1, v2; + +BOOST_AUTO_TEST_CASE(persistence_graph){ + // Random construction + for (int i = 0; i < n1; i++) { + double a = unif(re); + double b = unif(re); + v1.emplace_back(std::min(a,b), std::max(a,b)); + } + for (int i = 0; i < n2; i++) { + double a = unif(re); + double b = unif(re); + v2.emplace_back(std::min(a,b), std::max(a,b)); + } + Persistence_graph g(v1, v2, 0.); + std::vector<double> d(g.sorted_distances()); + // + BOOST_CHECK(!g.on_the_u_diagonal(n1-1)); + BOOST_CHECK(!g.on_the_u_diagonal(n1)); + BOOST_CHECK(!g.on_the_u_diagonal(n2-1)); + BOOST_CHECK(g.on_the_u_diagonal(n2)); + BOOST_CHECK(!g.on_the_v_diagonal(n1-1)); + BOOST_CHECK(g.on_the_v_diagonal(n1)); + BOOST_CHECK(g.on_the_v_diagonal(n2-1)); + BOOST_CHECK(g.on_the_v_diagonal(n2)); + // + BOOST_CHECK(g.corresponding_point_in_u(0)==n2); + BOOST_CHECK(g.corresponding_point_in_u(n1)==0); + BOOST_CHECK(g.corresponding_point_in_v(0)==n1); + BOOST_CHECK(g.corresponding_point_in_v(n2)==0); + // + BOOST_CHECK(g.size()==(n1+n2)); + // + BOOST_CHECK((int) d.size() == (n1+n2)*(n1+n2) + n1 + n2 + 1); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,0))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,(n1+n2)-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,0))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,(n1+n2)-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,0))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2-1))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2))>0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,(n1+n2)-1))>0); +} + +BOOST_AUTO_TEST_CASE(neighbors_finder) { + Persistence_graph g(v1, v2, 0.); + Neighbors_finder nf(g, 1.); + for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) + nf.add(v_point_index); + // + int v_point_index_1 = nf.pull_near(n2/2); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); + std::vector<int> l = nf.pull_all_near(n2/2); + bool v = true; + for(auto it = l.cbegin(); it != l.cend(); ++it) + v = v && (g.distance(n2/2,*it)>1.); + BOOST_CHECK(v); + int v_point_index_2 = nf.pull_near(n2/2); + BOOST_CHECK(v_point_index_2 == -1); +} + +BOOST_AUTO_TEST_CASE(layered_neighbors_finder) { + Persistence_graph g(v1, v2, 0.); + Layered_neighbors_finder lnf(g, 1.); + for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) + lnf.add(v_point_index, v_point_index % 7); + // + int v_point_index_1 = lnf.pull_near(n2/2,6); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); + int v_point_index_2 = lnf.pull_near(n2/2,6); + BOOST_CHECK(v_point_index_2 == -1); + v_point_index_1 = lnf.pull_near(n2/2,0); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); + v_point_index_2 = lnf.pull_near(n2/2,0); + BOOST_CHECK(v_point_index_2 == -1); +} + +BOOST_AUTO_TEST_CASE(graph_matching) { + Persistence_graph g(v1, v2, 0.); + Graph_matching m1(g); + m1.set_r(0.); + int e = 0; + while (m1.multi_augment()) + ++e; + BOOST_CHECK(e > 0); + BOOST_CHECK(e <= 2*sqrt(2*(n1+n2))); + Graph_matching m2 = m1; + BOOST_CHECK(!m2.multi_augment()); + m2.set_r(upper_bound); + e = 0; + while (m2.multi_augment()) + ++e; + BOOST_CHECK(e <= 2*sqrt(2*(n1+n2))); + BOOST_CHECK(m2.perfect()); + BOOST_CHECK(!m1.perfect()); +} + +BOOST_AUTO_TEST_CASE(global){ + std::uniform_real_distribution<double> unif1(0.,upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound/10000.,upper_bound/100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n1; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a,b), std::max(a,b)); + v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); + if(i%5==0) + v1.emplace_back(std::min(a,b),std::min(a,b)+x); + if(i%3==0) + v2.emplace_back(std::max(a,b),std::max(a,b)+y); + } + BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound/100.); + BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound/10000.) <= upper_bound/100. + upper_bound/10000.); + BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound/10000.)) <= upper_bound/10000.); +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e26b2d25..1fa62e35 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -115,6 +115,7 @@ else() add_subdirectory(example/Spatial_searching) add_subdirectory(example/Subsampling) add_subdirectory(example/Tangential_complex) + add_subdirectory(example/Bottleneck_distance) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Doxyfile b/src/Doxyfile index 943869ad..71aa038d 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -848,7 +848,8 @@ IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Bitmap_cubical_complex/ \ doc/Subsampling/ \ doc/Spatial_searching/ \ - doc/Tangential_complex/ + doc/Tangential_complex/ \ + doc/Bottleneck_distance/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/src/cmake/modules/GUDHI_user_version_target.txt b/src/cmake/modules/GUDHI_user_version_target.txt index 51553e7e..25db1c87 100644 --- a/src/cmake/modules/GUDHI_user_version_target.txt +++ b/src/cmake/modules/GUDHI_user_version_target.txt @@ -48,7 +48,7 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI) - set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Contraction;Hasse_complex;Persistent_cohomology;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex") + set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Bottleneck_distance;Contraction;Hasse_complex;Persistent_cohomology;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex") foreach(GUDHI_MODULE ${GUDHI_MODULES}) # doc files diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 1a2cb6ba..cf67b558 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -28,6 +28,7 @@ <b>Author:</b> Vincent Rouvreau<br> <b>Introduced in:</b> GUDHI 1.3.0<br> <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref cgal ≥ 4.7.0 and \ref eigen3 </td> <td width="75%"> Alpha_complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation.<br> @@ -130,6 +131,26 @@ </table> \section Toolbox Toolbox + \subsection BottleneckDistanceToolbox Bottleneck distance + \image html "perturb_pd.png" "Bottleneck distance is the length of the longest edge" +<table border="0"> + <tr> + <td width="25%"> + <b>Author:</b> François Godi<br> + <b>Introduced in:</b> GUDHI 1.4.0<br> + <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref cgal ≥ 4.8.0 and \ref eigen3 + </td> + <td width="75%"> + Bottleneck distance measures the similarity between two persistence diagrams. + It's the shortest distance b for which there exists a perfect matching between + the points of the two diagrams (+ all the diagonal points) such that + any couple of matched points are at distance at most b. + <br> + <b>User manual:</b> \ref bottleneck_distance + </td> + </tr> +</table> \subsection ContractionToolbox Contraction \image html "sphere_contraction_representation.png" "Sphere contraction example" <table border="0"> @@ -196,7 +217,7 @@ make \endverbatim * \verbatim make test \endverbatim * * \section optionallibrary Optional third-party library - * \subsection gmp GMP: + * \subsection gmp GMP * The multi-field persistent homology algorithm requires GMP which is a free library for arbitrary-precision * arithmetic, operating on signed integers, rational numbers, and floating point numbers. * @@ -209,11 +230,11 @@ make \endverbatim * * Having GMP version 4.2 or higher installed is recommended. * - * \subsection cgal CGAL: - * The \ref alpha_complex data structure and few examples requires CGAL, which is a C++ library which provides easy - * access to efficient and reliable geometric algorithms. + * \subsection cgal CGAL + * The \ref alpha_complex data structure, \ref bottleneck_distance, and few examples requires CGAL, which is a C++ + * library which provides easy access to efficient and reliable geometric algorithms. * - * Having CGAL version 4.4 or higher installed is recommended. The procedure to install this library according to + * Having CGAL version 4.4.0 or higher installed is recommended. The procedure to install this library according to * your operating system is detailed here http://doc.cgal.org/latest/Manual/installation.html * * The following examples require the <a target="_blank" href="http://www.cgal.org/">Computational Geometry Algorithms @@ -223,11 +244,11 @@ make \endverbatim * \li <a href="_simplex_tree_2example_alpha_shapes_3_simplex_tree_from_off_file_8cpp-example.html"> * Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp</a> * - * The following example requires CGAL version ≥ 4.6: + * The following example requires CGAL version ≥ 4.6.0: * \li <a href="_witness_complex_2witness_complex_sphere_8cpp-example.html"> * Witness_complex/witness_complex_sphere.cpp</a> * - * The following example requires CGAL version ≥ 4.7: + * The following example requires CGAL version ≥ 4.7.0: * \li <a href="_alpha_complex_2_alpha_complex_from_off_8cpp-example.html"> * Alpha_complex/Alpha_complex_from_off.cpp</a> * \li <a href="_alpha_complex_2_alpha_complex_from_points_8cpp-example.html"> @@ -239,7 +260,11 @@ make \endverbatim * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> * Persistent_cohomology/custom_persistence_sort.cpp</a> * - * \subsection eigen3 Eigen3: + * The following example requires CGAL version ≥ 4.8.0: + * \li <a href="_bottleneck_distance_2_bottleneck_example_8cpp-example.html"> + * + * Bottleneck_distance/bottleneck_example.cpp</a> + * \subsection eigen3 Eigen3 * The \ref alpha_complex data structure and few examples requires * <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen3</a> is a C++ template library for linear algebra: * matrices, vectors, numerical solvers, and related algorithms. @@ -247,9 +272,11 @@ make \endverbatim * The following example requires the <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen3</a> and will not be * built if Eigen3 is not installed: * \li <a href="_alpha_complex_2_alpha_complex_from_off_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_off.cpp</a> (requires also Eigen3) + * Alpha_complex/Alpha_complex_from_off.cpp</a> * \li <a href="_alpha_complex_2_alpha_complex_from_points_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_points.cpp</a> (requires also Eigen3) + * Alpha_complex/Alpha_complex_from_points.cpp</a> + * \li <a href="_bottleneck_distance_2_bottleneck_example_8cpp-example.html"> + * Bottleneck_distance/bottleneck_example.cpp</a> * \li <a href="_persistent_cohomology_2alpha_complex_persistence_8cpp-example.html"> * Persistent_cohomology/alpha_complex_persistence.cpp</a> * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html"> @@ -257,7 +284,7 @@ make \endverbatim * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> * Persistent_cohomology/custom_persistence_sort.cpp</a> * - * \subsection tbb Threading Building Blocks: + * \subsection tbb Threading Building Blocks * <a target="_blank" href="https://www.threadingbuildingblocks.org/">Intel® TBB</a> lets you easily write parallel * C++ programs that take full advantage of multicore performance, that are portable and composable, and that have * future-proof scalability. @@ -331,6 +358,7 @@ make \endverbatim /*! @file Examples * @example Alpha_complex/Alpha_complex_from_off.cpp * @example Alpha_complex/Alpha_complex_from_points.cpp + * @example Bottleneck_distance/bottleneck_example.cpp * @example Bitmap_cubical_complex/Bitmap_cubical_complex.cpp * @example Bitmap_cubical_complex/Bitmap_cubical_complex_periodic_boundary_conditions.cpp * @example Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp |
