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}