UnitSphereSampler.java
- /*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
- package org.apache.commons.rng.sampling;
- import org.apache.commons.rng.UniformRandomProvider;
- import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
- import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
- /**
- * Generate vectors <a href="http://mathworld.wolfram.com/SpherePointPicking.html">
- * isotropically located on the surface of a sphere</a>.
- *
- * <p>Sampling in 2 or more dimensions uses:</p>
- *
- * <ul>
- * <li>{@link UniformRandomProvider#nextLong()}
- * <li>{@link UniformRandomProvider#nextDouble()}
- * </ul>
- *
- * <p>Sampling in 1D uses the sign bit from {@link UniformRandomProvider#nextInt()} to set
- * the direction of the vector.
- *
- * @since 1.1
- */
- public class UnitSphereSampler implements SharedStateObjectSampler<double[]> {
- /** The dimension for 1D sampling. */
- private static final int ONE_D = 1;
- /** The dimension for 2D sampling. */
- private static final int TWO_D = 2;
- /** The dimension for 3D sampling. */
- private static final int THREE_D = 3;
- /**
- * The mask to extract the second bit from an integer
- * (naming starts at bit 1 for the least significant bit).
- * The masked integer will have a value 0 or 2.
- */
- private static final int MASK_BIT_2 = 0x2;
- /** The internal sampler optimised for the dimension. */
- private final UnitSphereSampler delegate;
- /**
- * Sample uniformly from the ends of a 1D unit line.
- */
- private static final class UnitSphereSampler1D extends UnitSphereSampler {
- /** The source of randomness. */
- private final UniformRandomProvider rng;
- /**
- * @param rng Source of randomness.
- */
- UnitSphereSampler1D(UniformRandomProvider rng) {
- this.rng = rng;
- }
- @Override
- public double[] sample() {
- // Either:
- // 1 - 0 = 1
- // 1 - 2 = -1
- // Use the sign bit
- return new double[] {1.0 - ((rng.nextInt() >>> 30) & MASK_BIT_2)};
- }
- @Override
- public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
- return new UnitSphereSampler1D(rng);
- }
- }
- /**
- * Sample uniformly from a 2D unit circle.
- * This is a 2D specialisation of the UnitSphereSamplerND.
- */
- private static final class UnitSphereSampler2D extends UnitSphereSampler {
- /** Sampler used for generating the individual components of the vectors. */
- private final NormalizedGaussianSampler sampler;
- /**
- * @param rng Source of randomness.
- */
- UnitSphereSampler2D(UniformRandomProvider rng) {
- sampler = ZigguratSampler.NormalizedGaussian.of(rng);
- }
- @Override
- public double[] sample() {
- final double x = sampler.sample();
- final double y = sampler.sample();
- final double sum = x * x + y * y;
- if (sum == 0) {
- // Zero-norm vector is discarded.
- return sample();
- }
- final double f = 1.0 / Math.sqrt(sum);
- return new double[] {x * f, y * f};
- }
- @Override
- public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
- return new UnitSphereSampler2D(rng);
- }
- }
- /**
- * Sample uniformly from a 3D unit sphere.
- * This is a 3D specialisation of the UnitSphereSamplerND.
- */
- private static final class UnitSphereSampler3D extends UnitSphereSampler {
- /** Sampler used for generating the individual components of the vectors. */
- private final NormalizedGaussianSampler sampler;
- /**
- * @param rng Source of randomness.
- */
- UnitSphereSampler3D(UniformRandomProvider rng) {
- sampler = ZigguratSampler.NormalizedGaussian.of(rng);
- }
- @Override
- public double[] sample() {
- final double x = sampler.sample();
- final double y = sampler.sample();
- final double z = sampler.sample();
- final double sum = x * x + y * y + z * z;
- if (sum == 0) {
- // Zero-norm vector is discarded.
- return sample();
- }
- final double f = 1.0 / Math.sqrt(sum);
- return new double[] {x * f, y * f, z * f};
- }
- @Override
- public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
- return new UnitSphereSampler3D(rng);
- }
- }
- /**
- * Sample uniformly from a ND unit sphere.
- */
- private static final class UnitSphereSamplerND extends UnitSphereSampler {
- /** Space dimension. */
- private final int dimension;
- /** Sampler used for generating the individual components of the vectors. */
- private final NormalizedGaussianSampler sampler;
- /**
- * @param rng Source of randomness.
- * @param dimension Space dimension.
- */
- UnitSphereSamplerND(UniformRandomProvider rng, int dimension) {
- this.dimension = dimension;
- sampler = ZigguratSampler.NormalizedGaussian.of(rng);
- }
- @Override
- public double[] sample() {
- final double[] v = new double[dimension];
- // Pick a point by choosing a standard Gaussian for each element,
- // and then normalize to unit length.
- double sum = 0;
- for (int i = 0; i < dimension; i++) {
- final double x = sampler.sample();
- v[i] = x;
- sum += x * x;
- }
- if (sum == 0) {
- // Zero-norm vector is discarded.
- // Using recursion as it is highly unlikely to generate more
- // than a few such vectors. It also protects against infinite
- // loop (in case a buggy generator is used), by eventually
- // raising a "StackOverflowError".
- return sample();
- }
- final double f = 1 / Math.sqrt(sum);
- for (int i = 0; i < dimension; i++) {
- v[i] *= f;
- }
- return v;
- }
- @Override
- public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
- return new UnitSphereSamplerND(rng, dimension);
- }
- }
- /**
- * This instance delegates sampling. Use the factory method
- * {@link #of(UniformRandomProvider, int)} to create an optimal sampler.
- *
- * @param dimension Space dimension.
- * @param rng Generator for the individual components of the vectors.
- * A shallow copy will be stored in this instance.
- * @throws IllegalArgumentException If {@code dimension <= 0}
- * @deprecated Use {@link #of(UniformRandomProvider, int)}.
- */
- @Deprecated
- public UnitSphereSampler(int dimension,
- UniformRandomProvider rng) {
- this(of(rng, dimension));
- }
- /**
- * Private constructor used by deprecated constructor used to prevent partially
- * initialized object if the construction of the delegate throws.
- * In future versions the public constructor should be removed and the class made abstract.
- *
- * @param delegate Delegate.
- */
- private UnitSphereSampler(UnitSphereSampler delegate) {
- this.delegate = delegate;
- }
- /**
- * Private constructor used by sub-class specialisations.
- * In future versions the public constructor should be removed and the class made abstract.
- */
- private UnitSphereSampler() {
- delegate = null;
- }
- /**
- * @return a random normalized Cartesian vector.
- * @since 1.4
- */
- @Override
- public double[] sample() {
- return delegate.sample();
- }
- /**
- * @return a random normalized Cartesian vector.
- * @deprecated Use {@link #sample()}.
- */
- @Deprecated
- public double[] nextVector() {
- return sample();
- }
- /**
- * {@inheritDoc}
- *
- * @since 1.3
- */
- @Override
- public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
- return delegate.withUniformRandomProvider(rng);
- }
- /**
- * Create a unit sphere sampler for the given dimension.
- *
- * @param rng Generator for the individual components of the vectors. A shallow
- * copy will be stored in the sampler.
- * @param dimension Space dimension.
- * @return the sampler
- * @throws IllegalArgumentException If {@code dimension <= 0}
- *
- * @since 1.4
- */
- public static UnitSphereSampler of(UniformRandomProvider rng,
- int dimension) {
- if (dimension <= 0) {
- throw new IllegalArgumentException("Dimension must be strictly positive");
- } else if (dimension == ONE_D) {
- return new UnitSphereSampler1D(rng);
- } else if (dimension == TWO_D) {
- return new UnitSphereSampler2D(rng);
- } else if (dimension == THREE_D) {
- return new UnitSphereSampler3D(rng);
- }
- return new UnitSphereSamplerND(rng, dimension);
- }
- }