diff options
Diffstat (limited to 'include/gudhi/Miniball.hpp')
| -rw-r--r-- | include/gudhi/Miniball.hpp | 523 |
1 files changed, 0 insertions, 523 deletions
diff --git a/include/gudhi/Miniball.hpp b/include/gudhi/Miniball.hpp deleted file mode 100644 index ce6cbb5b..00000000 --- a/include/gudhi/Miniball.hpp +++ /dev/null @@ -1,523 +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 <http://www.gnu.org/licenses/>. -// -// 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 <cassert> -#include <algorithm> -#include <list> -#include <ctime> -#include <limits> - -namespace Gudhi { - -namespace Miniball { - - // Global Functions - // ================ - template <typename NT> - 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<Pit_, Cit_*> { - typedef Pit_ Pit; - typedef Cit_* Cit; - inline Cit operator() (Pit it) const { return *it; } - }; - - // Class Declaration - // ================= - - template <typename CoordAccessor> - 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<Cit>::value_type NT; - // The iterator to go through the support points - typedef typename std::list<Pit>::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<Pit> 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<Pit>::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<NT>::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 <typename CoordAccessor> - Miniball<CoordAccessor>::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<d; ++j) c[0][j] = nt0; - current_c = c[0]; - - // compute miniball - pivot_mb (points_end); - - // update time - time = (clock() - time) / CLOCKS_PER_SEC; - } - - template <typename CoordAccessor> - Miniball<CoordAccessor>::~Miniball() - { - delete_arrays(); - } - - template <typename CoordAccessor> - void Miniball<CoordAccessor>::create_arrays() - { - c = new NT*[d+1]; - v = new NT*[d+1]; - a = new NT*[d+1]; - for (int i=0; i<d+1; ++i) { - c[i] = new NT[d]; - v[i] = new NT[d]; - a[i] = new NT[d]; - } - sqr_r = new NT[d+1]; - q0 = new NT[d]; - z = new NT[d+1]; - f = new NT[d+1]; - } - - template <typename CoordAccessor> - void Miniball<CoordAccessor>::delete_arrays() - { - delete[] f; - delete[] z; - delete[] q0; - delete[] sqr_r; - for (int i=0; i<d+1; ++i) { - delete[] a[i]; - delete[] v[i]; - delete[] c[i]; - } - delete[] a; - delete[] v; - delete[] c; - } - - template <typename CoordAccessor> - const typename Miniball<CoordAccessor>::NT* - Miniball<CoordAccessor>::center () const - { - return current_c; - } - - template <typename CoordAccessor> - typename Miniball<CoordAccessor>::NT - Miniball<CoordAccessor>::squared_radius () const - { - return current_sqr_r; - } - - template <typename CoordAccessor> - int Miniball<CoordAccessor>::nr_support_points () const - { - assert (ssize < d+2); - return ssize; - } - - template <typename CoordAccessor> - typename Miniball<CoordAccessor>::SupportPointIterator - Miniball<CoordAccessor>::support_points_begin () const - { - return L.begin(); - } - - template <typename CoordAccessor> - typename Miniball<CoordAccessor>::SupportPointIterator - Miniball<CoordAccessor>::support_points_end () const - { - return support_end; - } - - template <typename CoordAccessor> - typename Miniball<CoordAccessor>::NT - Miniball<CoordAccessor>::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 <typename CoordAccessor> - bool Miniball<CoordAccessor>::is_valid (NT tol) const - { - NT suboptimality; - return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); - } - - template <typename CoordAccessor> - double Miniball<CoordAccessor>::get_time() const - { - return time; - } - - template <typename CoordAccessor> - void Miniball<CoordAccessor>::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 <typename CoordAccessor> - void Miniball<CoordAccessor>::mtf_move_to_front (Sit j) - { - if (support_end == j) - support_end++; - L.splice (L.begin(), L, j); - } - - template <typename CoordAccessor> - void Miniball<CoordAccessor>::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<d; ++j) - e += mb_sqr<NT>(*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 <typename CoordAccessor> - void Miniball<CoordAccessor>::pivot_move_to_front (Pit j) - { - L.push_front(j); - if (std::distance(L.begin(), support_end) == d+2) - support_end--; - } - - template <typename CoordAccessor> - inline typename Miniball<CoordAccessor>::NT - Miniball<CoordAccessor>::excess (Pit pit) const - { - Cit p = coord_accessor(pit); - NT e = -current_sqr_r; - NT* c = current_c; - for (int k=0; k<d; ++k){ - e += mb_sqr<NT>(*p++-*c++); - } - return e; - } - - template <typename CoordAccessor> - void Miniball<CoordAccessor>::pop () - { - --fsize; - } - - template <typename CoordAccessor> - bool Miniball<CoordAccessor>::push (Pit pit) - { - int i, j; - NT eps = mb_sqr<NT>(std::numeric_limits<NT>::epsilon()); - - Cit cit = coord_accessor(pit); - Cit p = cit; - - if (fsize==0) { - for (i=0; i<d; ++i) - q0[i] = *p++; - for (i=0; i<d; ++i) - c[0][i] = q0[i]; - sqr_r[0] = nt0; - } - else { - // set v_fsize to Q_fsize - for (i=0; i<d; ++i) - //v[fsize][i] = p[i]-q0[i]; - v[fsize][i] = *p++-q0[i]; - - // compute the a_{fsize,i}, i< fsize - for (i=1; i<fsize; ++i) { - a[fsize][i] = nt0; - for (j=0; j<d; ++j) - a[fsize][i] += v[i][j] * v[fsize][j]; - a[fsize][i]*=(2/z[i]); - } - - // update v_fsize to Q_fsize-\bar{Q}_fsize - for (i=1; i<fsize; ++i) { - for (j=0; j<d; ++j) - v[fsize][j] -= a[fsize][i]*v[i][j]; - } - - // compute z_fsize - z[fsize]=nt0; - for (j=0; j<d; ++j) - z[fsize] += mb_sqr<NT>(v[fsize][j]); - z[fsize]*=2; - - // reject push if z_fsize too small - if (z[fsize]<eps*current_sqr_r) { - return false; - } - - // update c, sqr_r - p=cit; - NT e = -sqr_r[fsize-1]; - for (i=0; i<d; ++i) - e += mb_sqr<NT>(*p++-c[fsize-1][i]); - f[fsize]=e/z[fsize]; - - for (i=0; i<d; ++i) - c[fsize][i] = c[fsize-1][i]+f[fsize]*v[fsize][i]; - sqr_r[fsize] = sqr_r[fsize-1] + e*f[fsize]/2; - } - current_c = c[fsize]; - current_sqr_r = sqr_r[fsize]; - ssize = ++fsize; - return true; - } - - template <typename CoordAccessor> - typename Miniball<CoordAccessor>::NT - Miniball<CoordAccessor>::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_ |