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 */ 017 018package org.apache.commons.rng.sampling.shape; 019 020import org.apache.commons.rng.UniformRandomProvider; 021import org.apache.commons.rng.sampling.SharedStateObjectSampler; 022import org.apache.commons.rng.sampling.distribution.ContinuousSampler; 023import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler; 024import org.apache.commons.rng.sampling.distribution.ZigguratSampler; 025 026/** 027 * Generate coordinates <a href="http://mathworld.wolfram.com/BallPointPicking.html"> 028 * uniformly distributed within the unit n-ball</a>. 029 * 030 * <p>Sampling uses:</p> 031 * 032 * <ul> 033 * <li>{@link UniformRandomProvider#nextLong()} 034 * <li>{@link UniformRandomProvider#nextDouble()} (only for dimensions above 2) 035 * </ul> 036 * 037 * @since 1.4 038 */ 039public abstract class UnitBallSampler implements SharedStateObjectSampler<double[]> { 040 /** The dimension for 1D sampling. */ 041 private static final int ONE_D = 1; 042 /** The dimension for 2D sampling. */ 043 private static final int TWO_D = 2; 044 /** The dimension for 3D sampling. */ 045 private static final int THREE_D = 3; 046 /** 047 * The multiplier to convert the least significant 53-bits of a {@code long} to a {@code double}. 048 * Taken from o.a.c.rng.core.utils.NumberFactory. 049 * 050 * <p>This is equivalent to {@code 1.0 / (1L << 53)}. 051 */ 052 private static final double DOUBLE_MULTIPLIER = 0x1.0p-53d; 053 054 /** 055 * Sample uniformly from a 1D unit line. 056 */ 057 private static final class UnitBallSampler1D extends UnitBallSampler { 058 /** The source of randomness. */ 059 private final UniformRandomProvider rng; 060 061 /** 062 * @param rng Source of randomness. 063 */ 064 UnitBallSampler1D(UniformRandomProvider rng) { 065 this.rng = rng; 066 } 067 068 @Override 069 public double[] sample() { 070 return new double[] {makeSignedDouble(rng.nextLong())}; 071 } 072 073 @Override 074 public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) { 075 return new UnitBallSampler1D(rng); 076 } 077 } 078 079 /** 080 * Sample uniformly from a 2D unit disk. 081 */ 082 private static final class UnitBallSampler2D extends UnitBallSampler { 083 /** The source of randomness. */ 084 private final UniformRandomProvider rng; 085 086 /** 087 * @param rng Source of randomness. 088 */ 089 UnitBallSampler2D(UniformRandomProvider rng) { 090 this.rng = rng; 091 } 092 093 @Override 094 public double[] sample() { 095 // Generate via rejection method of a circle inside a square of edge length 2. 096 // This should compute approximately 2^2 / pi = 1.27 square positions per sample. 097 double x; 098 double y; 099 do { 100 x = makeSignedDouble(rng.nextLong()); 101 y = makeSignedDouble(rng.nextLong()); 102 } while (x * x + y * y > 1.0); 103 return new double[] {x, y}; 104 } 105 106 @Override 107 public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) { 108 return new UnitBallSampler2D(rng); 109 } 110 } 111 112 /** 113 * Sample uniformly from a 3D unit ball. This is an non-array based specialisation of 114 * {@link UnitBallSamplerND} for performance. 115 */ 116 private static final class UnitBallSampler3D extends UnitBallSampler { 117 /** The standard normal distribution. */ 118 private final NormalizedGaussianSampler normal; 119 /** The exponential distribution (mean=1). */ 120 private final ContinuousSampler exp; 121 122 /** 123 * @param rng Source of randomness. 124 */ 125 UnitBallSampler3D(UniformRandomProvider rng) { 126 normal = ZigguratSampler.NormalizedGaussian.of(rng); 127 // Require an Exponential(mean=2). 128 // Here we use mean = 1 and scale the output later. 129 exp = ZigguratSampler.Exponential.of(rng); 130 } 131 132 @Override 133 public double[] sample() { 134 final double x = normal.sample(); 135 final double y = normal.sample(); 136 final double z = normal.sample(); 137 // Include the exponential sample. It has mean 1 so multiply by 2. 138 final double sum = exp.sample() * 2 + x * x + y * y + z * z; 139 // Note: Handle the possibility of a zero sum and invalid inverse 140 if (sum == 0) { 141 return sample(); 142 } 143 final double f = 1.0 / Math.sqrt(sum); 144 return new double[] {x * f, y * f, z * f}; 145 } 146 147 @Override 148 public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) { 149 return new UnitBallSampler3D(rng); 150 } 151 } 152 153 /** 154 * Sample using ball point picking. 155 * @see <a href="https://mathworld.wolfram.com/BallPointPicking.html">Ball point picking</a> 156 */ 157 private static final class UnitBallSamplerND extends UnitBallSampler { 158 /** The dimension. */ 159 private final int dimension; 160 /** The standard normal distribution. */ 161 private final NormalizedGaussianSampler normal; 162 /** The exponential distribution (mean=1). */ 163 private final ContinuousSampler exp; 164 165 /** 166 * @param rng Source of randomness. 167 * @param dimension Space dimension. 168 */ 169 UnitBallSamplerND(UniformRandomProvider rng, int dimension) { 170 this.dimension = dimension; 171 normal = ZigguratSampler.NormalizedGaussian.of(rng); 172 // Require an Exponential(mean=2). 173 // Here we use mean = 1 and scale the output later. 174 exp = ZigguratSampler.Exponential.of(rng); 175 } 176 177 @Override 178 public double[] sample() { 179 final double[] sample = new double[dimension]; 180 // Include the exponential sample. It has mean 1 so multiply by 2. 181 double sum = exp.sample() * 2; 182 for (int i = 0; i < dimension; i++) { 183 final double x = normal.sample(); 184 sum += x * x; 185 sample[i] = x; 186 } 187 // Note: Handle the possibility of a zero sum and invalid inverse 188 if (sum == 0) { 189 return sample(); 190 } 191 final double f = 1.0 / Math.sqrt(sum); 192 for (int i = 0; i < dimension; i++) { 193 sample[i] *= f; 194 } 195 return sample; 196 } 197 198 @Override 199 public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) { 200 return new UnitBallSamplerND(rng, dimension); 201 } 202 } 203 204 /** 205 * Create an instance. 206 */ 207 public UnitBallSampler() {} 208 209 /** 210 * @return a random Cartesian coordinate within the unit n-ball. 211 */ 212 @Override 213 public abstract double[] sample(); 214 215 /** {@inheritDoc} */ 216 // Redeclare the signature to return a UnitBallSampler not a SharedStateObjectSampler<double[]> 217 @Override 218 public abstract UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng); 219 220 /** 221 * Create a unit n-ball sampler for the given dimension. 222 * Sampled points are uniformly distributed within the unit n-ball. 223 * 224 * <p>Sampling is supported in dimensions of 1 or above. 225 * 226 * @param rng Source of randomness. 227 * @param dimension Space dimension. 228 * @return the sampler 229 * @throws IllegalArgumentException If {@code dimension <= 0} 230 */ 231 public static UnitBallSampler of(UniformRandomProvider rng, 232 int dimension) { 233 if (dimension <= 0) { 234 throw new IllegalArgumentException("Dimension must be strictly positive"); 235 } else if (dimension == ONE_D) { 236 return new UnitBallSampler1D(rng); 237 } else if (dimension == TWO_D) { 238 return new UnitBallSampler2D(rng); 239 } else if (dimension == THREE_D) { 240 return new UnitBallSampler3D(rng); 241 } 242 return new UnitBallSamplerND(rng, dimension); 243 } 244 245 /** 246 * Creates a signed double in the range {@code [-1, 1)}. The magnitude is sampled evenly 247 * from the 2<sup>54</sup> dyadic rationals in the range. 248 * 249 * <p>Note: This method will not return samples for both -0.0 and 0.0. 250 * 251 * @param bits the bits 252 * @return the double 253 */ 254 private static double makeSignedDouble(long bits) { 255 // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a signed 256 // shift of 10 in place of an unsigned shift of 11. 257 // Use the upper 54 bits on the assumption they are more random. 258 // The sign bit is maintained by the signed shift. 259 // The next 53 bits generates a magnitude in the range [0, 2^53) or [-2^53, 0). 260 return (bits >> 10) * DOUBLE_MULTIPLIER; 261 } 262}