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 }