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    }