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}