From 0d0ca116e7fef77cc950b7e85380495661311d91 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 22 Mar 2018 08:36:00 +0000 Subject: Move Miniball in gudhi git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3302 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 736e5c33a9593e4dfb5b19ddc745ab9852d71b40 --- src/CMakeLists.txt | 1 - .../benchmark/cech_complex_benchmark.cpp | 5 +- src/Cech_complex/example/CMakeLists.txt | 2 + .../example/cech_complex_step_by_step.cpp | 4 +- .../include/Miniball/Miniball.COPYRIGHT | 4 - src/Cech_complex/include/Miniball/Miniball.README | 23 - src/Cech_complex/include/Miniball/Miniball.hpp | 520 -------------------- src/Cech_complex/include/gudhi/Miniball.COPYRIGHT | 4 + src/Cech_complex/include/gudhi/Miniball.README | 23 + src/Cech_complex/include/gudhi/Miniball.hpp | 523 +++++++++++++++++++++ src/Cech_complex/test/test_cech_complex.cpp | 5 +- src/common/include/gudhi/distance_functions.h | 2 +- 12 files changed, 559 insertions(+), 557 deletions(-) delete mode 100644 src/Cech_complex/include/Miniball/Miniball.COPYRIGHT delete mode 100644 src/Cech_complex/include/Miniball/Miniball.README delete mode 100644 src/Cech_complex/include/Miniball/Miniball.hpp create mode 100644 src/Cech_complex/include/gudhi/Miniball.COPYRIGHT create mode 100644 src/Cech_complex/include/gudhi/Miniball.README create mode 100644 src/Cech_complex/include/gudhi/Miniball.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ff7acafc..7cccb19f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -105,5 +105,4 @@ install(FILES # install the include file on "make install" install(DIRECTORY include/gudhi DESTINATION include) -install(DIRECTORY include/Miniball DESTINATION include) #--------------------------------------------------------------------------------------- diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 83ef9dca..6ba52419 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -27,8 +27,7 @@ #include #include #include - -#include +#include #include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations @@ -57,7 +56,7 @@ class Radius_distance { using Point_cloud = std::vector; using Point_iterator = typename Point_cloud::const_iterator; using Coordinate_iterator = typename Point::const_iterator; - using Min_sphere = typename Miniball::Miniball>; + using Min_sphere = typename Gudhi::Miniball::Miniball>; Point_cloud point_cloud; point_cloud.push_back(p1); diff --git a/src/Cech_complex/example/CMakeLists.txt b/src/Cech_complex/example/CMakeLists.txt index ac32ff95..ab391215 100644 --- a/src/Cech_complex/example/CMakeLists.txt +++ b/src/Cech_complex/example/CMakeLists.txt @@ -6,6 +6,8 @@ target_link_libraries(Cech_complex_example_step_by_step ${Boost_PROGRAM_OPTIONS_ if (TBB_FOUND) target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES}) endif() +add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3") add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp) if (TBB_FOUND) diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index 8705a3e5..0e7c4bbd 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -25,7 +25,7 @@ #include #include -#include +#include #include @@ -55,7 +55,7 @@ class Cech_blocker { using Point_cloud = std::vector; using Point_iterator = Point_cloud::const_iterator; using Coordinate_iterator = Point::const_iterator; - using Min_sphere = Miniball::Miniball >; + using Min_sphere = Gudhi::Miniball::Miniball >; public: bool operator()(Simplex_handle sh) { std::vector points; diff --git a/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT b/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT deleted file mode 100644 index dbe4c553..00000000 --- a/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT +++ /dev/null @@ -1,4 +0,0 @@ -The miniball software is available under the GNU General Public License (GPLv3 - https://www.gnu.org/copyleft/gpl.html). -If your intended use is not compliant with this license, please buy a commercial license (EUR 500 - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/license.html). -You need a license if the software that you develop using Miniball V3.0 is not open source. - diff --git a/src/Cech_complex/include/Miniball/Miniball.README b/src/Cech_complex/include/Miniball/Miniball.README deleted file mode 100644 index 86a96f08..00000000 --- a/src/Cech_complex/include/Miniball/Miniball.README +++ /dev/null @@ -1,23 +0,0 @@ -https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html - -Smallest Enclosing Balls of Points - Fast and Robust in C++. -(high-quality software for smallest enclosing balls of balls is available in the computational geometry algorithms library CGAL) - - -This is the miniball software (V3.0) for computing smallest enclosing balls of points in arbitrary dimensions. It consists of a C++ header file Miniball.hpp (around 500 lines of code) and two example programs miniball_example.cpp and miniball_example_containers.cpp that demonstrate the usage. The first example stores the coordinates of the input points in a two-dimensional array, the second example uses a list of vectors to show how generic containers can be used. - -Credits: Aditya Gupta and Alexandros Konstantinakis-Karmis have significantly contributed to this version of the software. - -Changes - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/changes.txt - from previous versions. - -The theory - https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf - behind the miniball software (Proc. 7th Annual European Symposium on Algorithms (ESA), Lecture Notes in Computer Science 1643, Springer-Verlag, pp.325-338, 1999). - -Main Features: - - Very fast in low dimensions. 1 million points in 5-space are processed within 0.05 seconds on any recent machine. - - High numerical stability. Almost all input degeneracies (cospherical points, multiple points, points very close together) are routinely handled. - - Easily integrates into your code. You can freely choose the coordinate type of your points and the container to store the points. If you still need to adapt the code, the header is small and readable and contains documentation for all major methods. - - diff --git a/src/Cech_complex/include/Miniball/Miniball.hpp b/src/Cech_complex/include/Miniball/Miniball.hpp deleted file mode 100644 index a42d62a7..00000000 --- a/src/Cech_complex/include/Miniball/Miniball.hpp +++ /dev/null @@ -1,520 +0,0 @@ -// Copright (C) 1999-2013, Bernd Gaertner -// $Rev: 3581 $ -// -// 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 . -// -// Contact: -// -------- -// Bernd Gaertner -// Institute of Theoretical Computer Science -// ETH Zuerich -// CAB G31.1 -// CH-8092 Zuerich, Switzerland -// http://www.inf.ethz.ch/personal/gaertner - -#ifndef MINIBALL_HPP_ -#define MINIBALL_HPP_ - -#include -#include -#include -#include -#include - -namespace Miniball { - - // Global Functions - // ================ - template - inline NT mb_sqr (NT r) {return r*r;} - - // Functors - // ======== - - // functor to map a point iterator to the corresponding coordinate iterator; - // generic version for points whose coordinate containers have begin() - template < typename Pit_, typename Cit_ > - struct CoordAccessor { - typedef Pit_ Pit; - typedef Cit_ Cit; - inline Cit operator() (Pit it) const { return (*it).begin(); } - }; - - // partial specialization for points whose coordinate containers are arrays - template < typename Pit_, typename Cit_ > - struct CoordAccessor { - typedef Pit_ Pit; - typedef Cit_* Cit; - inline Cit operator() (Pit it) const { return *it; } - }; - - // Class Declaration - // ================= - - template - class Miniball { - private: - // types - // The iterator type to go through the input points - typedef typename CoordAccessor::Pit Pit; - // The iterator type to go through the coordinates of a single point. - typedef typename CoordAccessor::Cit Cit; - // The coordinate type - typedef typename std::iterator_traits::value_type NT; - // The iterator to go through the support points - typedef typename std::list::iterator Sit; - - // data members... - const int d; // dimension - Pit points_begin; - Pit points_end; - CoordAccessor coord_accessor; - double time; - const NT nt0; // NT(0) - - //...for the algorithms - std::list L; - Sit support_end; - int fsize; // number of forced points - int ssize; // number of support points - - // ...for the ball updates - NT* current_c; - NT current_sqr_r; - NT** c; - NT* sqr_r; - - // helper arrays - NT* q0; - NT* z; - NT* f; - NT** v; - NT** a; - - public: - // The iterator type to go through the support points - typedef typename std::list::const_iterator SupportPointIterator; - - // PRE: [begin, end) is a nonempty range - // POST: computes the smallest enclosing ball of the points in the range - // [begin, end); the functor a maps a point iterator to an iterator - // through the d coordinates of the point - Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor()); - - // POST: returns a pointer to the first element of an array that holds - // the d coordinates of the center of the computed ball - const NT* center () const; - - // POST: returns the squared radius of the computed ball - NT squared_radius () const; - - // POST: returns the number of support points of the computed ball; - // the support points form a minimal set with the same smallest - // enclosing ball as the input set; in particular, the support - // points are on the boundary of the computed ball, and their - // number is at most d+1 - int nr_support_points () const; - - // POST: returns an iterator to the first support point - SupportPointIterator support_points_begin () const; - - // POST: returns a past-the-end iterator for the range of support points - SupportPointIterator support_points_end () const; - - // POST: returns the maximum excess of any input point w.r.t. the computed - // ball, divided by the squared radius of the computed ball. The - // excess of a point is the difference between its squared distance - // from the center and the squared radius; Ideally, the return value - // is 0. subopt is set to the absolute value of the most negative - // coefficient in the affine combination of the support points that - // yields the center. Ideally, this is a convex combination, and there - // is no negative coefficient in which case subopt is set to 0. - NT relative_error (NT& subopt) const; - - // POST: return true if the relative error is at most tol, and the - // suboptimality is 0; the default tolerance is 10 times the - // coordinate type's machine epsilon - bool is_valid (NT tol = NT(10) * std::numeric_limits::epsilon()) const; - - // POST: returns the time in seconds taken by the constructor call for - // computing the smallest enclosing ball - double get_time() const; - - // POST: deletes dynamically allocated arrays - ~Miniball(); - - private: - void mtf_mb (Sit n); - void mtf_move_to_front (Sit j); - void pivot_mb (Pit n); - void pivot_move_to_front (Pit j); - NT excess (Pit pit) const; - void pop (); - bool push (Pit pit); - NT suboptimality () const; - void create_arrays(); - void delete_arrays(); - }; - - // Class Definition - // ================ - template - Miniball::Miniball (int d_, Pit begin, Pit end, - CoordAccessor ca) - : d (d_), - points_begin (begin), - points_end (end), - coord_accessor (ca), - time (clock()), - nt0 (NT(0)), - L(), - support_end (L.begin()), - fsize(0), - ssize(0), - current_c (NULL), - current_sqr_r (NT(-1)), - c (NULL), - sqr_r (NULL), - q0 (NULL), - z (NULL), - f (NULL), - v (NULL), - a (NULL) - { - assert (points_begin != points_end); - create_arrays(); - - // set initial center - for (int j=0; j - Miniball::~Miniball() - { - delete_arrays(); - } - - template - void Miniball::create_arrays() - { - c = new NT*[d+1]; - v = new NT*[d+1]; - a = new NT*[d+1]; - for (int i=0; i - void Miniball::delete_arrays() - { - delete[] f; - delete[] z; - delete[] q0; - delete[] sqr_r; - for (int i=0; i - const typename Miniball::NT* - Miniball::center () const - { - return current_c; - } - - template - typename Miniball::NT - Miniball::squared_radius () const - { - return current_sqr_r; - } - - template - int Miniball::nr_support_points () const - { - assert (ssize < d+2); - return ssize; - } - - template - typename Miniball::SupportPointIterator - Miniball::support_points_begin () const - { - return L.begin(); - } - - template - typename Miniball::SupportPointIterator - Miniball::support_points_end () const - { - return support_end; - } - - template - typename Miniball::NT - Miniball::relative_error (NT& subopt) const - { - NT e, max_e = nt0; - // compute maximum absolute excess of support points - for (SupportPointIterator it = support_points_begin(); - it != support_points_end(); ++it) { - e = excess (*it); - if (e < nt0) e = -e; - if (e > max_e) { - max_e = e; - } - } - // compute maximum excess of any point - for (Pit i = points_begin; i != points_end; ++i) - if ((e = excess (i)) > max_e) - max_e = e; - - subopt = suboptimality(); - assert (current_sqr_r > nt0 || max_e == nt0); - return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r); - } - - template - bool Miniball::is_valid (NT tol) const - { - NT suboptimality; - return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); - } - - template - double Miniball::get_time() const - { - return time; - } - - template - void Miniball::mtf_mb (Sit n) - { - // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n) - // B: the set of forced points, defining the current ball - // S: the superset of support points computed by the algorithm - // -------------------------------------------------------------- - // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, - // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf - - // PRE: B = S - assert (fsize == ssize); - - support_end = L.begin(); - if ((fsize) == d+1) return; - - // incremental construction - for (Sit i = L.begin(); i != n;) - { - // INV: (support_end - L.begin() == |S|-|B|) - assert (std::distance (L.begin(), support_end) == ssize - fsize); - - Sit j = i++; - if (excess(*j) > nt0) - if (push(*j)) { // B := B + p_i - mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i) - pop(); // B := B - p_i - mtf_move_to_front(j); - } - } - // POST: the range [L.begin(), support_end) stores the set S\B - } - - template - void Miniball::mtf_move_to_front (Sit j) - { - if (support_end == j) - support_end++; - L.splice (L.begin(), L, j); - } - - template - void Miniball::pivot_mb (Pit n) - { - // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n) - // -------------------------------------------------------------- - // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, - // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf - NT old_sqr_r; - const NT* c; - Pit pivot, k; - NT e, max_e, sqr_r; - Cit p; - do { - old_sqr_r = current_sqr_r; - sqr_r = current_sqr_r; - - pivot = points_begin; - max_e = nt0; - for (k = points_begin; k != n; ++k) { - p = coord_accessor(k); - e = -sqr_r; - c = current_c; - for (int j=0; j(*p++-*c++); - if (e > max_e) { - max_e = e; - pivot = k; - } - } - - if (max_e > nt0) { - // check if the pivot is already contained in the support set - if (std::find(L.begin(), support_end, pivot) == support_end) { - assert (fsize == 0); - if (push (pivot)) { - mtf_mb(support_end); - pop(); - pivot_move_to_front(pivot); - } - } - } - } while (old_sqr_r < current_sqr_r); - } - - template - void Miniball::pivot_move_to_front (Pit j) - { - L.push_front(j); - if (std::distance(L.begin(), support_end) == d+2) - support_end--; - } - - template - inline typename Miniball::NT - Miniball::excess (Pit pit) const - { - Cit p = coord_accessor(pit); - NT e = -current_sqr_r; - NT* c = current_c; - for (int k=0; k(*p++-*c++); - } - return e; - } - - template - void Miniball::pop () - { - --fsize; - } - - template - bool Miniball::push (Pit pit) - { - int i, j; - NT eps = mb_sqr(std::numeric_limits::epsilon()); - - Cit cit = coord_accessor(pit); - Cit p = cit; - - if (fsize==0) { - for (i=0; i(v[fsize][j]); - z[fsize]*=2; - - // reject push if z_fsize too small - if (z[fsize](*p++-c[fsize-1][i]); - f[fsize]=e/z[fsize]; - - for (i=0; i - typename Miniball::NT - Miniball::suboptimality () const - { - NT* l = new NT[d+1]; - NT min_l = nt0; - l[0] = NT(1); - for (int i=ssize-1; i>0; --i) { - l[i] = f[i]; - for (int k=ssize-1; k>i; --k) - l[i]-=a[k][i]*l[k]; - if (l[i] < min_l) min_l = l[i]; - l[0] -= l[i]; - } - if (l[0] < min_l) min_l = l[0]; - delete[] l; - if (min_l < nt0) - return -min_l; - return nt0; - } - -} // end Namespace Miniball - -#endif // MINIBALL_HPP_ diff --git a/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT new file mode 100644 index 00000000..dbe4c553 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT @@ -0,0 +1,4 @@ +The miniball software is available under the GNU General Public License (GPLv3 - https://www.gnu.org/copyleft/gpl.html). +If your intended use is not compliant with this license, please buy a commercial license (EUR 500 - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/license.html). +You need a license if the software that you develop using Miniball V3.0 is not open source. + diff --git a/src/Cech_complex/include/gudhi/Miniball.README b/src/Cech_complex/include/gudhi/Miniball.README new file mode 100644 index 00000000..86a96f08 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.README @@ -0,0 +1,23 @@ +https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html + +Smallest Enclosing Balls of Points - Fast and Robust in C++. +(high-quality software for smallest enclosing balls of balls is available in the computational geometry algorithms library CGAL) + + +This is the miniball software (V3.0) for computing smallest enclosing balls of points in arbitrary dimensions. It consists of a C++ header file Miniball.hpp (around 500 lines of code) and two example programs miniball_example.cpp and miniball_example_containers.cpp that demonstrate the usage. The first example stores the coordinates of the input points in a two-dimensional array, the second example uses a list of vectors to show how generic containers can be used. + +Credits: Aditya Gupta and Alexandros Konstantinakis-Karmis have significantly contributed to this version of the software. + +Changes - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/changes.txt - from previous versions. + +The theory - https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf - behind the miniball software (Proc. 7th Annual European Symposium on Algorithms (ESA), Lecture Notes in Computer Science 1643, Springer-Verlag, pp.325-338, 1999). + +Main Features: + + Very fast in low dimensions. 1 million points in 5-space are processed within 0.05 seconds on any recent machine. + + High numerical stability. Almost all input degeneracies (cospherical points, multiple points, points very close together) are routinely handled. + + Easily integrates into your code. You can freely choose the coordinate type of your points and the container to store the points. If you still need to adapt the code, the header is small and readable and contains documentation for all major methods. + + diff --git a/src/Cech_complex/include/gudhi/Miniball.hpp b/src/Cech_complex/include/gudhi/Miniball.hpp new file mode 100644 index 00000000..ce6cbb5b --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.hpp @@ -0,0 +1,523 @@ +// Copright (C) 1999-2013, Bernd Gaertner +// $Rev: 3581 $ +// +// 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 . +// +// Contact: +// -------- +// Bernd Gaertner +// Institute of Theoretical Computer Science +// ETH Zuerich +// CAB G31.1 +// CH-8092 Zuerich, Switzerland +// http://www.inf.ethz.ch/personal/gaertner + +#ifndef MINIBALL_HPP_ +#define MINIBALL_HPP_ + +#include +#include +#include +#include +#include + +namespace Gudhi { + +namespace Miniball { + + // Global Functions + // ================ + template + inline NT mb_sqr (NT r) {return r*r;} + + // Functors + // ======== + + // functor to map a point iterator to the corresponding coordinate iterator; + // generic version for points whose coordinate containers have begin() + template < typename Pit_, typename Cit_ > + struct CoordAccessor { + typedef Pit_ Pit; + typedef Cit_ Cit; + inline Cit operator() (Pit it) const { return (*it).begin(); } + }; + + // partial specialization for points whose coordinate containers are arrays + template < typename Pit_, typename Cit_ > + struct CoordAccessor { + typedef Pit_ Pit; + typedef Cit_* Cit; + inline Cit operator() (Pit it) const { return *it; } + }; + + // Class Declaration + // ================= + + template + class Miniball { + private: + // types + // The iterator type to go through the input points + typedef typename CoordAccessor::Pit Pit; + // The iterator type to go through the coordinates of a single point. + typedef typename CoordAccessor::Cit Cit; + // The coordinate type + typedef typename std::iterator_traits::value_type NT; + // The iterator to go through the support points + typedef typename std::list::iterator Sit; + + // data members... + const int d; // dimension + Pit points_begin; + Pit points_end; + CoordAccessor coord_accessor; + double time; + const NT nt0; // NT(0) + + //...for the algorithms + std::list L; + Sit support_end; + int fsize; // number of forced points + int ssize; // number of support points + + // ...for the ball updates + NT* current_c; + NT current_sqr_r; + NT** c; + NT* sqr_r; + + // helper arrays + NT* q0; + NT* z; + NT* f; + NT** v; + NT** a; + + public: + // The iterator type to go through the support points + typedef typename std::list::const_iterator SupportPointIterator; + + // PRE: [begin, end) is a nonempty range + // POST: computes the smallest enclosing ball of the points in the range + // [begin, end); the functor a maps a point iterator to an iterator + // through the d coordinates of the point + Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor()); + + // POST: returns a pointer to the first element of an array that holds + // the d coordinates of the center of the computed ball + const NT* center () const; + + // POST: returns the squared radius of the computed ball + NT squared_radius () const; + + // POST: returns the number of support points of the computed ball; + // the support points form a minimal set with the same smallest + // enclosing ball as the input set; in particular, the support + // points are on the boundary of the computed ball, and their + // number is at most d+1 + int nr_support_points () const; + + // POST: returns an iterator to the first support point + SupportPointIterator support_points_begin () const; + + // POST: returns a past-the-end iterator for the range of support points + SupportPointIterator support_points_end () const; + + // POST: returns the maximum excess of any input point w.r.t. the computed + // ball, divided by the squared radius of the computed ball. The + // excess of a point is the difference between its squared distance + // from the center and the squared radius; Ideally, the return value + // is 0. subopt is set to the absolute value of the most negative + // coefficient in the affine combination of the support points that + // yields the center. Ideally, this is a convex combination, and there + // is no negative coefficient in which case subopt is set to 0. + NT relative_error (NT& subopt) const; + + // POST: return true if the relative error is at most tol, and the + // suboptimality is 0; the default tolerance is 10 times the + // coordinate type's machine epsilon + bool is_valid (NT tol = NT(10) * std::numeric_limits::epsilon()) const; + + // POST: returns the time in seconds taken by the constructor call for + // computing the smallest enclosing ball + double get_time() const; + + // POST: deletes dynamically allocated arrays + ~Miniball(); + + private: + void mtf_mb (Sit n); + void mtf_move_to_front (Sit j); + void pivot_mb (Pit n); + void pivot_move_to_front (Pit j); + NT excess (Pit pit) const; + void pop (); + bool push (Pit pit); + NT suboptimality () const; + void create_arrays(); + void delete_arrays(); + }; + + // Class Definition + // ================ + template + Miniball::Miniball (int d_, Pit begin, Pit end, + CoordAccessor ca) + : d (d_), + points_begin (begin), + points_end (end), + coord_accessor (ca), + time (clock()), + nt0 (NT(0)), + L(), + support_end (L.begin()), + fsize(0), + ssize(0), + current_c (NULL), + current_sqr_r (NT(-1)), + c (NULL), + sqr_r (NULL), + q0 (NULL), + z (NULL), + f (NULL), + v (NULL), + a (NULL) + { + assert (points_begin != points_end); + create_arrays(); + + // set initial center + for (int j=0; j + Miniball::~Miniball() + { + delete_arrays(); + } + + template + void Miniball::create_arrays() + { + c = new NT*[d+1]; + v = new NT*[d+1]; + a = new NT*[d+1]; + for (int i=0; i + void Miniball::delete_arrays() + { + delete[] f; + delete[] z; + delete[] q0; + delete[] sqr_r; + for (int i=0; i + const typename Miniball::NT* + Miniball::center () const + { + return current_c; + } + + template + typename Miniball::NT + Miniball::squared_radius () const + { + return current_sqr_r; + } + + template + int Miniball::nr_support_points () const + { + assert (ssize < d+2); + return ssize; + } + + template + typename Miniball::SupportPointIterator + Miniball::support_points_begin () const + { + return L.begin(); + } + + template + typename Miniball::SupportPointIterator + Miniball::support_points_end () const + { + return support_end; + } + + template + typename Miniball::NT + Miniball::relative_error (NT& subopt) const + { + NT e, max_e = nt0; + // compute maximum absolute excess of support points + for (SupportPointIterator it = support_points_begin(); + it != support_points_end(); ++it) { + e = excess (*it); + if (e < nt0) e = -e; + if (e > max_e) { + max_e = e; + } + } + // compute maximum excess of any point + for (Pit i = points_begin; i != points_end; ++i) + if ((e = excess (i)) > max_e) + max_e = e; + + subopt = suboptimality(); + assert (current_sqr_r > nt0 || max_e == nt0); + return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r); + } + + template + bool Miniball::is_valid (NT tol) const + { + NT suboptimality; + return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); + } + + template + double Miniball::get_time() const + { + return time; + } + + template + void Miniball::mtf_mb (Sit n) + { + // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n) + // B: the set of forced points, defining the current ball + // S: the superset of support points computed by the algorithm + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + + // PRE: B = S + assert (fsize == ssize); + + support_end = L.begin(); + if ((fsize) == d+1) return; + + // incremental construction + for (Sit i = L.begin(); i != n;) + { + // INV: (support_end - L.begin() == |S|-|B|) + assert (std::distance (L.begin(), support_end) == ssize - fsize); + + Sit j = i++; + if (excess(*j) > nt0) + if (push(*j)) { // B := B + p_i + mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i) + pop(); // B := B - p_i + mtf_move_to_front(j); + } + } + // POST: the range [L.begin(), support_end) stores the set S\B + } + + template + void Miniball::mtf_move_to_front (Sit j) + { + if (support_end == j) + support_end++; + L.splice (L.begin(), L, j); + } + + template + void Miniball::pivot_mb (Pit n) + { + // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n) + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + NT old_sqr_r; + const NT* c; + Pit pivot, k; + NT e, max_e, sqr_r; + Cit p; + do { + old_sqr_r = current_sqr_r; + sqr_r = current_sqr_r; + + pivot = points_begin; + max_e = nt0; + for (k = points_begin; k != n; ++k) { + p = coord_accessor(k); + e = -sqr_r; + c = current_c; + for (int j=0; j(*p++-*c++); + if (e > max_e) { + max_e = e; + pivot = k; + } + } + + if (max_e > nt0) { + // check if the pivot is already contained in the support set + if (std::find(L.begin(), support_end, pivot) == support_end) { + assert (fsize == 0); + if (push (pivot)) { + mtf_mb(support_end); + pop(); + pivot_move_to_front(pivot); + } + } + } + } while (old_sqr_r < current_sqr_r); + } + + template + void Miniball::pivot_move_to_front (Pit j) + { + L.push_front(j); + if (std::distance(L.begin(), support_end) == d+2) + support_end--; + } + + template + inline typename Miniball::NT + Miniball::excess (Pit pit) const + { + Cit p = coord_accessor(pit); + NT e = -current_sqr_r; + NT* c = current_c; + for (int k=0; k(*p++-*c++); + } + return e; + } + + template + void Miniball::pop () + { + --fsize; + } + + template + bool Miniball::push (Pit pit) + { + int i, j; + NT eps = mb_sqr(std::numeric_limits::epsilon()); + + Cit cit = coord_accessor(pit); + Cit p = cit; + + if (fsize==0) { + for (i=0; i(v[fsize][j]); + z[fsize]*=2; + + // reject push if z_fsize too small + if (z[fsize](*p++-c[fsize-1][i]); + f[fsize]=e/z[fsize]; + + for (i=0; i + typename Miniball::NT + Miniball::suboptimality () const + { + NT* l = new NT[d+1]; + NT min_l = nt0; + l[0] = NT(1); + for (int i=ssize-1; i>0; --i) { + l[i] = f[i]; + for (int k=ssize-1; k>i; --k) + l[i]-=a[k][i]*l[k]; + if (l[i] < min_l) min_l = l[i]; + l[0] -= l[i]; + } + if (l[0] < min_l) min_l = l[0]; + delete[] l; + if (min_l < nt0) + return -min_l; + return nt0; + } +} // namespace Miniball + +} // namespace Gudhi + +#endif // MINIBALL_HPP_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 8cbfe431..4aa85057 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -36,8 +36,7 @@ #include #include #include - -#include +#include // Type definitions using Simplex_tree = Gudhi::Simplex_tree<>; @@ -49,7 +48,7 @@ using Cech_complex = Gudhi::cech_complex::Cech_complex>; +using Min_sphere = Gudhi::Miniball::Miniball>; BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index 7749b98e..93956a69 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -25,7 +25,7 @@ #include -#include +#include #include -- cgit v1.2.3