From eaae83fc6b471bb05891addd35863fe00c351565 Mon Sep 17 00:00:00 2001 From: fgodi Date: Sun, 29 Mar 2015 22:37:12 +0000 Subject: Copyrigth Sophia + tests git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@512 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 353baa0c191df44a3bc420efd4424de2f5c6859e --- src/Bottleneck/include/gudhi/Graph_matching.h | 53 +++++----- .../include/gudhi/Layered_neighbors_finder.h | 2 +- src/Bottleneck/include/gudhi/Neighbors_finder.h | 2 +- .../include/gudhi/Persistence_diagrams_graph.h | 113 +++++++++++---------- .../include/gudhi/Planar_neighbors_finder.h | 3 +- 5 files changed, 87 insertions(+), 86 deletions(-) (limited to 'src/Bottleneck/include/gudhi') diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h index ea47e1d5..4cd3180e 100644 --- a/src/Bottleneck/include/gudhi/Graph_matching.h +++ b/src/Bottleneck/include/gudhi/Graph_matching.h @@ -4,7 +4,7 @@ * * Author(s): Francois Godi * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 INRIA Sophia-Antipolis (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 @@ -34,11 +34,11 @@ namespace Gudhi { namespace bottleneck { template -double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.); +double bottleneck_distance(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); class Graph_matching { public: - Graph_matching(const Persistence_diagrams_graph& g); + explicit Graph_matching(const Persistence_diagrams_graph& g); Graph_matching& operator=(const Graph_matching& m); bool perfect() const; bool multi_augment(); @@ -52,7 +52,7 @@ class Graph_matching { Layered_neighbors_finder* layering() const; bool augment(Layered_neighbors_finder* layered_nf, int u_start_index, int max_depth); - void update(std::deque& path); + void update(std::deque* path); }; Graph_matching::Graph_matching(const Persistence_diagrams_graph& g) @@ -128,36 +128,37 @@ Layered_neighbors_finder* Graph_matching::layering() const { } bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) { - std::deque path; - path.emplace_back(u_start_index); + std::deque* path = new std::deque(); + path->emplace_back(u_start_index); // start is a point from U do { - if (static_cast(path.size()) > max_depth) { - path.pop_back(); - path.pop_back(); + if (static_cast(path->size()) > max_depth) { + path->pop_back(); + path->pop_back(); } - if (path.empty()) + 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()) + 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->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(); + path->emplace_back(v_to_u.at(path->back())); + } while (path->back() != null_point_index()); + path->pop_back(); update(path); + delete path; return true; } -void Graph_matching::update(std::deque& path) { - unmatched_in_u.remove(path.front()); - for (auto it = path.cbegin(); it != path.cend(); ++it) { +void Graph_matching::update(std::deque *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; @@ -165,14 +166,14 @@ void Graph_matching::update(std::deque& path) { } template -double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) { +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { Persistence_diagrams_graph g(diag1, diag2, e); std::vector* 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; + Graph_matching biggest_unperfect(g); while (idmin != idmax) { int pas = static_cast((idmax - idmin) / alpha); m.set_r(sd->at(idmin + pas)); diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h index de36e00b..80e5b5ed 100644 --- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h @@ -4,7 +4,7 @@ * * Author(s): Francois Godi * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 INRIA Sophia-Antipolis (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 diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h index 98256571..6a05e1a0 100644 --- a/src/Bottleneck/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h @@ -4,7 +4,7 @@ * * Author(s): Francois Godi * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 INRIA Sophia-Antipolis (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 diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h index 7e278209..8e9eaaf6 100644 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h @@ -4,7 +4,7 @@ * * Author(s): Francois Godi * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 INRIA Sophia-Antipolis (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 @@ -44,100 +44,101 @@ int null_point_index(); 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 - 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* sorted_distances(); + // 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 + Persistence_diagrams_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); + 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* sorted_distances(); private: - std::vector u; - std::vector v; - Diagram_point get_u_point(int u_point_index) const; - Diagram_point get_v_point(int v_point_index) const; + std::vector u; + std::vector 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; + return -1; } template -Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) +Persistence_diagrams_graph::Persistence_diagrams_graph(const Persistence_diagram1 &diag1, + const 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); + 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::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 (u.size()); + return u_point_index >= static_cast (u.size()); } inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const { - return v_point_index >= static_cast (v.size()); + return v_point_index >= static_cast (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 (v.size()) : v_point_index + static_cast (u.size()); + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast (v.size()) : v_point_index + static_cast (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 (u.size()) : u_point_index + static_cast (v.size()); + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast (u.size()) : u_point_index + static_cast (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)); + // 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 (u.size() + v.size()); + return static_cast (u.size() + v.size()); } inline std::vector* Persistence_diagrams_graph::sorted_distances() { - // could be optimized - std::set 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(sorted_distances.cbegin(), sorted_distances.cend()); + // could be optimized + std::set 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(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); + 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); + 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 diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h index 4af672e4..3d04aa97 100644 --- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h @@ -4,7 +4,7 @@ * * Author(s): Francois Godi * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 INRIA Sophia-Antipolis (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 @@ -24,7 +24,6 @@ #define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ #include -#include #include #include "Persistence_diagrams_graph.h" -- cgit v1.2.3 From ef48b9447c7270ba63d45a76e4d2e6ecb1a7d199 Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 25 Jun 2015 10:19:49 +0000 Subject: tests unitaires terminés MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@647 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: aca61ea029d72ddfbf8eca3e4db00da3e7b268c3 --- src/Bottleneck/example/CMakeLists.txt | 23 +- src/Bottleneck/include/gudhi/Graph_matching.h | 284 +++++++++++---------- .../include/gudhi/Layered_neighbors_finder.h | 8 +- src/Bottleneck/include/gudhi/Neighbors_finder.h | 93 ++++--- .../include/gudhi/Persistence_diagrams_graph.h | 27 +- .../include/gudhi/Planar_neighbors_finder.h | 16 +- src/Bottleneck/test/CMakeLists.txt | 6 + src/Bottleneck/test/bottleneck_unit_test.cpp | 187 +++++++++----- 8 files changed, 385 insertions(+), 259 deletions(-) (limited to 'src/Bottleneck/include/gudhi') diff --git a/src/Bottleneck/example/CMakeLists.txt b/src/Bottleneck/example/CMakeLists.txt index 2ff009c4..8bf61c50 100644 --- a/src/Bottleneck/example/CMakeLists.txt +++ b/src/Bottleneck/example/CMakeLists.txt @@ -1,5 +1,24 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIBottleneckExample) -add_executable ( RandomDiagrams random_diagrams.cpp ) -add_test(RandomDiagrams ${CMAKE_CURRENT_BINARY_DIR}/RandomDiagrams) +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") +endif() + message("CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}") + message("CMAKE_CXX_FLAGS_DEBUG = ${CMAKE_CXX_FLAGS_DEBUG}") + message("CMAKE_CXX_FLAGS_RELEASE = ${CMAKE_CXX_FLAGS_RELEASE}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic -Wsign-compare") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") + +add_executable ( BottleneckExample example.cpp ) + diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h index 4cd3180e..9f36e936 100644 --- a/src/Bottleneck/include/gudhi/Graph_matching.h +++ b/src/Bottleneck/include/gudhi/Graph_matching.h @@ -22,12 +22,12 @@ #ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ #define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ - +//#define DEBUG #include #include #include -#include "gudhi/Layered_neighbors_finder.h" +#include "Layered_neighbors_finder.h" namespace Gudhi { @@ -37,158 +37,182 @@ template double bottleneck_distance(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); class Graph_matching { - public: - explicit 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 v_to_u; - std::list 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* path); +public: + explicit 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 v_to_u; + std::list unmatched_in_u; + + std::unique_ptr layering() const; + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + void update(std::deque & 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); + : 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); } 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; + 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::perfect() const { +#ifdef DEBUG + std::cout << " perfect? unmatched_in_u.size = " << unmatched_in_u.size() << std::endl << std::flush; +#endif + 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* tries = new std::list(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 */ bool Graph_matching::multi_augment() { + if (perfect()) + return false; +#ifdef DEBUG + std::cout << " multi augment" << std::endl << std::flush; +#endif + Layered_neighbors_finder layered_nf = *layering(); + double rn = sqrt(g.size()); + int max_depth = layered_nf.vlayers_number()*2 - 1; +#ifdef DEBUG + std::cout<< " nb_max_layer = " << max_depth << std::endl << std::flush; +#endif + // verification of a necessary criterion + if (max_depth <0 || (unmatched_in_u.size() > rn && max_depth >= rn)) + return false; + bool successful = false; + std::list tries(unmatched_in_u); + for (auto it = tries.cbegin(); it != tries.cend(); it++){ + const bool tmp = augment(layered_nf, *it, max_depth); //force augment evaluation even if successful is already true + successful = successful || tmp; + } + return successful; } -inline void Graph_matching::set_r(double r) { - this->r = r; +/* 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 u_vertices(unmatched_in_u); - std::list 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* 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) { +#ifdef DEBUG + std::cout << " augment" << std::endl << std::flush; +#endif +#ifdef DEBUG + std::cout << " u_start_index = " << u_start_index << std::endl << std::flush; +#endif + std::deque path; + path.emplace_back(u_start_index); + // u_start is a point from U + do { +#ifdef DEBUG + std::cout << " do" << std::endl << std::flush; + std::cout << " path.size = " << static_cast(path.size()) << std::endl << std::flush; +#endif + if (static_cast(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(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())); +#ifdef DEBUG + std::cout << " v_to_u = " << path.back() << std::endl << std::flush; +#endif + } while (path.back() != null_point_index()); + path.pop_back(); + update(path); + return true; } -bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) { - std::deque* path = new std::deque(); - path->emplace_back(u_start_index); - // start is a point from U - do { - if (static_cast(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)); +std::unique_ptr Graph_matching::layering() const { +#ifdef DEBUG + std::cout << " layering" << std::endl << std::flush; +#endif + bool end = false; + int layer = 0; + std::list u_vertices(unmatched_in_u); + std::list 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); + std::unique_ptr layered_nf(new Layered_neighbors_finder(g, r)); + while (!u_vertices.empty()) { + for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) { + std::unique_ptr< std::list > u_succ = std::move(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); + } + } + 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++; } - path->emplace_back(v_to_u.at(path->back())); - } while (path->back() != null_point_index()); - path->pop_back(); - update(path); - delete path; - return true; + return layered_nf; } -void Graph_matching::update(std::deque *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; - } +void Graph_matching::update(std::deque& path) { +#ifdef DEBUG + std::cout << " update" << std::endl << std::flush; +#endif + 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 double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { - Persistence_diagrams_graph g(diag1, diag2, e); - std::vector* 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(g); - while (idmin != idmax) { - int pas = static_cast((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; + Persistence_diagrams_graph g(diag1, diag2, e); + std::unique_ptr< std::vector > sd = std::move(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(g); + while (idmin != idmax) { + int pas = static_cast((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; + double b = sd->at(idmin); + return b; } } // namespace bottleneck diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h index 80e5b5ed..1cf83393 100644 --- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h @@ -51,20 +51,20 @@ class Layered_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) { +/* 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) { +/* inline */ int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { if (static_cast (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(); +/* inline */ int Layered_neighbors_finder::vlayers_number() const { + return static_cast(neighbors_finder.size()); } } // namespace bottleneck diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h index 6a05e1a0..0fad9889 100644 --- a/src/Bottleneck/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h @@ -26,7 +26,7 @@ #include #include -#include "gudhi/Planar_neighbors_finder.h" +#include "Planar_neighbors_finder.h" namespace Gudhi { @@ -37,56 +37,65 @@ namespace bottleneck { // 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* 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 projections_f; +public: + Neighbors_finder(const Persistence_diagrams_graph& g, double r); + void add(int v_point_index); + int pull_near(int u_point_index); + std::unique_ptr< std::list > 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 projections_f; + void remove(int v_point_index); + bool contains(int v_point_index); }; 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 */ 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 */ void Neighbors_finder::remove(int v_point_index) { + if(v_point_index == null_point_index()) + return; + if (g.on_the_v_diagonal(v_point_index)) + projections_f.erase(v_point_index); + else + planar_neighbors_f.remove(v_point_index); +} + +/* inline */ bool Neighbors_finder::contains(int v_point_index) { + return planar_neighbors_f.contains(v_point_index) || (projections_f.count(v_point_index)>0); } -inline std::list* Neighbors_finder::pull_all_near(int u_point_index) { - std::list* 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; +/* 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.cbegin()!=projections_f.cend()) + tmp = *projections_f.cbegin(); + else if (contains(c) && (g.distance(u_point_index, c) <= r)) + tmp = c; + else + tmp = planar_neighbors_f.pull_near(u_point_index); + remove(tmp); + return tmp; +} + +std::unique_ptr< std::list > Neighbors_finder::pull_all_near(int u_point_index) { + std::unique_ptr< std::list > all_pull = std::move(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 diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h index 105f9059..644c1cd8 100644 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h @@ -55,7 +55,7 @@ class Persistence_diagrams_graph { 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* sorted_distances(); + std::unique_ptr< std::vector > sorted_distances(); private: std::vector u; @@ -64,7 +64,7 @@ class Persistence_diagrams_graph { Diagram_point get_v_point(int v_point_index) const; }; -inline int null_point_index() { +/* inline */ int null_point_index() { return -1; } @@ -85,47 +85,48 @@ Persistence_diagrams_graph::Persistence_diagrams_graph(const Persistence_diagram Persistence_diagrams_graph::Persistence_diagrams_graph() : u(), v() { } -inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const { +/* inline */ bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const { return u_point_index >= static_cast (u.size()); } -inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const { +/* inline */ bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const { return v_point_index >= static_cast (v.size()); } -inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const { +/* 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 (v.size()) : v_point_index + static_cast (u.size()); } -inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const { +/* 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 (u.size()) : u_point_index + static_cast (v.size()); } -inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) const { +/* 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)); + 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 { +/* inline */ int Persistence_diagrams_graph::size() const { return static_cast (u.size() + v.size()); } -inline std::vector* Persistence_diagrams_graph::sorted_distances() { +/* inline */ std::unique_ptr< std::vector > Persistence_diagrams_graph::sorted_distances() { // could be optimized std::set 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(sorted_distances.cbegin(), sorted_distances.cend()); + std::unique_ptr< std::vector > sd_up(new std::vector(sorted_distances.cbegin(), sorted_distances.cend())); + return sd_up; } -inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const { +/* 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)); @@ -133,7 +134,7 @@ inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) return Diagram_point(x, x); } -inline Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const { +/* 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)); diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h index 3d04aa97..e558239d 100644 --- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h @@ -45,7 +45,7 @@ class Abstract_planar_neighbors_finder { 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* pull_all_near(int u_point_index); + virtual std::unique_ptr< std::list > pull_all_near(int u_point_index); protected: const Persistence_diagrams_graph& g; @@ -74,10 +74,10 @@ 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 */ Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() { } -inline std::list* Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { - std::list* all_pull = new std::list(); +/* inline */ std::unique_ptr< std::list > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { + std::unique_ptr< std::list > all_pull(new std::list); int last_pull = pull_near(u_point_index); while (last_pull != null_point_index()) { all_pull->emplace_back(last_pull); @@ -89,19 +89,19 @@ inline std::list* Abstract_planar_neighbors_finder::pull_all_near(int u_poi 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) { +/* inline */ void Naive_pnf::add(int v_point_index) { candidates.emplace(v_point_index); } -inline void Naive_pnf::remove(int 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 { +/* 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) { +/* 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; diff --git a/src/Bottleneck/test/CMakeLists.txt b/src/Bottleneck/test/CMakeLists.txt index 3dfd80cd..8e4371b8 100644 --- a/src/Bottleneck/test/CMakeLists.txt +++ b/src/Bottleneck/test/CMakeLists.txt @@ -13,6 +13,12 @@ if (GPROF_PATH) set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") endif() + message("CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}") + message("CMAKE_CXX_FLAGS_DEBUG = ${CMAKE_CXX_FLAGS_DEBUG}") + message("CMAKE_CXX_FLAGS_RELEASE = ${CMAKE_CXX_FLAGS_RELEASE}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic -Wsign-compare") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") add_executable ( BottleneckUT bottleneck_unit_test.cpp ) target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp index ea7f51f2..77f45867 100644 --- a/src/Bottleneck/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck/test/bottleneck_unit_test.cpp @@ -2,14 +2,15 @@ #include #include -#include "gudhi/Graph_matching.h" +#include "../include/gudhi/Graph_matching.h" using namespace Gudhi::bottleneck; -Persistence_diagrams_graph* random_graph_generator(){ - int n1 = 100; - int n2 = 120; - double upper_bound = 80; +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 400.5; // any real >0 + +std::unique_ptr random_graph_generator(){ // Random construction std::uniform_real_distribution unif(0.,upper_bound); std::default_random_engine re; @@ -24,73 +25,83 @@ Persistence_diagrams_graph* random_graph_generator(){ double b = unif(re); v2.emplace_back(std::min(a,b), std::max(a,b)); } - return new Persistence_diagrams_graph(v1, v2, 0.); + return std::unique_ptr(new Persistence_diagrams_graph(v1, v2, 0.)); } + BOOST_AUTO_TEST_CASE(global){ - int n = 100; - // Random construction + std::uniform_real_distribution unif1(0.,upper_bound); + std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); + std::default_random_engine re; std::vector< std::pair > 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)); + 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, 1.) == 98); + BOOST_CHECK(bottleneck_distance(v1, v2) <= upper_bound/100.); } BOOST_AUTO_TEST_CASE(persistence_diagrams_graph) { - Persistence_diagrams_graph* g = random_graph_generator(); - std::vector* d = g->sorted_distances(); + std::unique_ptr g = std::move(random_graph_generator()); + std::unique_ptr< std::vector > d = std::move(g->sorted_distances()); // - BOOST_CHECK(!g->on_the_u_diagonal(99)); - BOOST_CHECK(!g->on_the_u_diagonal(100)); - BOOST_CHECK(!g->on_the_u_diagonal(119)); - BOOST_CHECK(g->on_the_u_diagonal(120)); - BOOST_CHECK(!g->on_the_v_diagonal(99)); - BOOST_CHECK(g->on_the_v_diagonal(100)); - BOOST_CHECK(g->on_the_v_diagonal(119)); - BOOST_CHECK(g->on_the_v_diagonal(120)); + 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)==120); - BOOST_CHECK(g->corresponding_point_in_u(100)==0); - BOOST_CHECK(g->corresponding_point_in_v(0)==100); - BOOST_CHECK(g->corresponding_point_in_v(120)==0); + 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()==220); + BOOST_CHECK(g->size()==(n1+n2)); // - BOOST_CHECK((int) d->size() <= 220*220 - 100*120 + 1); // could be more strict + BOOST_CHECK((int) d->size() <= (n1+n2)*(n1+n2) - n1*n2 + 1); BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,0))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,99))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,100))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,119))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,120))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,219))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,0))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,99))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,100))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,119))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,120))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,219))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,0))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,99))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,100))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,119))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,120))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,219))==1); - // - delete g; - delete d; + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,(n1+n2)-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,0))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,(n1+n2)-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,0))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,(n1+n2)-1))==1); } -BOOST_AUTO_TEST_CASE(planar_nf) { - Persistence_diagrams_graph* g = random_graph_generator(); +BOOST_AUTO_TEST_CASE(planar_neighbors_finder) { + std::unique_ptr g = std::move(random_graph_generator()); Planar_neighbors_finder pnf = Planar_neighbors_finder(*g,1.); - for(int v_point_index=0; v_point_index<100; v_point_index+=2) + for(int v_point_index=0; v_point_indexdistance(120/2,v_point_index_1)<1.)); + int v_point_index_1 = pnf.pull_near(n2/2); + BOOST_CHECK((v_point_index_1 == -1) || (g->distance(n2/2,v_point_index_1)<=1.)); BOOST_CHECK(!pnf.contains(v_point_index_1)); - int v_point_index_2 = pnf.pull_near(120/2); - BOOST_CHECK((v_point_index_2 == -1) || (g->distance(120/2,v_point_index_2)<1.)); - BOOST_CHECK(!pnf.contains(v_point_index_2)); - BOOST_CHECK((v_point_index_2 != -1) || (v_point_index_1 == -1)); + std::list l = *pnf.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 = pnf.pull_near(n2/2); + BOOST_CHECK(v_point_index_2 == -1); pnf.add(v_point_index_1); BOOST_CHECK(pnf.contains(v_point_index_1)); } + +BOOST_AUTO_TEST_CASE(neighbors_finder) { + std::unique_ptr g = std::move(random_graph_generator()); + Neighbors_finder nf = Neighbors_finder(*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::list 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) { + std::unique_ptr g = std::move(random_graph_generator()); + Layered_neighbors_finder lnf = Layered_neighbors_finder(*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) { + std::unique_ptr g = std::move(random_graph_generator()); + Graph_matching m1(*g); + m1.set_r(0.); + int e = 0; + while (m1.multi_augment()) + ++e; + 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()); +} -- cgit v1.2.3 From 4ffac47ea18e63f40e6175055166d6bf3d08cff9 Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 29 Jun 2015 10:01:45 +0000 Subject: persistence_diagrams_graph static + others things git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@659 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 2562cbd49534373dae6d6d23ec52ba5984080729 --- src/Bottleneck/include/gudhi/Graph_matching.h | 126 +++++++++---------- .../include/gudhi/Layered_neighbors_finder.h | 56 +++++---- src/Bottleneck/include/gudhi/Neighbors_finder.h | 45 ++++--- .../include/gudhi/Persistence_diagrams_graph.h | 132 ++++++++++++-------- .../include/gudhi/Planar_neighbors_finder.h | 130 +++++++++++--------- src/Bottleneck/test/bottleneck_unit_test.cpp | 134 ++++++++++----------- 6 files changed, 332 insertions(+), 291 deletions(-) (limited to 'src/Bottleneck/include/gudhi') diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h index 9f36e936..17412e6c 100644 --- a/src/Bottleneck/include/gudhi/Graph_matching.h +++ b/src/Bottleneck/include/gudhi/Graph_matching.h @@ -22,7 +22,7 @@ #ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ #define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ -//#define DEBUG + #include #include #include @@ -33,31 +33,50 @@ namespace Gudhi { namespace bottleneck { +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. + * + * + * + * \ingroup bottleneck_distance + */ template double bottleneck_distance(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); +/** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph. + * + * \ingroup bottleneck_distance + */ class Graph_matching { public: - explicit Graph_matching(const Persistence_diagrams_graph& g); + /** \internal \brief Constructor constructing an empty matching. */ + explicit Graph_matching(); + /** \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: - const Persistence_diagrams_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 v_to_u; + /** \internal \brief All the unmatched points in U. */ std::list unmatched_in_u; + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ std::unique_ptr 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::deque & path); }; -Graph_matching::Graph_matching(const Persistence_diagrams_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) +Graph_matching::Graph_matching() + : 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); } @@ -68,56 +87,36 @@ Graph_matching& Graph_matching::operator=(const Graph_matching& m) { return *this; } -/* inline */ bool Graph_matching::perfect() const { -#ifdef DEBUG - std::cout << " perfect? unmatched_in_u.size = " << unmatched_in_u.size() << std::endl << std::flush; -#endif +inline bool Graph_matching::perfect() const { return unmatched_in_u.empty(); } -/* inline */ bool Graph_matching::multi_augment() { +inline bool Graph_matching::multi_augment() { if (perfect()) return false; -#ifdef DEBUG - std::cout << " multi augment" << std::endl << std::flush; -#endif Layered_neighbors_finder layered_nf = *layering(); - double rn = sqrt(g.size()); int max_depth = layered_nf.vlayers_number()*2 - 1; -#ifdef DEBUG - std::cout<< " nb_max_layer = " << max_depth << std::endl << std::flush; -#endif - // verification of a necessary criterion + 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 tries(unmatched_in_u); - for (auto it = tries.cbegin(); it != tries.cend(); it++){ - const bool tmp = augment(layered_nf, *it, max_depth); //force augment evaluation even if successful is already true - successful = successful || tmp; - } + 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) { +inline void Graph_matching::set_r(double r) { this->r = r; } -bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { -#ifdef DEBUG - std::cout << " augment" << std::endl << std::flush; -#endif -#ifdef DEBUG - std::cout << " u_start_index = " << u_start_index << std::endl << std::flush; -#endif +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::deque path; path.emplace_back(u_start_index); - // u_start is a point from U do { -#ifdef DEBUG - std::cout << " do" << std::endl << std::flush; - std::cout << " path.size = " << static_cast(path.size()) << std::endl << std::flush; -#endif if (static_cast(path.size()) > max_depth) { path.pop_back(); path.pop_back(); @@ -134,28 +133,22 @@ bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_ path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); } path.emplace_back(v_to_u.at(path.back())); -#ifdef DEBUG - std::cout << " v_to_u = " << path.back() << std::endl << std::flush; -#endif } 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; } -std::unique_ptr Graph_matching::layering() const { -#ifdef DEBUG - std::cout << " layering" << std::endl << std::flush; -#endif - bool end = false; - int layer = 0; +inline std::unique_ptr Graph_matching::layering() const { std::list u_vertices(unmatched_in_u); std::list v_vertices; - Neighbors_finder nf(g, r); - for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) + Neighbors_finder nf(r); + for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index) nf.add(v_point_index); - std::unique_ptr layered_nf(new Layered_neighbors_finder(g, r)); - while (!u_vertices.empty()) { + std::unique_ptr layered_nf(new Layered_neighbors_finder(r)); + for(int layer = 0; !u_vertices.empty(); layer++) { + // one layer is one step in the BFS for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) { std::unique_ptr< std::list > u_succ = std::move(nf.pull_all_near(*it)); for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) { @@ -163,27 +156,27 @@ std::unique_ptr Graph_matching::layering() const { v_vertices.emplace_back(*it); } } + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS u_vertices.clear(); - for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) { + 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(); - layer++; } return layered_nf; } -void Graph_matching::update(std::deque& path) { -#ifdef DEBUG - std::cout << " update" << std::endl << std::flush; -#endif +inline void Graph_matching::update(std::deque& 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; ++it; v_to_u[*it] = tmp; @@ -192,27 +185,28 @@ void Graph_matching::update(std::deque& path) { template double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { - Persistence_diagrams_graph g(diag1, diag2, e); - std::unique_ptr< std::vector > sd = std::move(g.sorted_distances()); + G::initialize(diag1, diag2, e); + std::unique_ptr< std::vector > sd = std::move(G::sorted_distances()); int idmin = 0; int idmax = sd->size() - 1; + // alpha can be modified, this will change the complexity double alpha = pow(sd->size(), 0.25); - Graph_matching m(g); - Graph_matching biggest_unperfect(g); + Graph_matching m; + Graph_matching biggest_unperfect; while (idmin != idmax) { - int pas = static_cast((idmax - idmin) / alpha); - m.set_r(sd->at(idmin + pas)); + int step = static_cast((idmax - idmin) / alpha); + m.set_r(sd->at(idmin + step)); while (m.multi_augment()); + // The above while compute a maximum matching (according to the r setted before) if (m.perfect()) { - idmax = idmin + pas; + idmax = idmin + step; m = biggest_unperfect; } else { biggest_unperfect = m; - idmin = idmin + pas + 1; + idmin = idmin + step + 1; } } - double b = sd->at(idmin); - return b; + return sd->at(idmin); } } // namespace bottleneck diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h index 1cf83393..e6b7f30d 100644 --- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h @@ -27,44 +27,50 @@ #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 { +/** \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: - 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; +public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Layered_neighbors_finder(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_diagrams_graph& g; - const double r; - std::vector neighbors_finder; +private: + const double r; + std::vector neighbors_finder; }; -Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r), neighbors_finder() { } +Layered_neighbors_finder::Layered_neighbors_finder(double r) : + 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 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(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 (neighbors_finder.size()) <= vlayer) - return null_point_index(); - return neighbors_finder.at(vlayer).pull_near(u_point_index); +inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { + if (static_cast (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(neighbors_finder.size()); +inline int Layered_neighbors_finder::vlayers_number() const { + return static_cast(neighbors_finder.size()); } } // namespace bottleneck diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h index 0fad9889..ac3f8600 100644 --- a/src/Bottleneck/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h @@ -32,19 +32,25 @@ 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. - +/** \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 { public: - Neighbors_finder(const Persistence_diagrams_graph& g, double r); + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(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::unique_ptr< std::list > 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 projections_f; @@ -52,43 +58,46 @@ private: bool contains(int v_point_index); }; -Neighbors_finder::Neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r), planar_neighbors_f(g, r), projections_f() { } +Neighbors_finder::Neighbors_finder(double r) : + r(r), planar_neighbors_f(r), projections_f() { } -/* inline */ void Neighbors_finder::add(int v_point_index) { - if (g.on_the_v_diagonal(v_point_index)) +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 */ void Neighbors_finder::remove(int v_point_index) { +inline void Neighbors_finder::remove(int v_point_index) { if(v_point_index == null_point_index()) return; - if (g.on_the_v_diagonal(v_point_index)) + if (G::on_the_v_diagonal(v_point_index)) projections_f.erase(v_point_index); else planar_neighbors_f.remove(v_point_index); } -/* inline */ bool Neighbors_finder::contains(int v_point_index) { +inline bool Neighbors_finder::contains(int v_point_index) { return planar_neighbors_f.contains(v_point_index) || (projections_f.count(v_point_index)>0); } -/* inline */ int Neighbors_finder::pull_near(int u_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.cbegin()!=projections_f.cend()) + int c = G::corresponding_point_in_v(u_point_index); + if (G::on_the_u_diagonal(u_point_index) && !projections_f.empty()) + //All projections are at distance 0 tmp = *projections_f.cbegin(); - else if (contains(c) && (g.distance(u_point_index, c) <= r)) + else if (contains(c) && (G::distance(u_point_index, c) <= r)) + //Is the query point near to its projection ? tmp = c; else + //Is the query point near to a V point in the plane ? tmp = planar_neighbors_f.pull_near(u_point_index); remove(tmp); return tmp; } -std::unique_ptr< std::list > Neighbors_finder::pull_all_near(int u_point_index) { +inline std::unique_ptr< std::list > Neighbors_finder::pull_all_near(int u_point_index) { std::unique_ptr< std::list > all_pull = std::move(planar_neighbors_f.pull_all_near(u_point_index)); int last_pull = pull_near(u_point_index); while (last_pull != null_point_index()) { diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h index 644c1cd8..66e43145 100644 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h @@ -26,52 +26,70 @@ #include #include #include -#include // for pair<> -#include // for max +#include +#include namespace Gudhi { namespace bottleneck { -// Diagram_point is the type of the persistence diagram's points -typedef std::pair Diagram_point; - -// Return the used index for encoding none of the points +/** \internal \brief Returns 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. - +/** \internal \brief Structure representing an euclidean bipartite graph containing + * the points from the two persistence diagrams (including the projections). + * + * \ingroup bottleneck_distance + */ 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). +public: + /** \internal \brief Initializer taking 2 Point (concept) ranges as parameters. */ template - Persistence_diagrams_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); - 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::unique_ptr< std::vector > sorted_distances(); - - private: - std::vector u; - std::vector v; - Diagram_point get_u_point(int u_point_index) const; - Diagram_point get_v_point(int v_point_index) const; + static void initialize(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 ? */ + static bool on_the_u_diagonal(int u_point_index); + /** \internal \brief Is the given point from V the projection of a point in U ? */ + static bool on_the_v_diagonal(int v_point_index); + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + static int corresponding_point_in_u(int v_point_index); + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + static int corresponding_point_in_v(int u_point_index); + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + static double distance(int u_point_index, int v_point_index); + /** \internal \brief Returns size = |U| = |V|. */ + static int size(); + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + static std::unique_ptr< std::vector > sorted_distances(); + /** \internal \brief Compare points regarding x coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ + struct Compare_x{bool operator() (const int point_index_1, const int point_index_2) const;}; + /** \internal \brief Compare points regarding y coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ + struct Compare_y{bool operator() (const int point_index_1, const int point_index_2) const;}; + +private: + /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ + typedef std::pair Internal_point; + static std::vector u; + static std::vector v; + static Internal_point get_u_point(int u_point_index); + static Internal_point get_v_point(int v_point_index); }; -/* inline */ int null_point_index() { +/** \internal \typedef \brief Alias used outside. */ +typedef Persistence_diagrams_graph G; + +// static initialization, seems to work but strange +std::vector Persistence_diagrams_graph::u = [] {return std::vector();}(); +std::vector Persistence_diagrams_graph::v = [] {return std::vector();}(); + +inline int null_point_index() { return -1; } template -Persistence_diagrams_graph::Persistence_diagrams_graph(const Persistence_diagram1 &diag1, - const Persistence_diagram2 &diag2, double e) - : u(), v() { +void Persistence_diagrams_graph::initialize(const Persistence_diagram1 &diag1, + const Persistence_diagram2 &diag2, double e){ + u.clear(); + v.clear(); for (auto it = diag1.cbegin(); it != diag1.cend(); ++it) if (it->second - it->first > e) u.emplace_back(*it); @@ -82,41 +100,37 @@ Persistence_diagrams_graph::Persistence_diagrams_graph(const Persistence_diagram 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 { +inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) { return u_point_index >= static_cast (u.size()); } -/* inline */ bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const { +inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) { return v_point_index >= static_cast (v.size()); } -/* inline */ int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const { +inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) { return on_the_v_diagonal(v_point_index) ? v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); } -/* inline */ int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const { +inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) { return on_the_u_diagonal(u_point_index) ? u_point_index - static_cast (u.size()) : u_point_index + static_cast (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 double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) { + 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.first - p_v.first), std::fabs(p_u.second - p_v.second)); } -/* inline */ int Persistence_diagrams_graph::size() const { +inline int Persistence_diagrams_graph::size() { return static_cast (u.size() + v.size()); } -/* inline */ std::unique_ptr< std::vector > Persistence_diagrams_graph::sorted_distances() { +inline std::unique_ptr< std::vector > Persistence_diagrams_graph::sorted_distances() { // could be optimized std::set sorted_distances; for (int u_point_index = 0; u_point_index < size(); ++u_point_index) @@ -126,20 +140,32 @@ Persistence_diagrams_graph::Persistence_diagrams_graph() return sd_up; } -/* inline */ Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const { +inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_u_point(int u_point_index) { 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)); + Internal_point projector = v.at(corresponding_point_in_v(u_point_index)); double x = (projector.first + projector.second) / 2; - return Diagram_point(x, x); + return Internal_point(x, x); } -/* inline */ Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const { +inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_v_point(int v_point_index) { 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)); + Internal_point projector = u.at(corresponding_point_in_u(v_point_index)); double x = (projector.first + projector.second) / 2; - return Diagram_point(x, x); + return Internal_point(x, x); +} + +inline bool Persistence_diagrams_graph::Compare_x::operator() (const int point_index_1, const int point_index_2) const{ + G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); + G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); + return p1.first < p2.first; +} + +inline bool Persistence_diagrams_graph::Compare_y::operator() (const int point_index_1, const int point_index_2) const{ + G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); + G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); + return p1.second < p2.second; } } // namespace bottleneck diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h index e558239d..4d820c7f 100644 --- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h @@ -27,88 +27,102 @@ #include #include "Persistence_diagrams_graph.h" +#include "Grid_cell.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. - +/** \internal \brief Structure used to find any point in V near (according to the planar distance) to a query point from U. + * + * V points have to be added manually using their index and before the first remove/pull. A neighbor pulled is automatically removed. but we can also + * remove points manually using their index. + * + * \ingroup bottleneck_distance + */ 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::unique_ptr< std::list > pull_all_near(int u_point_index); - - protected: - const Persistence_diagrams_graph& g; - const double r; +public: + /** \internal \brief Constructor TODO. */ + Abstract_planar_neighbors_finder(double r); + virtual ~Abstract_planar_neighbors_finder() = 0; + /** \internal \brief A point added will be possibly pulled. */ + virtual void add(int v_point_index) = 0; + /** \internal \brief A point manually removed will no longer be possibly pulled. */ + virtual void remove(int v_point_index) = 0; + /** \internal \brief Can the point given as parameter be returned ? */ + virtual bool contains(int v_point_index) const = 0; + /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ + virtual int pull_near(int u_point_index) = 0; + /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */ + virtual std::unique_ptr< std::list > pull_all_near(int u_point_index); + +protected: + const double r; }; - -// Naive_pnf is a nave implementation of Abstract_planar_neighbors_finder - +/** \internal \brief Naive_pnf is an naïve Abstract_planar_neighbors_finder implementation + * + * \ingroup bottleneck_distance + */ 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 candidates; +public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Naive_pnf(double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief A point manually removed will no longer be possibly pulled. */ + void remove(int v_point_index); + /** \internal \brief Can the point given as parameter be returned ? */ + bool contains(int v_point_index) const; + /** \internal \brief Provide 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); + +private: + std::set candidates; }; - -// Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation +/** \internal \typedef \brief 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() { } +Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) : + r(r) { } -/* inline */ std::unique_ptr< std::list > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { - std::unique_ptr< std::list > all_pull(new std::list); - 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 Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() {} + +inline std::unique_ptr< std::list > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { + std::unique_ptr< std::list > all_pull(new std::list); + 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() { } +Naive_pnf::Naive_pnf(double r) : + Abstract_planar_neighbors_finder(r), candidates() { } -/* inline */ void Naive_pnf::add(int v_point_index) { - candidates.emplace(v_point_index); +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 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 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(); +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 diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp index 77f45867..926c8430 100644 --- a/src/Bottleneck/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck/test/bottleneck_unit_test.cpp @@ -10,25 +10,6 @@ int n1 = 81; // a natural number >0 int n2 = 180; // a natural number >0 double upper_bound = 400.5; // any real >0 -std::unique_ptr random_graph_generator(){ - // Random construction - std::uniform_real_distribution unif(0.,upper_bound); - std::default_random_engine re; - std::vector< std::pair > v1, v2; - 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)); - } - return std::unique_ptr(new Persistence_diagrams_graph(v1, v2, 0.)); -} - - BOOST_AUTO_TEST_CASE(global){ std::uniform_real_distribution unif1(0.,upper_bound); std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); @@ -49,51 +30,64 @@ BOOST_AUTO_TEST_CASE(global){ BOOST_CHECK(bottleneck_distance(v1, v2) <= upper_bound/100.); } - -BOOST_AUTO_TEST_CASE(persistence_diagrams_graph) { - std::unique_ptr g = std::move(random_graph_generator()); - std::unique_ptr< std::vector > d = std::move(g->sorted_distances()); +BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ + // Random construction + std::uniform_real_distribution unif(0.,upper_bound); + std::default_random_engine re; + std::vector< std::pair > v1, v2; + 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)); + } + G::initialize(v1, v2, 0.); + std::unique_ptr< std::vector > d = std::move(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::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::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(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))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n1-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n2-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n2))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,(n1+n2)-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,0))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n1-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n2-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n2))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,(n1+n2)-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,0))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n1-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n2-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n2))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,(n1+n2)-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,0))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,(n1+n2)-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,0))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,(n1+n2)-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,0))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,(n1+n2)-1))==1); } + BOOST_AUTO_TEST_CASE(planar_neighbors_finder) { - std::unique_ptr g = std::move(random_graph_generator()); - Planar_neighbors_finder pnf = Planar_neighbors_finder(*g,1.); + Planar_neighbors_finder pnf = Planar_neighbors_finder(1.); for(int v_point_index=0; v_point_indexdistance(n2/2,v_point_index_1)<=1.)); + BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.)); BOOST_CHECK(!pnf.contains(v_point_index_1)); std::list l = *pnf.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.); + v = v && (G::distance(n2/2,*it)>1.); BOOST_CHECK(v); int v_point_index_2 = pnf.pull_near(n2/2); BOOST_CHECK(v_point_index_2 == -1); - pnf.add(v_point_index_1); - BOOST_CHECK(pnf.contains(v_point_index_1)); } BOOST_AUTO_TEST_CASE(neighbors_finder) { - std::unique_ptr g = std::move(random_graph_generator()); - Neighbors_finder nf = Neighbors_finder(*g,1.); + Neighbors_finder nf = Neighbors_finder(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.)); + BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.)); std::list 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.); + 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) { - std::unique_ptr g = std::move(random_graph_generator()); - Layered_neighbors_finder lnf = Layered_neighbors_finder(*g,1.); + Layered_neighbors_finder lnf = Layered_neighbors_finder(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.)); + 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.)); + 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) { - std::unique_ptr g = std::move(random_graph_generator()); - Graph_matching m1(*g); + Graph_matching m1; m1.set_r(0.); int e = 0; while (m1.multi_augment()) @@ -177,3 +165,7 @@ BOOST_AUTO_TEST_CASE(graph_matching) { BOOST_CHECK(m2.perfect()); BOOST_CHECK(!m1.perfect()); } + +BOOST_AUTO_TEST_CASE(grid_cell) { + +} -- cgit v1.2.3 From 6bc84d6e645cc4aada745fe779e985a2f62718e7 Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 29 Jun 2015 18:16:07 +0000 Subject: Ajout de Grid_cell.h git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@660 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: e1b4ccdcb482c992c444adc91d67531bfa3e90dc --- src/Bottleneck/include/gudhi/Graph_matching.h | 4 +- src/Bottleneck/include/gudhi/Grid_cell.h | 188 +++++++++++++++++++++ .../include/gudhi/Layered_neighbors_finder.h | 2 +- src/Bottleneck/include/gudhi/Neighbors_finder.h | 2 +- .../include/gudhi/Persistence_diagrams_graph.h | 59 ++++--- .../include/gudhi/Planar_neighbors_finder.h | 4 +- 6 files changed, 231 insertions(+), 28 deletions(-) create mode 100644 src/Bottleneck/include/gudhi/Grid_cell.h (limited to 'src/Bottleneck/include/gudhi') diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h index 17412e6c..2bcc6a61 100644 --- a/src/Bottleneck/include/gudhi/Graph_matching.h +++ b/src/Bottleneck/include/gudhi/Graph_matching.h @@ -74,13 +74,13 @@ private: void update(std::deque & path); }; -Graph_matching::Graph_matching() +inline Graph_matching::Graph_matching() : 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); } -Graph_matching& Graph_matching::operator=(const Graph_matching& m) { +inline 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; diff --git a/src/Bottleneck/include/gudhi/Grid_cell.h b/src/Bottleneck/include/gudhi/Grid_cell.h new file mode 100644 index 00000000..eee938cb --- /dev/null +++ b/src/Bottleneck/include/gudhi/Grid_cell.h @@ -0,0 +1,188 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ + +#include +#include +#include + +#include "Persistence_diagrams_graph.h" + +namespace Gudhi { + +namespace bottleneck { + +/** \internal \brief TODO + * + * \ingroup bottleneck_distance + */ +class Grid_cell { +public: + Grid_cell(double r); + void add(int v_point_index); + void remove(int v_point_index); + bool contains(int v_point_index) const; + int pull_center(); + int pull_xi(int u_point_index); + int pull_xd(int u_point_index); + int pull_yi(int u_point_index); + int pull_yd(int u_point_index); + int pull_xi_yi(int u_point_index); + int pull_xi_yd(int u_point_index); + int pull_xd_yi(int u_point_index); + int pull_xd_yd(int u_point_index); + +private: + double r; + std::set xi_order; + std::set yi_order; + struct Hidden_points_tree{int point; std::list hidden;}; + typedef std::map, G::Compare_x> Corner_tree; + Corner_tree xi_yi_order; + Corner_tree xi_yd_order; + Corner_tree xd_yi_order; + Corner_tree xd_yd_order; + void remove_aux(Corner_tree t, int v_point_index); + void build_xi_yi(); + void build_xi_yd(); + void build_xd_yi(); + void build_xd_yd(); +}; + + +inline Grid_cell::Grid_cell(double r) + : r(r), xi_order(G::Compare_x(r)), yi_order(G::Compare_y(r)), xi_yi_order(G::Compare_x(r)), + xi_yd_order(G::Compare_x(r)), xd_yi_order(G::Compare_x(r)), xd_yd_order(G::Compare_x(r)) {} + +inline void Grid_cell::add(int v_point_index){ + xi_order.emplace(v_point_index); +} + +inline bool Grid_cell::contains(int v_point_index) const{ + return (xi_order.count(v_point_index) > 0); +} + +inline void Grid_cell::remove_aux(Corner_tree t, int v_point_index){ + if(t.empty()) + return; + std::list hidden_points = t.at(v_point_index); + t.erase(v_point_index); + for(auto it = hidden_points.begin(); it != hidden_points.end(); ++it) + t.emplace(it->point,it->hidden); + +} + +inline void Grid_cell::remove(int v_point_index){ + xi_order.erase(v_point_index); + yi_order.erase(v_point_index); + remove_aux(xi_yi_order,v_point_index); + remove_aux(xi_yd_order,v_point_index); + remove_aux(xd_yi_order,v_point_index); + remove_aux(xd_yd_order,v_point_index); +} + +//factorization needed \/ \/ \/ + +inline int Grid_cell::pull_center(){ + if(xi_order.empty()) + return null_point_index(); + int v_point_index = *xi_order.begin(); + remove(v_point_index); + return v_point_index; +} + +inline int Grid_cell::pull_xi(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + int v_point_index = *xi_order.begin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_xd(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + int v_point_index = *xi_order.rbegin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_yi(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + if(yi_order.empty()) + for(auto it = xi_order.begin(); it!= xi_order.end(); ++it) + yi_order.emplace(*it); + int v_point_index = *yi_order.begin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_yd(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + if(yi_order.empty()) + for(auto it = xi_order.begin(); it!= xi_order.end(); ++it) + yi_order.emplace(*it); + int v_point_index = *yi_order.rbegin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_xi_yi(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + if(xi_yi_order.empty()) + build_xi_yi(); + auto it = xi_yi_order.upper_bound(u_point_index+G::size()); + if(it==xi_yi_order.cbegin()) //! + return null_point_index(); + it--; //! + int v_point_index = it->first; + for(auto it2=it->second.begin();it2!=it->second.end();it2++) + xi_yi_order.emplace(it2->point,it2->hidden); + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h index e6b7f30d..58805b86 100644 --- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h @@ -54,7 +54,7 @@ private: std::vector neighbors_finder; }; -Layered_neighbors_finder::Layered_neighbors_finder(double r) : +inline Layered_neighbors_finder::Layered_neighbors_finder(double r) : r(r), neighbors_finder() { } inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h index ac3f8600..be81877a 100644 --- a/src/Bottleneck/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h @@ -58,7 +58,7 @@ private: bool contains(int v_point_index); }; -Neighbors_finder::Neighbors_finder(double r) : +inline Neighbors_finder::Neighbors_finder(double r) : r(r), planar_neighbors_f(r), projections_f() { } inline void Neighbors_finder::add(int v_point_index) { diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h index 66e43145..aed328e2 100644 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h @@ -28,6 +28,7 @@ #include #include #include +#include namespace Gudhi { @@ -60,10 +61,10 @@ public: static int size(); /** \internal \brief Returns the O(n^2) sorted distances between the points. */ static std::unique_ptr< std::vector > sorted_distances(); - /** \internal \brief Compare points regarding x coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ - struct Compare_x{bool operator() (const int point_index_1, const int point_index_2) const;}; - /** \internal \brief Compare points regarding y coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ - struct Compare_y{bool operator() (const int point_index_1, const int point_index_2) const;}; + /** \internal \brief Compare points regarding x%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ + struct Compare_x{double r; Compare_x(double r); bool operator()(const int point_index_1, const int point_index_2) const;}; + /** \internal \brief Compare points regarding y%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ + struct Compare_y{double r; Compare_y(double r);bool operator()(const int point_index_1, const int point_index_2) const;}; private: /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ @@ -74,20 +75,20 @@ private: static Internal_point get_v_point(int v_point_index); }; -/** \internal \typedef \brief Alias used outside. */ +/** \internal \typedef \brief Shorter alias */ typedef Persistence_diagrams_graph G; // static initialization, seems to work but strange -std::vector Persistence_diagrams_graph::u = [] {return std::vector();}(); -std::vector Persistence_diagrams_graph::v = [] {return std::vector();}(); +std::vector G::u = [] {return std::vector();}(); +std::vector G::v = [] {return std::vector();}(); inline int null_point_index() { return -1; } template -void Persistence_diagrams_graph::initialize(const Persistence_diagram1 &diag1, - const Persistence_diagram2 &diag2, double e){ +inline void G::initialize(const Persistence_diagram1 &diag1, + const Persistence_diagram2 &diag2, double e){ u.clear(); v.clear(); for (auto it = diag1.cbegin(); it != diag1.cend(); ++it) @@ -100,25 +101,25 @@ void Persistence_diagrams_graph::initialize(const Persistence_diagram1 &diag1, swap(u, v); } -inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) { +inline bool G::on_the_u_diagonal(int u_point_index) { return u_point_index >= static_cast (u.size()); } -inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) { +inline bool G::on_the_v_diagonal(int v_point_index) { return v_point_index >= static_cast (v.size()); } -inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) { +inline int G::corresponding_point_in_u(int v_point_index) { return on_the_v_diagonal(v_point_index) ? v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); } -inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) { +inline int G::corresponding_point_in_v(int u_point_index) { return on_the_u_diagonal(u_point_index) ? u_point_index - static_cast (u.size()) : u_point_index + static_cast (v.size()); } -inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) { +inline double G::distance(int u_point_index, int v_point_index) { 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); @@ -126,11 +127,11 @@ inline double Persistence_diagrams_graph::distance(int u_point_index, int v_poin 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() { +inline int G::size() { return static_cast (u.size() + v.size()); } -inline std::unique_ptr< std::vector > Persistence_diagrams_graph::sorted_distances() { +inline std::unique_ptr< std::vector > G::sorted_distances() { // could be optimized std::set sorted_distances; for (int u_point_index = 0; u_point_index < size(); ++u_point_index) @@ -140,7 +141,7 @@ inline std::unique_ptr< std::vector > Persistence_diagrams_graph::sorted return sd_up; } -inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_u_point(int u_point_index) { +inline G::Internal_point G::get_u_point(int u_point_index) { 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)); @@ -148,7 +149,7 @@ inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::ge return Internal_point(x, x); } -inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_v_point(int v_point_index) { +inline G::Internal_point G::get_v_point(int v_point_index) { 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)); @@ -156,16 +157,30 @@ inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::ge return Internal_point(x, x); } -inline bool Persistence_diagrams_graph::Compare_x::operator() (const int point_index_1, const int point_index_2) const{ +G::Compare_x::Compare_x(double r) + : r(r){ } + +G::Compare_y::Compare_y(double r) + : r(r){ } + +inline bool G::Compare_x::operator()(const int point_index_1, const int point_index_2) const{ G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); - return p1.first < p2.first; + double x1 = fmod(p1.first,r); + double x2 = fmod(p2.first,r); + if(x1 == x2) + return point_index_1 > point_index_2; + return x1 < x2; } -inline bool Persistence_diagrams_graph::Compare_y::operator() (const int point_index_1, const int point_index_2) const{ +inline bool G::Compare_y::operator()(const int point_index_1, const int point_index_2) const{ G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); - return p1.second < p2.second; + double y1 = fmod(p1.second,r); + double y2 = fmod(p2.second,r); + if(y1 == y2) + return point_index_1 > point_index_2; + return y1 < y2; } } // namespace bottleneck diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h index 4d820c7f..e403735c 100644 --- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h @@ -85,7 +85,7 @@ private: typedef Naive_pnf Planar_neighbors_finder; -Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) : +inline Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) : r(r) { } inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() {} @@ -100,7 +100,7 @@ inline std::unique_ptr< std::list > Abstract_planar_neighbors_finder::pull_ return all_pull; } -Naive_pnf::Naive_pnf(double r) : +inline Naive_pnf::Naive_pnf(double r) : Abstract_planar_neighbors_finder(r), candidates() { } inline void Naive_pnf::add(int v_point_index) { -- cgit v1.2.3 From c0543ba5794592af4b0ea01d35710d44e206a8a3 Mon Sep 17 00:00:00 2001 From: fgodi Date: Tue, 30 Jun 2015 05:54:52 +0000 Subject: renommage git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@661 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 3316ed6b6f7d865605786547b62ee1e5d9118d19 --- .../concept/Persistence_diagram.h | 7 + src/Bipartite_graph_matching/example/example.cpp | 44 +++++ .../include/gudhi/Graph_matching.h | 216 +++++++++++++++++++++ .../include/gudhi/Grid_cell.h | 188 ++++++++++++++++++ .../include/gudhi/Layered_neighbors_finder.h | 80 ++++++++ .../include/gudhi/Neighbors_finder.h | 114 +++++++++++ .../include/gudhi/Persistence_diagrams_graph.h | 190 ++++++++++++++++++ .../include/gudhi/Planar_neighbors_finder.h | 132 +++++++++++++ src/Bipartite_graph_matching/test/CMakeLists.txt | 31 +++ src/Bipartite_graph_matching/test/README | 12 ++ .../test/bottleneck_unit_test.cpp | 171 ++++++++++++++++ src/Bottleneck/concept/Persistence_diagram.h | 7 - src/Bottleneck/example/example.cpp | 44 ----- src/Bottleneck/include/gudhi/Graph_matching.h | 216 --------------------- src/Bottleneck/include/gudhi/Grid_cell.h | 188 ------------------ .../include/gudhi/Layered_neighbors_finder.h | 80 -------- src/Bottleneck/include/gudhi/Neighbors_finder.h | 114 ----------- .../include/gudhi/Persistence_diagrams_graph.h | 190 ------------------ .../include/gudhi/Planar_neighbors_finder.h | 132 ------------- src/Bottleneck/test/CMakeLists.txt | 31 --- src/Bottleneck/test/README | 12 -- src/Bottleneck/test/bottleneck_unit_test.cpp | 171 ---------------- 22 files changed, 1185 insertions(+), 1185 deletions(-) create mode 100644 src/Bipartite_graph_matching/concept/Persistence_diagram.h create mode 100644 src/Bipartite_graph_matching/example/example.cpp create mode 100644 src/Bipartite_graph_matching/include/gudhi/Graph_matching.h create mode 100644 src/Bipartite_graph_matching/include/gudhi/Grid_cell.h create mode 100644 src/Bipartite_graph_matching/include/gudhi/Layered_neighbors_finder.h create mode 100644 src/Bipartite_graph_matching/include/gudhi/Neighbors_finder.h create mode 100644 src/Bipartite_graph_matching/include/gudhi/Persistence_diagrams_graph.h create mode 100644 src/Bipartite_graph_matching/include/gudhi/Planar_neighbors_finder.h create mode 100644 src/Bipartite_graph_matching/test/CMakeLists.txt create mode 100644 src/Bipartite_graph_matching/test/README create mode 100644 src/Bipartite_graph_matching/test/bottleneck_unit_test.cpp delete mode 100644 src/Bottleneck/concept/Persistence_diagram.h delete mode 100644 src/Bottleneck/example/example.cpp delete mode 100644 src/Bottleneck/include/gudhi/Graph_matching.h delete mode 100644 src/Bottleneck/include/gudhi/Grid_cell.h delete mode 100644 src/Bottleneck/include/gudhi/Layered_neighbors_finder.h delete mode 100644 src/Bottleneck/include/gudhi/Neighbors_finder.h delete mode 100644 src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h delete mode 100644 src/Bottleneck/include/gudhi/Planar_neighbors_finder.h delete mode 100644 src/Bottleneck/test/CMakeLists.txt delete mode 100644 src/Bottleneck/test/README delete mode 100644 src/Bottleneck/test/bottleneck_unit_test.cpp (limited to 'src/Bottleneck/include/gudhi') diff --git a/src/Bipartite_graph_matching/concept/Persistence_diagram.h b/src/Bipartite_graph_matching/concept/Persistence_diagram.h new file mode 100644 index 00000000..eaaf8bc5 --- /dev/null +++ b/src/Bipartite_graph_matching/concept/Persistence_diagram.h @@ -0,0 +1,7 @@ +typedef typename std::pair Diagram_point; + +struct Persistence_Diagram +{ + const_iterator cbegin() const; + const_iterator cend() const; +}; diff --git a/src/Bipartite_graph_matching/example/example.cpp b/src/Bipartite_graph_matching/example/example.cpp new file mode 100644 index 00000000..683f7723 --- /dev/null +++ b/src/Bipartite_graph_matching/example/example.cpp @@ -0,0 +1,44 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#include "../include/gudhi/Graph_matching.h" +#include + +using namespace Gudhi::bottleneck; + +int main() { + + std::vector< std::pair > v1, v2; + + v1.emplace_back(2.7,3.7); + v1.emplace_back(9.6,14); + v1.emplace_back(34.2,34.974); + + v2.emplace_back(2.8,4.45); + v2.emplace_back(9.5,14.1); + + + double b = bottleneck_distance(v1, v2); + + std::cout << "Bottleneck distance = " << b << std::endl; + +} diff --git a/src/Bipartite_graph_matching/include/gudhi/Graph_matching.h b/src/Bipartite_graph_matching/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..2bcc6a61 --- /dev/null +++ b/src/Bipartite_graph_matching/include/gudhi/Graph_matching.h @@ -0,0 +1,216 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ + +#include +#include +#include + +#include "Layered_neighbors_finder.h" + +namespace Gudhi { + +namespace bottleneck { + +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. + * + * + * + * \ingroup bottleneck_distance + */ +template +double bottleneck_distance(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); + +/** \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(); + /** \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: + double r; + /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ + std::vector v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::list unmatched_in_u; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + std::unique_ptr 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::deque & path); +}; + +inline Graph_matching::Graph_matching() + : 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) { + 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 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::deque path; + path.emplace_back(u_start_index); + do { + if (static_cast(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(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 std::unique_ptr Graph_matching::layering() const { + std::list u_vertices(unmatched_in_u); + std::list v_vertices; + Neighbors_finder nf(r); + for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index) + nf.add(v_point_index); + std::unique_ptr layered_nf(new Layered_neighbors_finder(r)); + for(int layer = 0; !u_vertices.empty(); layer++) { + // one layer is one step in the BFS + for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) { + std::unique_ptr< std::list > u_succ = std::move(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); + } + } + // 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::deque& 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; + ++it; + v_to_u[*it] = tmp; + } +} + +template +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { + G::initialize(diag1, diag2, e); + std::unique_ptr< std::vector > sd = std::move(G::sorted_distances()); + int idmin = 0; + int idmax = sd->size() - 1; + // alpha can be modified, this will change the complexity + double alpha = pow(sd->size(), 0.25); + Graph_matching m; + Graph_matching biggest_unperfect; + while (idmin != idmax) { + int step = static_cast((idmax - idmin) / alpha); + m.set_r(sd->at(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; + } + } + return sd->at(idmin); +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bipartite_graph_matching/include/gudhi/Grid_cell.h b/src/Bipartite_graph_matching/include/gudhi/Grid_cell.h new file mode 100644 index 00000000..eee938cb --- /dev/null +++ b/src/Bipartite_graph_matching/include/gudhi/Grid_cell.h @@ -0,0 +1,188 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ + +#include +#include +#include + +#include "Persistence_diagrams_graph.h" + +namespace Gudhi { + +namespace bottleneck { + +/** \internal \brief TODO + * + * \ingroup bottleneck_distance + */ +class Grid_cell { +public: + Grid_cell(double r); + void add(int v_point_index); + void remove(int v_point_index); + bool contains(int v_point_index) const; + int pull_center(); + int pull_xi(int u_point_index); + int pull_xd(int u_point_index); + int pull_yi(int u_point_index); + int pull_yd(int u_point_index); + int pull_xi_yi(int u_point_index); + int pull_xi_yd(int u_point_index); + int pull_xd_yi(int u_point_index); + int pull_xd_yd(int u_point_index); + +private: + double r; + std::set xi_order; + std::set yi_order; + struct Hidden_points_tree{int point; std::list hidden;}; + typedef std::map, G::Compare_x> Corner_tree; + Corner_tree xi_yi_order; + Corner_tree xi_yd_order; + Corner_tree xd_yi_order; + Corner_tree xd_yd_order; + void remove_aux(Corner_tree t, int v_point_index); + void build_xi_yi(); + void build_xi_yd(); + void build_xd_yi(); + void build_xd_yd(); +}; + + +inline Grid_cell::Grid_cell(double r) + : r(r), xi_order(G::Compare_x(r)), yi_order(G::Compare_y(r)), xi_yi_order(G::Compare_x(r)), + xi_yd_order(G::Compare_x(r)), xd_yi_order(G::Compare_x(r)), xd_yd_order(G::Compare_x(r)) {} + +inline void Grid_cell::add(int v_point_index){ + xi_order.emplace(v_point_index); +} + +inline bool Grid_cell::contains(int v_point_index) const{ + return (xi_order.count(v_point_index) > 0); +} + +inline void Grid_cell::remove_aux(Corner_tree t, int v_point_index){ + if(t.empty()) + return; + std::list hidden_points = t.at(v_point_index); + t.erase(v_point_index); + for(auto it = hidden_points.begin(); it != hidden_points.end(); ++it) + t.emplace(it->point,it->hidden); + +} + +inline void Grid_cell::remove(int v_point_index){ + xi_order.erase(v_point_index); + yi_order.erase(v_point_index); + remove_aux(xi_yi_order,v_point_index); + remove_aux(xi_yd_order,v_point_index); + remove_aux(xd_yi_order,v_point_index); + remove_aux(xd_yd_order,v_point_index); +} + +//factorization needed \/ \/ \/ + +inline int Grid_cell::pull_center(){ + if(xi_order.empty()) + return null_point_index(); + int v_point_index = *xi_order.begin(); + remove(v_point_index); + return v_point_index; +} + +inline int Grid_cell::pull_xi(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + int v_point_index = *xi_order.begin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_xd(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + int v_point_index = *xi_order.rbegin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_yi(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + if(yi_order.empty()) + for(auto it = xi_order.begin(); it!= xi_order.end(); ++it) + yi_order.emplace(*it); + int v_point_index = *yi_order.begin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_yd(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + if(yi_order.empty()) + for(auto it = xi_order.begin(); it!= xi_order.end(); ++it) + yi_order.emplace(*it); + int v_point_index = *yi_order.rbegin(); //! + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +inline int Grid_cell::pull_xi_yi(int u_point_index){ + if(xi_order.empty()) + return null_point_index(); + if(xi_yi_order.empty()) + build_xi_yi(); + auto it = xi_yi_order.upper_bound(u_point_index+G::size()); + if(it==xi_yi_order.cbegin()) //! + return null_point_index(); + it--; //! + int v_point_index = it->first; + for(auto it2=it->second.begin();it2!=it->second.end();it2++) + xi_yi_order.emplace(it2->point,it2->hidden); + if(G::distance(u_point_index,v_point_index)<=r){ + remove(v_point_index); + return v_point_index; + } + return null_point_index(); +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ diff --git a/src/Bipartite_graph_matching/include/gudhi/Layered_neighbors_finder.h b/src/Bipartite_graph_matching/include/gudhi/Layered_neighbors_finder.h new file mode 100644 index 00000000..58805b86 --- /dev/null +++ b/src/Bipartite_graph_matching/include/gudhi/Layered_neighbors_finder.h @@ -0,0 +1,80 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ + +#include + +#include "Neighbors_finder.h" + +namespace Gudhi { + +namespace bottleneck { + +/** \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(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 double r; + std::vector neighbors_finder; +}; + +inline Layered_neighbors_finder::Layered_neighbors_finder(double r) : + 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(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 (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(neighbors_finder.size()); +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ diff --git a/src/Bipartite_graph_matching/include/gudhi/Neighbors_finder.h b/src/Bipartite_graph_matching/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..be81877a --- /dev/null +++ b/src/Bipartite_graph_matching/include/gudhi/Neighbors_finder.h @@ -0,0 +1,114 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ + +#include +#include + +#include "Planar_neighbors_finder.h" + +namespace Gudhi { + +namespace bottleneck { + +/** \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 { +public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(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::unique_ptr< std::list > pull_all_near(int u_point_index); + +private: + const double r; + Planar_neighbors_finder planar_neighbors_f; + std::unordered_set projections_f; + void remove(int v_point_index); + bool contains(int v_point_index); +}; + +inline Neighbors_finder::Neighbors_finder(double r) : + r(r), planar_neighbors_f(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 void Neighbors_finder::remove(int v_point_index) { + if(v_point_index == null_point_index()) + return; + if (G::on_the_v_diagonal(v_point_index)) + projections_f.erase(v_point_index); + else + planar_neighbors_f.remove(v_point_index); +} + +inline bool Neighbors_finder::contains(int v_point_index) { + return planar_neighbors_f.contains(v_point_index) || (projections_f.count(v_point_index)>0); +} + +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()) + //All projections are at distance 0 + tmp = *projections_f.cbegin(); + else if (contains(c) && (G::distance(u_point_index, c) <= r)) + //Is the query point near to its projection ? + tmp = c; + else + //Is the query point near to a V point in the plane ? + tmp = planar_neighbors_f.pull_near(u_point_index); + remove(tmp); + return tmp; +} + +inline std::unique_ptr< std::list > Neighbors_finder::pull_all_near(int u_point_index) { + std::unique_ptr< std::list > all_pull = std::move(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/Bipartite_graph_matching/include/gudhi/Persistence_diagrams_graph.h b/src/Bipartite_graph_matching/include/gudhi/Persistence_diagrams_graph.h new file mode 100644 index 00000000..aed328e2 --- /dev/null +++ b/src/Bipartite_graph_matching/include/gudhi/Persistence_diagrams_graph.h @@ -0,0 +1,190 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ + +#include +#include +#include +#include +#include +#include + +namespace Gudhi { + +namespace bottleneck { + +/** \internal \brief Returns the used index for encoding none of the points */ +int null_point_index(); + +/** \internal \brief Structure representing an euclidean bipartite graph containing + * the points from the two persistence diagrams (including the projections). + * + * \ingroup bottleneck_distance + */ +class Persistence_diagrams_graph { +public: + /** \internal \brief Initializer taking 2 Point (concept) ranges as parameters. */ + template + static void initialize(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 ? */ + static bool on_the_u_diagonal(int u_point_index); + /** \internal \brief Is the given point from V the projection of a point in U ? */ + static bool on_the_v_diagonal(int v_point_index); + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + static int corresponding_point_in_u(int v_point_index); + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + static int corresponding_point_in_v(int u_point_index); + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + static double distance(int u_point_index, int v_point_index); + /** \internal \brief Returns size = |U| = |V|. */ + static int size(); + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + static std::unique_ptr< std::vector > sorted_distances(); + /** \internal \brief Compare points regarding x%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ + struct Compare_x{double r; Compare_x(double r); bool operator()(const int point_index_1, const int point_index_2) const;}; + /** \internal \brief Compare points regarding y%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ + struct Compare_y{double r; Compare_y(double r);bool operator()(const int point_index_1, const int point_index_2) const;}; + +private: + /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ + typedef std::pair Internal_point; + static std::vector u; + static std::vector v; + static Internal_point get_u_point(int u_point_index); + static Internal_point get_v_point(int v_point_index); +}; + +/** \internal \typedef \brief Shorter alias */ +typedef Persistence_diagrams_graph G; + +// static initialization, seems to work but strange +std::vector G::u = [] {return std::vector();}(); +std::vector G::v = [] {return std::vector();}(); + +inline int null_point_index() { + return -1; +} + +template +inline void G::initialize(const Persistence_diagram1 &diag1, + const Persistence_diagram2 &diag2, double e){ + u.clear(); + v.clear(); + 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); +} + +inline bool G::on_the_u_diagonal(int u_point_index) { + return u_point_index >= static_cast (u.size()); +} + +inline bool G::on_the_v_diagonal(int v_point_index) { + return v_point_index >= static_cast (v.size()); +} + +inline int G::corresponding_point_in_u(int v_point_index) { + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); +} + +inline int G::corresponding_point_in_v(int u_point_index) { + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast (u.size()) : u_point_index + static_cast (v.size()); +} + +inline double G::distance(int u_point_index, int v_point_index) { + 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.first - p_v.first), std::fabs(p_u.second - p_v.second)); +} + +inline int G::size() { + return static_cast (u.size() + v.size()); +} + +inline std::unique_ptr< std::vector > G::sorted_distances() { + // could be optimized + std::set 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)); + std::unique_ptr< std::vector > sd_up(new std::vector(sorted_distances.cbegin(), sorted_distances.cend())); + return sd_up; +} + +inline G::Internal_point G::get_u_point(int u_point_index) { + 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 x = (projector.first + projector.second) / 2; + return Internal_point(x, x); +} + +inline G::Internal_point G::get_v_point(int v_point_index) { + 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 x = (projector.first + projector.second) / 2; + return Internal_point(x, x); +} + +G::Compare_x::Compare_x(double r) + : r(r){ } + +G::Compare_y::Compare_y(double r) + : r(r){ } + +inline bool G::Compare_x::operator()(const int point_index_1, const int point_index_2) const{ + G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); + G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); + double x1 = fmod(p1.first,r); + double x2 = fmod(p2.first,r); + if(x1 == x2) + return point_index_1 > point_index_2; + return x1 < x2; +} + +inline bool G::Compare_y::operator()(const int point_index_1, const int point_index_2) const{ + G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); + G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); + double y1 = fmod(p1.second,r); + double y2 = fmod(p2.second,r); + if(y1 == y2) + return point_index_1 > point_index_2; + return y1 < y2; +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ diff --git a/src/Bipartite_graph_matching/include/gudhi/Planar_neighbors_finder.h b/src/Bipartite_graph_matching/include/gudhi/Planar_neighbors_finder.h new file mode 100644 index 00000000..e403735c --- /dev/null +++ b/src/Bipartite_graph_matching/include/gudhi/Planar_neighbors_finder.h @@ -0,0 +1,132 @@ +/* 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 Sophia-Antipolis (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 . + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ + +#include +#include + +#include "Persistence_diagrams_graph.h" +#include "Grid_cell.h" + +namespace Gudhi { + +namespace bottleneck { + +/** \internal \brief Structure used to find any point in V near (according to the planar distance) to a query point from U. + * + * V points have to be added manually using their index and before the first remove/pull. A neighbor pulled is automatically removed. but we can also + * remove points manually using their index. + * + * \ingroup bottleneck_distance + */ +class Abstract_planar_neighbors_finder { +public: + /** \internal \brief Constructor TODO. */ + Abstract_planar_neighbors_finder(double r); + virtual ~Abstract_planar_neighbors_finder() = 0; + /** \internal \brief A point added will be possibly pulled. */ + virtual void add(int v_point_index) = 0; + /** \internal \brief A point manually removed will no longer be possibly pulled. */ + virtual void remove(int v_point_index) = 0; + /** \internal \brief Can the point given as parameter be returned ? */ + virtual bool contains(int v_point_index) const = 0; + /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ + virtual int pull_near(int u_point_index) = 0; + /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */ + virtual std::unique_ptr< std::list > pull_all_near(int u_point_index); + +protected: + const double r; +}; + +/** \internal \brief Naive_pnf is an naïve Abstract_planar_neighbors_finder implementation + * + * \ingroup bottleneck_distance + */ +class Naive_pnf : public Abstract_planar_neighbors_finder { +public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Naive_pnf(double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief A point manually removed will no longer be possibly pulled. */ + void remove(int v_point_index); + /** \internal \brief Can the point given as parameter be returned ? */ + bool contains(int v_point_index) const; + /** \internal \brief Provide 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); + +private: + std::set candidates; +}; + +/** \internal \typedef \brief Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation. */ +typedef Naive_pnf Planar_neighbors_finder; + + +inline Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) : + r(r) { } + +inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() {} + +inline std::unique_ptr< std::list > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { + std::unique_ptr< std::list > all_pull(new std::list); + 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 Naive_pnf::Naive_pnf(double r) : + Abstract_planar_neighbors_finder(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/Bipartite_graph_matching/test/CMakeLists.txt b/src/Bipartite_graph_matching/test/CMakeLists.txt new file mode 100644 index 00000000..8e4371b8 --- /dev/null +++ b/src/Bipartite_graph_matching/test/CMakeLists.txt @@ -0,0 +1,31 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIBottleneckUT) + +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") +endif() + message("CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}") + message("CMAKE_CXX_FLAGS_DEBUG = ${CMAKE_CXX_FLAGS_DEBUG}") + message("CMAKE_CXX_FLAGS_RELEASE = ${CMAKE_CXX_FLAGS_RELEASE}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic -Wsign-compare") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") + +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/Bipartite_graph_matching/test/README b/src/Bipartite_graph_matching/test/README new file mode 100644 index 00000000..0e7b8673 --- /dev/null +++ b/src/Bipartite_graph_matching/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./BottleneckUnitTest --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Bipartite_graph_matching/test/bottleneck_unit_test.cpp b/src/Bipartite_graph_matching/test/bottleneck_unit_test.cpp new file mode 100644 index 00000000..926c8430 --- /dev/null +++ b/src/Bipartite_graph_matching/test/bottleneck_unit_test.cpp @@ -0,0 +1,171 @@ +#define BOOST_TEST_MODULE bottleneck test + +#include +#include +#include "../include/gudhi/Graph_matching.h" + +using namespace Gudhi::bottleneck; + +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 400.5; // any real >0 + +BOOST_AUTO_TEST_CASE(global){ + std::uniform_real_distribution unif1(0.,upper_bound); + std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); + std::default_random_engine re; + std::vector< std::pair > 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) <= upper_bound/100.); +} + +BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ + // Random construction + std::uniform_real_distribution unif(0.,upper_bound); + std::default_random_engine re; + std::vector< std::pair > v1, v2; + 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)); + } + G::initialize(v1, v2, 0.); + std::unique_ptr< std::vector > d = std::move(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))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,(n1+n2)-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,0))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,(n1+n2)-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,0))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2-1))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2))==1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,(n1+n2)-1))==1); +} + + +BOOST_AUTO_TEST_CASE(planar_neighbors_finder) { + Planar_neighbors_finder pnf = Planar_neighbors_finder(1.); + for(int v_point_index=0; v_point_index l = *pnf.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 = pnf.pull_near(n2/2); + BOOST_CHECK(v_point_index_2 == -1); +} + + +BOOST_AUTO_TEST_CASE(neighbors_finder) { + Neighbors_finder nf = Neighbors_finder(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::list 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) { + Layered_neighbors_finder lnf = Layered_neighbors_finder(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) { + Graph_matching m1; + m1.set_r(0.); + int e = 0; + while (m1.multi_augment()) + ++e; + 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(grid_cell) { + +} 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 Diagram_point; - -struct Persistence_Diagram -{ - const_iterator cbegin() const; - const_iterator cend() const; -}; diff --git a/src/Bottleneck/example/example.cpp b/src/Bottleneck/example/example.cpp deleted file mode 100644 index 683f7723..00000000 --- a/src/Bottleneck/example/example.cpp +++ /dev/null @@ -1,44 +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 Sophia-Antipolis (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 . - */ - -#include "../include/gudhi/Graph_matching.h" -#include - -using namespace Gudhi::bottleneck; - -int main() { - - std::vector< std::pair > v1, v2; - - v1.emplace_back(2.7,3.7); - v1.emplace_back(9.6,14); - v1.emplace_back(34.2,34.974); - - v2.emplace_back(2.8,4.45); - v2.emplace_back(9.5,14.1); - - - double b = bottleneck_distance(v1, v2); - - std::cout << "Bottleneck distance = " << b << std::endl; - -} diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h deleted file mode 100644 index 2bcc6a61..00000000 --- a/src/Bottleneck/include/gudhi/Graph_matching.h +++ /dev/null @@ -1,216 +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 Sophia-Antipolis (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 . - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ - -#include -#include -#include - -#include "Layered_neighbors_finder.h" - -namespace Gudhi { - -namespace bottleneck { - -/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. - * - * - * - * \ingroup bottleneck_distance - */ -template -double bottleneck_distance(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); - -/** \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(); - /** \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: - double r; - /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ - std::vector v_to_u; - /** \internal \brief All the unmatched points in U. */ - std::list unmatched_in_u; - - /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ - std::unique_ptr 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::deque & path); -}; - -inline Graph_matching::Graph_matching() - : 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) { - 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 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::deque path; - path.emplace_back(u_start_index); - do { - if (static_cast(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(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 std::unique_ptr Graph_matching::layering() const { - std::list u_vertices(unmatched_in_u); - std::list v_vertices; - Neighbors_finder nf(r); - for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index) - nf.add(v_point_index); - std::unique_ptr layered_nf(new Layered_neighbors_finder(r)); - for(int layer = 0; !u_vertices.empty(); layer++) { - // one layer is one step in the BFS - for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) { - std::unique_ptr< std::list > u_succ = std::move(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); - } - } - // 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::deque& 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; - ++it; - v_to_u[*it] = tmp; - } -} - -template -double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { - G::initialize(diag1, diag2, e); - std::unique_ptr< std::vector > sd = std::move(G::sorted_distances()); - int idmin = 0; - int idmax = sd->size() - 1; - // alpha can be modified, this will change the complexity - double alpha = pow(sd->size(), 0.25); - Graph_matching m; - Graph_matching biggest_unperfect; - while (idmin != idmax) { - int step = static_cast((idmax - idmin) / alpha); - m.set_r(sd->at(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; - } - } - return sd->at(idmin); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bottleneck/include/gudhi/Grid_cell.h b/src/Bottleneck/include/gudhi/Grid_cell.h deleted file mode 100644 index eee938cb..00000000 --- a/src/Bottleneck/include/gudhi/Grid_cell.h +++ /dev/null @@ -1,188 +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 Sophia-Antipolis (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 . - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_ - -#include -#include -#include - -#include "Persistence_diagrams_graph.h" - -namespace Gudhi { - -namespace bottleneck { - -/** \internal \brief TODO - * - * \ingroup bottleneck_distance - */ -class Grid_cell { -public: - Grid_cell(double r); - void add(int v_point_index); - void remove(int v_point_index); - bool contains(int v_point_index) const; - int pull_center(); - int pull_xi(int u_point_index); - int pull_xd(int u_point_index); - int pull_yi(int u_point_index); - int pull_yd(int u_point_index); - int pull_xi_yi(int u_point_index); - int pull_xi_yd(int u_point_index); - int pull_xd_yi(int u_point_index); - int pull_xd_yd(int u_point_index); - -private: - double r; - std::set xi_order; - std::set yi_order; - struct Hidden_points_tree{int point; std::list hidden;}; - typedef std::map, G::Compare_x> Corner_tree; - Corner_tree xi_yi_order; - Corner_tree xi_yd_order; - Corner_tree xd_yi_order; - Corner_tree xd_yd_order; - void remove_aux(Corner_tree t, int v_point_index); - void build_xi_yi(); - void build_xi_yd(); - void build_xd_yi(); - void build_xd_yd(); -}; - - -inline Grid_cell::Grid_cell(double r) - : r(r), xi_order(G::Compare_x(r)), yi_order(G::Compare_y(r)), xi_yi_order(G::Compare_x(r)), - xi_yd_order(G::Compare_x(r)), xd_yi_order(G::Compare_x(r)), xd_yd_order(G::Compare_x(r)) {} - -inline void Grid_cell::add(int v_point_index){ - xi_order.emplace(v_point_index); -} - -inline bool Grid_cell::contains(int v_point_index) const{ - return (xi_order.count(v_point_index) > 0); -} - -inline void Grid_cell::remove_aux(Corner_tree t, int v_point_index){ - if(t.empty()) - return; - std::list hidden_points = t.at(v_point_index); - t.erase(v_point_index); - for(auto it = hidden_points.begin(); it != hidden_points.end(); ++it) - t.emplace(it->point,it->hidden); - -} - -inline void Grid_cell::remove(int v_point_index){ - xi_order.erase(v_point_index); - yi_order.erase(v_point_index); - remove_aux(xi_yi_order,v_point_index); - remove_aux(xi_yd_order,v_point_index); - remove_aux(xd_yi_order,v_point_index); - remove_aux(xd_yd_order,v_point_index); -} - -//factorization needed \/ \/ \/ - -inline int Grid_cell::pull_center(){ - if(xi_order.empty()) - return null_point_index(); - int v_point_index = *xi_order.begin(); - remove(v_point_index); - return v_point_index; -} - -inline int Grid_cell::pull_xi(int u_point_index){ - if(xi_order.empty()) - return null_point_index(); - int v_point_index = *xi_order.begin(); //! - if(G::distance(u_point_index,v_point_index)<=r){ - remove(v_point_index); - return v_point_index; - } - return null_point_index(); -} - -inline int Grid_cell::pull_xd(int u_point_index){ - if(xi_order.empty()) - return null_point_index(); - int v_point_index = *xi_order.rbegin(); //! - if(G::distance(u_point_index,v_point_index)<=r){ - remove(v_point_index); - return v_point_index; - } - return null_point_index(); -} - -inline int Grid_cell::pull_yi(int u_point_index){ - if(xi_order.empty()) - return null_point_index(); - if(yi_order.empty()) - for(auto it = xi_order.begin(); it!= xi_order.end(); ++it) - yi_order.emplace(*it); - int v_point_index = *yi_order.begin(); //! - if(G::distance(u_point_index,v_point_index)<=r){ - remove(v_point_index); - return v_point_index; - } - return null_point_index(); -} - -inline int Grid_cell::pull_yd(int u_point_index){ - if(xi_order.empty()) - return null_point_index(); - if(yi_order.empty()) - for(auto it = xi_order.begin(); it!= xi_order.end(); ++it) - yi_order.emplace(*it); - int v_point_index = *yi_order.rbegin(); //! - if(G::distance(u_point_index,v_point_index)<=r){ - remove(v_point_index); - return v_point_index; - } - return null_point_index(); -} - -inline int Grid_cell::pull_xi_yi(int u_point_index){ - if(xi_order.empty()) - return null_point_index(); - if(xi_yi_order.empty()) - build_xi_yi(); - auto it = xi_yi_order.upper_bound(u_point_index+G::size()); - if(it==xi_yi_order.cbegin()) //! - return null_point_index(); - it--; //! - int v_point_index = it->first; - for(auto it2=it->second.begin();it2!=it->second.end();it2++) - xi_yi_order.emplace(it2->point,it2->hidden); - if(G::distance(u_point_index,v_point_index)<=r){ - remove(v_point_index); - return v_point_index; - } - return null_point_index(); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_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 58805b86..00000000 --- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h +++ /dev/null @@ -1,80 +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 Sophia-Antipolis (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 . - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ - -#include - -#include "Neighbors_finder.h" - -namespace Gudhi { - -namespace bottleneck { - -/** \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(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 double r; - std::vector neighbors_finder; -}; - -inline Layered_neighbors_finder::Layered_neighbors_finder(double r) : - 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(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 (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(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 be81877a..00000000 --- a/src/Bottleneck/include/gudhi/Neighbors_finder.h +++ /dev/null @@ -1,114 +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 Sophia-Antipolis (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 . - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ - -#include -#include - -#include "Planar_neighbors_finder.h" - -namespace Gudhi { - -namespace bottleneck { - -/** \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 { -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Neighbors_finder(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::unique_ptr< std::list > pull_all_near(int u_point_index); - -private: - const double r; - Planar_neighbors_finder planar_neighbors_f; - std::unordered_set projections_f; - void remove(int v_point_index); - bool contains(int v_point_index); -}; - -inline Neighbors_finder::Neighbors_finder(double r) : - r(r), planar_neighbors_f(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 void Neighbors_finder::remove(int v_point_index) { - if(v_point_index == null_point_index()) - return; - if (G::on_the_v_diagonal(v_point_index)) - projections_f.erase(v_point_index); - else - planar_neighbors_f.remove(v_point_index); -} - -inline bool Neighbors_finder::contains(int v_point_index) { - return planar_neighbors_f.contains(v_point_index) || (projections_f.count(v_point_index)>0); -} - -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()) - //All projections are at distance 0 - tmp = *projections_f.cbegin(); - else if (contains(c) && (G::distance(u_point_index, c) <= r)) - //Is the query point near to its projection ? - tmp = c; - else - //Is the query point near to a V point in the plane ? - tmp = planar_neighbors_f.pull_near(u_point_index); - remove(tmp); - return tmp; -} - -inline std::unique_ptr< std::list > Neighbors_finder::pull_all_near(int u_point_index) { - std::unique_ptr< std::list > all_pull = std::move(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 aed328e2..00000000 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ /dev/null @@ -1,190 +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 Sophia-Antipolis (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 . - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ - -#include -#include -#include -#include -#include -#include - -namespace Gudhi { - -namespace bottleneck { - -/** \internal \brief Returns the used index for encoding none of the points */ -int null_point_index(); - -/** \internal \brief Structure representing an euclidean bipartite graph containing - * the points from the two persistence diagrams (including the projections). - * - * \ingroup bottleneck_distance - */ -class Persistence_diagrams_graph { -public: - /** \internal \brief Initializer taking 2 Point (concept) ranges as parameters. */ - template - static void initialize(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 ? */ - static bool on_the_u_diagonal(int u_point_index); - /** \internal \brief Is the given point from V the projection of a point in U ? */ - static bool on_the_v_diagonal(int v_point_index); - /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ - static int corresponding_point_in_u(int v_point_index); - /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ - static int corresponding_point_in_v(int u_point_index); - /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ - static double distance(int u_point_index, int v_point_index); - /** \internal \brief Returns size = |U| = |V|. */ - static int size(); - /** \internal \brief Returns the O(n^2) sorted distances between the points. */ - static std::unique_ptr< std::vector > sorted_distances(); - /** \internal \brief Compare points regarding x%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ - struct Compare_x{double r; Compare_x(double r); bool operator()(const int point_index_1, const int point_index_2) const;}; - /** \internal \brief Compare points regarding y%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */ - struct Compare_y{double r; Compare_y(double r);bool operator()(const int point_index_1, const int point_index_2) const;}; - -private: - /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ - typedef std::pair Internal_point; - static std::vector u; - static std::vector v; - static Internal_point get_u_point(int u_point_index); - static Internal_point get_v_point(int v_point_index); -}; - -/** \internal \typedef \brief Shorter alias */ -typedef Persistence_diagrams_graph G; - -// static initialization, seems to work but strange -std::vector G::u = [] {return std::vector();}(); -std::vector G::v = [] {return std::vector();}(); - -inline int null_point_index() { - return -1; -} - -template -inline void G::initialize(const Persistence_diagram1 &diag1, - const Persistence_diagram2 &diag2, double e){ - u.clear(); - v.clear(); - 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); -} - -inline bool G::on_the_u_diagonal(int u_point_index) { - return u_point_index >= static_cast (u.size()); -} - -inline bool G::on_the_v_diagonal(int v_point_index) { - return v_point_index >= static_cast (v.size()); -} - -inline int G::corresponding_point_in_u(int v_point_index) { - return on_the_v_diagonal(v_point_index) ? - v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); -} - -inline int G::corresponding_point_in_v(int u_point_index) { - return on_the_u_diagonal(u_point_index) ? - u_point_index - static_cast (u.size()) : u_point_index + static_cast (v.size()); -} - -inline double G::distance(int u_point_index, int v_point_index) { - 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.first - p_v.first), std::fabs(p_u.second - p_v.second)); -} - -inline int G::size() { - return static_cast (u.size() + v.size()); -} - -inline std::unique_ptr< std::vector > G::sorted_distances() { - // could be optimized - std::set 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)); - std::unique_ptr< std::vector > sd_up(new std::vector(sorted_distances.cbegin(), sorted_distances.cend())); - return sd_up; -} - -inline G::Internal_point G::get_u_point(int u_point_index) { - 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 x = (projector.first + projector.second) / 2; - return Internal_point(x, x); -} - -inline G::Internal_point G::get_v_point(int v_point_index) { - 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 x = (projector.first + projector.second) / 2; - return Internal_point(x, x); -} - -G::Compare_x::Compare_x(double r) - : r(r){ } - -G::Compare_y::Compare_y(double r) - : r(r){ } - -inline bool G::Compare_x::operator()(const int point_index_1, const int point_index_2) const{ - G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); - G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); - double x1 = fmod(p1.first,r); - double x2 = fmod(p2.first,r); - if(x1 == x2) - return point_index_1 > point_index_2; - return x1 < x2; -} - -inline bool G::Compare_y::operator()(const int point_index_1, const int point_index_2) const{ - G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size()); - G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size()); - double y1 = fmod(p1.second,r); - double y2 = fmod(p2.second,r); - if(y1 == y2) - return point_index_1 > point_index_2; - return y1 < y2; -} - -} // 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 e403735c..00000000 --- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h +++ /dev/null @@ -1,132 +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 Sophia-Antipolis (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 . - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ - -#include -#include - -#include "Persistence_diagrams_graph.h" -#include "Grid_cell.h" - -namespace Gudhi { - -namespace bottleneck { - -/** \internal \brief Structure used to find any point in V near (according to the planar distance) to a query point from U. - * - * V points have to be added manually using their index and before the first remove/pull. A neighbor pulled is automatically removed. but we can also - * remove points manually using their index. - * - * \ingroup bottleneck_distance - */ -class Abstract_planar_neighbors_finder { -public: - /** \internal \brief Constructor TODO. */ - Abstract_planar_neighbors_finder(double r); - virtual ~Abstract_planar_neighbors_finder() = 0; - /** \internal \brief A point added will be possibly pulled. */ - virtual void add(int v_point_index) = 0; - /** \internal \brief A point manually removed will no longer be possibly pulled. */ - virtual void remove(int v_point_index) = 0; - /** \internal \brief Can the point given as parameter be returned ? */ - virtual bool contains(int v_point_index) const = 0; - /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ - virtual int pull_near(int u_point_index) = 0; - /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */ - virtual std::unique_ptr< std::list > pull_all_near(int u_point_index); - -protected: - const double r; -}; - -/** \internal \brief Naive_pnf is an naïve Abstract_planar_neighbors_finder implementation - * - * \ingroup bottleneck_distance - */ -class Naive_pnf : public Abstract_planar_neighbors_finder { -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Naive_pnf(double r); - /** \internal \brief A point added will be possibly pulled. */ - void add(int v_point_index); - /** \internal \brief A point manually removed will no longer be possibly pulled. */ - void remove(int v_point_index); - /** \internal \brief Can the point given as parameter be returned ? */ - bool contains(int v_point_index) const; - /** \internal \brief Provide 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); - -private: - std::set candidates; -}; - -/** \internal \typedef \brief Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation. */ -typedef Naive_pnf Planar_neighbors_finder; - - -inline Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) : - r(r) { } - -inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() {} - -inline std::unique_ptr< std::list > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { - std::unique_ptr< std::list > all_pull(new std::list); - 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 Naive_pnf::Naive_pnf(double r) : - Abstract_planar_neighbors_finder(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 8e4371b8..00000000 --- a/src/Bottleneck/test/CMakeLists.txt +++ /dev/null @@ -1,31 +0,0 @@ -cmake_minimum_required(VERSION 2.6) -project(GUDHIBottleneckUT) - -if (GCOVR_PATH) - # for gcovr to make coverage reports - Corbera Jenkins plugin - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") -endif() -if (GPROF_PATH) - # for gprof to make coverage reports - Jenkins - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") -endif() - message("CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}") - message("CMAKE_CXX_FLAGS_DEBUG = ${CMAKE_CXX_FLAGS_DEBUG}") - message("CMAKE_CXX_FLAGS_RELEASE = ${CMAKE_CXX_FLAGS_RELEASE}") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic -Wsign-compare") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") - -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/README b/src/Bottleneck/test/README deleted file mode 100644 index 0e7b8673..00000000 --- a/src/Bottleneck/test/README +++ /dev/null @@ -1,12 +0,0 @@ -To compile: -*********** - -cmake . -make - -To launch with details: -*********************** - -./BottleneckUnitTest --report_level=detailed --log_level=all - - ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp deleted file mode 100644 index 926c8430..00000000 --- a/src/Bottleneck/test/bottleneck_unit_test.cpp +++ /dev/null @@ -1,171 +0,0 @@ -#define BOOST_TEST_MODULE bottleneck test - -#include -#include -#include "../include/gudhi/Graph_matching.h" - -using namespace Gudhi::bottleneck; - -int n1 = 81; // a natural number >0 -int n2 = 180; // a natural number >0 -double upper_bound = 400.5; // any real >0 - -BOOST_AUTO_TEST_CASE(global){ - std::uniform_real_distribution unif1(0.,upper_bound); - std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); - std::default_random_engine re; - std::vector< std::pair > 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) <= upper_bound/100.); -} - -BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ - // Random construction - std::uniform_real_distribution unif(0.,upper_bound); - std::default_random_engine re; - std::vector< std::pair > v1, v2; - 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)); - } - G::initialize(v1, v2, 0.); - std::unique_ptr< std::vector > d = std::move(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))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,(n1+n2)-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,0))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,(n1+n2)-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,0))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2-1))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2))==1); - BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,(n1+n2)-1))==1); -} - - -BOOST_AUTO_TEST_CASE(planar_neighbors_finder) { - Planar_neighbors_finder pnf = Planar_neighbors_finder(1.); - for(int v_point_index=0; v_point_index l = *pnf.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 = pnf.pull_near(n2/2); - BOOST_CHECK(v_point_index_2 == -1); -} - - -BOOST_AUTO_TEST_CASE(neighbors_finder) { - Neighbors_finder nf = Neighbors_finder(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::list 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) { - Layered_neighbors_finder lnf = Layered_neighbors_finder(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) { - Graph_matching m1; - m1.set_r(0.); - int e = 0; - while (m1.multi_augment()) - ++e; - 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(grid_cell) { - -} -- cgit v1.2.3