You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/04/15 22:00:58 UTC
[commons-rng] branch master updated (c3f7dd5 -> a089841)
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a change to branch master
in repository https://gitbox.apache.org/repos/asf/commons-rng.git.
from c3f7dd5 [RNG-128] Add UnitBallSampler
new 5692fe1 [RNG-129:130] Improve performance of UnitSphereSampler
new a089841 UnitBallSampler: Use the sign bit for the sign of the range [-1, 1)
The 2 revisions listed above as "new" are entirely new to this
repository and will be described in separate emails. The revisions
listed as "add" were already present in the repository and have only
been added to this reference.
Summary of changes:
.../jmh/sampling/UnitBallSamplerBenchmark.java | 14 +-
.../jmh/sampling/UnitSphereSamplerBenchmark.java | 446 ++++++++++++++++++
.../commons/rng/sampling/UnitBallSampler.java | 11 +-
.../commons/rng/sampling/UnitSphereSampler.java | 255 ++++++++--
.../rng/sampling/UnitSphereSamplerTest.java | 521 ++++++++++++++++++---
src/changes/changes.xml | 7 +
6 files changed, 1144 insertions(+), 110 deletions(-)
create mode 100644 commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java
[commons-rng] 01/02: [RNG-129:130] Improve performance of
UnitSphereSampler
Posted by ah...@apache.org.
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 <ah...@apache.org>
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>
[commons-rng] 02/02: UnitBallSampler: Use the sign bit for the sign
of the range [-1, 1)
Posted by ah...@apache.org.
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 a089841b408de17f19a0dea43f11af2b18a8a304
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Apr 15 23:00:54 2021 +0100
UnitBallSampler: Use the sign bit for the sign of the range [-1, 1)
Changes the use of the upper 54 bits of the long from 53 bits for the
double then 1 bit for the sign to 1 bit for the sign then 53 bits for
the double.
---
.../examples/jmh/sampling/UnitBallSamplerBenchmark.java | 14 ++++++++------
.../org/apache/commons/rng/sampling/UnitBallSampler.java | 11 +++++++----
2 files changed, 15 insertions(+), 10 deletions(-)
diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitBallSamplerBenchmark.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitBallSamplerBenchmark.java
index 9a7f883..450f85c 100644
--- a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitBallSamplerBenchmark.java
+++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitBallSamplerBenchmark.java
@@ -46,7 +46,6 @@ import java.util.concurrent.TimeUnit;
@State(Scope.Benchmark)
@Fork(value = 1, jvmArgs = { "-server", "-Xms128M", "-Xmx128M" })
public class UnitBallSamplerBenchmark {
-
/** Name for the baseline method. */
private static final String BASELINE = "Baseline";
/** Name for the rejection method. */
@@ -61,6 +60,8 @@ public class UnitBallSamplerBenchmark {
private static final String HYPERSPHERE_DISCARD = "HypersphereDiscard";
/** Error message for an unknown sampler type. */
private static final String UNKNOWN_SAMPLER = "Unknown sampler type: ";
+ /** The mask to extract the lower 53-bits from a long. */
+ private static final long LOWER_53_BITS = -1L >>> 11;
/**
* The sampler.
@@ -738,10 +739,11 @@ public class UnitBallSamplerBenchmark {
* @return the double
*/
private static double makeSignedDouble(long bits) {
- // Upper 53-bits generates a positive number in the range [0, 1).
- // This has 1 optionally subtracted. Do not use the lower bits on the
- // assumption they are less random.
- return ((bits >>> 11) * 0x1.0p-53d) - ((bits >>> 10) & 0x1L);
+ // Use the upper 54 bits on the assumption they are more random.
+ // The sign bit generates a value of 0 or 1 for subtraction.
+ // The next 53 bits generates a positive number in the range [0, 1).
+ // [0, 1) - (0 or 1) => [-1, 1)
+ return (((bits >>> 10) & LOWER_53_BITS) * 0x1.0p-53d) - (bits >>> 63);
}
/**
@@ -763,7 +765,7 @@ public class UnitBallSamplerBenchmark {
* @param bh Data sink
* @param data Input data.
*/
- //@Benchmark
+ @Benchmark
public void create1D(Blackhole bh, Sampler1D data) {
runSampler(bh, data);
}
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitBallSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitBallSampler.java
index e0c2b69..4c359b5 100644
--- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitBallSampler.java
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitBallSampler.java
@@ -41,6 +41,8 @@ public abstract class UnitBallSampler implements SharedStateSampler<UnitBallSamp
private static final int TWO_D = 2;
/** The dimension for 3D sampling. */
private static final int THREE_D = 3;
+ /** The mask to extract the lower 53-bits from a long. */
+ private static final long LOWER_53_BITS = -1L >>> 11;
/**
* Sample uniformly from a 1D unit line.
@@ -228,9 +230,10 @@ public abstract class UnitBallSampler implements SharedStateSampler<UnitBallSamp
* @return the double
*/
private static double makeSignedDouble(long bits) {
- // Upper 53-bits generates a positive number in the range [0, 1).
- // This has 1 optionally subtracted. Do not use the lower bits on the
- // assumption they are less random.
- return ((bits >>> 11) * 0x1.0p-53d) - ((bits >>> 10) & 0x1L);
+ // Use the upper 54 bits on the assumption they are more random.
+ // The sign bit generates a value of 0 or 1 for subtraction.
+ // The next 53 bits generates a positive number in the range [0, 1).
+ // [0, 1) - (0 or 1) => [-1, 1)
+ return (((bits >>> 10) & LOWER_53_BITS) * 0x1.0p-53d) - (bits >>> 63);
}
}