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/example/random_diagrams.cpp | 2 +- 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 +- src/Bottleneck/test/bottleneck_unit_test.cpp | 122 ++++++++++++++++++--- 7 files changed, 192 insertions(+), 105 deletions(-) diff --git a/src/Bottleneck/example/random_diagrams.cpp b/src/Bottleneck/example/random_diagrams.cpp index 71f152a6..a9d0968d 100644 --- a/src/Bottleneck/example/random_diagrams.cpp +++ b/src/Bottleneck/example/random_diagrams.cpp @@ -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/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" diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp index 068b8690..ea7f51f2 100644 --- a/src/Bottleneck/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck/test/bottleneck_unit_test.cpp @@ -1,26 +1,112 @@ #define BOOST_TEST_MODULE bottleneck test #include - +#include #include "gudhi/Graph_matching.h" -#include using namespace Gudhi::bottleneck; -BOOST_AUTO_TEST_CASE(random_diagrams) { - int n = 100; - // Random construction - 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)); - } - // v1 and v2 are persistence diagrams containing each 100 randoms points. - double b = bottleneck_distance(v1, v2, 0); - // - std::cout << b << std::endl; - const double EXPECTED_DISTANCE = 98.5; - BOOST_CHECK(b == EXPECTED_DISTANCE); +Persistence_diagrams_graph* random_graph_generator(){ + int n1 = 100; + int n2 = 120; + double upper_bound = 80; + // 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 new Persistence_diagrams_graph(v1, v2, 0.); +} + +BOOST_AUTO_TEST_CASE(global){ + int n = 100; + // Random construction + 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)); + } + // + BOOST_CHECK(bottleneck_distance(v1, v2, 1.) == 98); } + + +BOOST_AUTO_TEST_CASE(persistence_diagrams_graph) { + Persistence_diagrams_graph* g = random_graph_generator(); + std::vector* d = 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->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->size()==220); + // + BOOST_CHECK((int) d->size() <= 220*220 - 100*120 + 1); // could be more strict + 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_AUTO_TEST_CASE(planar_nf) { + Persistence_diagrams_graph* g = 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) + pnf.add(v_point_index); + pnf.remove(0); + pnf.remove(1); + // + BOOST_CHECK(!pnf.contains(0)); + BOOST_CHECK(!pnf.contains(1)); + BOOST_CHECK(pnf.contains(2)); + BOOST_CHECK(!pnf.contains(3)); + // + int v_point_index_1 = pnf.pull_near(120/2); + BOOST_CHECK((v_point_index_1 == -1) || (g->distance(120/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)); + pnf.add(v_point_index_1); + BOOST_CHECK(pnf.contains(v_point_index_1)); +} + -- cgit v1.2.3