Miniball.hpp
1 // Copright (C) 1999-2013, Bernd Gaertner
2 // $Rev: 3581 $
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 
14 // You should have received a copy of the GNU General Public License
15 // along with this program. If not, see <http://www.gnu.org/licenses/>.
16 //
17 // Contact:
18 // --------
19 // Bernd Gaertner
20 // Institute of Theoretical Computer Science
21 // ETH Zuerich
22 // CAB G31.1
23 // CH-8092 Zuerich, Switzerland
24 // http://www.inf.ethz.ch/personal/gaertner
25 
26 #ifndef MINIBALL_HPP_
27 #define MINIBALL_HPP_
28 
29 #include <cassert>
30 #include <algorithm>
31 #include <list>
32 #include <ctime>
33 #include <limits>
34 
35 namespace Gudhi {
36 
37 namespace Miniball {
38 
39  // Global Functions
40  // ================
41  template <typename NT>
42  inline NT mb_sqr (NT r) {return r*r;}
43 
44  // Functors
45  // ========
46 
47  // functor to map a point iterator to the corresponding coordinate iterator;
48  // generic version for points whose coordinate containers have begin()
49  template < typename Pit_, typename Cit_ >
50  struct CoordAccessor {
51  typedef Pit_ Pit;
52  typedef Cit_ Cit;
53  inline Cit operator() (Pit it) const { return (*it).begin(); }
54  };
55 
56  // partial specialization for points whose coordinate containers are arrays
57  template < typename Pit_, typename Cit_ >
58  struct CoordAccessor<Pit_, Cit_*> {
59  typedef Pit_ Pit;
60  typedef Cit_* Cit;
61  inline Cit operator() (Pit it) const { return *it; }
62  };
63 
64  // Class Declaration
65  // =================
66 
67  template <typename CoordAccessor>
68  class Miniball {
69  private:
70  // types
71  // The iterator type to go through the input points
72  typedef typename CoordAccessor::Pit Pit;
73  // The iterator type to go through the coordinates of a single point.
74  typedef typename CoordAccessor::Cit Cit;
75  // The coordinate type
76  typedef typename std::iterator_traits<Cit>::value_type NT;
77  // The iterator to go through the support points
78  typedef typename std::list<Pit>::iterator Sit;
79 
80  // data members...
81  const int d; // dimension
82  Pit points_begin;
83  Pit points_end;
84  CoordAccessor coord_accessor;
85  double time;
86  const NT nt0; // NT(0)
87 
88  //...for the algorithms
89  std::list<Pit> L;
90  Sit support_end;
91  int fsize; // number of forced points
92  int ssize; // number of support points
93 
94  // ...for the ball updates
95  NT* current_c;
96  NT current_sqr_r;
97  NT** c;
98  NT* sqr_r;
99 
100  // helper arrays
101  NT* q0;
102  NT* z;
103  NT* f;
104  NT** v;
105  NT** a;
106 
107  public:
108  // The iterator type to go through the support points
109  typedef typename std::list<Pit>::const_iterator SupportPointIterator;
110 
111  // PRE: [begin, end) is a nonempty range
112  // POST: computes the smallest enclosing ball of the points in the range
113  // [begin, end); the functor a maps a point iterator to an iterator
114  // through the d coordinates of the point
115  Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor());
116 
117  // POST: returns a pointer to the first element of an array that holds
118  // the d coordinates of the center of the computed ball
119  const NT* center () const;
120 
121  // POST: returns the squared radius of the computed ball
122  NT squared_radius () const;
123 
124  // POST: returns the number of support points of the computed ball;
125  // the support points form a minimal set with the same smallest
126  // enclosing ball as the input set; in particular, the support
127  // points are on the boundary of the computed ball, and their
128  // number is at most d+1
129  int nr_support_points () const;
130 
131  // POST: returns an iterator to the first support point
132  SupportPointIterator support_points_begin () const;
133 
134  // POST: returns a past-the-end iterator for the range of support points
135  SupportPointIterator support_points_end () const;
136 
137  // POST: returns the maximum excess of any input point w.r.t. the computed
138  // ball, divided by the squared radius of the computed ball. The
139  // excess of a point is the difference between its squared distance
140  // from the center and the squared radius; Ideally, the return value
141  // is 0. subopt is set to the absolute value of the most negative
142  // coefficient in the affine combination of the support points that
143  // yields the center. Ideally, this is a convex combination, and there
144  // is no negative coefficient in which case subopt is set to 0.
145  NT relative_error (NT& subopt) const;
146 
147  // POST: return true if the relative error is at most tol, and the
148  // suboptimality is 0; the default tolerance is 10 times the
149  // coordinate type's machine epsilon
150  bool is_valid (NT tol = NT(10) * std::numeric_limits<NT>::epsilon()) const;
151 
152  // POST: returns the time in seconds taken by the constructor call for
153  // computing the smallest enclosing ball
154  double get_time() const;
155 
156  // POST: deletes dynamically allocated arrays
157  ~Miniball();
158 
159  private:
160  void mtf_mb (Sit n);
161  void mtf_move_to_front (Sit j);
162  void pivot_mb (Pit n);
163  void pivot_move_to_front (Pit j);
164  NT excess (Pit pit) const;
165  void pop ();
166  bool push (Pit pit);
167  NT suboptimality () const;
168  void create_arrays();
169  void delete_arrays();
170  };
171 
172  // Class Definition
173  // ================
174  template <typename CoordAccessor>
175  Miniball<CoordAccessor>::Miniball (int d_, Pit begin, Pit end,
176  CoordAccessor ca)
177  : d (d_),
178  points_begin (begin),
179  points_end (end),
180  coord_accessor (ca),
181  time (clock()),
182  nt0 (NT(0)),
183  L(),
184  support_end (L.begin()),
185  fsize(0),
186  ssize(0),
187  current_c (NULL),
188  current_sqr_r (NT(-1)),
189  c (NULL),
190  sqr_r (NULL),
191  q0 (NULL),
192  z (NULL),
193  f (NULL),
194  v (NULL),
195  a (NULL)
196  {
197  assert (points_begin != points_end);
198  create_arrays();
199 
200  // set initial center
201  for (int j=0; j<d; ++j) c[0][j] = nt0;
202  current_c = c[0];
203 
204  // compute miniball
205  pivot_mb (points_end);
206 
207  // update time
208  time = (clock() - time) / CLOCKS_PER_SEC;
209  }
210 
211  template <typename CoordAccessor>
212  Miniball<CoordAccessor>::~Miniball()
213  {
214  delete_arrays();
215  }
216 
217  template <typename CoordAccessor>
218  void Miniball<CoordAccessor>::create_arrays()
219  {
220  c = new NT*[d+1];
221  v = new NT*[d+1];
222  a = new NT*[d+1];
223  for (int i=0; i<d+1; ++i) {
224  c[i] = new NT[d];
225  v[i] = new NT[d];
226  a[i] = new NT[d];
227  }
228  sqr_r = new NT[d+1];
229  q0 = new NT[d];
230  z = new NT[d+1];
231  f = new NT[d+1];
232  }
233 
234  template <typename CoordAccessor>
235  void Miniball<CoordAccessor>::delete_arrays()
236  {
237  delete[] f;
238  delete[] z;
239  delete[] q0;
240  delete[] sqr_r;
241  for (int i=0; i<d+1; ++i) {
242  delete[] a[i];
243  delete[] v[i];
244  delete[] c[i];
245  }
246  delete[] a;
247  delete[] v;
248  delete[] c;
249  }
250 
251  template <typename CoordAccessor>
252  const typename Miniball<CoordAccessor>::NT*
253  Miniball<CoordAccessor>::center () const
254  {
255  return current_c;
256  }
257 
258  template <typename CoordAccessor>
259  typename Miniball<CoordAccessor>::NT
260  Miniball<CoordAccessor>::squared_radius () const
261  {
262  return current_sqr_r;
263  }
264 
265  template <typename CoordAccessor>
266  int Miniball<CoordAccessor>::nr_support_points () const
267  {
268  assert (ssize < d+2);
269  return ssize;
270  }
271 
272  template <typename CoordAccessor>
273  typename Miniball<CoordAccessor>::SupportPointIterator
274  Miniball<CoordAccessor>::support_points_begin () const
275  {
276  return L.begin();
277  }
278 
279  template <typename CoordAccessor>
280  typename Miniball<CoordAccessor>::SupportPointIterator
281  Miniball<CoordAccessor>::support_points_end () const
282  {
283  return support_end;
284  }
285 
286  template <typename CoordAccessor>
287  typename Miniball<CoordAccessor>::NT
288  Miniball<CoordAccessor>::relative_error (NT& subopt) const
289  {
290  NT e, max_e = nt0;
291  // compute maximum absolute excess of support points
292  for (SupportPointIterator it = support_points_begin();
293  it != support_points_end(); ++it) {
294  e = excess (*it);
295  if (e < nt0) e = -e;
296  if (e > max_e) {
297  max_e = e;
298  }
299  }
300  // compute maximum excess of any point
301  for (Pit i = points_begin; i != points_end; ++i)
302  if ((e = excess (i)) > max_e)
303  max_e = e;
304 
305  subopt = suboptimality();
306  assert (current_sqr_r > nt0 || max_e == nt0);
307  return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r);
308  }
309 
310  template <typename CoordAccessor>
311  bool Miniball<CoordAccessor>::is_valid (NT tol) const
312  {
313  NT suboptimality;
314  return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) );
315  }
316 
317  template <typename CoordAccessor>
318  double Miniball<CoordAccessor>::get_time() const
319  {
320  return time;
321  }
322 
323  template <typename CoordAccessor>
324  void Miniball<CoordAccessor>::mtf_mb (Sit n)
325  {
326  // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n)
327  // B: the set of forced points, defining the current ball
328  // S: the superset of support points computed by the algorithm
329  // --------------------------------------------------------------
330  // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999,
331  // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf
332 
333  // PRE: B = S
334  assert (fsize == ssize);
335 
336  support_end = L.begin();
337  if ((fsize) == d+1) return;
338 
339  // incremental construction
340  for (Sit i = L.begin(); i != n;)
341  {
342  // INV: (support_end - L.begin() == |S|-|B|)
343  assert (std::distance (L.begin(), support_end) == ssize - fsize);
344 
345  Sit j = i++;
346  if (excess(*j) > nt0)
347  if (push(*j)) { // B := B + p_i
348  mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i)
349  pop(); // B := B - p_i
350  mtf_move_to_front(j);
351  }
352  }
353  // POST: the range [L.begin(), support_end) stores the set S\B
354  }
355 
356  template <typename CoordAccessor>
357  void Miniball<CoordAccessor>::mtf_move_to_front (Sit j)
358  {
359  if (support_end == j)
360  support_end++;
361  L.splice (L.begin(), L, j);
362  }
363 
364  template <typename CoordAccessor>
365  void Miniball<CoordAccessor>::pivot_mb (Pit n)
366  {
367  // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n)
368  // --------------------------------------------------------------
369  // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999,
370  // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf
371  NT old_sqr_r;
372  const NT* c;
373  Pit pivot, k;
374  NT e, max_e, sqr_r;
375  Cit p;
376  do {
377  old_sqr_r = current_sqr_r;
378  sqr_r = current_sqr_r;
379 
380  pivot = points_begin;
381  max_e = nt0;
382  for (k = points_begin; k != n; ++k) {
383  p = coord_accessor(k);
384  e = -sqr_r;
385  c = current_c;
386  for (int j=0; j<d; ++j)
387  e += mb_sqr<NT>(*p++-*c++);
388  if (e > max_e) {
389  max_e = e;
390  pivot = k;
391  }
392  }
393 
394  if (max_e > nt0) {
395  // check if the pivot is already contained in the support set
396  if (std::find(L.begin(), support_end, pivot) == support_end) {
397  assert (fsize == 0);
398  if (push (pivot)) {
399  mtf_mb(support_end);
400  pop();
401  pivot_move_to_front(pivot);
402  }
403  }
404  }
405  } while (old_sqr_r < current_sqr_r);
406  }
407 
408  template <typename CoordAccessor>
409  void Miniball<CoordAccessor>::pivot_move_to_front (Pit j)
410  {
411  L.push_front(j);
412  if (std::distance(L.begin(), support_end) == d+2)
413  support_end--;
414  }
415 
416  template <typename CoordAccessor>
417  inline typename Miniball<CoordAccessor>::NT
418  Miniball<CoordAccessor>::excess (Pit pit) const
419  {
420  Cit p = coord_accessor(pit);
421  NT e = -current_sqr_r;
422  NT* c = current_c;
423  for (int k=0; k<d; ++k){
424  e += mb_sqr<NT>(*p++-*c++);
425  }
426  return e;
427  }
428 
429  template <typename CoordAccessor>
430  void Miniball<CoordAccessor>::pop ()
431  {
432  --fsize;
433  }
434 
435  template <typename CoordAccessor>
436  bool Miniball<CoordAccessor>::push (Pit pit)
437  {
438  int i, j;
439  NT eps = mb_sqr<NT>(std::numeric_limits<NT>::epsilon());
440 
441  Cit cit = coord_accessor(pit);
442  Cit p = cit;
443 
444  if (fsize==0) {
445  for (i=0; i<d; ++i)
446  q0[i] = *p++;
447  for (i=0; i<d; ++i)
448  c[0][i] = q0[i];
449  sqr_r[0] = nt0;
450  }
451  else {
452  // set v_fsize to Q_fsize
453  for (i=0; i<d; ++i)
454  //v[fsize][i] = p[i]-q0[i];
455  v[fsize][i] = *p++-q0[i];
456 
457  // compute the a_{fsize,i}, i< fsize
458  for (i=1; i<fsize; ++i) {
459  a[fsize][i] = nt0;
460  for (j=0; j<d; ++j)
461  a[fsize][i] += v[i][j] * v[fsize][j];
462  a[fsize][i]*=(2/z[i]);
463  }
464 
465  // update v_fsize to Q_fsize-\bar{Q}_fsize
466  for (i=1; i<fsize; ++i) {
467  for (j=0; j<d; ++j)
468  v[fsize][j] -= a[fsize][i]*v[i][j];
469  }
470 
471  // compute z_fsize
472  z[fsize]=nt0;
473  for (j=0; j<d; ++j)
474  z[fsize] += mb_sqr<NT>(v[fsize][j]);
475  z[fsize]*=2;
476 
477  // reject push if z_fsize too small
478  if (z[fsize]<eps*current_sqr_r) {
479  return false;
480  }
481 
482  // update c, sqr_r
483  p=cit;
484  NT e = -sqr_r[fsize-1];
485  for (i=0; i<d; ++i)
486  e += mb_sqr<NT>(*p++-c[fsize-1][i]);
487  f[fsize]=e/z[fsize];
488 
489  for (i=0; i<d; ++i)
490  c[fsize][i] = c[fsize-1][i]+f[fsize]*v[fsize][i];
491  sqr_r[fsize] = sqr_r[fsize-1] + e*f[fsize]/2;
492  }
493  current_c = c[fsize];
494  current_sqr_r = sqr_r[fsize];
495  ssize = ++fsize;
496  return true;
497  }
498 
499  template <typename CoordAccessor>
500  typename Miniball<CoordAccessor>::NT
501  Miniball<CoordAccessor>::suboptimality () const
502  {
503  NT* l = new NT[d+1];
504  NT min_l = nt0;
505  l[0] = NT(1);
506  for (int i=ssize-1; i>0; --i) {
507  l[i] = f[i];
508  for (int k=ssize-1; k>i; --k)
509  l[i]-=a[k][i]*l[k];
510  if (l[i] < min_l) min_l = l[i];
511  l[0] -= l[i];
512  }
513  if (l[0] < min_l) min_l = l[0];
514  delete[] l;
515  if (min_l < nt0)
516  return -min_l;
517  return nt0;
518  }
519 } // namespace Miniball
520 
521 } // namespace Gudhi
522 
523 #endif // MINIBALL_HPP_
Definition: SimplicialComplexForAlpha.h:26
GUDHI  Version 2.2.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : GPL v3 Generated on Tue Aug 14 2018 11:45:43 for GUDHI by Doxygen 1.8.13