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