$search
00001 package edu.tum.cs.vis.model.util.algorithm; 00002 00003 import java.util.ArrayList; 00004 00061 public class Miniball { 00062 00063 class Miniball_b { 00064 00065 int m, s; // size and number of support points 00066 double[] q0 = new double[d]; 00067 double[] z = new double[d + 1]; 00068 double[] f = new double[d + 1]; 00069 double[][] v = new double[d + 1][d]; 00070 double[][] a = new double[d + 1][d]; 00071 double[][] c = new double[d + 1][d]; 00072 double[] sqr_r = new double[d + 1]; 00073 double[] current_c = new double[d]; // refers to some c[j] 00074 double current_sqr_r; 00075 00076 double excess(double[] p) { 00077 double e = -current_sqr_r; 00078 for (int k = 0; k < d; ++k) { 00079 e += mb_sqr(p[k] - current_c[k]); 00080 } 00081 return e; 00082 } 00083 00084 double[] getCenter() { 00085 return (current_c); 00086 } 00087 00088 void pop() { 00089 --m; 00090 } 00091 00092 boolean push(double[] p) { 00093 // System.out.println("Miniball_b:push"); 00094 int i, j; 00095 double eps = 1e-32; 00096 00097 if (m == 0) { 00098 for (i = 0; i < d; ++i) { 00099 q0[i] = p[i]; 00100 } 00101 for (i = 0; i < d; ++i) { 00102 c[0][i] = q0[i]; 00103 } 00104 sqr_r[0] = 0; 00105 } else { 00106 // set v_m to Q_m 00107 for (i = 0; i < d; ++i) { 00108 v[m][i] = p[i] - q0[i]; 00109 } 00110 00111 // compute the a_{m,i}, i< m 00112 for (i = 1; i < m; ++i) { 00113 a[m][i] = 0; 00114 for (j = 0; j < d; ++j) { 00115 a[m][i] += v[i][j] * v[m][j]; 00116 } 00117 a[m][i] *= (2 / z[i]); 00118 } 00119 00120 // update v_m to Q_m-\bar{Q}_m 00121 for (i = 1; i < m; ++i) { 00122 for (j = 0; j < d; ++j) { 00123 v[m][j] -= a[m][i] * v[i][j]; 00124 } 00125 } 00126 00127 // compute z_m 00128 z[m] = 0; 00129 for (j = 0; j < d; ++j) { 00130 z[m] += mb_sqr(v[m][j]); 00131 } 00132 z[m] *= 2; 00133 00134 // reject push if z_m too small 00135 if (z[m] < eps * current_sqr_r) { 00136 return false; 00137 } 00138 00139 // update c, sqr_r 00140 double e = -sqr_r[m - 1]; 00141 for (i = 0; i < d; ++i) { 00142 e += mb_sqr(p[i] - c[m - 1][i]); 00143 } 00144 f[m] = e / z[m]; 00145 00146 for (i = 0; i < d; ++i) { 00147 c[m][i] = c[m - 1][i] + f[m] * v[m][i]; 00148 } 00149 sqr_r[m] = sqr_r[m - 1] + e * f[m] / 2; 00150 } 00151 current_c = c[m]; 00152 current_sqr_r = sqr_r[m]; 00153 s = ++m; 00154 return true; 00155 } 00156 00157 void reset() { 00158 m = 0; 00159 s = 0; 00160 // we misuse c[0] for the center of the empty sphere 00161 for (int j = 0; j < d; j++) { 00162 c[0][j] = 0; 00163 } 00164 current_c = c[0]; 00165 current_sqr_r = -1; 00166 } 00167 00168 int size() { 00169 return m; 00170 } 00171 00172 double slack() { 00173 double min_l = 0; 00174 double[] l = new double[d + 1]; 00175 l[0] = 1; 00176 for (int i = s - 1; i > 0; --i) { 00177 l[i] = f[i]; 00178 for (int k = s - 1; k > i; --k) { 00179 l[i] -= a[k][i] * l[k]; 00180 } 00181 if (l[i] < min_l) { 00182 min_l = l[i]; 00183 } 00184 l[0] -= l[i]; 00185 } 00186 if (l[0] < min_l) { 00187 min_l = l[0]; 00188 } 00189 return ((min_l < 0) ? -min_l : 0); 00190 } 00191 00192 double squared_radius() { 00193 return current_sqr_r; 00194 } 00195 00196 int support_size() { 00197 return s; 00198 } 00199 } 00200 00201 class pvt { 00202 int val; 00203 00204 pvt() { 00205 val = 0; 00206 } 00207 00208 int getVal() { 00209 return (val); 00210 } 00211 00212 void setVal(int i) { 00213 val = i; 00214 } 00215 } 00216 00217 static double mb_sqr(double r) { 00218 return r * r; 00219 } 00220 00221 static int points_begin() { 00222 return (0); 00223 } 00224 00225 static int support_points_begin() { 00226 return (0); 00227 } 00228 00229 int d; 00230 00231 @SuppressWarnings("rawtypes") 00232 ArrayList L; 00233 00234 Miniball_b B; 00235 00236 int support_end = 0; 00237 00238 @SuppressWarnings("rawtypes") 00239 public Miniball(int dim) { 00240 d = dim; 00241 L = new ArrayList(); 00242 B = new Miniball_b(); 00243 } 00244 00245 double accuracy(double slack) { 00246 double e, max_e = 0; 00247 int n_supp = 0; 00248 00249 int i; 00250 for (i = points_begin(); i != support_end; ++i, ++n_supp) { 00251 double[] sp = (double[]) L.get(i); 00252 if ((e = Math.abs(B.excess(sp))) > max_e) { 00253 max_e = e; 00254 } 00255 } 00256 if (n_supp == nr_support_points()) { 00257 System.out.println("Miniball.accuracy WARNING: STRANGE PROBLEM HERE!"); 00258 } 00259 for (i = support_end; i != points_end(); ++i) { 00260 double[] sp = (double[]) L.get(i); 00261 if ((e = B.excess(sp)) > max_e) { 00262 max_e = e; 00263 } 00264 } 00265 // slack = B.slack(); 00266 return (max_e / squared_radius()); 00267 } 00268 00273 public void build() { 00274 B.reset(); 00275 support_end = 0; 00276 pivot_mb(points_end()); 00277 } 00278 00284 public double[] center() { 00285 return B.getCenter(); 00286 } 00287 00295 @SuppressWarnings("unchecked") 00296 public void check_in(double[] p) { 00297 if (p != null) { 00298 L.add(p); 00299 } else { 00300 System.out.println("Miniball.check_in WARNING: Skipping null point"); 00301 } 00302 } 00303 00311 public void clear() { 00312 L.clear(); 00313 } 00314 00315 boolean is_valid(double tolerance) { 00316 double slack = 0.0; 00317 return ((accuracy(slack) < tolerance) && (slack == 0)); 00318 } 00319 00320 double max_excess(int t, int i, pvt pivot) { 00321 double[] c = B.getCenter(); 00322 double sqr_r = B.squared_radius(); 00323 double e, max_e = 0; 00324 for (int k = t; k != i; ++k) { 00325 double[] p = (double[]) L.get(k); 00326 00327 e = -sqr_r; 00328 for (int j = 0; j < d; ++j) { 00329 e += mb_sqr(p[j] - c[j]); 00330 } 00331 if (e > max_e) { 00332 max_e = e; 00333 pivot.setVal(k); 00334 } 00335 } 00336 return max_e; 00337 } 00338 00339 @SuppressWarnings("unchecked") 00340 void move_to_front(int j) { 00341 00342 if (support_end <= j) { 00343 support_end++; 00344 } 00345 // L.splice (L.begin(), L, j); 00346 double[] sp = (double[]) L.get(j); 00347 L.remove(j); 00348 L.add(0, sp); 00349 } 00350 00351 void mtf_mb(int i) { 00352 int pj = 0; 00353 support_end = points_begin(); 00354 if ((B.size()) == d + 1) { 00355 return; 00356 } 00357 for (int k = points_begin(); k != i;) { 00358 pj = pj + 1; 00359 int j = k++; 00360 double[] sp = (double[]) L.get(j); 00361 if (B.excess(sp) > 0) { 00362 if (B.push(sp)) { 00363 mtf_mb(j); 00364 B.pop(); 00365 move_to_front(j); 00366 } 00367 } 00368 } 00369 } 00370 00376 public int nr_points() { 00377 return L.size(); 00378 } 00379 00386 public int nr_support_points() { 00387 return B.support_size(); 00388 } 00389 00390 void pivot_mb(int i) { 00391 int t = 1; 00392 mtf_mb(t); 00393 double max_e = 0.0, old_sqr_r = -1; 00394 pvt pivot = new pvt(); 00395 do { 00396 max_e = max_excess(t, i, pivot); 00397 if (max_e > 0) { 00398 t = support_end; 00399 if (t == pivot.getVal()) { 00400 ++t; 00401 } 00402 old_sqr_r = B.squared_radius(); 00403 double[] sp = (double[]) L.get(pivot.getVal()); 00404 B.push(sp); 00405 mtf_mb(support_end); 00406 B.pop(); 00407 move_to_front(pivot.getVal()); 00408 } 00409 } while ((max_e > 0) && (B.squared_radius() > old_sqr_r)); 00410 } 00411 00412 int points_end() { 00413 return (L.size()); 00414 } 00415 00421 public double radius() { 00422 return ((1 + 0.00001) * Math.sqrt(B.squared_radius())); 00423 } 00424 00430 public double squared_radius() { 00431 return B.squared_radius(); 00432 } 00433 00434 int support_points_end() { 00435 return support_end; 00436 } 00437 }