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 */
017package org.apache.commons.math3.geometry.partitioning.utilities;
018
019import java.util.Arrays;
020
021import 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 * @since 3.0
092 * @deprecated as of 3.4, this class is not used anymore and considered
093 * to be out of scope of Apache Commons Math
094 */
095@Deprecated
096public class OrderedTuple implements Comparable<OrderedTuple> {
097
098    /** Sign bit mask. */
099    private static final long SIGN_MASK     = 0x8000000000000000L;
100
101    /** Exponent bits mask. */
102    private static final long EXPONENT_MASK = 0x7ff0000000000000L;
103
104    /** Mantissa bits mask. */
105    private static final long MANTISSA_MASK = 0x000fffffffffffffL;
106
107    /** Implicit MSB for normalized numbers. */
108    private static final long IMPLICIT_ONE  = 0x0010000000000000L;
109
110    /** Double components of the T-uple. */
111    private double[] components;
112
113    /** Offset scale. */
114    private int offset;
115
116    /** Least Significant Bit scale. */
117    private int lsb;
118
119    /** Ordering encoding of the double components. */
120    private long[] encoding;
121
122    /** Positive infinity marker. */
123    private boolean posInf;
124
125    /** Negative infinity marker. */
126    private boolean negInf;
127
128    /** Not A Number marker. */
129    private boolean nan;
130
131    /** Build an ordered T-uple from its components.
132     * @param components double components of the T-uple
133     */
134    public OrderedTuple(final double ... components) {
135        this.components = components.clone();
136        int msb = Integer.MIN_VALUE;
137        lsb     = Integer.MAX_VALUE;
138        posInf  = false;
139        negInf  = false;
140        nan     = false;
141        for (int i = 0; i < components.length; ++i) {
142            if (Double.isInfinite(components[i])) {
143                if (components[i] < 0) {
144                    negInf = true;
145                } else {
146                    posInf = true;
147                }
148            } else if (Double.isNaN(components[i])) {
149                nan = true;
150            } else {
151                final long b = Double.doubleToLongBits(components[i]);
152                final long m = mantissa(b);
153                if (m != 0) {
154                    final int e = exponent(b);
155                    msb = FastMath.max(msb, e + computeMSB(m));
156                    lsb = FastMath.min(lsb, e + computeLSB(m));
157                }
158            }
159        }
160
161        if (posInf && negInf) {
162            // instance cannot be sorted logically
163            posInf = false;
164            negInf = false;
165            nan    = true;
166        }
167
168        if (lsb <= msb) {
169            // encode the T-upple with the specified offset
170            encode(msb + 16);
171        } else {
172            encoding = new long[] {
173                0x0L
174            };
175        }
176
177    }
178
179    /** Encode the T-uple with a given offset.
180     * @param minOffset minimal scale of the offset to add to all
181     * components (must be greater than the MSBs of all components)
182     */
183    private void encode(final int minOffset) {
184
185        // choose an offset with some margins
186        offset  = minOffset + 31;
187        offset -= offset % 32;
188
189        if ((encoding != null) && (encoding.length == 1) && (encoding[0] == 0x0L)) {
190            // the components are all zeroes
191            return;
192        }
193
194        // allocate an integer array to encode the components (we use only
195        // 63 bits per element because there is no unsigned long in Java)
196        final int neededBits  = offset + 1 - lsb;
197        final int neededLongs = (neededBits + 62) / 63;
198        encoding = new long[components.length * neededLongs];
199
200        // mix the bits from all components
201        int  eIndex = 0;
202        int  shift  = 62;
203        long word   = 0x0L;
204        for (int k = offset; eIndex < encoding.length; --k) {
205            for (int vIndex = 0; vIndex < components.length; ++vIndex) {
206                if (getBit(vIndex, k) != 0) {
207                    word |= 0x1L << shift;
208                }
209                if (shift-- == 0) {
210                    encoding[eIndex++] = word;
211                    word  = 0x0L;
212                    shift = 62;
213                }
214            }
215        }
216
217    }
218
219    /** Compares this ordered T-uple with the specified object.
220
221     * <p>The ordering method is detailed in the general description of
222     * the class. Its main property is to be consistent with distance:
223     * geometrically close T-uples stay close to each other when stored
224     * in a sorted collection using this comparison method.</p>
225
226     * <p>T-uples with negative infinite, positive infinite are sorted
227     * logically.</p>
228
229     * <p>Some arbitrary choices have been made to handle specific
230     * cases. The rationale for these choices is to keep
231     * <em>normal</em> and consistent T-uples together.</p>
232     * <ul>
233     * <li>instances with different dimensions are sorted according to
234     * their dimension regardless of their components values</li>
235     * <li>instances with {@code Double.NaN} components are sorted
236     * after all other ones (evan after instances with positive infinite
237     * components</li>
238     * <li>instances with both positive and negative infinite components
239     * are considered as if they had {@code Double.NaN}
240     * components</li>
241     * </ul>
242
243     * @param ot T-uple to compare instance with
244     * @return a negative integer if the instance is less than the
245     * object, zero if they are equal, or a positive integer if the
246     * instance is greater than the object
247
248     */
249    public int compareTo(final OrderedTuple ot) {
250        if (components.length == ot.components.length) {
251            if (nan) {
252                return +1;
253            } else if (ot.nan) {
254                return -1;
255            } else if (negInf || ot.posInf) {
256                return -1;
257            } else if (posInf || ot.negInf) {
258                return +1;
259            } else {
260
261                if (offset < ot.offset) {
262                    encode(ot.offset);
263                } else if (offset > ot.offset) {
264                    ot.encode(offset);
265                }
266
267                final int limit = FastMath.min(encoding.length, ot.encoding.length);
268                for (int i = 0; i < limit; ++i) {
269                    if (encoding[i] < ot.encoding[i]) {
270                        return -1;
271                    } else if (encoding[i] > ot.encoding[i]) {
272                        return +1;
273                    }
274                }
275
276                if (encoding.length < ot.encoding.length) {
277                    return -1;
278                } else if (encoding.length > ot.encoding.length) {
279                    return +1;
280                } else {
281                    return 0;
282                }
283
284            }
285        }
286
287        return components.length - ot.components.length;
288
289    }
290
291    /** {@inheritDoc} */
292    @Override
293    public boolean equals(final Object other) {
294        if (this == other) {
295            return true;
296        } else if (other instanceof OrderedTuple) {
297            return compareTo((OrderedTuple) other) == 0;
298        } else {
299            return false;
300        }
301    }
302
303    /** {@inheritDoc} */
304    @Override
305    public int hashCode() {
306        // the following constants are arbitrary small primes
307        final int multiplier = 37;
308        final int trueHash   = 97;
309        final int falseHash  = 71;
310
311        // hash fields and combine them
312        // (we rely on the multiplier to have different combined weights
313        //  for all int fields and all boolean fields)
314        int hash = Arrays.hashCode(components);
315        hash = hash * multiplier + offset;
316        hash = hash * multiplier + lsb;
317        hash = hash * multiplier + (posInf ? trueHash : falseHash);
318        hash = hash * multiplier + (negInf ? trueHash : falseHash);
319        hash = hash * multiplier + (nan    ? trueHash : falseHash);
320
321        return hash;
322
323    }
324
325    /** Get the components array.
326     * @return array containing the T-uple components
327     */
328    public double[] getComponents() {
329        return components.clone();
330    }
331
332    /** Extract the sign from the bits of a double.
333     * @param bits binary representation of the double
334     * @return sign bit (zero if positive, non zero if negative)
335     */
336    private static long sign(final long bits) {
337        return bits & SIGN_MASK;
338    }
339
340    /** Extract the exponent from the bits of a double.
341     * @param bits binary representation of the double
342     * @return exponent
343     */
344    private static int exponent(final long bits) {
345        return ((int) ((bits & EXPONENT_MASK) >> 52)) - 1075;
346    }
347
348    /** Extract the mantissa from the bits of a double.
349     * @param bits binary representation of the double
350     * @return mantissa
351     */
352    private static long mantissa(final long bits) {
353        return ((bits & EXPONENT_MASK) == 0) ?
354               ((bits & MANTISSA_MASK) << 1) :          // subnormal number
355               (IMPLICIT_ONE | (bits & MANTISSA_MASK)); // normal number
356    }
357
358    /** Compute the most significant bit of a long.
359     * @param l long from which the most significant bit is requested
360     * @return scale of the most significant bit of {@code l},
361     * or 0 if {@code l} is zero
362     * @see #computeLSB
363     */
364    private static int computeMSB(final long l) {
365
366        long ll = l;
367        long mask  = 0xffffffffL;
368        int  scale = 32;
369        int  msb   = 0;
370
371        while (scale != 0) {
372            if ((ll & mask) != ll) {
373                msb |= scale;
374                ll >>= scale;
375            }
376            scale >>= 1;
377            mask >>= scale;
378        }
379
380        return msb;
381
382    }
383
384    /** Compute the least significant bit of a long.
385     * @param l long from which the least significant bit is requested
386     * @return scale of the least significant bit of {@code l},
387     * or 63 if {@code l} is zero
388     * @see #computeMSB
389     */
390    private static int computeLSB(final long l) {
391
392        long ll = l;
393        long mask  = 0xffffffff00000000L;
394        int  scale = 32;
395        int  lsb   = 0;
396
397        while (scale != 0) {
398            if ((ll & mask) == ll) {
399                lsb |= scale;
400                ll >>= scale;
401            }
402            scale >>= 1;
403            mask >>= scale;
404        }
405
406        return lsb;
407
408    }
409
410    /** Get a bit from the mantissa of a double.
411     * @param i index of the component
412     * @param k scale of the requested bit
413     * @return the specified bit (either 0 or 1), after the offset has
414     * been added to the double
415     */
416    private int getBit(final int i, final int k) {
417        final long bits = Double.doubleToLongBits(components[i]);
418        final int e = exponent(bits);
419        if ((k < e) || (k > offset)) {
420            return 0;
421        } else if (k == offset) {
422            return (sign(bits) == 0L) ? 1 : 0;
423        } else if (k > (e + 52)) {
424            return (sign(bits) == 0L) ? 0 : 1;
425        } else {
426            final long m = (sign(bits) == 0L) ? mantissa(bits) : -mantissa(bits);
427            return (int) ((m >> (k - e)) & 0x1L);
428        }
429    }
430
431}