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.ZigguratNormalizedGaussianSampler;
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 uses:</p>
029 *
030 * <ul>
031 *   <li>{@link UniformRandomProvider#nextLong()}
032 *   <li>{@link UniformRandomProvider#nextDouble()}
033 * </ul>
034 *
035 * @since 1.1
036 */
037public class UnitSphereSampler implements SharedStateSampler<UnitSphereSampler> {
038    /** Sampler used for generating the individual components of the vectors. */
039    private final NormalizedGaussianSampler sampler;
040    /** Space dimension. */
041    private final int dimension;
042
043    /**
044     * @param dimension Space dimension.
045     * @param rng Generator for the individual components of the vectors.
046     * A shallow copy will be stored in this instance.
047     * @throws IllegalArgumentException If {@code dimension <= 0}
048     */
049    public UnitSphereSampler(int dimension,
050                             UniformRandomProvider rng) {
051        if (dimension <= 0) {
052            throw new IllegalArgumentException("Dimension must be strictly positive");
053        }
054
055        this.dimension = dimension;
056        sampler = new ZigguratNormalizedGaussianSampler(rng);
057    }
058
059    /**
060     * @param rng Generator for the individual components of the vectors.
061     * @param source Source to copy.
062     */
063    private UnitSphereSampler(UniformRandomProvider rng,
064                              UnitSphereSampler source) {
065        // The Gaussian sampler has no shared state so create a new instance
066        sampler = new ZigguratNormalizedGaussianSampler(rng);
067        dimension = source.dimension;
068    }
069
070    /**
071     * @return a random normalized Cartesian vector.
072     */
073    public double[] nextVector() {
074        final double[] v = new double[dimension];
075
076        // Pick a point by choosing a standard Gaussian for each element,
077        // and then normalize to unit length.
078        double normSq = 0;
079        for (int i = 0; i < dimension; i++) {
080            final double comp = sampler.sample();
081            v[i] = comp;
082            normSq += comp * comp;
083        }
084
085        if (normSq == 0) {
086            // Zero-norm vector is discarded.
087            // Using recursion as it is highly unlikely to generate more
088            // than a few such vectors. It also protects against infinite
089            // loop (in case a buggy generator is used), by eventually
090            // raising a "StackOverflowError".
091            return nextVector();
092        }
093
094        final double f = 1 / Math.sqrt(normSq);
095        for (int i = 0; i < dimension; i++) {
096            v[i] *= f;
097        }
098
099        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}