36 #ifndef WFMATH_MINIBALL_FUNCS_H
37 #define WFMATH_MINIBALL_FUNCS_H
39 namespace WFMath {
namespace _miniball {
45 void Miniball<d>::check_in (
const Point& p)
52 void Miniball<d>::build (
bool pivoting)
55 support_end = L.begin();
64 void Miniball<d>::mtf_mb (It i)
66 support_end = L.begin();
67 if ((B.size())==d+1)
return;
68 for (It k=L.begin(); k!=i;) {
70 if (B.excess(*j) > 0) {
81 void Miniball<d>::move_to_front (It j)
85 L.splice (L.begin(), L, j);
90 void Miniball<d>::pivot_mb (It i)
94 double max_e, old_sqr_r;
97 max_e = max_excess (t, i, pivot);
102 old_sqr_r = B.squared_radius();
104 mtf_mb (support_end);
106 move_to_front (pivot);
107 }
while (B.squared_radius() > old_sqr_r);
112 double Miniball<d>::max_excess (It t, It i, It& pivot)
const
114 const double *c = B.center(), sqr_r = B.squared_radius();
116 for (It k=t; k!=i; ++k) {
117 const double *p = (*k).begin();
119 for (
int j=0; j<d; ++j)
132 typename Miniball<d>::Point Miniball<d>::center ()
const
134 return Point(B.center());
138 double Miniball<d>::squared_radius ()
const
140 return B.squared_radius();
145 int Miniball<d>::nr_points ()
const
151 typename Miniball<d>::Cit Miniball<d>::points_begin ()
const
157 typename Miniball<d>::Cit Miniball<d>::points_end ()
const
164 int Miniball<d>::nr_support_points ()
const
166 return B.support_size();
170 typename Miniball<d>::Cit Miniball<d>::support_points_begin ()
const
176 typename Miniball<d>::Cit Miniball<d>::support_points_end ()
const
184 double Miniball<d>::accuracy (
double& slack)
const
189 for (i=L.begin(); i!=support_end; ++i,++n_supp)
190 if ((e = abs (B.excess (*i))) > max_e)
194 assert (n_supp == nr_support_points());
196 for (i=support_end; i!=L.end(); ++i)
197 if ((e = B.excess (*i)) > max_e)
201 return (max_e/squared_radius());
206 bool Miniball<d>::is_valid (
double tolerance)
const
209 return ( (accuracy (slack) < tolerance) && (slack == 0) );
218 const double* Basis<d>::center ()
const
224 double Basis<d>::squared_radius()
const
226 return current_sqr_r;
230 int Basis<d>::size()
const
236 int Basis<d>::support_size()
const
242 double Basis<d>::excess (
const Point& p)
const
244 double e = -current_sqr_r;
245 for (
int k=0; k<d; ++k)
246 e += sqr(p[k]-current_c[k]);
253 void Basis<d>::reset ()
257 for (
int j=0; j<d; ++j)
265 Basis<d>::Basis () : m(0), s(0), current_c(0), current_sqr_r(-1.)
272 void Basis<d>::pop ()
279 bool Basis<d>::push (
const Point& p)
292 v[m][i] = p[i]-q0[i];
295 for (i=1; i<m; ++i) {
298 a[m][i] += v[i][j] * v[m][j];
303 for (i=1; i<m; ++i) {
305 v[m][j] -= a[m][i]*v[i][j];
311 z[m] += sqr(v[m][j]);
315 if (z[m]<eps*current_sqr_r) {
320 double e = -sqr_r[m-1];
322 e += sqr(p[i]-c[m-1][i]);
326 c[m][i] = c[m-1][i]+f[m]*v[m][i];
327 sqr_r[m] = sqr_r[m-1] + e*f[m]/2;
330 current_sqr_r = sqr_r[m];
338 double Basis<d>::slack ()
const
340 double l[d+1], min_l=0;
342 for (
int i=s-1; i>0; --i) {
344 for (
int k=s-1; k>i; --k)
346 if (l[i] < min_l) min_l = l[i];
349 if (l[0] < min_l) min_l = l[0];
350 return ( (min_l < 0) ? -min_l : 0);
Generic library namespace.