UnitSphereSampler.java

  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. package org.apache.commons.rng.sampling;

  18. import org.apache.commons.rng.UniformRandomProvider;
  19. import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
  20. import org.apache.commons.rng.sampling.distribution.ZigguratSampler;

  21. /**
  22.  * Generate vectors <a href="http://mathworld.wolfram.com/SpherePointPicking.html">
  23.  * isotropically located on the surface of a sphere</a>.
  24.  *
  25.  * <p>Sampling in 2 or more dimensions uses:</p>
  26.  *
  27.  * <ul>
  28.  *   <li>{@link UniformRandomProvider#nextLong()}
  29.  *   <li>{@link UniformRandomProvider#nextDouble()}
  30.  * </ul>
  31.  *
  32.  * <p>Sampling in 1D uses the sign bit from {@link UniformRandomProvider#nextInt()} to set
  33.  * the direction of the vector.
  34.  *
  35.  * @since 1.1
  36.  */
  37. public class UnitSphereSampler implements SharedStateObjectSampler<double[]> {
  38.     /** The dimension for 1D sampling. */
  39.     private static final int ONE_D = 1;
  40.     /** The dimension for 2D sampling. */
  41.     private static final int TWO_D = 2;
  42.     /** The dimension for 3D sampling. */
  43.     private static final int THREE_D = 3;
  44.     /**
  45.      * The mask to extract the second bit from an integer
  46.      * (naming starts at bit 1 for the least significant bit).
  47.      * The masked integer will have a value 0 or 2.
  48.      */
  49.     private static final int MASK_BIT_2 = 0x2;

  50.     /** The internal sampler optimised for the dimension. */
  51.     private final UnitSphereSampler delegate;

  52.     /**
  53.      * Sample uniformly from the ends of a 1D unit line.
  54.      */
  55.     private static final class UnitSphereSampler1D extends UnitSphereSampler {
  56.         /** The source of randomness. */
  57.         private final UniformRandomProvider rng;

  58.         /**
  59.          * @param rng Source of randomness.
  60.          */
  61.         UnitSphereSampler1D(UniformRandomProvider rng) {
  62.             this.rng = rng;
  63.         }

  64.         @Override
  65.         public double[] sample() {
  66.             // Either:
  67.             // 1 - 0 = 1
  68.             // 1 - 2 = -1
  69.             // Use the sign bit
  70.             return new double[] {1.0 - ((rng.nextInt() >>> 30) & MASK_BIT_2)};
  71.         }

  72.         @Override
  73.         public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
  74.             return new UnitSphereSampler1D(rng);
  75.         }
  76.     }

  77.     /**
  78.      * Sample uniformly from a 2D unit circle.
  79.      * This is a 2D specialisation of the UnitSphereSamplerND.
  80.      */
  81.     private static final class UnitSphereSampler2D extends UnitSphereSampler {
  82.         /** Sampler used for generating the individual components of the vectors. */
  83.         private final NormalizedGaussianSampler sampler;

  84.         /**
  85.          * @param rng Source of randomness.
  86.          */
  87.         UnitSphereSampler2D(UniformRandomProvider rng) {
  88.             sampler = ZigguratSampler.NormalizedGaussian.of(rng);
  89.         }

  90.         @Override
  91.         public double[] sample() {
  92.             final double x = sampler.sample();
  93.             final double y = sampler.sample();
  94.             final double sum = x * x + y * y;

  95.             if (sum == 0) {
  96.                 // Zero-norm vector is discarded.
  97.                 return sample();
  98.             }

  99.             final double f = 1.0 / Math.sqrt(sum);
  100.             return new double[] {x * f, y * f};
  101.         }

  102.         @Override
  103.         public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
  104.             return new UnitSphereSampler2D(rng);
  105.         }
  106.     }

  107.     /**
  108.      * Sample uniformly from a 3D unit sphere.
  109.      * This is a 3D specialisation of the UnitSphereSamplerND.
  110.      */
  111.     private static final class UnitSphereSampler3D extends UnitSphereSampler {
  112.         /** Sampler used for generating the individual components of the vectors. */
  113.         private final NormalizedGaussianSampler sampler;

  114.         /**
  115.          * @param rng Source of randomness.
  116.          */
  117.         UnitSphereSampler3D(UniformRandomProvider rng) {
  118.             sampler = ZigguratSampler.NormalizedGaussian.of(rng);
  119.         }

  120.         @Override
  121.         public double[] sample() {
  122.             final double x = sampler.sample();
  123.             final double y = sampler.sample();
  124.             final double z = sampler.sample();
  125.             final double sum = x * x + y * y + z * z;

  126.             if (sum == 0) {
  127.                 // Zero-norm vector is discarded.
  128.                 return sample();
  129.             }

  130.             final double f = 1.0 / Math.sqrt(sum);
  131.             return new double[] {x * f, y * f, z * f};
  132.         }

  133.         @Override
  134.         public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
  135.             return new UnitSphereSampler3D(rng);
  136.         }
  137.     }

  138.     /**
  139.      * Sample uniformly from a ND unit sphere.
  140.      */
  141.     private static final class UnitSphereSamplerND extends UnitSphereSampler {
  142.         /** Space dimension. */
  143.         private final int dimension;
  144.         /** Sampler used for generating the individual components of the vectors. */
  145.         private final NormalizedGaussianSampler sampler;

  146.         /**
  147.          * @param rng Source of randomness.
  148.          * @param dimension Space dimension.
  149.          */
  150.         UnitSphereSamplerND(UniformRandomProvider rng, int dimension) {
  151.             this.dimension = dimension;
  152.             sampler = ZigguratSampler.NormalizedGaussian.of(rng);
  153.         }

  154.         @Override
  155.         public double[] sample() {
  156.             final double[] v = new double[dimension];

  157.             // Pick a point by choosing a standard Gaussian for each element,
  158.             // and then normalize to unit length.
  159.             double sum = 0;
  160.             for (int i = 0; i < dimension; i++) {
  161.                 final double x = sampler.sample();
  162.                 v[i] = x;
  163.                 sum += x * x;
  164.             }

  165.             if (sum == 0) {
  166.                 // Zero-norm vector is discarded.
  167.                 // Using recursion as it is highly unlikely to generate more
  168.                 // than a few such vectors. It also protects against infinite
  169.                 // loop (in case a buggy generator is used), by eventually
  170.                 // raising a "StackOverflowError".
  171.                 return sample();
  172.             }

  173.             final double f = 1 / Math.sqrt(sum);
  174.             for (int i = 0; i < dimension; i++) {
  175.                 v[i] *= f;
  176.             }

  177.             return v;
  178.         }

  179.         @Override
  180.         public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
  181.             return new UnitSphereSamplerND(rng, dimension);
  182.         }
  183.     }

  184.     /**
  185.      * This instance delegates sampling. Use the factory method
  186.      * {@link #of(UniformRandomProvider, int)} to create an optimal sampler.
  187.      *
  188.      * @param dimension Space dimension.
  189.      * @param rng Generator for the individual components of the vectors.
  190.      * A shallow copy will be stored in this instance.
  191.      * @throws IllegalArgumentException If {@code dimension <= 0}
  192.      * @deprecated Use {@link #of(UniformRandomProvider, int)}.
  193.      */
  194.     @Deprecated
  195.     public UnitSphereSampler(int dimension,
  196.                              UniformRandomProvider rng) {
  197.         this(of(rng, dimension));
  198.     }

  199.     /**
  200.      * Private constructor used by deprecated constructor used to prevent partially
  201.      * initialized object if the construction of the delegate throws.
  202.      * In future versions the public constructor should be removed and the class made abstract.
  203.      *
  204.      * @param delegate Delegate.
  205.      */
  206.     private UnitSphereSampler(UnitSphereSampler delegate) {
  207.         this.delegate = delegate;
  208.     }

  209.     /**
  210.      * Private constructor used by sub-class specialisations.
  211.      * In future versions the public constructor should be removed and the class made abstract.
  212.      */
  213.     private UnitSphereSampler() {
  214.         delegate = null;
  215.     }

  216.     /**
  217.      * @return a random normalized Cartesian vector.
  218.      * @since 1.4
  219.      */
  220.     @Override
  221.     public double[] sample() {
  222.         return delegate.sample();
  223.     }

  224.     /**
  225.      * @return a random normalized Cartesian vector.
  226.      * @deprecated Use {@link #sample()}.
  227.      */
  228.     @Deprecated
  229.     public double[] nextVector() {
  230.         return sample();
  231.     }

  232.     /**
  233.      * {@inheritDoc}
  234.      *
  235.      * @since 1.3
  236.      */
  237.     @Override
  238.     public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) {
  239.         return delegate.withUniformRandomProvider(rng);
  240.     }

  241.     /**
  242.      * Create a unit sphere sampler for the given dimension.
  243.      *
  244.      * @param rng Generator for the individual components of the vectors. A shallow
  245.      * copy will be stored in the sampler.
  246.      * @param dimension Space dimension.
  247.      * @return the sampler
  248.      * @throws IllegalArgumentException If {@code dimension <= 0}
  249.      *
  250.      * @since 1.4
  251.      */
  252.     public static UnitSphereSampler of(UniformRandomProvider rng,
  253.                                        int dimension) {
  254.         if (dimension <= 0) {
  255.             throw new IllegalArgumentException("Dimension must be strictly positive");
  256.         } else if (dimension == ONE_D) {
  257.             return new UnitSphereSampler1D(rng);
  258.         } else if (dimension == TWO_D) {
  259.             return new UnitSphereSampler2D(rng);
  260.         } else if (dimension == THREE_D) {
  261.             return new UnitSphereSampler3D(rng);
  262.         }
  263.         return new UnitSphereSamplerND(rng, dimension);
  264.     }
  265. }