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
35namespace Gudhi {
36
37namespace 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_
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Tue Aug 16 2022 14:01:50 for GUDHIdev by Doxygen 1.9.4