From 3b5b94c90c81ae4938a89b2b4df8b1173f100ee3 Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 2 Jun 2016 17:09:48 +0000 Subject: renaming git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1240 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 53bcbf20e176a21dc41ef3c4d4295b1e9aa959b0 --- src/Bottleneck_distance/concept/Diagram_point.h | 6 + .../concept/Persistence_diagram.h | 7 + src/Bottleneck_distance/example/CMakeLists.txt | 36 ++++ .../example/bottleneck_example.cpp | 43 ++++ .../include/CGAL/Miscellaneous.h | 51 +++++ .../include/gudhi/Graph_matching.h | 233 +++++++++++++++++++++ .../include/gudhi/Internal_point.h | 74 +++++++ .../include/gudhi/Neighbors_finder.h | 154 ++++++++++++++ .../include/gudhi/Persistence_diagrams_graph.h | 168 +++++++++++++++ .../include/gudhi/Planar_neighbors_finder.h | 218 +++++++++++++++++++ src/Bottleneck_distance/test/CMakeLists.txt | 63 ++++++ src/Bottleneck_distance/test/README | 12 ++ src/Bottleneck_distance/test/bottleneck_chrono.cpp | 61 ++++++ .../test/bottleneck_unit_test.cpp | 188 +++++++++++++++++ 14 files changed, 1314 insertions(+) create mode 100644 src/Bottleneck_distance/concept/Diagram_point.h create mode 100644 src/Bottleneck_distance/concept/Persistence_diagram.h create mode 100644 src/Bottleneck_distance/example/CMakeLists.txt create mode 100644 src/Bottleneck_distance/example/bottleneck_example.cpp create mode 100644 src/Bottleneck_distance/include/CGAL/Miscellaneous.h create mode 100644 src/Bottleneck_distance/include/gudhi/Graph_matching.h create mode 100644 src/Bottleneck_distance/include/gudhi/Internal_point.h create mode 100644 src/Bottleneck_distance/include/gudhi/Neighbors_finder.h create mode 100644 src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h create mode 100644 src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h create mode 100644 src/Bottleneck_distance/test/CMakeLists.txt create mode 100644 src/Bottleneck_distance/test/README create mode 100644 src/Bottleneck_distance/test/bottleneck_chrono.cpp create mode 100644 src/Bottleneck_distance/test/bottleneck_unit_test.cpp (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/concept/Diagram_point.h b/src/Bottleneck_distance/concept/Diagram_point.h new file mode 100644 index 00000000..ce9c6d5b --- /dev/null +++ b/src/Bottleneck_distance/concept/Diagram_point.h @@ -0,0 +1,6 @@ + + +struct Diagram_point{ + double first; + double second; +}; diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h new file mode 100644 index 00000000..4951af32 --- /dev/null +++ b/src/Bottleneck_distance/concept/Persistence_diagram.h @@ -0,0 +1,7 @@ + + +struct Persistence_Diagram +{ + const_iterator cbegin() const; + const_iterator cend() const; +}; diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt new file mode 100644 index 00000000..4e8766cc --- /dev/null +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -0,0 +1,36 @@ +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} ) + + 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 + + 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/Bottleneck_distance/example/bottleneck_example.cpp b/src/Bottleneck_distance/example/bottleneck_example.cpp new file mode 100644 index 00000000..bbd83547 --- /dev/null +++ b/src/Bottleneck_distance/example/bottleneck_example.cpp @@ -0,0 +1,43 @@ +/* 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 + +using namespace Gudhi::Bottleneck_distance; + +int main() { + std::vector< std::pair > v1, v2; + + v1.push_back(std::pair(2.7,3.7)); + v1.push_back(std::pair(9.6,14)); + v1.push_back(std::pair(34.2,34.974)); + + v2.push_back(std::pair(2.8,4.45)); + v2.push_back(std::pair(9.5,14.1)); + + + double b = bottleneck_distance(v1, v2); + + std::cout << "Bottleneck distance = " << b << std::endl; + +} diff --git a/src/Bottleneck_distance/include/CGAL/Miscellaneous.h b/src/Bottleneck_distance/include/CGAL/Miscellaneous.h new file mode 100644 index 00000000..4de787fb --- /dev/null +++ b/src/Bottleneck_distance/include/CGAL/Miscellaneous.h @@ -0,0 +1,51 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(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_CGAL_MISCELLANEOUS_H_ +#define SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_ + +#include + +namespace CGAL { + +typedef Gudhi::bipartite_graph_matching::Internal_point Internal_point; + +template <> +struct Kernel_traits { + struct Kernel { + typedef double FT; + typedef double RT; + }; +}; + + +struct Construct_coord_iterator { + typedef const double* result_type; + const double* operator()(const Internal_point& p) const + { return static_cast(p.vec); } + const double* operator()(const Internal_point& p, int) const + { return static_cast(p.vec+2); } +}; + +} //namespace CGAL + +#endif // SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..d8860841 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -0,0 +1,233 @@ +/* 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 + +namespace Gudhi { + +namespace Bottleneck_distance { + +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. + * + * + * + * \ingroup bottleneck_distance + */ +template +double compute(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::shared_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::shared_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::shared_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 it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { + std::shared_ptr> u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ->begin(); it2 != u_succ->end(); ++it2) { + layered_nf->add(*it2, layer); + v_vertices.emplace_back(*it2); + } + } + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS + u_vertices.clear(); + bool end = false; + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) + if (v_to_u.at(*it) == null_point_index()) + // we stop when a nearest exposed V vertex (from U exposed vertices) has been found + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + // When the above for finishes, we have progress of one half-step (from V to U) in the BFS + if (end) + return layered_nf; + v_vertices.clear(); + } + return layered_nf; +} + +inline void Graph_matching::update(std::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; + v_to_u[*(++it)] = tmp; + } +} + +template +double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2) { + G::initialize(diag1, diag2, 0.); + std::shared_ptr< std::vector > sd(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); +} + +template +double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { + if(e< std::numeric_limits::min()) + return compute_exactly(diag1, diag2); + G::initialize(diag1, diag2, e); + double d = 0.; + double f = G::diameter(); + while (f-d > e){ + Graph_matching m; + m.set_r((d+f)/2.); + while (m.multi_augment()); + //The above while compute a maximum matching (according to the r setted before) + if (m.perfect()) + f = (d+f)/2.; + else + d= (d+f)/2.; + } + return d; +} + +} // namespace Bottleneck_distance + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h new file mode 100644 index 00000000..2080900d --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -0,0 +1,74 @@ +/* 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_INTERNAL_POINT_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_ + +namespace Gudhi { + +namespace Bottleneck_distance { + +/** \internal \brief Returns the used index for encoding none of the points */ +int null_point_index(); + +/** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ +struct Internal_point { + double vec[2]; + int point_index; + Internal_point() {} + Internal_point(double x, double y, int p_i = null_point_index()) { vec[0]=x; vec[1]=y; point_index = p_i; } + double x() const { return vec[ 0 ]; } + double y() const { return vec[ 1 ]; } + double& x() { return vec[ 0 ]; } + double& y() { return vec[ 1 ]; } + bool operator==(const Internal_point& p) const + { + return point_index==p.point_index; + } + bool operator!=(const Internal_point& p) const { return ! (*this == p); } +/* +Useless + friend void swap(Internal_point& a, Internal_point& b) + { + using std::swap; + double x_tmp = a.vec[0]; + double y_tmp = a.vec[1]; + int pi_tmp = a.point_index; + a.vec[0] = b.vec[0]; + a.vec[1] = b.vec[1]; + a.point_index = b.point_index; + b.vec[0] = x_tmp; + b.vec[1] = y_tmp; + b.point_index = pi_tmp; + } +*/ +}; + +inline int null_point_index() { + return -1; +} + +} // namespace Bottleneck_distance + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..7ad79126 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -0,0 +1,154 @@ +/* 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 + +namespace Gudhi { + +namespace Bottleneck_distance { + +/** \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::shared_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); +}; + +/** \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 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()) + //Any pair of projection is 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::shared_ptr< std::list > Neighbors_finder::pull_all_near(int u_point_index) { + std::shared_ptr< 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 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(std::shared_ptr(new 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_distance + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h new file mode 100644 index 00000000..e8a46f1c --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h @@ -0,0 +1,168 @@ +/* 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 +#include +#include + +namespace Gudhi { + +namespace Bottleneck_distance { + + +/** \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::shared_ptr< std::vector > sorted_distances(); + /** \internal \brief Returns an upper bound of the diameter of the convex hull */ + static double diameter(); + +private: + 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); + + friend class Naive_pnf; + friend class Cgal_pnf; +}; + +/** \internal \typedef \brief Shorter alias */ +typedef Persistence_diagrams_graph G; + +// static initialization +std::vector G::u = [] {return std::vector();}(); +std::vector G::v = [] {return std::vector();}(); + +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.push_back(Internal_point(it->first, it->second, u.size())); + for (auto it = diag2.cbegin(); it != diag2.cend(); ++it) + if (it->second - it->first > e) + v.push_back(Internal_point(it->first, it->second, v.size())); + 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.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); +} + +inline int G::size() { + return static_cast (u.size() + v.size()); +} + +inline std::shared_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::shared_ptr< std::vector > sd_up(new std::vector(sorted_distances.cbegin(), sorted_distances.cend())); + return sd_up; +} + +inline 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 m = (projector.x() + projector.y()) / 2; + return Internal_point(m,m,u_point_index); +} + +inline 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 m = (projector.x() + projector.y()) / 2; + return Internal_point(m,m,v_point_index); +} + +inline double G::diameter() { + double max = 0.; + for(it = u.cbegin(); it != u.cend(); it++) + max = std::max(max,it->y()); + for(it = v.cbegin(); it != v.cend(); it++) + max = std::max(max,it->y()); + return max; +} + +} // namespace Bottleneck_distance + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h new file mode 100644 index 00000000..c0b6383f --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h @@ -0,0 +1,218 @@ +/* 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 +#include +#include +#include +#include + + +namespace Gudhi { + +namespace Bottleneck_distance { + +/** \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 Naive_pnf { +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::shared_ptr< std::list > pull_all_near(int u_point_index); + +private: + double r; + std::pair get_v_key(int v_point_index) const; + std::multimap,int> grid; +}; + +class Cgal_pnf { + + typedef CGAL::Dimension_tag<2> D; + typedef CGAL::Search_traits Traits; + typedef CGAL::Weighted_Minkowski_distance Distance; + typedef CGAL::Orthogonal_incremental_neighbor_search K_neighbor_search; + typedef K_neighbor_search::Tree Kd_tree; + + +public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Cgal_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 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::shared_ptr< std::list > pull_all_near(int u_point_index); + +private: + double r; + std::set contents; + Kd_tree kd_t; +}; + +/** \internal \typedef \brief Planar_neighbors_finder is the used implementation. */ +typedef Naive_pnf Planar_neighbors_finder; + +inline Naive_pnf::Naive_pnf(double r_) : + r(r_), grid() { } + + +inline std::pair Naive_pnf::get_v_key(int v_point_index) const{ + Internal_point v_point = G::get_v_point(v_point_index); + return std::make_pair(static_cast(v_point.x()/r), static_cast(v_point.y()/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) { + if(v_point_index != null_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) { + Internal_point u_point = G::get_u_point(u_point_index); + int i0 = static_cast(u_point.x()/r); + int j0 = static_cast(u_point.y()/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::shared_ptr< std::list > Naive_pnf::pull_all_near(int u_point_index) { + std::shared_ptr< std::list > all_pull(new std::list); + Internal_point u_point = G::get_u_point(u_point_index); + int i0 = static_cast(u_point.x()/r); + int j0 = static_cast(u_point.y()/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();) + if (G::distance(u_point_index, it->second) <= r) { + all_pull->emplace_back(it->second); + it = grid.erase(it); + } + else it++; + return all_pull; +} + + +/** \internal \brief Constructor taking the near distance definition as parameter. */ +inline Cgal_pnf::Cgal_pnf(double r_) + : r(r_), contents(), kd_t() {} + + +/** \internal \brief A point added will be possibly pulled. */ +inline void Cgal_pnf::add(int v_point_index){ + if(v_point_index == null_point_index()) + return; + contents.insert(v_point_index); + kd_t.insert(G::get_v_point(v_point_index)); +} + +/** \internal \brief A point manually removed will no longer be possibly pulled. */ +inline void Cgal_pnf::remove(int v_point_index){ + if(contains(v_point_index)){ + contents.erase(v_point_index); + kd_t.remove(G::get_v_point(v_point_index)); + } +} + +/** \internal \brief Can the point given as parameter be returned ? */ +inline bool Cgal_pnf::contains(int v_point_index) const{ + if(v_point_index == null_point_index()) + return false; + return contents.count(v_point_index)>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. */ +inline int Cgal_pnf::pull_near(int u_point_index){ + Internal_point u_point = G::get_u_point(u_point_index); + std::vector w = {1., 1.}; + K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w)); + auto it = search.begin(); + if(it==search.end() || G::distance(u_point_index, it->first.point_index) > r) + return null_point_index(); + int tmp = it->first.point_index; + remove(tmp); + return tmp; +} + +inline std::shared_ptr< std::list > Cgal_pnf::pull_all_near(int u_point_index) { + std::shared_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; +} + + +} // namespace Bottleneck_distance + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt new file mode 100644 index 00000000..ab0461cc --- /dev/null +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -0,0 +1,63 @@ +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} ) + + 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 + + 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/Bottleneck_distance/test/README b/src/Bottleneck_distance/test/README new file mode 100644 index 00000000..0e7b8673 --- /dev/null +++ b/src/Bottleneck_distance/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/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp new file mode 100644 index 00000000..8f187a91 --- /dev/null +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -0,0 +1,61 @@ +/* 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 +#include + +using namespace Gudhi::Bottleneck_distance; + + +double upper_bound = 400.; // any real >0 + +int main(){ + std::ofstream objetfichier; + objetfichier.open("results.csv", std::ios::out); + + for(int n=0; n<=2500; n+=250){ + 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++) { + 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_approx(v1,v2, 0.); + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + typedef std::chrono::duration millisecs_t; + millisecs_t duration(std::chrono::duration_cast(end-start)); + objetfichier << n << ";" << duration.count() << ";" << b << std::endl; + } + objetfichier.close(); +} diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp new file mode 100644 index 00000000..33ac3f1b --- /dev/null +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -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 . + */ + + +#define BOOST_TEST_MODULE bottleneck test + +#include +#include +#include + +using namespace Gudhi::Bottleneck_distance; + +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 406.43; // any real >0 + +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::shared_ptr< std::vector > d(G::sorted_distances()); + // + BOOST_CHECK(!G::on_the_u_diagonal(n1-1)); + BOOST_CHECK(!G::on_the_u_diagonal(n1)); + BOOST_CHECK(!G::on_the_u_diagonal(n2-1)); + BOOST_CHECK(G::on_the_u_diagonal(n2)); + BOOST_CHECK(!G::on_the_v_diagonal(n1-1)); + BOOST_CHECK(G::on_the_v_diagonal(n1)); + BOOST_CHECK(G::on_the_v_diagonal(n2-1)); + BOOST_CHECK(G::on_the_v_diagonal(n2)); + // + BOOST_CHECK(G::corresponding_point_in_u(0)==n2); + BOOST_CHECK(G::corresponding_point_in_u(n1)==0); + BOOST_CHECK(G::corresponding_point_in_v(0)==n1); + BOOST_CHECK(G::corresponding_point_in_v(n2)==0); + // + BOOST_CHECK(G::size()==(n1+n2)); + // + BOOST_CHECK((int) d->size() <= (n1+n2)*(n1+n2) - n1*n2 + 1); + BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,0))==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(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(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(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(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.); +} -- cgit v1.2.3 From bd9b3b687b715a968216032f20f08fbaea7f3baa Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 2 Jun 2016 17:16:04 +0000 Subject: lists git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1242 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 134d7d5c1f940ac342eb1721de6ee9cb41f1102b --- CMakeLists.txt | 6 +++--- src/Bottleneck_distance/example/CMakeLists.txt | 2 +- src/Bottleneck_distance/test/CMakeLists.txt | 3 ++- 3 files changed, 6 insertions(+), 5 deletions(-) (limited to 'src/Bottleneck_distance') diff --git a/CMakeLists.txt b/CMakeLists.txt index 066a49cf..dccd59d8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -57,7 +57,7 @@ else() include_directories(src/common/include/) include_directories(src/Alpha_shapes/include/) - include_directories(src/Bipartite_graphs_matching/include/) + include_directories(src/Bottleneck_distance/include/) include_directories(src/Bottleneck/include/) include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) @@ -75,8 +75,8 @@ else() add_subdirectory(src/Hasse_complex/example) add_subdirectory(src/Alpha_shapes/example) add_subdirectory(src/Alpha_shapes/test) - add_subdirectory(src/Bipartite_graphs_matching/example) - add_subdirectory(src/Bipartite_graphs_matching/test) + add_subdirectory(src/Bottleneck_distance/example) + add_subdirectory(src/Bottleneck_distance/test) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt index 4e8766cc..cc1b9091 100644 --- a/src/Bottleneck_distance/example/CMakeLists.txt +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -24,7 +24,7 @@ if(CGAL_FOUND) include( ${EIGEN3_USE_FILE} ) include_directories (BEFORE "../../include") - add_executable (basic_example basic.cpp) + add_executable (bottleneck_example bottleneck_example.cpp) #add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3) else() diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt index ab0461cc..1338ed6b 100644 --- a/src/Bottleneck_distance/test/CMakeLists.txt +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -44,7 +44,8 @@ endif() 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 ) +add_executable ( bottleneckUT bottleneck_unit_test.cpp ) +add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests -- cgit v1.2.3 From ab8c111e6a9abc4ba718083266c3dbb238ab94e3 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 2 Jun 2016 18:23:10 +0000 Subject: CMake fix after executable renaming git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1244 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 5238ff5ecdc70b6e97539625b687ac0626c49ef2 --- src/Bottleneck_distance/test/CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt index 1338ed6b..54a2fe95 100644 --- a/src/Bottleneck_distance/test/CMakeLists.txt +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -46,13 +46,13 @@ endif() add_executable ( bottleneckUT bottleneck_unit_test.cpp ) add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) -target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +target_link_libraries(bottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests -add_test(NAME BottleneckUT - COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUT +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) + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/bottleneckUT.xml --log_level=test_suite --report_level=no) else() -- cgit v1.2.3 From 714a97d8b9c64a0a5ccc837184def8122ec8b4cd Mon Sep 17 00:00:00 2001 From: fgodi Date: Fri, 3 Jun 2016 09:32:46 +0000 Subject: compile git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1245 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: a26ddc43a1f71052b594d5e0d65d28d952a0917e --- src/Bottleneck_distance/example/bottleneck_example.cpp | 3 +-- src/Bottleneck_distance/include/CGAL/Miscellaneous.h | 2 +- src/Bottleneck_distance/include/gudhi/Graph_matching.h | 3 +++ src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h | 4 ++-- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 2 +- src/Bottleneck_distance/test/bottleneck_unit_test.cpp | 2 +- 6 files changed, 9 insertions(+), 7 deletions(-) (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/example/bottleneck_example.cpp b/src/Bottleneck_distance/example/bottleneck_example.cpp index bbd83547..c6a06235 100644 --- a/src/Bottleneck_distance/example/bottleneck_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_example.cpp @@ -23,7 +23,6 @@ #include #include -using namespace Gudhi::Bottleneck_distance; int main() { std::vector< std::pair > v1, v2; @@ -36,7 +35,7 @@ int main() { v2.push_back(std::pair(9.5,14.1)); - double b = bottleneck_distance(v1, v2); + double b = Gudhi::Bottleneck_distance::compute(v1, v2); std::cout << "Bottleneck distance = " << b << std::endl; diff --git a/src/Bottleneck_distance/include/CGAL/Miscellaneous.h b/src/Bottleneck_distance/include/CGAL/Miscellaneous.h index 4de787fb..21be0bc7 100644 --- a/src/Bottleneck_distance/include/CGAL/Miscellaneous.h +++ b/src/Bottleneck_distance/include/CGAL/Miscellaneous.h @@ -27,7 +27,7 @@ namespace CGAL { -typedef Gudhi::bipartite_graph_matching::Internal_point Internal_point; +typedef Gudhi::Bottleneck_distance::Internal_point Internal_point; template <> struct Kernel_traits { diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index d8860841..f0ce633f 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -40,6 +40,9 @@ namespace Bottleneck_distance { template double compute(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); +template +double compute_exactly(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2); + /** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph. * * \ingroup bottleneck_distance diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h index e8a46f1c..52e7c583 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h @@ -154,9 +154,9 @@ inline Internal_point G::get_v_point(int v_point_index) { inline double G::diameter() { double max = 0.; - for(it = u.cbegin(); it != u.cend(); it++) + for(auto it = u.cbegin(); it != u.cend(); it++) max = std::max(max,it->y()); - for(it = v.cbegin(); it != v.cend(); it++) + for(auto it = v.cbegin(); it != v.cend(); it++) max = std::max(max,it->y()); return max; } diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index 8f187a91..cde8afb1 100644 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -51,7 +51,7 @@ int main(){ 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_approx(v1,v2, 0.); + double b = compute(v1,v2, 0.01); std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); typedef std::chrono::duration millisecs_t; millisecs_t duration(std::chrono::duration_cast(end-start)); diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp index 33ac3f1b..d5b71156 100644 --- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -184,5 +184,5 @@ BOOST_AUTO_TEST_CASE(global){ 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_CHECK(compute(v1, v2) <= upper_bound/100.); } -- cgit v1.2.3 From d85314eb2cd5b8487ee5cc5c9069e1c8ff311931 Mon Sep 17 00:00:00 2001 From: fgodi Date: Fri, 3 Jun 2016 10:22:53 +0000 Subject: CGAL+ VERSION git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1246 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: b9bde0a97e56042b84d1850dea34e559dc196b91 --- .../example/bottleneck_example.cpp | 6 +-- .../include/gudhi/Graph_matching.h | 48 ++++++++++++++++++---- .../include/gudhi/Persistence_diagrams_graph.h | 1 + .../include/gudhi/Planar_neighbors_finder.h | 2 +- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 4 +- 5 files changed, 47 insertions(+), 14 deletions(-) (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/example/bottleneck_example.cpp b/src/Bottleneck_distance/example/bottleneck_example.cpp index c6a06235..b3d26f3c 100644 --- a/src/Bottleneck_distance/example/bottleneck_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_example.cpp @@ -34,9 +34,9 @@ int main() { v2.push_back(std::pair(2.8,4.45)); v2.push_back(std::pair(9.5,14.1)); + double b = Gudhi::Bottleneck_distance::compute(v1, v2, 0.0001); + std::cout << "Approx. bottleneck distance = " << b << std::endl; - double b = Gudhi::Bottleneck_distance::compute(v1, v2); - + b = Gudhi::Bottleneck_distance::compute(v1, v2); std::cout << "Bottleneck distance = " << b << std::endl; - } diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index f0ce633f..69470067 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -209,6 +209,43 @@ double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diag return sd->at(idmin); } +template +double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { + if(e< std::numeric_limits::min()) + return compute_exactly(diag1, diag2); + G::initialize(diag1, diag2, e); + int in = G::diameter()/e; + int idmin = 0; + int idmax = in; + // alpha can be modified, this will change the complexity + double alpha = pow(in, 0.25); + Graph_matching m; + Graph_matching biggest_unperfect; + while (idmin != idmax) { + int step = static_cast((idmax - idmin) / alpha); + m.set_r(e*(idmin + step)); + while (m.multi_augment()); + //The above while compute a maximum matching (according to the r setted before) + if (m.perfect()) { + idmax = idmin + step; + m = biggest_unperfect; + } else { + biggest_unperfect = m; + idmin = idmin + step + 1; + } + } + return e*(idmin); +} + + + +} // namespace Bottleneck_distance + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ + +/* Dichotomic version template double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { if(e< std::numeric_limits::min()) @@ -221,16 +258,11 @@ double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &di m.set_r((d+f)/2.); while (m.multi_augment()); //The above while compute a maximum matching (according to the r setted before) - if (m.perfect()) + if (m.perfect()) f = (d+f)/2.; - else + else d= (d+f)/2.; } return d; -} +} */ -} // namespace Bottleneck_distance - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h index 52e7c583..e31b2dae 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h @@ -129,6 +129,7 @@ inline int G::size() { inline std::shared_ptr< std::vector > G::sorted_distances() { // could be optimized std::set sorted_distances; + sorted_distances.emplace(0.); 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)); diff --git a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h index c0b6383f..3f3723c3 100644 --- a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h @@ -94,7 +94,7 @@ private: }; /** \internal \typedef \brief Planar_neighbors_finder is the used implementation. */ -typedef Naive_pnf Planar_neighbors_finder; +typedef Cgal_pnf Planar_neighbors_finder; inline Naive_pnf::Naive_pnf(double r_) : r(r_), grid() { } diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index cde8afb1..8b49b6be 100644 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -33,7 +33,7 @@ int main(){ std::ofstream objetfichier; objetfichier.open("results.csv", std::ios::out); - for(int n=0; n<=2500; n+=250){ + for(int n=0; n<=4000; n+=400){ std::uniform_real_distribution unif1(0.,upper_bound); std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); std::default_random_engine re; @@ -51,7 +51,7 @@ int main(){ 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 = compute(v1,v2, 0.01); + double b = compute(v1,v2, 0.0001); std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); typedef std::chrono::duration millisecs_t; millisecs_t duration(std::chrono::duration_cast(end-start)); -- cgit v1.2.3 From 284a45fc6ffcf049e667065acda71c17165b4f5c Mon Sep 17 00:00:00 2001 From: fgodi Date: Fri, 3 Jun 2016 12:28:19 +0000 Subject: concept git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1247 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 7bfce006ccae0784f3e0f3747f670eddaa9aba48 --- src/Bottleneck_distance/concept/Diagram_point.h | 6 ------ src/Bottleneck_distance/concept/Persistence_diagram.h | 15 +++++++++++++-- 2 files changed, 13 insertions(+), 8 deletions(-) delete mode 100644 src/Bottleneck_distance/concept/Diagram_point.h (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/concept/Diagram_point.h b/src/Bottleneck_distance/concept/Diagram_point.h deleted file mode 100644 index ce9c6d5b..00000000 --- a/src/Bottleneck_distance/concept/Diagram_point.h +++ /dev/null @@ -1,6 +0,0 @@ - - -struct Diagram_point{ - double first; - double second; -}; diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h index 4951af32..27a258d8 100644 --- a/src/Bottleneck_distance/concept/Persistence_diagram.h +++ b/src/Bottleneck_distance/concept/Persistence_diagram.h @@ -1,7 +1,18 @@ +namespace Gudhi { +namespace Bottleneck_distance { +namespace Concept { +struct Diagram_point{ + double first; + double second; +}; struct Persistence_Diagram { - const_iterator cbegin() const; - const_iterator cend() const; + const_iterator cbegin() const; + const_iterator cend() const; }; + +} //namespace Concept +} //namespace Bottleneck_distance +} //namespace Gudhi -- cgit v1.2.3 From e1d3c44120619ace7d64b1bb5151e873867c9639 Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 6 Jun 2016 08:27:19 +0000 Subject: kd_tree modified added git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1249 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: c2c617e0f808bb2c51042174c5d2b4bbdd3998c9 --- src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h index 3f3723c3..37c6a122 100644 --- a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h @@ -25,9 +25,7 @@ #include #include -#include -#include -#include +#include "../CGAL/Kd_tree.h" #include #include -- cgit v1.2.3 From 6c611b11c08bc314d8075ae6a038d4af549a7f2e Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 6 Jun 2016 08:27:57 +0000 Subject: kd_tree modified added git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1250 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 07254aea5d1921f8bb43f73ac1fe9d76b63c5c3e --- src/Bottleneck_distance/include/CGAL/Kd_tree.h | 517 +++++++++++++++++++++++++ 1 file changed, 517 insertions(+) create mode 100644 src/Bottleneck_distance/include/CGAL/Kd_tree.h (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/include/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/CGAL/Kd_tree.h new file mode 100644 index 00000000..4f874a8b --- /dev/null +++ b/src/Bottleneck_distance/include/CGAL/Kd_tree.h @@ -0,0 +1,517 @@ +// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Hans Tangelder (), +// : Waqar Khan + +#ifndef CGAL_KD_TREE_H +#define CGAL_KD_TREE_H + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include + +#ifdef CGAL_HAS_THREADS +#include +#endif + +namespace CGAL { + +//template , class UseExtendedNode = Tag_true > +template , class UseExtendedNode = Tag_true > +class Kd_tree { + +public: + typedef SearchTraits Traits; + typedef Splitter_ Splitter; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Splitter::Container Point_container; + + typedef typename SearchTraits::FT FT; + typedef Kd_tree_node Node; + typedef Kd_tree_leaf_node Leaf_node; + typedef Kd_tree_internal_node Internal_node; + typedef Kd_tree Tree; + typedef Kd_tree Self; + + typedef Node* Node_handle; + typedef const Node* Node_const_handle; + typedef Leaf_node* Leaf_node_handle; + typedef const Leaf_node* Leaf_node_const_handle; + typedef Internal_node* Internal_node_handle; + typedef const Internal_node* Internal_node_const_handle; + typedef typename std::vector::const_iterator Point_d_iterator; + typedef typename std::vector::const_iterator Point_d_const_iterator; + typedef typename Splitter::Separator Separator; + typedef typename std::vector::const_iterator iterator; + typedef typename std::vector::const_iterator const_iterator; + + typedef typename std::vector::size_type size_type; + + typedef typename internal::Get_dimension_tag::Dimension D; + +private: + SearchTraits traits_; + Splitter split; + + + // wokaround for https://svn.boost.org/trac/boost/ticket/9332 +#if (_MSC_VER == 1800) && (BOOST_VERSION == 105500) + std::deque internal_nodes; + std::deque leaf_nodes; +#else + boost::container::deque internal_nodes; + boost::container::deque leaf_nodes; +#endif + + Node_handle tree_root; + + Kd_tree_rectangle* bbox; + std::vector pts; + + // Instead of storing the points in arrays in the Kd_tree_node + // we put all the data in a vector in the Kd_tree. + // and we only store an iterator range in the Kd_tree_node. + // + std::vector data; + + +#ifdef CGAL_HAS_THREADS + mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build() +#endif + bool built_; + bool removed_; + + // protected copy constructor + Kd_tree(const Tree& tree) + : traits_(tree.traits_),built_(tree.built_) + {}; + + + // Instead of the recursive construction of the tree in the class Kd_tree_node + // we do this in the tree class. The advantage is that we then can optimize + // the allocation of the nodes. + + // The leaf node + Node_handle + create_leaf_node(Point_container& c) + { + Leaf_node node(true , static_cast(c.size())); + std::ptrdiff_t tmp = c.begin() - data.begin(); + node.data = pts.begin() + tmp; + + leaf_nodes.push_back(node); + Leaf_node_handle nh = &leaf_nodes.back(); + + + return nh; + } + + + // The internal node + + Node_handle + create_internal_node(Point_container& c, const Tag_true&) + { + return create_internal_node_use_extension(c); + } + + Node_handle + create_internal_node(Point_container& c, const Tag_false&) + { + return create_internal_node(c); + } + + + + // TODO: Similiar to the leaf_init function above, a part of the code should be + // moved to a the class Kd_tree_node. + // It is not proper yet, but the goal was to see if there is + // a potential performance gain through the Compact_container + Node_handle + create_internal_node_use_extension(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + + Separator sep; + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + int cd = nh->cutting_dimension(); + if(!c_low.empty()) + nh->low_val = c_low.tight_bounding_box().max_coord(cd); + else + nh->low_val = c_low.bounding_box().min_coord(cd); + if(!c.empty()) + nh->high_val = c.tight_bounding_box().min_coord(cd); + else + nh->high_val = c.bounding_box().max_coord(cd); + + CGAL_assertion(nh->cutting_value() >= nh->low_val); + CGAL_assertion(nh->cutting_value() <= nh->high_val); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node_use_extension(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node_use_extension(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + + return nh; + } + + + // Note also that I duplicated the code to get rid if the if's for + // the boolean use_extension which was constant over the construction + Node_handle + create_internal_node(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + Separator sep; + + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + return nh; + } + + + +public: + + Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + {} + + template + Kd_tree(InputIterator first, InputIterator beyond, + Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + { + pts.insert(pts.end(), first, beyond); + } + + bool empty() const { + return pts.empty(); + } + + void + build() + { + const Point_d& p = *pts.begin(); + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object(); + int dim = static_cast(std::distance(ccci(p), ccci(p,0))); + + data.reserve(pts.size()); + for(unsigned int i = 0; i < pts.size(); i++){ + data.push_back(&pts[i]); + } + Point_container c(dim, data.begin(), data.end(),traits_); + bbox = new Kd_tree_rectangle(c.bounding_box()); + if (c.size() <= split.bucket_size()){ + tree_root = create_leaf_node(c); + }else { + tree_root = create_internal_node(c, UseExtendedNode()); + } + + //Reorder vector for spatial locality + std::vector ptstmp; + ptstmp.resize(pts.size()); + for (std::size_t i = 0; i < pts.size(); ++i){ + ptstmp[i] = *data[i]; + } + for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ + std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); + leaf_nodes[i].data = ptstmp.begin() + tmp; + } + pts.swap(ptstmp); + + data.clear(); + + built_ = true; + } + +private: + //any call to this function is for the moment not threadsafe + void const_build() const { +#ifdef CGAL_HAS_THREADS + //this ensure that build() will be called once + CGAL_SCOPED_LOCK(building_mutex); + if(!is_built()) +#endif + const_cast(this)->build(); //THIS IS NOT THREADSAFE + } +public: + + bool is_built() const + { + return built_; + } + + void invalidate_built() + { + if(is_built()){ + internal_nodes.clear(); + leaf_nodes.clear(); + data.clear(); + delete bbox; + built_ = false; + } + } + + void clear() + { + invalidate_built(); + pts.clear(); + removed_ = false; + } + + void + insert(const Point_d& p) + { + if (removed_) throw std::logic_error("Once you start removing points, you cannot insert anymore, you need to start again from scratch."); + invalidate_built(); + pts.push_back(p); + } + + template + void + insert(InputIterator first, InputIterator beyond) + { + if (removed_ && first != beyond) throw std::logic_error("Once you start removing points, you cannot insert anymore, you need to start again from scratch."); + invalidate_built(); + pts.insert(pts.end(),first, beyond); + } + + void + remove(const Point_d& p) + { + // This does not actually remove points, and further insertions + // would make the points reappear, so we disallow it. + removed_ = true; + // Locate the point + Internal_node_handle grandparent = 0; + Internal_node_handle parent = 0; + bool islower = false, islower2; + Node_handle node = root(); // Calls build() if needed. + while (!node->is_leaf()) { + grandparent = parent; islower2 = islower; + parent = static_cast(node); + islower = traits().construct_cartesian_const_iterator_d_object()(p)[parent->cutting_dimension()] < parent->cutting_value(); + if (islower) { + node = parent->lower(); + } else { + node = parent->upper(); + } + } + Leaf_node_handle lnode = static_cast(node); + if (lnode->size() > 1) { + iterator pi = std::find(lnode->begin(), lnode->end(), p); + CGAL_assertion (pi != lnode->end()); + iterator lasti = lnode->end() - 1; + if (pi != lasti) { + // Hack to get a non-const iterator + std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin())); + } + lnode->drop_last_point(); + } else if (grandparent) { + CGAL_assertion (p == *lnode->begin()); + Node_handle brother = islower ? parent->upper() : parent->lower(); + if (islower2) + grandparent->set_lower(brother); + else + grandparent->set_upper(brother); + } else if (parent) { + tree_root = islower ? parent->upper() : parent->lower(); + } else { + clear(); + } + } + + //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known). + void reserve(size_t size) + { + pts.reserve(size); + } + + //Get the capacity of the underlying points vector. + size_t capacity() + { + return pts.capacity(); + } + + + template + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle b(*bbox); + return tree_root->search(it,q,b); + } + return it; + } + + + template + boost::optional + search_any_point(const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle b(*bbox); + return tree_root->search_any_point(q,b); + } + return boost::none; + } + + + ~Kd_tree() { + if(is_built()){ + delete bbox; + } + } + + + const SearchTraits& + traits() const + { + return traits_; + } + + Node_const_handle + root() const + { + if(! is_built()){ + const_build(); + } + return tree_root; + } + + Node_handle + root() + { + if(! is_built()){ + build(); + } + return tree_root; + } + + void + print() const + { + if(! is_built()){ + const_build(); + } + root()->print(); + } + + const Kd_tree_rectangle& + bounding_box() const + { + if(! is_built()){ + const_build(); + } + return *bbox; + } + + const_iterator + begin() const + { + return pts.begin(); + } + + const_iterator + end() const + { + return pts.end(); + } + + size_type + size() const + { + return pts.size(); + } + + // Print statistics of the tree. + std::ostream& + statistics(std::ostream& s) const + { + if(! is_built()){ + const_build(); + } + s << "Tree statistics:" << std::endl; + s << "Number of items stored: " + << root()->num_items() << std::endl; + s << "Number of nodes: " + << root()->num_nodes() << std::endl; + s << " Tree depth: " << root()->depth() << std::endl; + return s; + } + + +}; + +} // namespace CGAL + +#endif // CGAL_KD_TREE_H -- cgit v1.2.3 From aa9705ab541323e1d7aa0129c8faf0416b0489e0 Mon Sep 17 00:00:00 2001 From: fgodi Date: Sat, 23 Jul 2016 11:07:58 +0000 Subject: read files git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1393 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: d7f1b1166b65f7d7abc007abb47d86f8d4ba4157 --- .../example/bottleneck_example.cpp | 64 +++++++++++++++++----- 1 file changed, 49 insertions(+), 15 deletions(-) (limited to 'src/Bottleneck_distance') diff --git a/src/Bottleneck_distance/example/bottleneck_example.cpp b/src/Bottleneck_distance/example/bottleneck_example.cpp index b3d26f3c..b581a3be 100644 --- a/src/Bottleneck_distance/example/bottleneck_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_example.cpp @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author(s): Francois Godi, small modifications by Pawel Dlotko * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA Saclay (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -22,21 +22,55 @@ #include #include +#include +#include +#include +using namespace std; -int main() { - std::vector< std::pair > v1, v2; +//PD: I am using this type of procedures a lot in Gudhi stat. We should make one for all at some point. +std::vector< std::pair > read_diagram_from_file( const char* filename ) +{ + ifstream in; + in.open( filename ); + std::vector< std::pair > result; + if ( !in.is_open() ) + { + std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; + throw "File do not exist \n"; + } - v1.push_back(std::pair(2.7,3.7)); - v1.push_back(std::pair(9.6,14)); - v1.push_back(std::pair(34.2,34.974)); + std::string line; + while (!in.eof()) + { + getline(in,line); + if ( line.length() != 0 ) + { + std::stringstream lineSS; + lineSS << line; + double beginn, endd; + lineSS >> beginn; + lineSS >> endd; + result.push_back( std::make_pair( beginn , endd ) ); + } + } + in.close(); + return result; +}//read_diagram_from_file - v2.push_back(std::pair(2.8,4.45)); - v2.push_back(std::pair(9.5,14.1)); - - double b = Gudhi::Bottleneck_distance::compute(v1, v2, 0.0001); - std::cout << "Approx. bottleneck distance = " << b << std::endl; - - b = Gudhi::Bottleneck_distance::compute(v1, v2); - std::cout << "Bottleneck distance = " << b << std::endl; +int main( int argc , char** argv ) +{ + if ( argc < 3 ) + { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck distance (set by default to zero). The program will now terminate \n"; + } + std::vector< std::pair< double , double > > diag1 = read_diagram_from_file( argv[1] ); + std::vector< std::pair< double , double > > diag2 = read_diagram_from_file( argv[2] ); + double tolerance = 0; + if ( argc == 4 ) + { + tolerance = atof( argv[3] ); + } + double b = Gudhi::Bottleneck_distance::compute(diag1, diag2, tolerance); + std::cout << "The distance between the diagrams is : " << b << ". The tolerace is : " << tolerance << std::endl; } -- cgit v1.2.3