41 template <
typename NT>
42 inline NT mb_sqr (NT r) {
return r*r;}
49 template <
typename Pit_,
typename Cit_ >
50 struct CoordAccessor {
53 inline Cit operator() (Pit it)
const {
return (*it).begin(); }
57 template <
typename Pit_,
typename Cit_ >
58 struct CoordAccessor<Pit_, Cit_*> {
61 inline Cit operator() (Pit it)
const {
return *it; }
67 template <
typename CoordAccessor>
72 typedef typename CoordAccessor::Pit Pit;
74 typedef typename CoordAccessor::Cit Cit;
76 typedef typename std::iterator_traits<Cit>::value_type NT;
78 typedef typename std::list<Pit>::iterator Sit;
84 CoordAccessor coord_accessor;
109 typedef typename std::list<Pit>::const_iterator SupportPointIterator;
115 Miniball (
int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor());
119 const NT* center ()
const;
122 NT squared_radius ()
const;
129 int nr_support_points ()
const;
132 SupportPointIterator support_points_begin ()
const;
135 SupportPointIterator support_points_end ()
const;
145 NT relative_error (NT& subopt)
const;
150 bool is_valid (NT tol = NT(10) * std::numeric_limits<NT>::epsilon())
const;
154 double get_time()
const;
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;
167 NT suboptimality ()
const;
168 void create_arrays();
169 void delete_arrays();
174 template <
typename CoordAccessor>
175 Miniball<CoordAccessor>::Miniball (
int d_, Pit begin, Pit end,
178 points_begin (begin),
184 support_end (L.begin()),
188 current_sqr_r (NT(-1)),
197 assert (points_begin != points_end);
201 for (
int j=0; j<d; ++j) c[0][j] = nt0;
205 pivot_mb (points_end);
208 time = (clock() - time) / CLOCKS_PER_SEC;
211 template <
typename CoordAccessor>
212 Miniball<CoordAccessor>::~Miniball()
217 template <
typename CoordAccessor>
218 void Miniball<CoordAccessor>::create_arrays()
223 for (
int i=0; i<d+1; ++i) {
234 template <
typename CoordAccessor>
235 void Miniball<CoordAccessor>::delete_arrays()
241 for (
int i=0; i<d+1; ++i) {
251 template <
typename CoordAccessor>
252 const typename Miniball<CoordAccessor>::NT*
253 Miniball<CoordAccessor>::center ()
const 258 template <
typename CoordAccessor>
259 typename Miniball<CoordAccessor>::NT
260 Miniball<CoordAccessor>::squared_radius ()
const 262 return current_sqr_r;
265 template <
typename CoordAccessor>
266 int Miniball<CoordAccessor>::nr_support_points ()
const 268 assert (ssize < d+2);
272 template <
typename CoordAccessor>
273 typename Miniball<CoordAccessor>::SupportPointIterator
274 Miniball<CoordAccessor>::support_points_begin ()
const 279 template <
typename CoordAccessor>
280 typename Miniball<CoordAccessor>::SupportPointIterator
281 Miniball<CoordAccessor>::support_points_end ()
const 286 template <
typename CoordAccessor>
287 typename Miniball<CoordAccessor>::NT
288 Miniball<CoordAccessor>::relative_error (NT& subopt)
const 292 for (SupportPointIterator it = support_points_begin();
293 it != support_points_end(); ++it) {
301 for (Pit i = points_begin; i != points_end; ++i)
302 if ((e = excess (i)) > max_e)
305 subopt = suboptimality();
306 assert (current_sqr_r > nt0 || max_e == nt0);
307 return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r);
310 template <
typename CoordAccessor>
311 bool Miniball<CoordAccessor>::is_valid (NT tol)
const 314 return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) );
317 template <
typename CoordAccessor>
318 double Miniball<CoordAccessor>::get_time()
const 323 template <
typename CoordAccessor>
324 void Miniball<CoordAccessor>::mtf_mb (Sit n)
334 assert (fsize == ssize);
336 support_end = L.begin();
337 if ((fsize) == d+1)
return;
340 for (Sit i = L.begin(); i != n;)
343 assert (std::distance (L.begin(), support_end) == ssize - fsize);
346 if (excess(*j) > nt0)
350 mtf_move_to_front(j);
356 template <
typename CoordAccessor>
357 void Miniball<CoordAccessor>::mtf_move_to_front (Sit j)
359 if (support_end == j)
361 L.splice (L.begin(), L, j);
364 template <
typename CoordAccessor>
365 void Miniball<CoordAccessor>::pivot_mb (Pit n)
377 old_sqr_r = current_sqr_r;
378 sqr_r = current_sqr_r;
380 pivot = points_begin;
382 for (k = points_begin; k != n; ++k) {
383 p = coord_accessor(k);
386 for (
int j=0; j<d; ++j)
387 e += mb_sqr<NT>(*p++-*c++);
396 if (std::find(L.begin(), support_end, pivot) == support_end) {
401 pivot_move_to_front(pivot);
405 }
while (old_sqr_r < current_sqr_r);
408 template <
typename CoordAccessor>
409 void Miniball<CoordAccessor>::pivot_move_to_front (Pit j)
412 if (std::distance(L.begin(), support_end) == d+2)
416 template <
typename CoordAccessor>
417 inline typename Miniball<CoordAccessor>::NT
418 Miniball<CoordAccessor>::excess (Pit pit)
const 420 Cit p = coord_accessor(pit);
421 NT e = -current_sqr_r;
423 for (
int k=0; k<d; ++k){
424 e += mb_sqr<NT>(*p++-*c++);
429 template <
typename CoordAccessor>
430 void Miniball<CoordAccessor>::pop ()
435 template <
typename CoordAccessor>
436 bool Miniball<CoordAccessor>::push (Pit pit)
439 NT eps = mb_sqr<NT>(std::numeric_limits<NT>::epsilon());
441 Cit cit = coord_accessor(pit);
455 v[fsize][i] = *p++-q0[i];
458 for (i=1; i<fsize; ++i) {
461 a[fsize][i] += v[i][j] * v[fsize][j];
462 a[fsize][i]*=(2/z[i]);
466 for (i=1; i<fsize; ++i) {
468 v[fsize][j] -= a[fsize][i]*v[i][j];
474 z[fsize] += mb_sqr<NT>(v[fsize][j]);
478 if (z[fsize]<eps*current_sqr_r) {
484 NT e = -sqr_r[fsize-1];
486 e += mb_sqr<NT>(*p++-c[fsize-1][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;
493 current_c = c[fsize];
494 current_sqr_r = sqr_r[fsize];
499 template <
typename CoordAccessor>
500 typename Miniball<CoordAccessor>::NT
501 Miniball<CoordAccessor>::suboptimality ()
const 506 for (
int i=ssize-1; i>0; --i) {
508 for (
int k=ssize-1; k>i; --k)
510 if (l[i] < min_l) min_l = l[i];
513 if (l[0] < min_l) min_l = l[0];
523 #endif // MINIBALL_HPP_ Definition: SimplicialComplexForAlpha.h:14