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 }