    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    }