001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math3.geometry.partitioning.utilities; 018 019 import java.util.Arrays; 020 021 import org.apache.commons.math3.util.FastMath; 022 023 /** This class implements an ordering operation for T-uples. 024 * 025 * <p>Ordering is done by encoding all components of the T-uple into a 026 * single scalar value and using this value as the sorting 027 * key. Encoding is performed using the method invented by Georg 028 * Cantor in 1877 when he proved it was possible to establish a 029 * bijection between a line and a plane. The binary representations of 030 * the components of the T-uple are mixed together to form a single 031 * scalar. This means that the 2<sup>k</sup> bit of component 0 is 032 * followed by the 2<sup>k</sup> bit of component 1, then by the 033 * 2<sup>k</sup> bit of component 2 up to the 2<sup>k</sup> bit of 034 * component {@code t}, which is followed by the 2<sup>k-1</sup> 035 * bit of component 0, followed by the 2<sup>k-1</sup> bit of 036 * component 1 ... The binary representations are extended as needed 037 * to handle numbers with different scales and a suitable 038 * 2<sup>p</sup> offset is added to the components in order to avoid 039 * negative numbers (this offset is adjusted as needed during the 040 * comparison operations).</p> 041 * 042 * <p>The more interesting property of the encoding method for our 043 * purpose is that it allows to select all the points that are in a 044 * given range. This is depicted in dimension 2 by the following 045 * picture:</p> 046 * 047 * <img src="doc-files/OrderedTuple.png" /> 048 * 049 * <p>This picture shows a set of 100000 random 2-D pairs having their 050 * first component between -50 and +150 and their second component 051 * between -350 and +50. We wanted to extract all pairs having their 052 * first component between +30 and +70 and their second component 053 * between -120 and -30. We built the lower left point at coordinates 054 * (30, -120) and the upper right point at coordinates (70, -30). All 055 * points smaller than the lower left point are drawn in red and all 056 * points larger than the upper right point are drawn in blue. The 057 * green points are between the two limits. This picture shows that 058 * all the desired points are selected, along with spurious points. In 059 * this case, we get 15790 points, 4420 of which really belonging to 060 * the desired rectangle. It is possible to extract very small 061 * subsets. As an example extracting from the same 100000 points set 062 * the points having their first component between +30 and +31 and 063 * their second component between -91 and -90, we get a subset of 11 064 * points, 2 of which really belonging to the desired rectangle.</p> 065 * 066 * <p>the previous selection technique can be applied in all 067 * dimensions, still using two points to define the interval. The 068 * first point will have all its components set to their lower bounds 069 * while the second point will have all its components set to their 070 * upper bounds.</p> 071 * 072 * <p>T-uples with negative infinite or positive infinite components 073 * are sorted logically.</p> 074 * 075 * <p>Since the specification of the {@code Comparator} interface 076 * allows only {@code ClassCastException} errors, some arbitrary 077 * choices have been made to handle specific cases. The rationale for 078 * these choices is to keep <em>regular</em> and consistent T-uples 079 * together.</p> 080 * <ul> 081 * <li>instances with different dimensions are sorted according to 082 * their dimension regardless of their components values</li> 083 * <li>instances with {@code Double.NaN} components are sorted 084 * after all other ones (even after instances with positive infinite 085 * components</li> 086 * <li>instances with both positive and negative infinite components 087 * are considered as if they had {@code Double.NaN} 088 * components</li> 089 * </ul> 090 * 091 * @version $Id: OrderedTuple.java 1416643 2012-12-03 19:37:14Z tn $ 092 * @since 3.0 093 */ 094 public class OrderedTuple implements Comparable<OrderedTuple> { 095 096 /** Sign bit mask. */ 097 private static final long SIGN_MASK = 0x8000000000000000L; 098 099 /** Exponent bits mask. */ 100 private static final long EXPONENT_MASK = 0x7ff0000000000000L; 101 102 /** Mantissa bits mask. */ 103 private static final long MANTISSA_MASK = 0x000fffffffffffffL; 104 105 /** Implicit MSB for normalized numbers. */ 106 private static final long IMPLICIT_ONE = 0x0010000000000000L; 107 108 /** Double components of the T-uple. */ 109 private double[] components; 110 111 /** Offset scale. */ 112 private int offset; 113 114 /** Least Significant Bit scale. */ 115 private int lsb; 116 117 /** Ordering encoding of the double components. */ 118 private long[] encoding; 119 120 /** Positive infinity marker. */ 121 private boolean posInf; 122 123 /** Negative infinity marker. */ 124 private boolean negInf; 125 126 /** Not A Number marker. */ 127 private boolean nan; 128 129 /** Build an ordered T-uple from its components. 130 * @param components double components of the T-uple 131 */ 132 public OrderedTuple(final double ... components) { 133 this.components = components.clone(); 134 int msb = Integer.MIN_VALUE; 135 lsb = Integer.MAX_VALUE; 136 posInf = false; 137 negInf = false; 138 nan = false; 139 for (int i = 0; i < components.length; ++i) { 140 if (Double.isInfinite(components[i])) { 141 if (components[i] < 0) { 142 negInf = true; 143 } else { 144 posInf = true; 145 } 146 } else if (Double.isNaN(components[i])) { 147 nan = true; 148 } else { 149 final long b = Double.doubleToLongBits(components[i]); 150 final long m = mantissa(b); 151 if (m != 0) { 152 final int e = exponent(b); 153 msb = FastMath.max(msb, e + computeMSB(m)); 154 lsb = FastMath.min(lsb, e + computeLSB(m)); 155 } 156 } 157 } 158 159 if (posInf && negInf) { 160 // instance cannot be sorted logically 161 posInf = false; 162 negInf = false; 163 nan = true; 164 } 165 166 if (lsb <= msb) { 167 // encode the T-upple with the specified offset 168 encode(msb + 16); 169 } else { 170 encoding = new long[] { 171 0x0L 172 }; 173 } 174 175 } 176 177 /** Encode the T-uple with a given offset. 178 * @param minOffset minimal scale of the offset to add to all 179 * components (must be greater than the MSBs of all components) 180 */ 181 private void encode(final int minOffset) { 182 183 // choose an offset with some margins 184 offset = minOffset + 31; 185 offset -= offset % 32; 186 187 if ((encoding != null) && (encoding.length == 1) && (encoding[0] == 0x0L)) { 188 // the components are all zeroes 189 return; 190 } 191 192 // allocate an integer array to encode the components (we use only 193 // 63 bits per element because there is no unsigned long in Java) 194 final int neededBits = offset + 1 - lsb; 195 final int neededLongs = (neededBits + 62) / 63; 196 encoding = new long[components.length * neededLongs]; 197 198 // mix the bits from all components 199 int eIndex = 0; 200 int shift = 62; 201 long word = 0x0L; 202 for (int k = offset; eIndex < encoding.length; --k) { 203 for (int vIndex = 0; vIndex < components.length; ++vIndex) { 204 if (getBit(vIndex, k) != 0) { 205 word |= 0x1L << shift; 206 } 207 if (shift-- == 0) { 208 encoding[eIndex++] = word; 209 word = 0x0L; 210 shift = 62; 211 } 212 } 213 } 214 215 } 216 217 /** Compares this ordered T-uple with the specified object. 218 219 * <p>The ordering method is detailed in the general description of 220 * the class. Its main property is to be consistent with distance: 221 * geometrically close T-uples stay close to each other when stored 222 * in a sorted collection using this comparison method.</p> 223 224 * <p>T-uples with negative infinite, positive infinite are sorted 225 * logically.</p> 226 227 * <p>Some arbitrary choices have been made to handle specific 228 * cases. The rationale for these choices is to keep 229 * <em>normal</em> and consistent T-uples together.</p> 230 * <ul> 231 * <li>instances with different dimensions are sorted according to 232 * their dimension regardless of their components values</li> 233 * <li>instances with {@code Double.NaN} components are sorted 234 * after all other ones (evan after instances with positive infinite 235 * components</li> 236 * <li>instances with both positive and negative infinite components 237 * are considered as if they had {@code Double.NaN} 238 * components</li> 239 * </ul> 240 241 * @param ot T-uple to compare instance with 242 * @return a negative integer if the instance is less than the 243 * object, zero if they are equal, or a positive integer if the 244 * instance is greater than the object 245 246 */ 247 public int compareTo(final OrderedTuple ot) { 248 if (components.length == ot.components.length) { 249 if (nan) { 250 return +1; 251 } else if (ot.nan) { 252 return -1; 253 } else if (negInf || ot.posInf) { 254 return -1; 255 } else if (posInf || ot.negInf) { 256 return +1; 257 } else { 258 259 if (offset < ot.offset) { 260 encode(ot.offset); 261 } else if (offset > ot.offset) { 262 ot.encode(offset); 263 } 264 265 final int limit = FastMath.min(encoding.length, ot.encoding.length); 266 for (int i = 0; i < limit; ++i) { 267 if (encoding[i] < ot.encoding[i]) { 268 return -1; 269 } else if (encoding[i] > ot.encoding[i]) { 270 return +1; 271 } 272 } 273 274 if (encoding.length < ot.encoding.length) { 275 return -1; 276 } else if (encoding.length > ot.encoding.length) { 277 return +1; 278 } else { 279 return 0; 280 } 281 282 } 283 } 284 285 return components.length - ot.components.length; 286 287 } 288 289 /** {@inheritDoc} */ 290 @Override 291 public boolean equals(final Object other) { 292 if (this == other) { 293 return true; 294 } else if (other instanceof OrderedTuple) { 295 return compareTo((OrderedTuple) other) == 0; 296 } else { 297 return false; 298 } 299 } 300 301 /** {@inheritDoc} */ 302 @Override 303 public int hashCode() { 304 // the following constants are arbitrary small primes 305 final int multiplier = 37; 306 final int trueHash = 97; 307 final int falseHash = 71; 308 309 // hash fields and combine them 310 // (we rely on the multiplier to have different combined weights 311 // for all int fields and all boolean fields) 312 int hash = Arrays.hashCode(components); 313 hash = hash * multiplier + offset; 314 hash = hash * multiplier + lsb; 315 hash = hash * multiplier + (posInf ? trueHash : falseHash); 316 hash = hash * multiplier + (negInf ? trueHash : falseHash); 317 hash = hash * multiplier + (nan ? trueHash : falseHash); 318 319 return hash; 320 321 } 322 323 /** Get the components array. 324 * @return array containing the T-uple components 325 */ 326 public double[] getComponents() { 327 return components.clone(); 328 } 329 330 /** Extract the sign from the bits of a double. 331 * @param bits binary representation of the double 332 * @return sign bit (zero if positive, non zero if negative) 333 */ 334 private static long sign(final long bits) { 335 return bits & SIGN_MASK; 336 } 337 338 /** Extract the exponent from the bits of a double. 339 * @param bits binary representation of the double 340 * @return exponent 341 */ 342 private static int exponent(final long bits) { 343 return ((int) ((bits & EXPONENT_MASK) >> 52)) - 1075; 344 } 345 346 /** Extract the mantissa from the bits of a double. 347 * @param bits binary representation of the double 348 * @return mantissa 349 */ 350 private static long mantissa(final long bits) { 351 return ((bits & EXPONENT_MASK) == 0) ? 352 ((bits & MANTISSA_MASK) << 1) : // subnormal number 353 (IMPLICIT_ONE | (bits & MANTISSA_MASK)); // normal number 354 } 355 356 /** Compute the most significant bit of a long. 357 * @param l long from which the most significant bit is requested 358 * @return scale of the most significant bit of {@code l}, 359 * or 0 if {@code l} is zero 360 * @see #computeLSB 361 */ 362 private static int computeMSB(final long l) { 363 364 long ll = l; 365 long mask = 0xffffffffL; 366 int scale = 32; 367 int msb = 0; 368 369 while (scale != 0) { 370 if ((ll & mask) != ll) { 371 msb |= scale; 372 ll = ll >> scale; 373 } 374 scale = scale >> 1; 375 mask = mask >> scale; 376 } 377 378 return msb; 379 380 } 381 382 /** Compute the least significant bit of a long. 383 * @param l long from which the least significant bit is requested 384 * @return scale of the least significant bit of {@code l}, 385 * or 63 if {@code l} is zero 386 * @see #computeMSB 387 */ 388 private static int computeLSB(final long l) { 389 390 long ll = l; 391 long mask = 0xffffffff00000000L; 392 int scale = 32; 393 int lsb = 0; 394 395 while (scale != 0) { 396 if ((ll & mask) == ll) { 397 lsb |= scale; 398 ll = ll >> scale; 399 } 400 scale = scale >> 1; 401 mask = mask >> scale; 402 } 403 404 return lsb; 405 406 } 407 408 /** Get a bit from the mantissa of a double. 409 * @param i index of the component 410 * @param k scale of the requested bit 411 * @return the specified bit (either 0 or 1), after the offset has 412 * been added to the double 413 */ 414 private int getBit(final int i, final int k) { 415 final long bits = Double.doubleToLongBits(components[i]); 416 final int e = exponent(bits); 417 if ((k < e) || (k > offset)) { 418 return 0; 419 } else if (k == offset) { 420 return (sign(bits) == 0L) ? 1 : 0; 421 } else if (k > (e + 52)) { 422 return (sign(bits) == 0L) ? 0 : 1; 423 } else { 424 final long m = (sign(bits) == 0L) ? mantissa(bits) : -mantissa(bits); 425 return (int) ((m >> (k - e)) & 0x1L); 426 } 427 } 428 429 }