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.rng.sampling.distribution;
018
019import org.apache.commons.rng.UniformRandomProvider;
020
021/**
022 * Sampling from a uniform distribution.
023 *
024 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
025 *
026 * @since 1.0
027 */
028public class ContinuousUniformSampler
029    extends SamplerBase
030    implements SharedStateContinuousSampler {
031    /** The minimum ULP gap for the open interval when the doubles have the same sign. */
032    private static final int MIN_ULP_SAME_SIGN = 2;
033    /** The minimum ULP gap for the open interval when the doubles have the opposite sign. */
034    private static final int MIN_ULP_OPPOSITE_SIGN = 3;
035
036    /** Underlying source of randomness. */
037    protected final UniformRandomProvider rng;
038    /** Lower bound. */
039    private final double lo;
040    /** Higher bound. */
041    private final double hi;
042
043    /**
044     * Specialization to sample from an open interval {@code (lo, hi)} (see RNG-145).
045     *
046     * @since 1.4
047     */
048    private static final class OpenIntervalContinuousUniformSampler extends ContinuousUniformSampler {
049        /**
050         * @param rng Generator of uniformly distributed random numbers.
051         * @param lo Lower bound.
052         * @param hi Higher bound.
053         */
054        OpenIntervalContinuousUniformSampler(UniformRandomProvider rng, double lo, double hi) {
055            super(rng, lo, hi);
056        }
057
058        @Override
059        public double sample() {
060            final double x = super.sample();
061            // Due to rounding using a variate u in the open interval (lo, hi) with the original
062            // algorithm may generate a value at the bound. Thus the bound is explicitly tested
063            // and the sample repeated if necessary.
064            if (x == getHi() || x == getLo()) {
065                return sample();
066            }
067            return x;
068        }
069
070        @Override
071        public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
072            return new OpenIntervalContinuousUniformSampler(rng, getLo(), getHi());
073        }
074    }
075
076    /**
077     * Specialization to sample from an open interval {@code (0, 1)} (see RNG-190).
078     *
079     * @since 1.7
080     */
081    private static final class OpenIntervalContinuousUniformSampler01 extends ContinuousUniformSampler {
082        /**
083         * @param rng Generator of uniformly distributed random numbers.
084         */
085        OpenIntervalContinuousUniformSampler01(UniformRandomProvider rng) {
086            super(rng, 0, 1);
087        }
088
089        @Override
090        public double sample() {
091            // Adding a least significant bit before conversion to float.
092            // Output is 2^52 possible rationals.
093            return ((rng.nextLong() >>> 11) | 1L) * 0x1.0p-53d;
094        }
095
096        @Override
097        public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
098            return new OpenIntervalContinuousUniformSampler01(rng);
099        }
100    }
101
102    /**
103     * Create an instance.
104     *
105     * @param rng Generator of uniformly distributed random numbers.
106     * @param lo Lower bound.
107     * @param hi Higher bound.
108     */
109    public ContinuousUniformSampler(UniformRandomProvider rng,
110                                    double lo,
111                                    double hi) {
112        super(null);
113        this.rng = rng;
114        this.lo = lo;
115        this.hi = hi;
116    }
117
118    /** {@inheritDoc} */
119    @Override
120    public double sample() {
121        final double u = rng.nextDouble();
122        return u * hi + (1 - u) * lo;
123    }
124
125    /**
126     * Gets the lower bound. This is deliberately scoped as package private.
127     *
128     * @return the lower bound
129     */
130    double getLo() {
131        return lo;
132    }
133
134    /**
135     * Gets the higher bound. This is deliberately scoped as package private.
136     *
137     * @return the higher bound
138     */
139    double getHi() {
140        return hi;
141    }
142
143    /** {@inheritDoc} */
144    @Override
145    public String toString() {
146        return "Uniform deviate [" + rng.toString() + "]";
147    }
148
149    /**
150     * {@inheritDoc}
151     *
152     * @since 1.3
153     */
154    @Override
155    public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
156        return new ContinuousUniformSampler(rng, lo, hi);
157    }
158
159    /**
160     * Creates a new continuous uniform distribution sampler.
161     *
162     * @param rng Generator of uniformly distributed random numbers.
163     * @param lo Lower bound.
164     * @param hi Higher bound.
165     * @return the sampler
166     * @since 1.3
167     */
168    public static SharedStateContinuousSampler of(UniformRandomProvider rng,
169                                                  double lo,
170                                                  double hi) {
171        return new ContinuousUniformSampler(rng, lo, hi);
172    }
173
174    /**
175     * Creates a new continuous uniform distribution sampler.
176     *
177     * <p>The bounds can be optionally excluded to sample from the open interval
178     * {@code (lower, upper)}. In this case if the bounds have the same sign the
179     * open interval must contain at least 1 double value between the limits; if the
180     * bounds have opposite signs the open interval must contain at least 2 double values
181     * between the limits excluding {@code -0.0}. Thus the interval {@code (-x,x)} will
182     * raise an exception when {@code x} is {@link Double#MIN_VALUE}.
183     *
184     * @param rng Generator of uniformly distributed random numbers.
185     * @param lo Lower bound.
186     * @param hi Higher bound.
187     * @param excludeBounds Set to {@code true} to use the open interval
188     * {@code (lower, upper)}.
189     * @return the sampler
190     * @throws IllegalArgumentException If the open interval is invalid.
191     * @since 1.4
192     */
193    public static SharedStateContinuousSampler of(UniformRandomProvider rng,
194                                                  double lo,
195                                                  double hi,
196                                                  boolean excludeBounds) {
197        if (excludeBounds) {
198            if (!validateOpenInterval(lo, hi)) {
199                throw new IllegalArgumentException("Invalid open interval (" +
200                                                    lo + "," + hi + ")");
201            }
202            // Special case for (0, 1)
203            if (lo == 0 && hi == 1) {
204                return new OpenIntervalContinuousUniformSampler01(rng);
205            }
206            return new OpenIntervalContinuousUniformSampler(rng, lo, hi);
207        }
208        return new ContinuousUniformSampler(rng, lo, hi);
209    }
210
211    /**
212     * Check that the open interval is valid. It must contain at least one double value
213     * between the limits if the signs are the same, or two double values between the limits
214     * if the signs are different (excluding {@code -0.0}).
215     *
216     * @param lo Lower bound.
217     * @param hi Higher bound.
218     * @return false is the interval is invalid
219     */
220    private static boolean validateOpenInterval(double lo, double hi) {
221        // Get the raw bit representation.
222        long bitsx = Double.doubleToRawLongBits(lo);
223        long bitsy = Double.doubleToRawLongBits(hi);
224        // xor will set the sign bit if the signs are different
225        if ((bitsx ^ bitsy) < 0) {
226            // Opposite signs. Drop the sign bit to represent the count of
227            // bits above +0.0. When combined this should be above 2.
228            // This prevents the range (-Double.MIN_VALUE, Double.MIN_VALUE)
229            // which cannot be sampled unless the uniform deviate u=0.5.
230            // (MAX_VALUE has all bits set except the most significant sign bit.)
231            bitsx &= Long.MAX_VALUE;
232            bitsy &= Long.MAX_VALUE;
233            return !lessThanUnsigned(bitsx + bitsy, MIN_ULP_OPPOSITE_SIGN);
234        }
235        // Same signs, subtraction will count the ULP difference.
236        // This should be above 1.
237        return Math.abs(bitsx - bitsy) >= MIN_ULP_SAME_SIGN;
238    }
239
240    /**
241     * Compares two {@code long} values numerically treating the values as unsigned
242     * to test if the first value is less than the second value.
243     *
244     * <p>See Long.compareUnsigned(long, long) in JDK 1.8.
245     *
246     * @param x the first value
247     * @param y the second value
248     * @return true if {@code x < y}
249     */
250    private static boolean lessThanUnsigned(long x, long y) {
251        return x + Long.MIN_VALUE < y + Long.MIN_VALUE;
252    }
253}