commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DiscreteUniformSampler.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.rng.sampling.distribution;

import org.apache.commons.rng.UniformRandomProvider;

/**
 * Discrete uniform distribution sampler.
 *
 * <p>Sampling uses {@link UniformRandomProvider#nextInt}.</p>
 *
 * <p>When the range is a power of two the number of calls is 1 per sample.
 * Otherwise a rejection algorithm is used to ensure uniformity. In the worst
 * case scenario where the range spans half the range of an {@code int}
 * (2<sup>31</sup> + 1) the expected number of calls is 2 per sample.</p>
 *
 * <p>This sampler can be used as a replacement for {@link UniformRandomProvider#nextInt}
 * with appropriate adjustment of the upper bound to be inclusive and will outperform that
 * method when the range is not a power of two. The advantage is gained by pre-computation
 * of the rejection threshold.</p>
 *
 * <p>The sampling algorithm is described in:</p>
 *
 * <blockquote>
 *  Lemire, D (2019). <i>Fast Random Integer Generation in an Interval.</i>
 *  <strong>ACM Transactions on Modeling and Computer Simulation</strong> 29 (1).
 * </blockquote>
 *
 * <p>The number of {@code int} values required per sample follows a geometric distribution with
 * a probability of success p of {@code 1 - ((2^32 % n) / 2^32)}. This requires on average 1/p random
 * {@code int} values per sample.</p>
 *
 * @see <a href="https://arxiv.org/abs/1805.10941">Fast Random Integer Generation in an Interval</a>
 *
 * @since 1.0
 */
