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