View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.rng.sampling.distribution;
18  
19  import org.apache.commons.rng.UniformRandomProvider;
20  
21  /**
22   * Sampling from a uniform distribution.
23   *
24   * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
25   *
26   * @since 1.0
27   */
28  public class ContinuousUniformSampler
29      extends SamplerBase
30      implements SharedStateContinuousSampler {
31      /** The minimum ULP gap for the open interval when the doubles have the same sign. */
32      private static final int MIN_ULP_SAME_SIGN = 2;
33      /** The minimum ULP gap for the open interval when the doubles have the opposite sign. */
34      private static final int MIN_ULP_OPPOSITE_SIGN = 3;
35  
36      /** Lower bound. */
37      private final double lo;
38      /** Higher bound. */
39      private final double hi;
40      /** Underlying source of randomness. */
41      private final UniformRandomProvider rng;
42  
43      /**
44       * Specialization to sample from an open interval {@code (lo, hi)}.
45       */
46      private static final class OpenIntervalContinuousUniformSampler extends ContinuousUniformSampler {
47          /**
48           * @param rng Generator of uniformly distributed random numbers.
49           * @param lo Lower bound.
50           * @param hi Higher bound.
51           */
52          OpenIntervalContinuousUniformSampler(UniformRandomProvider rng, double lo, double hi) {
53              super(rng, lo, hi);
54          }
55  
56          @Override
57          public double sample() {
58              final double x = super.sample();
59              // Due to rounding using a variate u in the open interval (0,1) with the original
60              // algorithm may generate a value at the bound. Thus the bound is explicitly tested
61              // and the sample repeated if necessary.
62              if (x == getHi() || x == getLo()) {
63                  return sample();
64              }
65              return x;
66          }
67  
68          @Override
69          public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
70              return new OpenIntervalContinuousUniformSampler(rng, getLo(), getHi());
71          }
72      }
73  
74      /**
75       * Create an instance.
76       *
77       * @param rng Generator of uniformly distributed random numbers.
78       * @param lo Lower bound.
79       * @param hi Higher bound.
80       */
81      public ContinuousUniformSampler(UniformRandomProvider rng,
82                                      double lo,
83                                      double hi) {
84          super(null);
85          this.rng = rng;
86          this.lo = lo;
87          this.hi = hi;
88      }
89  
90      /** {@inheritDoc} */
91      @Override
92      public double sample() {
93          final double u = rng.nextDouble();
94          return u * hi + (1 - u) * lo;
95      }
96  
97      /**
98       * Gets the lower bound. This is deliberately scoped as package private.
99       *
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 }