diff options
Diffstat (limited to 'src')
11 files changed, 1070 insertions, 1 deletions
diff --git a/src/Bipartite_graphs_matching/concept/Persistence_diagram.h b/src/Bipartite_graphs_matching/concept/Persistence_diagram.h new file mode 100644 index 00000000..eaaf8bc5 --- /dev/null +++ b/src/Bipartite_graphs_matching/concept/Persistence_diagram.h @@ -0,0 +1,7 @@ +typedef typename std::pair<double,double> Diagram_point; + +struct Persistence_Diagram +{ + const_iterator<Diagram_point> cbegin() const; + const_iterator<Diagram_point> cend() const; +}; diff --git a/src/Bipartite_graphs_matching/example/CMakeLists.txt b/src/Bipartite_graphs_matching/example/CMakeLists.txt new file mode 100644 index 00000000..aaf32abd --- /dev/null +++ b/src/Bipartite_graphs_matching/example/CMakeLists.txt @@ -0,0 +1,27 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIBottleneckExample) + +# need CGAL 4.6 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + + add_executable (basic_example basic.cpp) + #add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3) + + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Bipartite_graphs_matching/example/basic.cpp b/src/Bipartite_graphs_matching/example/basic.cpp new file mode 100644 index 00000000..d190ab48 --- /dev/null +++ b/src/Bipartite_graphs_matching/example/basic.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 <http://www.gnu.org/licenses/>. + */ + +#include "../include/gudhi/Graph_matching.h" +#include <iostream> + +using namespace Gudhi::bipartite_graph_matching; + +int main() { + + std::vector< std::pair<double, double> > v1, v2; + + v1.emplace_back(2.7,3.7); + v1.emplace_back(9.6,14); + v1.emplace_back(34.2,34.974); + + 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_graphs_matching/include/gudhi/Graph_matching.h b/src/Bipartite_graphs_matching/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..a2754333 --- /dev/null +++ b/src/Bipartite_graphs_matching/include/gudhi/Graph_matching.h @@ -0,0 +1,213 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ + +#include <deque> + +#include "Neighbors_finder.h" + +namespace Gudhi { + +namespace bipartite_graph_matching { + +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. + * + * + * + * \ingroup bottleneck_distance + */ +template<typename Persistence_diagram1, typename Persistence_diagram2> +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<int> v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::list<int> unmatched_in_u; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + std::unique_ptr<Layered_neighbors_finder> layering() const; + /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ + void update(std::deque<int> & 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<int> tries(unmatched_in_u); + for (auto it = tries.cbegin(); it != tries.cend(); it++) + // 'augment' has side-effects which have to be always executed, don't change order + successful = augment(layered_nf, *it, max_depth) || successful; + return successful; +} + +inline void Graph_matching::set_r(double r) { + this->r = r; +} + +inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { + //V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::deque<int> path; + path.emplace_back(u_start_index); + do { + if (static_cast<int>(path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + path.emplace_back(layered_nf.pull_near(path.back(), static_cast<int>(path.size())/2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + //if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex + path.pop_back(); + update(path); + return true; +} + +inline std::unique_ptr<Layered_neighbors_finder> Graph_matching::layering() const { + std::list<int> u_vertices(unmatched_in_u); + std::list<int> 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_neighbors_finder> 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<int> > 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<int>& path) { + unmatched_in_u.remove(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + // Be careful, the iterator is incremented twice each time + int tmp = *it; + v_to_u[*(++it)] = tmp; + } +} + +template<typename Persistence_diagram1, typename Persistence_diagram2> +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { + G::initialize(diag1, diag2, e); + std::unique_ptr< std::vector<double> > 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<int>((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 bipartite_graph_matching + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h b/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..78b9debc --- /dev/null +++ b/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h @@ -0,0 +1,155 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ + +#include <unordered_set> + +#include "Planar_neighbors_finder.h" + +namespace Gudhi { + +namespace bipartite_graph_matching { + +/** \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<int> > pull_all_near(int u_point_index); + +private: + const double r; + Planar_neighbors_finder planar_neighbors_f; + std::unordered_set<int> projections_f; + void remove(int v_point_index); + bool contains(int v_point_index); +}; + +/** \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> neighbors_finder; +}; + +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<int> > Neighbors_finder::pull_all_near(int u_point_index) { + std::unique_ptr< std::list<int> > 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; +} + +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<int> (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer).pull_near(u_point_index); +} + +inline int Layered_neighbors_finder::vlayers_number() const { + return static_cast<int>(neighbors_finder.size()); +} + +} // namespace bipartite_graph_matching + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ diff --git a/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h b/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h new file mode 100644 index 00000000..6420b772 --- /dev/null +++ b/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h @@ -0,0 +1,164 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ + +#include <vector> +#include <set> +#include <cmath> +#include <utility> +#include <algorithm> +#include <math.h> +#include <memory> + +namespace Gudhi { + +namespace bipartite_graph_matching { + +/** \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<typename Persistence_diagram1, typename Persistence_diagram2> + 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<double> > sorted_distances(); + +private: + /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ + typedef std::pair<double, double> Internal_point; + static std::vector<Internal_point> u; + static std::vector<Internal_point> v; + static Internal_point get_u_point(int u_point_index); + static Internal_point get_v_point(int v_point_index); + + friend class Naive_pnf; + friend class Upper_envelope_tree; +}; + +/** \internal \typedef \brief Shorter alias */ +typedef Persistence_diagrams_graph G; + +// static initialization, seems to work but strange +std::vector<G::Internal_point> G::u = [] {return std::vector<G::Internal_point>();}(); +std::vector<G::Internal_point> G::v = [] {return std::vector<G::Internal_point>();}(); + +inline int null_point_index() { + return -1; +} + +template<typename Persistence_diagram1, typename Persistence_diagram2> +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<int> (u.size()); +} + +inline bool G::on_the_v_diagonal(int v_point_index) { + return v_point_index >= static_cast<int> (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<int> (v.size()) : v_point_index + static_cast<int> (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<int> (u.size()) : u_point_index + static_cast<int> (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<int> (u.size() + v.size()); +} + +inline std::unique_ptr< std::vector<double> > G::sorted_distances() { + // could be optimized + std::set<double> sorted_distances; + for (int u_point_index = 0; u_point_index < size(); ++u_point_index) + for (int v_point_index = 0; v_point_index < size(); ++v_point_index) + sorted_distances.emplace(distance(u_point_index, v_point_index)); + std::unique_ptr< std::vector<double> > sd_up(new std::vector<double>(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); +} + +} // namespace bipartite_graph_matching + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ diff --git a/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h new file mode 100644 index 00000000..3564dcf9 --- /dev/null +++ b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h @@ -0,0 +1,170 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ + +#include <list> +#include <map> + +#include "Persistence_diagrams_graph.h" + +namespace Gudhi { + +namespace bipartite_graph_matching { + +/** \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<int> > 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); + /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */ + virtual std::unique_ptr< std::list<int> > pull_all_near(int u_point_index); + +private: + std::pair<int,int> get_v_key(int v_point_index) const; + std::multimap<std::pair<int,int>,int> grid; +}; + +/** \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<int> > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { + std::unique_ptr< std::list<int> > all_pull(new std::list<int>); + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull->emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; +} + +inline Naive_pnf::Naive_pnf(double r) : + Abstract_planar_neighbors_finder(r), grid() { } + + +inline std::pair<int,int> Naive_pnf::get_v_key(int v_point_index) const{ + G::Internal_point v_point = G::get_v_point(v_point_index); + return std::make_pair(static_cast<int>(v_point.first/r), static_cast<int>(v_point.second/r)); +} + +inline void Naive_pnf::add(int v_point_index) { + grid.emplace(get_v_key(v_point_index),v_point_index); +} + +inline void Naive_pnf::remove(int v_point_index) { + for(auto it = grid.find(get_v_key(v_point_index)); it!=grid.end(); it++) + if(it->second==v_point_index){ + grid.erase(it); + return; + } +} + +inline bool Naive_pnf::contains(int v_point_index) const { + if(v_point_index == null_point_index()) + return false; + for(auto it = grid.find(get_v_key(v_point_index)); it!=grid.end(); it++) + if(it->second==v_point_index) + return true; + return false; +} + +inline int Naive_pnf::pull_near(int u_point_index) { + G::Internal_point u_point = G::get_u_point(u_point_index); + int i0 = static_cast<int>(u_point.first/r); + int j0 = static_cast<int>(u_point.second/r); + for(int i = 1; i<= 3; i++) + for(int j = 1; j<= 3; j++) + for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end(); it++) + if (G::distance(u_point_index, it->second) <= r) { + int tmp = it->second; + grid.erase(it); + return tmp; + } + return null_point_index(); +} + +inline std::unique_ptr< std::list<int> > Naive_pnf::pull_all_near(int u_point_index) { + std::unique_ptr< std::list<int> > all_pull(new std::list<int>); + G::Internal_point u_point = G::get_u_point(u_point_index); + int i0 = static_cast<int>(u_point.first/r); + int j0 = static_cast<int>(u_point.second/r); + for(int i = 1; i<= 3; i++) + for(int j = 1; j<= 3; j++) + for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end(); it++) + if (G::distance(u_point_index, it->second) <= r) { + int tmp = it->second; + grid.erase(it); + all_pull->emplace_back(tmp); + } + return all_pull; +} + +} // namespace bipartite_graph_matching + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ diff --git a/src/Bipartite_graphs_matching/test/CMakeLists.txt b/src/Bipartite_graphs_matching/test/CMakeLists.txt new file mode 100644 index 00000000..abf18b1b --- /dev/null +++ b/src/Bipartite_graphs_matching/test/CMakeLists.txt @@ -0,0 +1,54 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIBottleneckUT) + + + + +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + + 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) + + + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile bipartite graphs matching feature. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Bipartite_graphs_matching/test/README b/src/Bipartite_graphs_matching/test/README new file mode 100644 index 00000000..0e7b8673 --- /dev/null +++ b/src/Bipartite_graphs_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_graphs_matching/test/bottleneck_unit_test.cpp b/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp new file mode 100644 index 00000000..dc33f76b --- /dev/null +++ b/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp @@ -0,0 +1,207 @@ +#define BOOST_TEST_MODULE bottleneck test + +#include <boost/test/included/unit_test.hpp> +#include <random> +#include "../include/gudhi/Graph_matching.h" + +#include <chrono> +#include <fstream> + +using namespace Gudhi::bipartite_graph_matching; + +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<double> unif1(0.,upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound/1000.,upper_bound/100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < 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); + } + //99.5 and not 100 to avoid float errors. + BOOST_CHECK(bottleneck_distance(v1, v2) <= upper_bound/99.5); +} + +BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ + // Random construction + std::uniform_real_distribution<double> unif(0.,upper_bound); + std::default_random_engine re; + std::vector< std::pair<double, double> > 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<double> > 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<n1; v_point_index+=2) + pnf.add(v_point_index); + // + BOOST_CHECK(pnf.contains(0)); + BOOST_CHECK(!pnf.contains(1)); + BOOST_CHECK(pnf.contains(2)); + BOOST_CHECK(!pnf.contains(3)); + // + 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(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)); + std::list<int> 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<int> l = *nf.pull_all_near(n2/2); + bool v = true; + for(auto it = l.cbegin(); it != l.cend(); ++it) + v = v && (G::distance(n2/2,*it)>1.); + BOOST_CHECK(v); + int v_point_index_2 = nf.pull_near(n2/2); + BOOST_CHECK(v_point_index_2 == -1); +} + +BOOST_AUTO_TEST_CASE(layered_neighbors_finder) { + 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(chrono) { + std::ofstream objetfichier; + objetfichier.open("results.csv", std::ios::out); + + for(int n =50; n<=1000; n+=100){ + std::uniform_real_distribution<double> unif1(0.,upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound/1000.,upper_bound/100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a,b), std::max(a,b)); + v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); + if(i%5==0) + v1.emplace_back(std::min(a,b),std::min(a,b)+x); + if(i%3==0) + v2.emplace_back(std::max(a,b),std::max(a,b)+y); + } + + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + double b = bottleneck_distance(v1,v2); + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + + typedef std::chrono::duration<int,std::milli> millisecs_t; + millisecs_t duration(std::chrono::duration_cast<millisecs_t>(end-start)); + objetfichier << n << ";" << duration.count() << ";" << b << std::endl; + } + objetfichier.close(); +} +*/ + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ceb993fa..1436d65a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -40,4 +40,20 @@ else() add_subdirectory(example/Alpha_shapes) add_subdirectory(example/Bottleneck) -endif() +endif() + + include( ${CGAL_USE_FILE} ) + # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CXX_FLAGS are overwritten. + # cf. http://doc.cgal.org/latest/Manual/installation.html#title40 + # A workaround is to add "-std=c++11" again. + # A fix would be to use https://cmake.org/cmake/help/v3.1/prop_gbl/CMAKE_CXX_KNOWN_FEATURES.html + # or even better https://cmake.org/cmake/help/v3.1/variable/CMAKE_CXX_STANDARD.html + # but it implies to use cmake version 3.1 at least. + if(NOT MSVC) + include(CheckCXXCompilerFlag) + CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11) + if(COMPILER_SUPPORTS_CXX11) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") + endif() + endif() + # - End of workaround |