View Javadoc
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.ZigguratNormalizedGaussianSampler;
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 uses:</p>
29   *
30   * <ul>
31   *   <li>{@link UniformRandomProvider#nextLong()}
32   *   <li>{@link UniformRandomProvider#nextDouble()}
33   * </ul>
34   *
35   * @since 1.1
36   */
37  public class UnitSphereSampler implements SharedStateSampler<UnitSphereSampler> {
38      /** Sampler used for generating the individual components of the vectors. */
39      private final NormalizedGaussianSampler sampler;
40      /** Space dimension. */
41      private final int dimension;
42  
43      /**
44       * @param dimension Space dimension.
45       * @param rng Generator for the individual components of the vectors.
46       * A shallow copy will be stored in this instance.
47       * @throws IllegalArgumentException If {@code dimension <= 0}
48       */
49      public UnitSphereSampler(int dimension,
50                               UniformRandomProvider rng) {
51          if (dimension <= 0) {
52              throw new IllegalArgumentException("Dimension must be strictly positive");
53          }
54  
55          this.dimension = dimension;
56          sampler = new ZigguratNormalizedGaussianSampler(rng);
57      }
58  
59      /**
60       * @param rng Generator for the individual components of the vectors.
61       * @param source Source to copy.
62       */
63      private UnitSphereSampler(UniformRandomProvider rng,
64                                UnitSphereSampler source) {
65          // The Gaussian sampler has no shared state so create a new instance
66          sampler = new ZigguratNormalizedGaussianSampler(rng);
67          dimension = source.dimension;
68      }
69  
70      /**
71       * @return a random normalized Cartesian vector.
72       */
73      public double[] nextVector() {
74          final double[] v = new double[dimension];
75  
76          // Pick a point by choosing a standard Gaussian for each element,
77          // and then normalize to unit length.
78          double normSq = 0;
79          for (int i = 0; i < dimension; i++) {
80              final double comp = sampler.sample();
81              v[i] = comp;
82              normSq += comp * comp;
83          }
84  
85          if (normSq == 0) {
86              // Zero-norm vector is discarded.
87              // Using recursion as it is highly unlikely to generate more
88              // than a few such vectors. It also protects against infinite
89              // loop (in case a buggy generator is used), by eventually
90              // raising a "StackOverflowError".
91              return nextVector();
92          }
93  
94          final double f = 1 / Math.sqrt(normSq);
95          for (int i = 0; i < dimension; i++) {
96              v[i] *= f;
97          }
98  
99          return v;
100     }
101 
102     /**
103      * {@inheritDoc}
104      *
105      * @since 1.3
106      */
107     @Override
108     public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
109         return new UnitSphereSampler(rng, this);
110     }
111 }