// Copright (C) 1999-2021, Bernd Gaertner // November 12, 2021 // // 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 #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; // by how much do we allow points outside? NT default_tol; 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 () 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), default_tol (NT(10) * std::numeric_limits::epsilon()) { 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 () const { NT suboptimality; return ( (relative_error (suboptimality) <= default_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 const NT* c; Pit pivot, k; NT e, max_e, sqr_r; Cit p; unsigned int loops_without_progress = 0; NT best_sqr_r = current_sqr_r; do { 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 (sqr_r < nt0 || 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); } } } if (best_sqr_r < current_sqr_r) { best_sqr_r = current_sqr_r; loops_without_progress = 0; } else ++loops_without_progress; } while (loops_without_progress < 2); } 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