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; 019 020import org.apache.commons.rng.UniformRandomProvider; 021import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler; 022import org.apache.commons.rng.sampling.distribution.ZigguratSampler; 023 024/** 025 * Generate vectors <a href="http://mathworld.wolfram.com/SpherePointPicking.html"> 026 * isotropically located on the surface of a sphere</a>. 027 * 028 * <p>Sampling in 2 or more dimensions uses:</p> 029 * 030 * <ul> 031 * <li>{@link UniformRandomProvider#nextLong()} 032 * <li>{@link UniformRandomProvider#nextDouble()} 033 * </ul> 034 * 035 * <p>Sampling in 1D uses the sign bit from {@link UniformRandomProvider#nextInt()} to set 036 * the direction of the vector. 037 * 038 * @since 1.1 039 */ 040public class UnitSphereSampler implements SharedStateObjectSampler<double[]> { 041 /** The dimension for 1D sampling. */ 042 private static final int ONE_D = 1; 043 /** The dimension for 2D sampling. */ 044 private static final int TWO_D = 2; 045 /** The dimension for 3D sampling. */ 046 private static final int THREE_D = 3; 047 /** 048 * The mask to extract the second bit from an integer 049 * (naming starts at bit 1 for the least significant bit). 050 * The masked integer will have a value 0 or 2. 051 */ 052 private static final int MASK_BIT_2 = 0x2; 053 054 /** The internal sampler optimised for the dimension. */ 055 private final UnitSphereSampler delegate; 056 057 /** 058 * Sample uniformly from the ends of a 1D unit line. 059 */ 060 private static final class UnitSphereSampler1D extends UnitSphereSampler { 061 /** The source of randomness. */ 062 private final UniformRandomProvider rng; 063 064 /** 065 * @param rng Source of randomness. 066 */ 067 UnitSphereSampler1D(UniformRandomProvider rng) { 068 this.rng = rng; 069 } 070 071 @Override 072 public double[] sample() { 073 // Either: 074 // 1 - 0 = 1 075 // 1 - 2 = -1 076 // Use the sign bit 077 return new double[] {1.0 - ((rng.nextInt() >>> 30) & MASK_BIT_2)}; 078 } 079 080 @Override 081 public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { 082 return new UnitSphereSampler1D(rng); 083 } 084 } 085 086 /** 087 * Sample uniformly from a 2D unit circle. 088 * This is a 2D specialisation of the UnitSphereSamplerND. 089 */ 090 private static final class UnitSphereSampler2D extends UnitSphereSampler { 091 /** Sampler used for generating the individual components of the vectors. */ 092 private final NormalizedGaussianSampler sampler; 093 094 /** 095 * @param rng Source of randomness. 096 */ 097 UnitSphereSampler2D(UniformRandomProvider rng) { 098 sampler = ZigguratSampler.NormalizedGaussian.of(rng); 099 } 100 101 @Override 102 public double[] sample() { 103 final double x = sampler.sample(); 104 final double y = sampler.sample(); 105 final double sum = x * x + y * y; 106 107 if (sum == 0) { 108 // Zero-norm vector is discarded. 109 return sample(); 110 } 111 112 final double f = 1.0 / Math.sqrt(sum); 113 return new double[] {x * f, y * f}; 114 } 115 116 @Override 117 public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { 118 return new UnitSphereSampler2D(rng); 119 } 120 } 121 122 /** 123 * Sample uniformly from a 3D unit sphere. 124 * This is a 3D specialisation of the UnitSphereSamplerND. 125 */ 126 private static final class UnitSphereSampler3D extends UnitSphereSampler { 127 /** Sampler used for generating the individual components of the vectors. */ 128 private final NormalizedGaussianSampler sampler; 129 130 /** 131 * @param rng Source of randomness. 132 */ 133 UnitSphereSampler3D(UniformRandomProvider rng) { 134 sampler = ZigguratSampler.NormalizedGaussian.of(rng); 135 } 136 137 @Override 138 public double[] sample() { 139 final double x = sampler.sample(); 140 final double y = sampler.sample(); 141 final double z = sampler.sample(); 142 final double sum = x * x + y * y + z * z; 143 144 if (sum == 0) { 145 // Zero-norm vector is discarded. 146 return sample(); 147 } 148 149 final double f = 1.0 / Math.sqrt(sum); 150 return new double[] {x * f, y * f, z * f}; 151 } 152 153 @Override 154 public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { 155 return new UnitSphereSampler3D(rng); 156 } 157 } 158 159 /** 160 * Sample uniformly from a ND unit sphere. 161 */ 162 private static final class UnitSphereSamplerND extends UnitSphereSampler { 163 /** Space dimension. */ 164 private final int dimension; 165 /** Sampler used for generating the individual components of the vectors. */ 166 private final NormalizedGaussianSampler sampler; 167 168 /** 169 * @param rng Source of randomness. 170 * @param dimension Space dimension. 171 */ 172 UnitSphereSamplerND(UniformRandomProvider rng, int dimension) { 173 this.dimension = dimension; 174 sampler = ZigguratSampler.NormalizedGaussian.of(rng); 175 } 176 177 @Override 178 public double[] sample() { 179 final double[] v = new double[dimension]; 180 181 // Pick a point by choosing a standard Gaussian for each element, 182 // and then normalize to unit length. 183 double sum = 0; 184 for (int i = 0; i < dimension; i++) { 185 final double x = sampler.sample(); 186 v[i] = x; 187 sum += x * x; 188 } 189 190 if (sum == 0) { 191 // Zero-norm vector is discarded. 192 // Using recursion as it is highly unlikely to generate more 193 // than a few such vectors. It also protects against infinite 194 // loop (in case a buggy generator is used), by eventually 195 // raising a "StackOverflowError". 196 return sample(); 197 } 198 199 final double f = 1 / Math.sqrt(sum); 200 for (int i = 0; i < dimension; i++) { 201 v[i] *= f; 202 } 203 204 return v; 205 } 206 207 @Override 208 public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { 209 return new UnitSphereSamplerND(rng, dimension); 210 } 211 } 212 213 /** 214 * This instance delegates sampling. Use the factory method 215 * {@link #of(UniformRandomProvider, int)} to create an optimal sampler. 216 * 217 * @param dimension Space dimension. 218 * @param rng Generator for the individual components of the vectors. 219 * A shallow copy will be stored in this instance. 220 * @throws IllegalArgumentException If {@code dimension <= 0} 221 * @deprecated Use {@link #of(UniformRandomProvider, int)}. 222 */ 223 @Deprecated 224 public UnitSphereSampler(int dimension, 225 UniformRandomProvider rng) { 226 this(of(rng, dimension)); 227 } 228 229 /** 230 * Private constructor used by deprecated constructor used to prevent partially 231 * initialized object if the construction of the delegate throws. 232 * In future versions the public constructor should be removed and the class made abstract. 233 * 234 * @param delegate Delegate. 235 */ 236 private UnitSphereSampler(UnitSphereSampler delegate) { 237 this.delegate = delegate; 238 } 239 240 /** 241 * Private constructor used by sub-class specialisations. 242 * In future versions the public constructor should be removed and the class made abstract. 243 */ 244 private UnitSphereSampler() { 245 delegate = null; 246 } 247 248 /** 249 * @return a random normalized Cartesian vector. 250 * @since 1.4 251 */ 252 @Override 253 public double[] sample() { 254 return delegate.sample(); 255 } 256 257 /** 258 * @return a random normalized Cartesian vector. 259 * @deprecated Use {@link #sample()}. 260 */ 261 @Deprecated 262 public double[] nextVector() { 263 return sample(); 264 } 265 266 /** 267 * {@inheritDoc} 268 * 269 * @since 1.3 270 */ 271 @Override 272 public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { 273 return delegate.withUniformRandomProvider(rng); 274 } 275 276 /** 277 * Create a unit sphere sampler for the given dimension. 278 * 279 * @param rng Generator for the individual components of the vectors. A shallow 280 * copy will be stored in the sampler. 281 * @param dimension Space dimension. 282 * @return the sampler 283 * @throws IllegalArgumentException If {@code dimension <= 0} 284 * 285 * @since 1.4 286 */ 287 public static UnitSphereSampler of(UniformRandomProvider rng, 288 int dimension) { 289 if (dimension <= 0) { 290 throw new IllegalArgumentException("Dimension must be strictly positive"); 291 } else if (dimension == ONE_D) { 292 return new UnitSphereSampler1D(rng); 293 } else if (dimension == TWO_D) { 294 return new UnitSphereSampler2D(rng); 295 } else if (dimension == THREE_D) { 296 return new UnitSphereSampler3D(rng); 297 } 298 return new UnitSphereSamplerND(rng, dimension); 299 } 300}