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
18 package org.apache.commons.rng.sampling;
19
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
22 import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
23
24 /**
25 * Generate vectors <a href="http://mathworld.wolfram.com/SpherePointPicking.html">
26 * isotropically located on the surface of a sphere</a>.
27 *
28 * <p>Sampling in 2 or more dimensions uses:</p>
29 *
30 * <ul>
31 * <li>{@link UniformRandomProvider#nextLong()}
32 * <li>{@link UniformRandomProvider#nextDouble()}
33 * </ul>
34 *
35 * <p>Sampling in 1D uses the sign bit from {@link UniformRandomProvider#nextInt()} to set
36 * the direction of the vector.
37 *
38 * @since 1.1
39 */
40 public class UnitSphereSampler implements SharedStateObjectSampler<double[]> {
41 /** The dimension for 1D sampling. */
42 private static final int ONE_D = 1;
43 /** The dimension for 2D sampling. */
44 private static final int TWO_D = 2;
45 /** The dimension for 3D sampling. */
46 private static final int THREE_D = 3;
47 /**
48 * The mask to extract the second bit from an integer
49 * (naming starts at bit 1 for the least significant bit).
50 * The masked integer will have a value 0 or 2.
51 */
52 private static final int MASK_BIT_2 = 0x2;
53
54 /** The internal sampler optimised for the dimension. */
55 private final UnitSphereSampler delegate;
56
57 /**
58 * Sample uniformly from the ends of a 1D unit line.
59 */
60 private static final class UnitSphereSampler1D extends UnitSphereSampler {
61 /** The source of randomness. */
62 private final UniformRandomProvider rng;
63
64 /**
65 * @param rng Source of randomness.
66 */
67 UnitSphereSampler1D(UniformRandomProvider rng) {
68 this.rng = rng;
69 }
70
71 @Override
72 public double[] sample() {
73 // Either:
74 // 1 - 0 = 1
75 // 1 - 2 = -1
76 // Use the sign bit
77 return new double[] {1.0 - ((rng.nextInt() >>> 30) & MASK_BIT_2)};
78 }
79
80 @Override
81 public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
82 return new UnitSphereSampler1D(rng);
83 }
84 }
85
86 /**
87 * Sample uniformly from a 2D unit circle.
88 * This is a 2D specialisation of the UnitSphereSamplerND.
89 */
90 private static final class UnitSphereSampler2D extends UnitSphereSampler {
91 /** Sampler used for generating the individual components of the vectors. */
92 private final NormalizedGaussianSampler sampler;
93
94 /**
95 * @param rng Source of randomness.
96 */
97 UnitSphereSampler2D(UniformRandomProvider rng) {
98 sampler = ZigguratSampler.NormalizedGaussian.of(rng);
99 }
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 }