public class DiscreteUniformSampler
    extends SamplerBase
    implements SharedStateDiscreteSampler {

    /** The appropriate uniform sampler for the parameters. */
    private final SharedStateDiscreteSampler delegate;

    /**
     * Base class for a sampler from a discrete uniform distribution. This contains the
     * source of randomness.
     */
    private abstract static class AbstractDiscreteUniformSampler
            implements SharedStateDiscreteSampler {
        /** Underlying source of randomness. */
        protected final UniformRandomProvider rng;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         */
        AbstractDiscreteUniformSampler(UniformRandomProvider rng) {
            this.rng = rng;
        }

        @Override
        public String toString() {
            return "Uniform deviate [" + rng.toString() + "]";
        }
    }

    /**
     * Discrete uniform distribution sampler when the sample value is fixed.
     */
    private static class FixedDiscreteUniformSampler
            extends AbstractDiscreteUniformSampler {
        /** The value. */
        private final int value;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param value The value.
         */
        FixedDiscreteUniformSampler(UniformRandomProvider rng,
                                    int value) {
            super(rng);
            this.value = value;
        }

        @Override
        public int sample() {
            return value;
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            // No requirement for the RNG
            return this;
        }
    }

    /**
     * Discrete uniform distribution sampler when the range is a power of 2 and greater than 1.
     * This sampler assumes the lower bound of the range is 0.
     *
     * <p>Note: This cannot be used when the range is 1 (2^0) as the shift would be 32-bits
     * which is ignored by the shift operator.</p>
     */
    private static class PowerOf2RangeDiscreteUniformSampler
            extends AbstractDiscreteUniformSampler {
        /** Bit shift to apply to the integer sample. */
        private final int shift;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param range Maximum range of the sample (exclusive).
         * Must be a power of 2 greater than 2^0.
         */
        PowerOf2RangeDiscreteUniformSampler(UniformRandomProvider rng,
                                            int range) {
            super(rng);
            this.shift = Integer.numberOfLeadingZeros(range) + 1;
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        PowerOf2RangeDiscreteUniformSampler(UniformRandomProvider rng,
                                            PowerOf2RangeDiscreteUniformSampler source) {
            super(rng);
            this.shift = source.shift;
        }

        @Override
        public int sample() {
            // Use a bit shift to favour the most significant bits.
            // Note: The result would be the same as the rejection method used in the
            // small range sampler when there is no rejection threshold.
            return rng.nextInt() >>> shift;
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new PowerOf2RangeDiscreteUniformSampler(rng, this);
        }
    }

    /**
     * Discrete uniform distribution sampler when the range is small
     * enough to fit in a positive integer.
     * This sampler assumes the lower bound of the range is 0.
     *
     * <p>Implements the algorithm of Lemire (2019).</p>
     *
     * @see <a href="https://arxiv.org/abs/1805.10941">Fast Random Integer Generation in an Interval</a>
     */
    private static class SmallRangeDiscreteUniformSampler
            extends AbstractDiscreteUniformSampler {
        /** Maximum range of the sample (exclusive). */
        private final long n;

        /**
         * The level below which samples are rejected based on the fraction remainder.
         *
         * <p>Any remainder below this denotes that there are still floor(2^32 / n) more
         * observations of this sample from the interval [0, 2^32), where n is the range.</p>
         */
        private final long threshold;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param range Maximum range of the sample (exclusive).
         */
        SmallRangeDiscreteUniformSampler(UniformRandomProvider rng,
                                         int range) {
            super(rng);
            // Handle range as an unsigned 32-bit integer
            this.n = range & 0xffffffffL;
            // Compute 2^32 % n
            threshold = (1L << 32) % n;
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        SmallRangeDiscreteUniformSampler(UniformRandomProvider rng,
                                         SmallRangeDiscreteUniformSampler source) {
            super(rng);
            this.n = source.n;
            this.threshold = source.threshold;
        }

        @Override
        public int sample() {
            // Rejection method using multiply by a fraction:
            // n * [0, 2^32 - 1)
            //     -------------
            //         2^32
            // The result is mapped back to an integer and will be in the range [0, n).
            // Note this is comparable to range * rng.nextDouble() but with compensation for
            // non-uniformity due to round-off.
            long result;
            do {
                // Compute 64-bit unsigned product of n * [0, 2^32 - 1).
                // The upper 32-bits contains the sample value in the range [0, n), i.e. result / 2^32.
                // The lower 32-bits contains the remainder, i.e. result % 2^32.
                result = n * (rng.nextInt() & 0xffffffffL);

                // Test the sample uniformity.
                // Samples are observed on average (2^32 / n) times at a frequency of either
                // floor(2^32 / n) or ceil(2^32 / n).
                // To ensure all samples have a frequency of floor(2^32 / n) reject any results with
                // a remainder < (2^32 % n), i.e. the level below which denotes that there are still
                // floor(2^32 / n) more observations of this sample.
            } while ((result & 0xffffffffL) < threshold);
            // Divide by 2^32 to get the sample.
            return (int)(result >>> 32);
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new SmallRangeDiscreteUniformSampler(rng, this);
        }
    }

    /**
     * Discrete uniform distribution sampler when the range between lower and upper is too large
     * to fit in a positive integer.
     */
    private static class LargeRangeDiscreteUniformSampler
            extends AbstractDiscreteUniformSampler {
        /** Lower bound. */
        private final int lower;
        /** Upper bound. */
        private final int upper;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param lower Lower bound (inclusive) of the distribution.
         * @param upper Upper bound (inclusive) of the distribution.
         */
        LargeRangeDiscreteUniformSampler(UniformRandomProvider rng,
                                         int lower,
                                         int upper) {
            super(rng);
            this.lower = lower;
            this.upper = upper;
        }

        @Override
        public int sample() {
            // Use a simple rejection method.
            // This is used when (upper-lower) >= Integer.MAX_VALUE.
            // This will loop on average 2 times in the worst case scenario
            // when (upper-lower) == Integer.MAX_VALUE.
            while (true) {
                final int r = rng.nextInt();
                if (r >= lower &&
                    r <= upper) {
                    return r;
                }
            }
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new LargeRangeDiscreteUniformSampler(rng, lower, upper);
        }
    }

    /**
     * Adds an offset to an underlying discrete sampler.
     */
    private static class OffsetDiscreteUniformSampler
            extends AbstractDiscreteUniformSampler {
        /** The offset. */
        private final int offset;
        /** The discrete sampler. */
        private final SharedStateDiscreteSampler sampler;

        /**
         * @param offset The offset for the sample.
         * @param sampler The discrete sampler.
         */
        OffsetDiscreteUniformSampler(int offset,
                                     SharedStateDiscreteSampler sampler) {
            super(null);
            this.offset = offset;
            this.sampler = sampler;
        }

        @Override
        public int sample() {
            return offset + sampler.sample();
        }

        /** {@inheritDoc} */
        @Override
        public String toString() {
            return sampler.toString();
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new OffsetDiscreteUniformSampler(offset, sampler.withUniformRandomProvider(rng));
        }
    }

    /**
     * This instance delegates sampling. Use the factory method
     * {@link #of(UniformRandomProvider, int, int)} to create an optimal sampler.
     *
     * @param rng Generator of uniformly distributed random numbers.
     * @param lower Lower bound (inclusive) of the distribution.
     * @param upper Upper bound (inclusive) of the distribution.
     * @throws IllegalArgumentException if {@code lower > upper}.
     */
    public DiscreteUniformSampler(UniformRandomProvider rng,
                                  int lower,
                                  int upper) {
        super(null);
        delegate = of(rng, lower, upper);
    }

    /** {@inheritDoc} */
    @Override
    public int sample() {
        return delegate.sample();
    }

    /** {@inheritDoc} */
    @Override
    public String toString() {
        return delegate.toString();
    }

    /**
     * {@inheritDoc}
     *
     * @since 1.3
     */
    @Override
    public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
        // Direct return of the optimised sampler
        return delegate.withUniformRandomProvider(rng);
    }

    /**
     * Creates a new discrete uniform distribution sampler.
     *
     * @param rng Generator of uniformly distributed random numbers.
     * @param lower Lower bound (inclusive) of the distribution.
     * @param upper Upper bound (inclusive) of the distribution.
     * @return the sampler
     * @throws IllegalArgumentException if {@code lower > upper}.
     * @since 1.3
     */
    public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
                                                int lower,
                                                int upper) {
        if (lower > upper) {
            throw new IllegalArgumentException(lower  + " > " + upper);
        }

        // Choose the algorithm depending on the range

        // Edge case for no range.
        // This must be done first as the methods to handle lower == 0
        // do not handle upper == 0.
        if (upper == lower) {
            return new FixedDiscreteUniformSampler(rng, lower);
        }

        // Algorithms to ignore the lower bound if it is zero.
        if (lower == 0) {
            return createZeroBoundedSampler(rng, upper);
        }

        final int range = (upper - lower) + 1;
        // Check power of 2 first to handle range == 2^31.
        if (isPowerOf2(range)) {
            return new OffsetDiscreteUniformSampler(lower,
                                                    new PowerOf2RangeDiscreteUniformSampler(rng, range));
        }
        if (range <= 0) {
            // The range is too wide to fit in a positive int (larger
            // than 2^31); use a simple rejection method.
            // Note: if range == 0 then the input is [Integer.MIN_VALUE, Integer.MAX_VALUE].
            // No specialisation exists for this and it is handled as a large range.
            return new LargeRangeDiscreteUniformSampler(rng, lower, upper);
        }
        // Use a sample from the range added to the lower bound.
        return new OffsetDiscreteUniformSampler(lower,
                                                new SmallRangeDiscreteUniformSampler(rng, range));
    }

    /**
     * Create a new sampler for the range {@code 0} inclusive to {@code upper} inclusive.
     *
     * <p>This can handle any positive {@code upper}.
     *
     * @param rng Generator of uniformly distributed random numbers.
     * @param upper Upper bound (inclusive) of the distribution. Must be positive.
     * @return the sampler
     */
    private static AbstractDiscreteUniformSampler createZeroBoundedSampler(UniformRandomProvider rng,
                                                                           int upper) {
        // Note: Handle any range up to 2^31 (which is negative as a signed
        // 32-bit integer but handled as a power of 2)
        final int range = upper + 1;
        return isPowerOf2(range) ?
            new PowerOf2RangeDiscreteUniformSampler(rng, range) :
            new SmallRangeDiscreteUniformSampler(rng, range);
    }

    /**
     * Checks if the value is a power of 2.
     *
     * <p>This returns {@code true} for the value {@code Integer.MIN_VALUE} which can be
     * handled as an unsigned integer of 2^31.</p>
     *
     * @param value Value.
     * @return {@code true} if a power of 2
     */
    private static boolean isPowerOf2(final int value) {
        return value != 0 && (value & (value - 1)) == 0;
    }
}