36 #ifndef WFMATH_MINIBALL_FUNCS_H
37 #define WFMATH_MINIBALL_FUNCS_H
41 namespace WFMath {
namespace _miniball {
47 void Miniball<d>::check_in (
const Point& p)
54 void Miniball<d>::build (
bool pivoting)
57 support_end = L.begin();
66 void Miniball<d>::mtf_mb (It i)
68 support_end = L.begin();
69 if ((B.size())==d+1)
return;
70 for (It k=L.begin(); k!=i;) {
72 if (B.excess(*j) > 0) {
83 void Miniball<d>::move_to_front (It j)
87 L.splice (L.begin(), L, j);
92 void Miniball<d>::pivot_mb (It i)
96 double max_e, old_sqr_r;
99 max_e = max_excess (t, i, pivot);
104 old_sqr_r = B.squared_radius();
106 mtf_mb (support_end);
108 move_to_front (pivot);
109 }
while (B.squared_radius() > old_sqr_r);
114 double Miniball<d>::max_excess (It t, It i, It& pivot)
const
116 const double *c = B.center(), sqr_r = B.squared_radius();
118 for (It k=t; k!=i; ++k) {
119 const double *p = (*k).begin();
121 for (
int j=0; j<d; ++j)
134 typename Miniball<d>::Point Miniball<d>::center ()
const
136 return Point(B.center());
140 double Miniball<d>::squared_radius ()
const
142 return B.squared_radius();
147 int Miniball<d>::nr_points ()
const
153 typename Miniball<d>::Cit Miniball<d>::points_begin ()
const
159 typename Miniball<d>::Cit Miniball<d>::points_end ()
const
166 int Miniball<d>::nr_support_points ()
const
168 return B.support_size();
172 typename Miniball<d>::Cit Miniball<d>::support_points_begin ()
const
178 typename Miniball<d>::Cit Miniball<d>::support_points_end ()
const
186 double Miniball<d>::accuracy (
double& slack)
const
191 for (i=L.begin(); i!=support_end; ++i,++n_supp)
192 if ((e = abs (B.excess (*i))) > max_e)
196 assert (n_supp == nr_support_points());
198 for (i=support_end; i!=L.end(); ++i)
199 if ((e = B.excess (*i)) > max_e)
203 return (max_e/squared_radius());
208 bool Miniball<d>::is_valid (
double tolerance)
const
211 return ( (accuracy (slack) < tolerance) && (slack == 0) );
220 const double* Basis<d>::center ()
const
226 double Basis<d>::squared_radius()
const
228 return current_sqr_r;
232 int Basis<d>::size()
const
238 int Basis<d>::support_size()
const
244 double Basis<d>::excess (
const Point& p)
const
246 double e = -current_sqr_r;
247 for (
int k=0; k<d; ++k)
248 e += sqr(p[k]-current_c[k]);
255 void Basis<d>::reset ()
259 for (
int j=0; j<d; ++j)
267 Basis<d>::Basis () : m(0), s(0), current_c(0), current_sqr_r(-1.)
274 void Basis<d>::pop ()
281 bool Basis<d>::push (
const Point& p)
294 v[m][i] = p[i]-q0[i];
297 for (i=1; i<m; ++i) {
300 a[m][i] += v[i][j] * v[m][j];
305 for (i=1; i<m; ++i) {
307 v[m][j] -= a[m][i]*v[i][j];
313 z[m] += sqr(v[m][j]);
317 if (z[m]<eps*current_sqr_r) {
322 double e = -sqr_r[m-1];
324 e += sqr(p[i]-c[m-1][i]);
328 c[m][i] = c[m-1][i]+f[m]*v[m][i];
329 sqr_r[m] = sqr_r[m-1] + e*f[m]/2;
332 current_sqr_r = sqr_r[m];
340 double Basis<d>::slack ()
const
342 double l[d+1], min_l=0;
344 for (
int i=s-1; i>0; --i) {
346 for (
int k=s-1; k>i; --k)
348 if (l[i] < min_l) min_l = l[i];
351 if (l[0] < min_l) min_l = l[0];
352 return ( (min_l < 0) ? -min_l : 0);
Generic library namespace.