UnitBallSampler.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.  */

  18. package org.apache.commons.rng.sampling.shape;

  20. import org.apache.commons.rng.UniformRandomProvider;
  21. import org.apache.commons.rng.sampling.SharedStateObjectSampler;
  22. import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
  23. import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
  24. import org.apache.commons.rng.sampling.distribution.ZigguratSampler;

  26. /**
  27.  * Generate coordinates <a href="http://mathworld.wolfram.com/BallPointPicking.html">
  28.  * uniformly distributed within the unit n-ball</a>.
  29.  *
  30.  * <p>Sampling uses:</p>
  31.  *
  32.  * <ul>
  33.  *   <li>{@link UniformRandomProvider#nextLong()}
  34.  *   <li>{@link UniformRandomProvider#nextDouble()} (only for dimensions above 2)
  35.  * </ul>
  36.  *
  37.  * @since 1.4
  38.  */
  39. public abstract class UnitBallSampler implements SharedStateObjectSampler<double[]> {
  40.     /** The dimension for 1D sampling. */
  41.     private static final int ONE_D = 1;
  42.     /** The dimension for 2D sampling. */
  43.     private static final int TWO_D = 2;
  44.     /** The dimension for 3D sampling. */
  45.     private static final int THREE_D = 3;
  46.     /**
  47.      * The multiplier to convert the least significant 53-bits of a {@code long} to a {@code double}.
  48.      * Taken from o.a.c.rng.core.utils.NumberFactory.
  49.      *
  50.      * <p>This is equivalent to {@code 1.0 / (1L << 53)}.
  51.      */
  52.     private static final double DOUBLE_MULTIPLIER = 0x1.0p-53d;

  54.     /**
  55.      * Sample uniformly from a 1D unit line.
  56.      */
  57.     private static final class UnitBallSampler1D extends UnitBallSampler {
  58.         /** The source of randomness. */
  59.         private final UniformRandomProvider rng;

  61.         /**
  62.          * @param rng Source of randomness.
  63.          */
  64.         UnitBallSampler1D(UniformRandomProvider rng) {
  65.             this.rng = rng;
  66.         }

  68.         @Override
  69.         public double[] sample() {
  70.             return new double[] {makeSignedDouble(rng.nextLong())};
  71.         }

  73.         @Override
  74.         public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
  75.             return new UnitBallSampler1D(rng);
  76.         }
  77.     }

  79.     /**
  80.      * Sample uniformly from a 2D unit disk.
  81.      */
  82.     private static final class UnitBallSampler2D extends UnitBallSampler {
  83.         /** The source of randomness. */
  84.         private final UniformRandomProvider rng;

  86.         /**
  87.          * @param rng Source of randomness.
  88.          */
  89.         UnitBallSampler2D(UniformRandomProvider rng) {
  90.             this.rng = rng;
  91.         }

  93.         @Override
  94.         public double[] sample() {
  95.             // Generate via rejection method of a circle inside a square of edge length 2.
  96.             // This should compute approximately 2^2 / pi = 1.27 square positions per sample.
  97.             double x;
  98.             double y;
  99.             do {
  100.                 x = makeSignedDouble(rng.nextLong());
  101.                 y = makeSignedDouble(rng.nextLong());
  102.             } while (x * x + y * y > 1.0);
  103.             return new double[] {x, y};
  104.         }

  106.         @Override
  107.         public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
  108.             return new UnitBallSampler2D(rng);
  109.         }
  110.     }

  112.     /**
  113.      * Sample uniformly from a 3D unit ball. This is an non-array based specialisation of
  114.      * {@link UnitBallSamplerND} for performance.
  115.      */
  116.     private static final class UnitBallSampler3D extends UnitBallSampler {
  117.         /** The standard normal distribution. */
  118.         private final NormalizedGaussianSampler normal;
  119.         /** The exponential distribution (mean=1). */
  120.         private final ContinuousSampler exp;

  122.         /**
  123.          * @param rng Source of randomness.
  124.          */
  125.         UnitBallSampler3D(UniformRandomProvider rng) {
  126.             normal = ZigguratSampler.NormalizedGaussian.of(rng);
  127.             // Require an Exponential(mean=2).
  128.             // Here we use mean = 1 and scale the output later.
  129.             exp = ZigguratSampler.Exponential.of(rng);
  130.         }

  132.         @Override
  133.         public double[] sample() {
  134.             final double x = normal.sample();
  135.             final double y = normal.sample();
  136.             final double z = normal.sample();
  137.             // Include the exponential sample. It has mean 1 so multiply by 2.
  138.             final double sum = exp.sample() * 2 + x * x + y * y + z * z;
  139.             // Note: Handle the possibility of a zero sum and invalid inverse
  140.             if (sum == 0) {
  141.                 return sample();
  142.             }
  143.             final double f = 1.0 / Math.sqrt(sum);
  144.             return new double[] {x * f, y * f, z * f};
  145.         }

  147.         @Override
  148.         public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
  149.             return new UnitBallSampler3D(rng);
  150.         }
  151.     }

  153.     /**
  154.      * Sample using ball point picking.
  155.      * @see <a href="https://mathworld.wolfram.com/BallPointPicking.html">Ball point picking</a>
  156.      */
  157.     private static final class UnitBallSamplerND extends UnitBallSampler {
  158.         /** The dimension. */
  159.         private final int dimension;
  160.         /** The standard normal distribution. */
  161.         private final NormalizedGaussianSampler normal;
  162.         /** The exponential distribution (mean=1). */
  163.         private final ContinuousSampler exp;

  165.         /**
  166.          * @param rng Source of randomness.
  167.          * @param dimension Space dimension.
  168.          */
  169.         UnitBallSamplerND(UniformRandomProvider rng, int dimension) {
  170.             this.dimension  = dimension;
  171.             normal = ZigguratSampler.NormalizedGaussian.of(rng);
  172.             // Require an Exponential(mean=2).
  173.             // Here we use mean = 1 and scale the output later.
  174.             exp = ZigguratSampler.Exponential.of(rng);
  175.         }

  177.         @Override
  178.         public double[] sample() {
  179.             final double[] sample = new double[dimension];
  180.             // Include the exponential sample. It has mean 1 so multiply by 2.
  181.             double sum = exp.sample() * 2;
  182.             for (int i = 0; i < dimension; i++) {
  183.                 final double x = normal.sample();
  184.                 sum += x * x;
  185.                 sample[i] = x;
  186.             }
  187.             // Note: Handle the possibility of a zero sum and invalid inverse
  188.             if (sum == 0) {
  189.                 return sample();
  190.             }
  191.             final double f = 1.0 / Math.sqrt(sum);
  192.             for (int i = 0; i < dimension; i++) {
  193.                 sample[i] *= f;
  194.             }
  195.             return sample;
  196.         }

  198.         @Override
  199.         public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
  200.             return new UnitBallSamplerND(rng, dimension);
  201.         }
  202.     }

  204.     /**
  205.      * Create an instance.
  206.      */
  207.     public UnitBallSampler() {}

  209.     /**
  210.      * @return a random Cartesian coordinate within the unit n-ball.
  211.      */
  212.     @Override
  213.     public abstract double[] sample();

  215.     /** {@inheritDoc} */
  216.     // Redeclare the signature to return a UnitBallSampler not a SharedStateObjectSampler<double[]>
  217.     @Override
  218.     public abstract UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng);

  220.     /**
  221.      * Create a unit n-ball sampler for the given dimension.
  222.      * Sampled points are uniformly distributed within the unit n-ball.
  223.      *
  224.      * <p>Sampling is supported in dimensions of 1 or above.
  225.      *
  226.      * @param rng Source of randomness.
  227.      * @param dimension Space dimension.
  228.      * @return the sampler
  229.      * @throws IllegalArgumentException If {@code dimension <= 0}
  230.      */
  231.     public static UnitBallSampler of(UniformRandomProvider rng,
  232.                                      int dimension) {
  233.         if (dimension <= 0) {
  234.             throw new IllegalArgumentException("Dimension must be strictly positive");
  235.         } else if (dimension == ONE_D) {
  236.             return new UnitBallSampler1D(rng);
  237.         } else if (dimension == TWO_D) {
  238.             return new UnitBallSampler2D(rng);
  239.         } else if (dimension == THREE_D) {
  240.             return new UnitBallSampler3D(rng);
  241.         }
  242.         return new UnitBallSamplerND(rng, dimension);
  243.     }

  245.     /**
  246.      * Creates a signed double in the range {@code [-1, 1)}. The magnitude is sampled evenly
  247.      * from the 2<sup>54</sup> dyadic rationals in the range.
  248.      *
  249.      * <p>Note: This method will not return samples for both -0.0 and 0.0.
  250.      *
  251.      * @param bits the bits
  252.      * @return the double
  253.      */
  254.     private static double makeSignedDouble(long bits) {
  255.         // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a signed
  256.         // shift of 10 in place of an unsigned shift of 11.
  257.         // Use the upper 54 bits on the assumption they are more random.
  258.         // The sign bit is maintained by the signed shift.
  259.         // The next 53 bits generates a magnitude in the range [0, 2^53) or [-2^53, 0).
  260.         return (bits >> 10) * DOUBLE_MULTIPLIER;
  261.     }
  262. }
Created with JaCoCo 0.8.12.202403310830