This is an automated email from the ASF dual-hosted git repository. aherbert pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-rng.git
commit 5692fe119d9f5866e4448a47472ba20a4ea7a813 Author: Alex Herbert <[email protected]> AuthorDate: Mon Apr 12 13:39:08 2021 +0100 [RNG-129:130] Improve performance of UnitSphereSampler RNG-129: - Add a benchmark for specialisations for low dimensions. - Implement specialistions for 1D, 2D and 3D. - Use a delegate in a publicly constructed UnitSphereSampler. RNG-130: - Fix 1 dimension sampling to only return vectors containing 1 or -1. --- .../jmh/sampling/UnitSphereSamplerBenchmark.java | 446 ++++++++++++++++++ .../commons/rng/sampling/UnitSphereSampler.java | 255 ++++++++-- .../rng/sampling/UnitSphereSamplerTest.java | 521 ++++++++++++++++++--- src/changes/changes.xml | 7 + 4 files changed, 1129 insertions(+), 100 deletions(-) diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java new file mode 100644 index 0000000..7aae519 --- /dev/null +++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java @@ -0,0 +1,446 @@ +/* + * 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.examples.jmh.sampling; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler; +import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; +import org.apache.commons.rng.simple.RandomSource; +import org.openjdk.jmh.annotations.Benchmark; +import org.openjdk.jmh.annotations.BenchmarkMode; +import org.openjdk.jmh.annotations.Fork; +import org.openjdk.jmh.annotations.Measurement; +import org.openjdk.jmh.annotations.Mode; +import org.openjdk.jmh.annotations.OutputTimeUnit; +import org.openjdk.jmh.annotations.Param; +import org.openjdk.jmh.annotations.Scope; +import org.openjdk.jmh.annotations.Setup; +import org.openjdk.jmh.annotations.State; +import org.openjdk.jmh.annotations.Warmup; +import org.openjdk.jmh.infra.Blackhole; +import java.util.concurrent.TimeUnit; + +/** + * Executes benchmark to compare the speed of generating samples on the surface of an + * N-dimension unit sphere. + */ +@BenchmarkMode(Mode.AverageTime) +@OutputTimeUnit(TimeUnit.NANOSECONDS) +@Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS) +@Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS) +@State(Scope.Benchmark) +@Fork(value = 1, jvmArgs = { "-server", "-Xms128M", "-Xmx128M" }) +public class UnitSphereSamplerBenchmark { + /** Name for the baseline method. */ + private static final String BASELINE = "Baseline"; + /** Name for the non-array based method. */ + private static final String NON_ARRAY = "NonArray"; + /** Name for the array based method. */ + private static final String ARRAY = "Array"; + /** Error message for an unknown sampler type. */ + private static final String UNKNOWN_SAMPLER = "Unknown sampler type: "; + + /** + * The sampler. + */ + private interface Sampler { + /** + * Gets the next sample. + * + * @return the sample + */ + double[] sample(); + } + + /** + * Base class for the sampler data. + * Contains the source of randomness and the number of samples. + * The sampler should be created by a sub-class of the data. + */ + @State(Scope.Benchmark) + public abstract static class SamplerData { + /** The sampler. */ + private Sampler sampler; + + /** The number of samples. */ + @Param({"100"}) + private int size; + + /** + * Gets the size. + * + * @return the size + */ + public int getSize() { + return size; + } + + /** + * Gets the sampler. + * + * @return the sampler + */ + public Sampler getSampler() { + return sampler; + } + + /** + * Create the source of randomness. + */ + @Setup + public void setup() { + // This could be configured using @Param + final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP); + sampler = createSampler(rng); + } + + /** + * Creates the sampler. + * + * @param rng the source of randomness + * @return the sampler + */ + protected abstract Sampler createSampler(UniformRandomProvider rng); + } + + /** + * The 1D unit line sampler. + */ + @State(Scope.Benchmark) + public static class Sampler1D extends SamplerData { + /** Name for the signed double method. */ + private static final String SIGNED_DOUBLE = "signedDouble"; + /** Name for the boolean method. */ + private static final String BOOLEAN = "boolean"; + + /** The sampler type. */ + @Param({BASELINE, SIGNED_DOUBLE, BOOLEAN, ARRAY}) + private String type; + + /** {@inheritDoc} */ + @Override + protected Sampler createSampler(final UniformRandomProvider rng) { + if (BASELINE.equals(type)) { + return new Sampler() { + @Override + public double[] sample() { + return new double[] {1.0}; + } + }; + } else if (SIGNED_DOUBLE.equals(type)) { + return new Sampler() { + @Override + public double[] sample() { + // (1 - 0) or (1 - 2) + // Use the sign bit + return new double[] {1.0 - ((rng.nextInt() >>> 30) & 0x2)}; + } + }; + } else if (BOOLEAN.equals(type)) { + return new Sampler() { + @Override + public double[] sample() { + return new double[] {rng.nextBoolean() ? -1.0 : 1.0}; + } + }; + } else if (ARRAY.equals(type)) { + return new ArrayBasedUnitSphereSampler(1, rng); + } + throw new IllegalStateException(UNKNOWN_SAMPLER + type); + } + } + + /** + * The 2D unit circle sampler. + */ + @State(Scope.Benchmark) + public static class Sampler2D extends SamplerData { + /** The sampler type. */ + @Param({BASELINE, ARRAY, NON_ARRAY}) + private String type; + + /** {@inheritDoc} */ + @Override + protected Sampler createSampler(final UniformRandomProvider rng) { + if (BASELINE.equals(type)) { + return new Sampler() { + @Override + public double[] sample() { + return new double[] {1.0, 0.0}; + } + }; + } else if (ARRAY.equals(type)) { + return new ArrayBasedUnitSphereSampler(2, rng); + } else if (NON_ARRAY.equals(type)) { + return new UnitSphereSampler2D(rng); + } + throw new IllegalStateException(UNKNOWN_SAMPLER + type); + } + + /** + * Sample from a 2D unit sphere. + */ + private static class UnitSphereSampler2D implements Sampler { + /** Sampler used for generating the individual components of the vectors. */ + private final NormalizedGaussianSampler sampler; + + /** + * @param rng the source of randomness + */ + UnitSphereSampler2D(UniformRandomProvider rng) { + sampler = new ZigguratNormalizedGaussianSampler(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}; + } + } + } + + /** + * The 3D unit sphere sampler. + */ + @State(Scope.Benchmark) + public static class Sampler3D extends SamplerData { + /** The sampler type. */ + @Param({BASELINE, ARRAY, NON_ARRAY}) + private String type; + + /** {@inheritDoc} */ + @Override + protected Sampler createSampler(final UniformRandomProvider rng) { + if (BASELINE.equals(type)) { + return new Sampler() { + @Override + public double[] sample() { + return new double[] {1.0, 0.0, 0.0}; + } + }; + } else if (ARRAY.equals(type)) { + return new ArrayBasedUnitSphereSampler(3, rng); + } else if (NON_ARRAY.equals(type)) { + return new UnitSphereSampler3D(rng); + } + throw new IllegalStateException(UNKNOWN_SAMPLER + type); + } + + /** + * Sample from a 3D unit sphere. + */ + private static class UnitSphereSampler3D implements Sampler { + /** Sampler used for generating the individual components of the vectors. */ + private final NormalizedGaussianSampler sampler; + + /** + * @param rng the source of randomness + */ + UnitSphereSampler3D(UniformRandomProvider rng) { + sampler = new ZigguratNormalizedGaussianSampler(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}; + } + } + } + + /** + * The 4D unit hypersphere sampler. + */ + @State(Scope.Benchmark) + public static class Sampler4D extends SamplerData { + /** The sampler type. */ + @Param({BASELINE, ARRAY, NON_ARRAY}) + private String type; + + /** {@inheritDoc} */ + @Override + protected Sampler createSampler(final UniformRandomProvider rng) { + if (BASELINE.equals(type)) { + return new Sampler() { + @Override + public double[] sample() { + return new double[] {1.0, 0.0, 0.0, 0.0}; + } + }; + } else if (ARRAY.equals(type)) { + return new ArrayBasedUnitSphereSampler(4, rng); + } else if (NON_ARRAY.equals(type)) { + return new UnitSphereSampler4D(rng); + } + throw new IllegalStateException(UNKNOWN_SAMPLER + type); + } + + /** + * Sample from a 4D unit hypersphere. + */ + private static class UnitSphereSampler4D implements Sampler { + /** Sampler used for generating the individual components of the vectors. */ + private final NormalizedGaussianSampler sampler; + + /** + * @param rng the source of randomness + */ + UnitSphereSampler4D(UniformRandomProvider rng) { + sampler = new ZigguratNormalizedGaussianSampler(rng); + } + + @Override + public double[] sample() { + final double x = sampler.sample(); + final double y = sampler.sample(); + final double z = sampler.sample(); + final double a = sampler.sample(); + final double sum = x * x + y * y + z * z + a * a; + + 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, a * f}; + } + } + } + + /** + * Sample from a unit sphere using an array based method. + */ + private static class ArrayBasedUnitSphereSampler implements Sampler { + /** Space dimension. */ + private final int dimension; + /** Sampler used for generating the individual components of the vectors. */ + private final NormalizedGaussianSampler sampler; + + /** + * @param dimension space dimension + * @param rng the source of randomness + */ + ArrayBasedUnitSphereSampler(int dimension, UniformRandomProvider rng) { + this.dimension = dimension; + sampler = new ZigguratNormalizedGaussianSampler(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. + return sample(); + } + + final double f = 1 / Math.sqrt(sum); + for (int i = 0; i < dimension; i++) { + v[i] *= f; + } + + return v; + } + } + + /** + * Run the sampler for the configured number of samples. + * + * @param bh Data sink + * @param data Input data. + */ + private static void runSampler(Blackhole bh, SamplerData data) { + final Sampler sampler = data.getSampler(); + for (int i = data.getSize() - 1; i >= 0; i--) { + bh.consume(sampler.sample()); + } + } + + /** + * Generation of uniform samples on a 1D unit line. + * + * @param bh Data sink + * @param data Input data. + */ + @Benchmark + public void create1D(Blackhole bh, Sampler1D data) { + runSampler(bh, data); + } + + /** + * Generation of uniform samples from a 2D unit circle. + * + * @param bh Data sink + * @param data Input data. + */ + @Benchmark + public void create2D(Blackhole bh, Sampler2D data) { + runSampler(bh, data); + } + + /** + * Generation of uniform samples from a 3D unit sphere. + * + * @param bh Data sink + * @param data Input data. + */ + @Benchmark + public void create3D(Blackhole bh, Sampler3D data) { + runSampler(bh, data); + } + + /** + * Generation of uniform samples from a 4D unit sphere. + * + * @param bh Data sink + * @param data Input data. + */ + @Benchmark + public void create4D(Blackhole bh, Sampler4D data) { + runSampler(bh, data); + } +} diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java index 6ce8b74..eadaf9b 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java @@ -25,78 +25,218 @@ import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSa * Generate vectors <a href="http://mathworld.wolfram.com/SpherePointPicking.html";> * isotropically located on the surface of a sphere</a>. * - * <p>Sampling uses:</p> + * <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 SharedStateSampler<UnitSphereSampler> { - /** Sampler used for generating the individual components of the vectors. */ - private final NormalizedGaussianSampler sampler; - /** Space dimension. */ - private final int dimension; + /** 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; /** - * @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} + * Sample uniformly from the ends of a 1D unit line. */ - public UnitSphereSampler(int dimension, - UniformRandomProvider rng) { - if (dimension <= 0) { - throw new IllegalArgumentException("Dimension must be strictly positive"); + private static 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[] nextVector() { + // Either: + // 1 - 0 = 1 + // 1 - 2 = -1 + // Use the sign bit + return new double[] {1.0 - ((rng.nextInt() >>> 30) & MASK_BIT_2)}; } - this.dimension = dimension; - sampler = new ZigguratNormalizedGaussianSampler(rng); + @Override + public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { + return new UnitSphereSampler1D(rng); + } } /** - * @param rng Generator for the individual components of the vectors. - * @param source Source to copy. + * Sample uniformly from a 2D unit circle. + * This is a 2D specialisation of the UnitSphereSamplerND. */ - private UnitSphereSampler(UniformRandomProvider rng, - UnitSphereSampler source) { - // The Gaussian sampler has no shared state so create a new instance - sampler = new ZigguratNormalizedGaussianSampler(rng); - dimension = source.dimension; + private static 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 = new ZigguratNormalizedGaussianSampler(rng); + } + + @Override + public double[] nextVector() { + 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 nextVector(); + } + + 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); + } } /** - * @return a random normalized Cartesian vector. + * Sample uniformly from a 3D unit sphere. + * This is a 3D specialisation of the UnitSphereSamplerND. */ - public double[] nextVector() { - final double[] v = new double[dimension]; - - // Pick a point by choosing a standard Gaussian for each element, - // and then normalize to unit length. - double normSq = 0; - for (int i = 0; i < dimension; i++) { - final double comp = sampler.sample(); - v[i] = comp; - normSq += comp * comp; + private static 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 = new ZigguratNormalizedGaussianSampler(rng); } - if (normSq == 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 nextVector(); + @Override + public double[] nextVector() { + 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 nextVector(); + } + + final double f = 1.0 / Math.sqrt(sum); + return new double[] {x * f, y * f, z * f}; } - final double f = 1 / Math.sqrt(normSq); - for (int i = 0; i < dimension; i++) { - v[i] *= f; + @Override + public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { + return new UnitSphereSampler3D(rng); } + } + + /** + * Sample uniformly from a ND unit sphere. + */ + private static 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 dimension Space dimension. + * @param rng Source of randomness. + */ + UnitSphereSamplerND(int dimension, UniformRandomProvider rng) { + this.dimension = dimension; + sampler = new ZigguratNormalizedGaussianSampler(rng); + } + + @Override + public double[] nextVector() { + 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; + } - return v; + 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 nextVector(); + } + + 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(dimension, rng); + } + } + + /** + * This instance delegates sampling. Use the factory method + * {@link #of(int, UniformRandomProvider)} 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} + */ + public UnitSphereSampler(int dimension, + UniformRandomProvider rng) { + delegate = of(dimension, rng); + } + + /** + * 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. + */ + public double[] nextVector() { + return delegate.nextVector(); } /** @@ -106,6 +246,31 @@ public class UnitSphereSampler implements SharedStateSampler<UnitSphereSampler> */ @Override public UnitSphereSampler withUniformRandomProvider(UniformRandomProvider rng) { - return new UnitSphereSampler(rng, this); + return delegate.withUniformRandomProvider(rng); + } + + /** + * Create a unit sphere sampler for the given dimension. + * + * @param dimension Space dimension. + * @param rng Generator for the individual components of the vectors. A shallow + * copy will be stored in the sampler. + * @return the sampler + * @throws IllegalArgumentException If {@code dimension <= 0} + * + * @since 1.4 + */ + public static UnitSphereSampler of(int dimension, + UniformRandomProvider rng) { + 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(dimension, rng); } } diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/UnitSphereSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/UnitSphereSamplerTest.java index 3c9df0b..efd1a6d 100644 --- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/UnitSphereSamplerTest.java +++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/UnitSphereSamplerTest.java @@ -19,6 +19,8 @@ package org.apache.commons.rng.sampling; import org.junit.Assert; import org.junit.Test; import org.apache.commons.rng.simple.RandomSource; +import java.util.Arrays; +import org.apache.commons.math3.stat.inference.ChiSquareTest; import org.apache.commons.rng.UniformRandomProvider; import org.apache.commons.rng.core.source64.SplitMix64; @@ -26,77 +28,408 @@ import org.apache.commons.rng.core.source64.SplitMix64; * Test for {@link UnitSphereSampler}. */ public class UnitSphereSamplerTest { + /** 2 pi */ + private static final double TWO_PI = 2 * Math.PI; + + /** + * Test a non-positive dimension. + */ @Test(expected = IllegalArgumentException.class) - public void testPrecondition() { + public void testInvalidDimensionThrows() { new UnitSphereSampler(0, null); } /** + * Test a non-positive dimension. + */ + @Test(expected = IllegalArgumentException.class) + public void testInvalidDimensionThrowsWithFactoryConstructor() { + UnitSphereSampler.of(0, null); + } + + /** + * Test the distribution of points in one dimension. + */ + @Test + public void testDistribution1D() { + testDistribution1D(false); + } + + /** + * Test the distribution of points in one dimension with the factory constructor. + */ + @Test + public void testDistribution1DWithFactoryConstructor() { + testDistribution1D(true); + } + + /** + * Test the distribution of points in one dimension. + * RNG-130: All samples should be 1 or -1. + */ + private static void testDistribution1D(boolean factoryConstructor) { + final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP, 0x1a2b3cL); + final UnitSphereSampler generator = createUnitSphereSampler(1, rng, factoryConstructor); + final int samples = 10000; + // Count the negatives. + int count = 0; + for (int i = 0; i < samples; i++) { + final double[] v = generator.nextVector(); + Assert.assertEquals(1, v.length); + final double d = v[0]; + if (d == -1.0) { + count++; + } else if (d != 1.0) { + // RNG-130: All samples should be 1 or -1. + Assert.fail("Invalid unit length: " + d); + } + } + // Test the number of negatives is approximately 50% + assertMonobit(count, samples); + } + + /** + * Assert that the number of 1 bits is approximately 50%. This is based upon a fixed-step "random + * walk" of +1/-1 from zero. + * + * <p>The test is equivalent to the NIST Monobit test with a fixed p-value of 0.01. The number of + * bits is recommended to be above 100.</p> + * + * @see <a href="https://csrc.nist.gov/publications/detail/sp/800-22/rev-1a/final";>Bassham, et al + * (2010) NIST SP 800-22: A Statistical Test Suite for Random and Pseudorandom Number + * Generators for Cryptographic Applications. Section 2.1.</a> + * + * @param bitCount The bit count. + * @param numberOfBits Number of bits. + */ + private static void assertMonobit(int bitCount, int numberOfBits) { + // Convert the bit count into a number of +1/-1 steps. + final double sum = 2.0 * bitCount - numberOfBits; + // The reference distribution is Normal with a standard deviation of sqrt(n). + // Check the absolute position is not too far from the mean of 0 with a fixed + // p-value of 0.01 taken from a 2-tailed Normal distribution. Computation of + // the p-value requires the complimentary error function. + final double absSum = Math.abs(sum); + final double max = Math.sqrt(numberOfBits) * 2.576; + Assert.assertTrue( + "Walked too far astray: " + absSum + " > " + max + + " (test will fail randomly about 1 in 100 times)", + absSum <= max); + } + + /** * Test the distribution of points in two dimensions. */ @Test public void testDistribution2D() { - UniformRandomProvider rng = RandomSource.create(RandomSource.XOR_SHIFT_1024_S, 17399225432L); - UnitSphereSampler generator = new UnitSphereSampler(2, rng); + testDistribution2D(false); + } + + /** + * Test the distribution of points in two dimensions with the factory constructor. + */ + @Test + public void testDistribution2DWithFactoryConstructor() { + testDistribution2D(true); + } + + /** + * Test the distribution of points in two dimensions. + * Obtains polar coordinates and checks the angle distribution is uniform. + */ + private static void testDistribution2D(boolean factoryConstructor) { + final UniformRandomProvider rng = RandomSource.create(RandomSource.XOR_SHIFT_1024_S, 17399225432L); + final UnitSphereSampler generator = createUnitSphereSampler(2, rng, factoryConstructor); // In 2D, angles with a given vector should be uniformly distributed. - final int[] angleBuckets = new int[100]; - final int steps = 1000000; + final int angleBins = 200; + final long[] observed = new long[angleBins]; + final int steps = 100000; for (int i = 0; i < steps; ++i) { final double[] v = generator.nextVector(); Assert.assertEquals(2, v.length); - Assert.assertEquals(1, length(v), 1e-10); - // Compute angle formed with vector (1, 0)? - // Cosine of angle is their dot product, because both are unit length. - // Dot product here is just the first element of the vector by construction. - final double angle = Math.acos(v[0]); - final int bucket = (int) (angle / Math.PI * angleBuckets.length); - ++angleBuckets[bucket]; + Assert.assertEquals(1.0, length(v), 1e-10); + // Get the polar angle bin from xy + final int angleBin = angleBin(angleBins, v[0], v[1]); + observed[angleBin]++; } - // Simplistic test for roughly even distribution. - final int expectedBucketSize = steps / angleBuckets.length; - for (int bucket : angleBuckets) { - Assert.assertTrue("Bucket count " + bucket + " vs expected " + expectedBucketSize, - Math.abs(expectedBucketSize - bucket) < 350); + final double[] expected = new double[observed.length]; + Arrays.fill(expected, (double) steps / observed.length); + final double p = new ChiSquareTest().chiSquareTest(expected, observed); + Assert.assertFalse("p-value too small: " + p, p < 0.01); + } + + /** + * Test the distribution of points in three dimensions. + */ + @Test + public void testDistribution3D() { + testDistribution3D(false); + } + + /** + * Test the distribution of points in three dimensions with the factory constructor. + */ + @Test + public void testDistribution3DWithFactoryConstructor() { + testDistribution3D(true); + } + + /** + * Test the distribution of points in three dimensions. + * Obtains spherical coordinates and checks the distribution is uniform. + */ + private static void testDistribution3D(boolean factoryConstructor) { + final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_SHI_RO_256_PP, 0xabcdefL); + final UnitSphereSampler generator = createUnitSphereSampler(3, rng, factoryConstructor); + + // Get 3D spherical coordinates. Assign to a bin. + // + // polar angle (theta) in [0, 2pi) + // azimuthal angle (phi) in [0, pi) + // + // theta = arctan(y/x) and is uniformly distributed + // phi = arccos(z); and cos(phi) is uniformly distributed + final int angleBins = 20; + final int depthBins = 10; + final long[] observed = new long[angleBins * depthBins]; + final int steps = 1000000; + for (int i = 0; i < steps; ++i) { + final double[] v = generator.nextVector(); + Assert.assertEquals(3, v.length); + Assert.assertEquals(1.0, length(v), 1e-10); + // Get the polar angle bin from xy + final int angleBin = angleBin(angleBins, v[0], v[1]); + // Map cos(phi) = z from [-1, 1) to [0, 1) then assign a bin + final int depthBin = (int) (depthBins * (v[2] + 1) / 2); + observed[depthBin * angleBins + angleBin]++; + } + + final double[] expected = new double[observed.length]; + Arrays.fill(expected, (double) steps / observed.length); + final double p = new ChiSquareTest().chiSquareTest(expected, observed); + Assert.assertFalse("p-value too small: " + p, p < 0.01); + } + + /** + * Test the distribution of points in four dimensions. + */ + @Test + public void testDistribution4D() { + testDistribution4D(false); + } + + /** + * Test the distribution of points in four dimensions with the factory constructor. + */ + @Test + public void testDistribution4DWithFactoryConstructor() { + testDistribution4D(true); + } + + /** + * Test the distribution of points in four dimensions. + * Checks the surface of the 3-sphere can be used to generate uniform samples within a circle. + */ + private static void testDistribution4D(boolean factoryConstructor) { + final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_SHI_RO_512_PP, 0x9876543210L); + final UnitSphereSampler generator = createUnitSphereSampler(4, rng, factoryConstructor); + + // No uniform distribution of spherical coordinates for a 3-sphere. + // https://en.wikipedia.org/wiki/N-sphere#Spherical_coordinates + // Here we exploit the fact that the uniform distribution of a (n+1)-sphere + // when discarding two coordinates is uniform within a n-ball. + // Thus any two coordinates of the 4 are uniform within a circle. + // Here we separately test pairs (0, 1) and (2, 3). + // Note: We cannot create a single bin from the two bins in each circle as they are + // not independent. A point close to the edge in one circle requires a point close to + // the centre in the other circle to create a unit radius from all 4 coordinates. + // This test exercises the N-dimension sampler and demonstrates the vectors obey + // properties of the (n+1)-sphere. + + // Divide the circle into layers of concentric rings and an angle. + final int layers = 10; + final int angleBins = 20; + + // Compute the radius for each layer to have the same area + // (i.e. incrementally larger concentric circles must increase area by a constant). + // r = sqrt(fraction * maxArea / pi) + // Unit circle has area pi so we just use sqrt(fraction). + final double[] r = new double[layers]; + for (int i = 1; i < layers; i++) { + r[i - 1] = Math.sqrt((double) i / layers); + } + // The final radius should be 1.0. + r[layers - 1] = 1.0; + + final long[] observed1 = new long[layers * angleBins]; + final long[] observed2 = new long[observed1.length]; + final int steps = 1000000; + for (int i = 0; i < steps; ++i) { + final double[] v = generator.nextVector(); + Assert.assertEquals(4, v.length); + Assert.assertEquals(1.0, length(v), 1e-10); + // Circle 1 + observed1[circleBin(angleBins, r, v[0], v[1])]++; + // Circle 2 + observed2[circleBin(angleBins, r, v[2], v[3])]++; + } + + final double[] expected = new double[observed1.length]; + Arrays.fill(expected, (double) steps / observed1.length); + final ChiSquareTest chi = new ChiSquareTest(); + final double p1 = chi.chiSquareTest(expected, observed1); + Assert.assertFalse("Circle 1 p-value too small: " + p1, p1 < 0.01); + final double p2 = chi.chiSquareTest(expected, observed2); + Assert.assertFalse("Circle 2 p-value too small: " + p2, p2 < 0.01); + } + + /** + * Compute a bin inside the circle using the polar angle theta and the radius thresholds. + * + * @param angleBins the number of angle bins + * @param r the radius bin thresholds (ascending order) + * @param x the x coordinate + * @param y the y coordinate + * @return the bin + */ + private static int circleBin(int angleBins, double[] r, double x, double y) { + final int angleBin = angleBin(angleBins, x, y); + final int radiusBin = radiusBin(r, x, y); + return radiusBin * angleBins + angleBin; + } + + /** + * Compute an angle bin from the xy vector. The bin will represent the range [0, 2pi). + * + * @param angleBins the number of angle bins + * @param x the x coordinate + * @param y the y coordinate + * @return the bin + */ + private static int angleBin(int angleBins, double x, double y) { + final double angle = Math.atan2(y, x); + // Map [-pi, pi) to [0, 1) then assign a bin + return (int) (angleBins * (angle + Math.PI) / TWO_PI); + } + + /** + * Compute a radius bin from the xy vector. The bin is assigned if the length of the vector + * is above the threshold of the bin. + * + * @param r the radius bin thresholds (ascending order) + * @param x the x coordinate + * @param y the y coordinate + * @return the bin + */ + private static int radiusBin(double[] r, double x, double y) { + final double length = Math.sqrt(x * x + y * y); + + // Note: The bin should be uniformly distributed. + // A test using a custom binary search (avoiding double NaN checks) + // shows the simple loop over a small number of bins is comparable in speed. + // The loop is preferred for simplicity. A binary search may be better + // if the number of bins is increased. + for (int layer = 0; layer < r.length; layer++) { + if (length <= r[layer]) { + return layer; + } } + // Unreachable if the xy component is from a vector of length <= 1 + throw new AssertionError("Invalid sample length: " + length); } - /** Cf. RNG-55. */ + /** + * Test infinite recursion occurs with a bad provider in 2D. + */ @Test(expected = StackOverflowError.class) - public void testBadProvider1() { - final UniformRandomProvider bad = new UniformRandomProvider() { - // CHECKSTYLE: stop all - public long nextLong(long n) { return 0; } - public long nextLong() { return 0; } - public int nextInt(int n) { return 0; } - public int nextInt() { return 0; } - public float nextFloat() { return 0; } - public double nextDouble() { return 0;} - public void nextBytes(byte[] bytes, int start, int len) {} - public void nextBytes(byte[] bytes) {} - public boolean nextBoolean() { return false; } - // CHECKSTYLE: resume all - }; - - new UnitSphereSampler(1, bad).nextVector(); - } - - /** Cf. RNG-55. */ - @Test - public void testBadProvider1ThenGoodProvider() { + public void testBadProvider2D() { + testBadProvider(2); + } + + /** + * Test infinite recursion occurs with a bad provider in 3D. + */ + @Test(expected = StackOverflowError.class) + public void testBadProvider3D() { + testBadProvider(3); + } + + /** + * Test infinite recursion occurs with a bad provider in 4D. + */ + @Test(expected = StackOverflowError.class) + public void testBadProvider4D() { + testBadProvider(4); + } + + /** + * Test the edge case where the normalisation sum to divide by is always zero. + * This test requires generation of Gaussian samples with the value 0. + * The sample eventually fails due to infinite recursion. + * See RNG-55. + * + * @param dimension the dimension + */ + private static void testBadProvider(final int dimension) { + // A provider that will create zero valued Gaussian samples + // from the ZigguratNormalizedGaussianSampler. + final UniformRandomProvider bad = new SplitMix64(0L) { + @Override + public long nextLong() { + return 0; + } + }; + + UnitSphereSampler.of(dimension, bad).nextVector(); + } + + /** + * Test the edge case where the normalisation sum to divide by is zero for 2D. + */ + @Test + public void testInvalidInverseNormalisation2D() { + testInvalidInverseNormalisationND(2); + } + + /** + * Test the edge case where the normalisation sum to divide by is zero for 3D. + */ + @Test + public void testInvalidInverseNormalisation3D() { + testInvalidInverseNormalisationND(3); + } + + /** + * Test the edge case where the normalisation sum to divide by is zero for 4D. + */ + @Test + public void testInvalidInverseNormalisation4D() { + testInvalidInverseNormalisationND(4); + } + + /** + * Test the edge case where the normalisation sum to divide by is zero. + * This test requires generation of Gaussian samples with the value 0. + * See RNG-55. + */ + private static void testInvalidInverseNormalisationND(final int dimension) { // Create a provider that will create a bad first sample but then recover. // This checks recursion will return a good value. - final UniformRandomProvider bad = new SplitMix64(0L) { - private int count; - // CHECKSTYLE: stop all - public long nextLong() { return (count++ == 0) ? 0 : super.nextLong(); } - public double nextDouble() { return (count++ == 0) ? 0 : super.nextDouble(); } - // CHECKSTYLE: resume all - }; + final UniformRandomProvider bad = new SplitMix64(0x1a2b3cL) { + private int count = -2 * dimension; - final double[] vector = new UnitSphereSampler(1, bad).nextVector(); - Assert.assertEquals(1, vector.length); + @Override + public long nextLong() { + // Return enough zeros to create Gaussian samples of zero for all coordinates. + return count++ < 0 ? 0 : super.nextLong(); + } + }; + + final double[] vector = UnitSphereSampler.of(dimension, bad).nextVector(); + Assert.assertEquals(dimension, vector.length); + Assert.assertEquals(1.0, length(vector), 1e-10); } /** @@ -107,7 +440,7 @@ public class UnitSphereSamplerTest { public void testNextNormSquaredAfterZeroIsValid() { // The sampler explicitly handles length == 0 using recursion. // Anything above zero should be valid. - final double normSq = Math.nextAfter(0, 1); + final double normSq = Math.nextAfter(0.0, 1); // Map to the scaling factor final double f = 1 / Math.sqrt(normSq); // As long as this is finite positive then the sampler is valid @@ -115,15 +448,79 @@ public class UnitSphereSamplerTest { } /** - * Test the SharedStateSampler implementation. + * Test the SharedStateSampler implementation for 1D. + */ + @Test + public void testSharedStateSampler1D() { + testSharedStateSampler(1, false); + } + + /** + * Test the SharedStateSampler implementation for 2D. + */ + @Test + public void testSharedStateSampler2D() { + testSharedStateSampler(2, false); + } + + /** + * Test the SharedStateSampler implementation for 3D. + */ + @Test + public void testSharedStateSampler3D() { + testSharedStateSampler(3, false); + } + + /** + * Test the SharedStateSampler implementation for 4D. + */ + @Test + public void testSharedStateSampler4D() { + testSharedStateSampler(4, false); + } + + /** + * Test the SharedStateSampler implementation for 1D using the factory constructor. + */ + @Test + public void testSharedStateSampler1DWithFactoryConstructor() { + testSharedStateSampler(1, true); + } + + /** + * Test the SharedStateSampler implementation for 2D using the factory constructor. + */ + @Test + public void testSharedStateSampler2DWithFactoryConstructor() { + testSharedStateSampler(2, true); + } + + /** + * Test the SharedStateSampler implementation for 3D using the factory constructor. + */ + @Test + public void testSharedStateSampler3DWithFactoryConstructor() { + testSharedStateSampler(3, true); + } + + /** + * Test the SharedStateSampler implementation for 4D using the factory constructor. */ @Test - public void testSharedStateSampler() { + public void testSharedStateSampler4DWithFactoryConstructor() { + testSharedStateSampler(4, true); + } + + /** + * Test the SharedStateSampler implementation for the given dimension. + * + * @param dimension the dimension + * @param factoryConstructor true to use the factory constructor + */ + private static void testSharedStateSampler(int dimension, boolean factoryConstructor) { final UniformRandomProvider rng1 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); final UniformRandomProvider rng2 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); - final int n = 3; - final UnitSphereSampler sampler1 = - new UnitSphereSampler(n, rng1); + final UnitSphereSampler sampler1 = createUnitSphereSampler(dimension, rng1, factoryConstructor); final UnitSphereSampler sampler2 = sampler1.withUniformRandomProvider(rng2); RandomAssert.assertProduceSameSequence( new RandomAssert.Sampler<double[]>() { @@ -141,11 +538,25 @@ public class UnitSphereSamplerTest { } /** + * Creates a UnitSphereSampler. + * + * @param dimension the dimension + * @param rng the source of randomness + * @param factoryConstructor true to use the factory constructor + * @return the sampler + */ + private static UnitSphereSampler createUnitSphereSampler(int dimension, UniformRandomProvider rng, + boolean factoryConstructor) { + return factoryConstructor ? + UnitSphereSampler.of(dimension, rng) : new UnitSphereSampler(dimension, rng); + } + + /** * @return the length (L2-norm) of given vector. */ private static double length(double[] vector) { double total = 0; - for (double d : vector) { + for (final double d : vector) { total += d * d; } return Math.sqrt(total); diff --git a/src/changes/changes.xml b/src/changes/changes.xml index a8bbd88..2d24b13 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -75,6 +75,13 @@ re-run tests that fail, and pass the build if they succeed within the allotted number of reruns (the test will be marked as 'flaky' in the report). "> + <action dev="aherbert" type="fix" issue="130"> + "UnitSphereSampler": Fix 1 dimension sampling to only return vectors containing 1 or -1. + </action> + <action dev="aherbert" type="update" issue="129"> + "UnitSphereSampler": Improve performance with specialisations for low order dimensions. + Added a factory constructor to create the sampler. + </action> <action dev="aherbert" type="add" issue="128"> New "UnitBallSampler" to generate coordinates uniformly within an n-unit ball. </action>
