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 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       * @param rng Generator of uniformly distributed random numbers.
76       * @param lo Lower bound.
77       * @param hi Higher bound.
78       */
79      public ContinuousUniformSampler(UniformRandomProvider rng,
80                                      double lo,
81                                      double hi) {
82          super(null);
83          this.rng = rng;
84          this.lo = lo;
85          this.hi = hi;
86      }
87  
88      /** {@inheritDoc} */
89      @Override
90      public double sample() {
91          final double u = rng.nextDouble();
92          return u * hi + (1 - u) * lo;
93      }
94  
95      /**
96       * Gets the lower bound. This is deliberately scoped as package private.
97       *
98       * @return the lower bound
99       */
100     double getLo() {
101         return lo;
102     }
103 
104     /**
105      * Gets the higher bound. This is deliberately scoped as package private.
106      *
107      * @return the higher bound
108      */
109     double getHi() {
110         return hi;
111     }
112 
113     /** {@inheritDoc} */
114     @Override
115     public String toString() {
116         return "Uniform deviate [" + rng.toString() + "]";
117     }
118 
119     /**
120      * {@inheritDoc}
121      *
122      * @since 1.3
123      */
124     @Override
125     public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
126         return new ContinuousUniformSampler(rng, lo, hi);
127     }
128 
129     /**
130      * Creates a new continuous uniform distribution sampler.
131      *
132      * @param rng Generator of uniformly distributed random numbers.
133      * @param lo Lower bound.
134      * @param hi Higher bound.
135      * @return the sampler
136      * @since 1.3
137      */
138     public static SharedStateContinuousSampler of(UniformRandomProvider rng,
139                                                   double lo,
140                                                   double hi) {
141         return new ContinuousUniformSampler(rng, lo, hi);
142     }
143 
144     /**
145      * Creates a new continuous uniform distribution sampler.
146      *
147      * <p>The bounds can be optionally excluded to sample from the open interval
148      * {@code (lower, upper)}. In this case if the bounds have the same sign the
149      * open interval must contain at least 1 double value between the limits; if the
150      * bounds have opposite signs the open interval must contain at least 2 double values
151      * between the limits excluding {@code -0.0}. Thus the interval {@code (-x,x)} will
152      * raise an exception when {@code x} is {@link Double#MIN_VALUE}.
153      *
154      * @param rng Generator of uniformly distributed random numbers.
155      * @param lo Lower bound.
156      * @param hi Higher bound.
157      * @param excludeBounds Set to {@code true} to use the open interval
158      * {@code (lower, upper)}.
159      * @return the sampler
160      * @throws IllegalArgumentException If the open interval is invalid.
161      * @since 1.4
162      */
163     public static SharedStateContinuousSampler of(UniformRandomProvider rng,
164                                                   double lo,
165                                                   double hi,
166                                                   boolean excludeBounds) {
167         if (excludeBounds) {
168             if (!validateOpenInterval(lo, hi)) {
169                 throw new IllegalArgumentException("Invalid open interval (" +
170                                                     lo + "," + hi + ")");
171             }
172             return new OpenIntervalContinuousUniformSampler(rng, lo, hi);
173         }
174         return new ContinuousUniformSampler(rng, lo, hi);
175     }
176 
177     /**
178      * Check that the open interval is valid. It must contain at least one double value
179      * between the limits if the signs are the same, or two double values between the limits
180      * if the signs are different (excluding {@code -0.0}).
181      *
182      * @param lo Lower bound.
183      * @param hi Higher bound.
184      * @return false is the interval is invalid
185      */
186     private static boolean validateOpenInterval(double lo, double hi) {
187         // Get the raw bit representation.
188         long bitsx = Double.doubleToRawLongBits(lo);
189         long bitsy = Double.doubleToRawLongBits(hi);
190         // xor will set the sign bit if the signs are different
191         if ((bitsx ^ bitsy) < 0) {
192             // Opposite signs. Drop the sign bit to represent the count of
193             // bits above +0.0. When combined this should be above 2.
194             // This prevents the range (-Double.MIN_VALUE, Double.MIN_VALUE)
195             // which cannot be sampled unless the uniform deviate u=0.5.
196             // (MAX_VALUE has all bits set except the most significant sign bit.)
197             bitsx &= Long.MAX_VALUE;
198             bitsy &= Long.MAX_VALUE;
199             return !lessThanUnsigned(bitsx + bitsy, MIN_ULP_OPPOSITE_SIGN);
200         }
201         // Same signs, subtraction will count the ULP difference.
202         // This should be above 1.
203         return Math.abs(bitsx - bitsy) >= MIN_ULP_SAME_SIGN;
204     }
205 
206     /**
207      * Compares two {@code long} values numerically treating the values as unsigned
208      * to test if the first value is less than the second value.
209      *
210      * <p>See Long.compareUnsigned(long, long) in JDK 1.8.
211      *
212      * @param x the first value
213      * @param y the second value
214      * @return true if {@code x < y}
215      */
216     private static boolean lessThanUnsigned(long x, long y) {
217         return x + Long.MIN_VALUE < y + Long.MIN_VALUE;
218     }
219 }