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/13 07:36:44 UTC
[commons-rng] branch master updated: [RNG-128] Add UnitBallSampler
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
The following commit(s) were added to refs/heads/master by this push:
new c3f7dd5 [RNG-128] Add UnitBallSampler
c3f7dd5 is described below
commit c3f7dd513bff7c22b104c186044182e883ac27f7
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Apr 8 23:50:45 2021 +0100
[RNG-128] Add UnitBallSampler
---
.../jmh/sampling/UnitBallSamplerBenchmark.java | 803 +++++++++++++++++++++
.../commons/rng/sampling/UnitBallSampler.java | 236 ++++++
.../commons/rng/sampling/UnitBallSamplerTest.java | 423 +++++++++++
src/changes/changes.xml | 3 +
4 files changed, 1465 insertions(+)
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
new file mode 100644
index 0000000..9a7f883
--- /dev/null
+++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitBallSamplerBenchmark.java
@@ -0,0 +1,803 @@
+/*
+ * 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.AhrensDieterExponentialSampler;
+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 within an N-dimension unit ball.
+ */
+@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 UnitBallSamplerBenchmark {
+
+ /** Name for the baseline method. */
+ private static final String BASELINE = "Baseline";
+ /** Name for the rejection method. */
+ private static final String REJECTION = "Rejection";
+ /** Name for the disk-point picking method. */
+ private static final String DISK_POINT = "DiskPoint";
+ /** Name for the ball-point picking method. */
+ private static final String BALL_POINT = "BallPoint";
+ /** Name for the method moving from the surface of the hypersphere to an internal point. */
+ private static final String HYPERSPHERE_INTERNAL = "HypersphereInternal";
+ /** Name for the picking from the surface of a greater dimension hypersphere and discarding 2 points. */
+ private static final String HYPERSPHERE_DISCARD = "HypersphereDiscard";
+ /** 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 a sampler using a provided source of randomness.
+ */
+ private abstract static class BaseSampler implements Sampler {
+ /** The source of randomness. */
+ protected UniformRandomProvider rng;
+
+ /**
+ * Create an instance.
+ *
+ * @param rng the source of randomness
+ */
+ BaseSampler(UniformRandomProvider rng) {
+ this.rng = rng;
+ }
+ }
+
+ /**
+ * 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({//"1",
+ "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 two doubles method. */
+ private static final String TWO_DOUBLES = "twoDoubles";
+ /** Name for the boolean double method. */
+ private static final String BOOLEAN_DOUBLE = "booleanDouble";
+
+ /** The sampler type. */
+ @Param({BASELINE, SIGNED_DOUBLE, TWO_DOUBLES, BOOLEAN_DOUBLE})
+ 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[] {0.5};
+ }
+ };
+ } else if (SIGNED_DOUBLE.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ // Sample [-1, 1) uniformly
+ return new double[] {makeSignedDouble(rng.nextLong())};
+ }
+ };
+ } else if (TWO_DOUBLES.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ // Sample [-1, 1) excluding -0.0 but also missing the final 1.0 - 2^-53.
+ // The 1.0 could be adjusted to 1.0 - 2^-53 to create the interval (-1, 1).
+ return new double[] {rng.nextDouble() + rng.nextDouble() - 1.0};
+ }
+ };
+ } else if (BOOLEAN_DOUBLE.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ // This will sample (-1, 1) including -0.0 and 0.0
+ return new double[] {rng.nextBoolean() ? -rng.nextDouble() : rng.nextDouble()};
+ }
+ };
+ }
+ throw new IllegalStateException(UNKNOWN_SAMPLER + type);
+ }
+ }
+
+ /**
+ * The 2D unit disk sampler.
+ */
+ @State(Scope.Benchmark)
+ public static class Sampler2D extends SamplerData {
+ /** The sampler type. */
+ @Param({BASELINE, REJECTION, DISK_POINT, HYPERSPHERE_INTERNAL, HYPERSPHERE_DISCARD})
+ 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[] {0.5, 0};
+ }
+ };
+ } else if (REJECTION.equals(type)) {
+ return new RejectionSampler(rng);
+ } else if (DISK_POINT.equals(type)) {
+ return new DiskPointSampler(rng);
+ } else if (HYPERSPHERE_INTERNAL.equals(type)) {
+ return new HypersphereInternalSampler(rng);
+ } else if (HYPERSPHERE_DISCARD.equals(type)) {
+ return new HypersphereDiscardSampler(rng);
+ }
+ throw new IllegalStateException(UNKNOWN_SAMPLER + type);
+ }
+
+ /**
+ * Sample using a simple rejection method.
+ */
+ private static class RejectionSampler extends BaseSampler {
+ /**
+ * @param rng the source of randomness
+ */
+ RejectionSampler(UniformRandomProvider rng) {
+ super(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ // Generate via rejection method of a circle inside a square of edge length 2.
+ // This should compute approximately 2^2 / pi = 1.27 square positions per sample.
+ double x;
+ double y;
+ do {
+ x = makeSignedDouble(rng.nextLong());
+ y = makeSignedDouble(rng.nextLong());
+ } while (x * x + y * y > 1.0);
+ return new double[] {x, y};
+ }
+ }
+
+ /**
+ * Sample using disk point picking.
+ * @see <a href="https://mathworld.wolfram.com/DiskPointPicking.html">Disk point picking</a>
+ */
+ private static class DiskPointSampler extends BaseSampler {
+ /** 2 pi. */
+ private static final double TWO_PI = 2 * Math.PI;
+
+ /**
+ * @param rng the source of randomness
+ */
+ DiskPointSampler(UniformRandomProvider rng) {
+ super(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ final double t = TWO_PI * rng.nextDouble();
+ final double r = Math.sqrt(rng.nextDouble());
+ final double x = r * Math.cos(t);
+ final double y = r * Math.sin(t);
+ return new double[] {x, y};
+ }
+ }
+
+ /**
+ * Choose a uniform point X on the unit hypersphere, then multiply it by U<sup>1/n</sup>
+ * where U in [0, 1].
+ * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
+ * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
+ */
+ private static class HypersphereInternalSampler extends BaseSampler {
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+
+ /**
+ * @param rng the source of randomness
+ */
+ HypersphereInternalSampler(UniformRandomProvider rng) {
+ super(rng);
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ final double x = normal.sample();
+ final double y = normal.sample();
+ final double sum = x * x + y * y;
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ // Take a point on the unit hypersphere then multiply it by U^1/n
+ final double f = Math.sqrt(rng.nextDouble()) / Math.sqrt(sum);
+ return new double[] {x * f, y * f};
+ }
+ }
+
+ /**
+ * Take a random point on the (n+1)-dimensional hypersphere and drop two coordinates.
+ * Remember that the (n+1)-hypersphere is the unit sphere of R^(n+2).
+ * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
+ * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
+ */
+ private static class HypersphereDiscardSampler extends BaseSampler {
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+
+ /**
+ * @param rng the source of randomness
+ */
+ HypersphereDiscardSampler(UniformRandomProvider rng) {
+ super(rng);
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ // Discard 2 samples from the coordinate but include in the sum
+ final double x0 = normal.sample();
+ final double x1 = normal.sample();
+ final double x = normal.sample();
+ final double y = normal.sample();
+ final double sum = x0 * x0 + x1 * x1 + x * x + y * y;
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ final double f = 1.0 / Math.sqrt(sum);
+ return new double[] {x * f, y * f};
+ }
+ }
+ }
+
+ /**
+ * The 3D unit ball sampler.
+ */
+ @State(Scope.Benchmark)
+ public static class Sampler3D extends SamplerData {
+ /** The sampler type. */
+ @Param({BASELINE, REJECTION, BALL_POINT, HYPERSPHERE_INTERNAL, HYPERSPHERE_DISCARD})
+ 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[] {0.5, 0, 0};
+ }
+ };
+ } else if (REJECTION.equals(type)) {
+ return new RejectionSampler(rng);
+ } else if (BALL_POINT.equals(type)) {
+ return new BallPointSampler(rng);
+ } else if (HYPERSPHERE_INTERNAL.equals(type)) {
+ return new HypersphereInternalSampler(rng);
+ } else if (HYPERSPHERE_DISCARD.equals(type)) {
+ return new HypersphereDiscardSampler(rng);
+ }
+ throw new IllegalStateException(UNKNOWN_SAMPLER + type);
+ }
+
+ /**
+ * Sample using a simple rejection method.
+ */
+ private static class RejectionSampler extends BaseSampler {
+ /**
+ * @param rng the source of randomness
+ */
+ RejectionSampler(UniformRandomProvider rng) {
+ super(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ // Generate via rejection method of a ball inside a cube of edge length 2.
+ // This should compute approximately 2^3 / (4pi/3) = 1.91 cube positions per sample.
+ double x;
+ double y;
+ double z;
+ do {
+ x = makeSignedDouble(rng.nextLong());
+ y = makeSignedDouble(rng.nextLong());
+ z = makeSignedDouble(rng.nextLong());
+ } while (x * x + y * y + z * z > 1.0);
+ return new double[] {x, y, z};
+ }
+ }
+
+ /**
+ * Sample using ball point picking.
+ * @see <a href="https://mathworld.wolfram.com/BallPointPicking.html">Ball point picking</a>
+ */
+ private static class BallPointSampler extends BaseSampler {
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+ /** The exponential distribution. */
+ private final AhrensDieterExponentialSampler exp;
+
+ /**
+ * @param rng the source of randomness
+ */
+ BallPointSampler(UniformRandomProvider rng) {
+ super(rng);
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ // Exponential(mean=2) == Chi-squared distribution(degrees freedom=2)
+ // thus is the equivalent of the HypersphereDiscardSampler.
+ exp = new AhrensDieterExponentialSampler(rng, 2.0);
+ }
+
+ @Override
+ public double[] sample() {
+ final double x = normal.sample();
+ final double y = normal.sample();
+ final double z = normal.sample();
+ // Include the exponential sample
+ final double sum = exp.sample() + x * x + y * y + z * z;
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ final double f = 1.0 / Math.sqrt(sum);
+ return new double[] {x * f, y * f, z * f};
+ }
+ }
+
+ /**
+ * Choose a uniform point X on the unit hypersphere, then multiply it by U<sup>1/n</sup>
+ * where U in [0, 1].
+ * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
+ * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
+ */
+ private static class HypersphereInternalSampler extends BaseSampler {
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+
+ /**
+ * @param rng the source of randomness
+ */
+ HypersphereInternalSampler(UniformRandomProvider rng) {
+ super(rng);
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ final double x = normal.sample();
+ final double y = normal.sample();
+ final double z = normal.sample();
+ final double sum = x * x + y * y + z * z;
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ // Take a point on the unit hypersphere then multiply it by U^1/n
+ final double f = Math.cbrt(rng.nextDouble()) / Math.sqrt(sum);
+ return new double[] {x * f, y * f, z * f};
+ }
+ }
+
+ /**
+ * Take a random point on the (n+1)-dimensional hypersphere and drop two coordinates.
+ * Remember that the (n+1)-hypersphere is the unit sphere of R^(n+2).
+ * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
+ * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
+ */
+ private static class HypersphereDiscardSampler extends BaseSampler {
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+
+ /**
+ * @param rng the source of randomness
+ */
+ HypersphereDiscardSampler(UniformRandomProvider rng) {
+ super(rng);
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ // Discard 2 samples from the coordinate but include in the sum
+ final double x0 = normal.sample();
+ final double x1 = normal.sample();
+ final double x = normal.sample();
+ final double y = normal.sample();
+ final double z = normal.sample();
+ final double sum = x0 * x0 + x1 * x1 + x * x + y * y + z * z;
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ final double f = 1.0 / Math.sqrt(sum);
+ return new double[] {x * f, y * f, z * f};
+ }
+ }
+ }
+
+ /**
+ * The ND unit ball sampler.
+ */
+ @State(Scope.Benchmark)
+ public static class SamplerND extends SamplerData {
+ /** The sampler type. */
+ @Param({BASELINE, REJECTION, BALL_POINT, HYPERSPHERE_INTERNAL, HYPERSPHERE_DISCARD})
+ private String type;
+ /** The number of dimensions. */
+ @Param({"3", "4", "5"})
+ private int dimension;
+
+ /** {@inheritDoc} */
+ @Override
+ protected Sampler createSampler(final UniformRandomProvider rng) {
+ if (BASELINE.equals(type)) {
+ return new Sampler() {
+ private final int dim = dimension;
+ @Override
+ public double[] sample() {
+ final double[] sample = new double[dim];
+ for (int i = 0; i < dim; i++) {
+ sample[i] = 0.01;
+ }
+ return sample;
+ }
+ };
+ } else if (REJECTION.equals(type)) {
+ return new RejectionSampler(rng, dimension);
+ } else if (BALL_POINT.equals(type)) {
+ return new BallPointSampler(rng, dimension);
+ } else if (HYPERSPHERE_INTERNAL.equals(type)) {
+ return new HypersphereInternalSampler(rng, dimension);
+ } else if (HYPERSPHERE_DISCARD.equals(type)) {
+ return new HypersphereDiscardSampler(rng, dimension);
+ }
+ throw new IllegalStateException(UNKNOWN_SAMPLER + type);
+ }
+
+ /**
+ * Sample using a simple rejection method.
+ */
+ private static class RejectionSampler extends BaseSampler {
+ /** The dimension. */
+ private final int dimension;
+
+ /**
+ * @param rng the source of randomness
+ * @param dimension the dimension
+ */
+ RejectionSampler(UniformRandomProvider rng, int dimension) {
+ super(rng);
+ this.dimension = dimension;
+ }
+
+ @Override
+ public double[] sample() {
+ // Generate via rejection method of a ball inside a hypercube of edge length 2.
+ final double[] sample = new double[dimension];
+ double sum;
+ do {
+ sum = 0;
+ for (int i = 0; i < dimension; i++) {
+ final double x = makeSignedDouble(rng.nextLong());
+ sum += x * x;
+ sample[i] = x;
+ }
+ } while (sum > 1.0);
+ return sample;
+ }
+ }
+
+ /**
+ * Sample using ball point picking.
+ * @see <a href="https://mathworld.wolfram.com/BallPointPicking.html">Ball point picking</a>
+ */
+ private static class BallPointSampler extends BaseSampler {
+ /** The dimension. */
+ private final int dimension;
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+ /** The exponential distribution. */
+ private final AhrensDieterExponentialSampler exp;
+
+ /**
+ * @param rng the source of randomness
+ * @param dimension the dimension
+ */
+ BallPointSampler(UniformRandomProvider rng, int dimension) {
+ super(rng);
+ this.dimension = dimension;
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ // Exponential(mean=2) == Chi-squared distribution(degrees freedom=2)
+ // thus is the equivalent of the HypersphereDiscardSampler.
+ exp = new AhrensDieterExponentialSampler(rng, 2.0);
+ }
+
+ @Override
+ public double[] sample() {
+ final double[] sample = new double[dimension];
+ // Include the exponential sample
+ double sum = exp.sample();
+ for (int i = 0; i < dimension; i++) {
+ final double x = normal.sample();
+ sum += x * x;
+ sample[i] = x;
+ }
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ final double f = 1.0 / Math.sqrt(sum);
+ for (int i = 0; i < dimension; i++) {
+ sample[i] *= f;
+ }
+ return sample;
+ }
+ }
+
+ /**
+ * Choose a uniform point X on the unit hypersphere, then multiply it by U<sup>1/n</sup>
+ * where U in [0, 1].
+ * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
+ * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
+ */
+ private static class HypersphereInternalSampler extends BaseSampler {
+ /** The dimension. */
+ private final int dimension;
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+ /** Reciprocal of the dimension. */
+ private final double power;
+
+ /**
+ * @param rng the source of randomness
+ * @param dimension the dimension
+ */
+ HypersphereInternalSampler(UniformRandomProvider rng, int dimension) {
+ super(rng);
+ this.dimension = dimension;
+ power = 1.0 / dimension;
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ final double[] sample = new double[dimension];
+ double sum = 0;
+ for (int i = 0; i < dimension; i++) {
+ final double x = normal.sample();
+ sum += x * x;
+ sample[i] = x;
+ }
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ // Take a point on the unit hypersphere then multiply it by U^1/n
+ final double f = Math.pow(rng.nextDouble(), power) / Math.sqrt(sum);
+ for (int i = 0; i < dimension; i++) {
+ sample[i] *= f;
+ }
+ return sample;
+ }
+ }
+
+ /**
+ * Take a random point on the (n+1)-dimensional hypersphere and drop two coordinates.
+ * Remember that the (n+1)-hypersphere is the unit sphere of R^(n+2).
+ * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
+ * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
+ */
+ private static class HypersphereDiscardSampler extends BaseSampler {
+ /** The dimension. */
+ private final int dimension;
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+
+ /**
+ * @param rng the source of randomness
+ * @param dimension the dimension
+ */
+ HypersphereDiscardSampler(UniformRandomProvider rng, int dimension) {
+ super(rng);
+ this.dimension = dimension;
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ final double[] sample = new double[dimension];
+ // Discard 2 samples from the coordinate but include in the sum
+ final double x0 = normal.sample();
+ final double x1 = normal.sample();
+ double sum = x0 * x0 + x1 * x1;
+ for (int i = 0; i < dimension; i++) {
+ final double x = normal.sample();
+ sum += x * x;
+ sample[i] = x;
+ }
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ final double f = 1.0 / Math.sqrt(sum);
+ for (int i = 0; i < dimension; i++) {
+ sample[i] *= f;
+ }
+ return sample;
+ }
+ }
+ }
+
+ /**
+ * Creates a signed double in the range {@code [-1, 1)}. The magnitude is sampled evenly
+ * from the 2<sup>54</sup> dyadic rationals in the range.
+ *
+ * <p>Note: This method will not return samples for both -0.0 and 0.0.
+ *
+ * @param bits the bits
+ * @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);
+ }
+
+ /**
+ * 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 disk.
+ *
+ * @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 ball.
+ *
+ * @param bh Data sink
+ * @param data Input data.
+ */
+ @Benchmark
+ public void create3D(Blackhole bh, Sampler3D data) {
+ runSampler(bh, data);
+ }
+
+ /**
+ * Generation of uniform samples from an ND unit ball.
+ *
+ * @param bh Data sink
+ * @param data Input data.
+ */
+ @Benchmark
+ public void createND(Blackhole bh, SamplerND 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
new file mode 100644
index 0000000..e0c2b69
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitBallSampler.java
@@ -0,0 +1,236 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.rng.sampling;
+
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
+import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
+
+/**
+ * Generate coordinates <a href="http://mathworld.wolfram.com/BallPointPicking.html">
+ * uniformly distributed within the unit n-ball</a>.
+ *
+ * <p>Sampling uses:</p>
+ *
+ * <ul>
+ * <li>{@link UniformRandomProvider#nextLong()}
+ * <li>{@link UniformRandomProvider#nextDouble()} (only for dimensions above 2)
+ * </ul>
+ *
+ * @since 1.4
+ */
+public abstract class UnitBallSampler implements SharedStateSampler<UnitBallSampler> {
+ /** 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;
+
+ /**
+ * Sample uniformly from a 1D unit line.
+ */
+ private static class UnitBallSampler1D extends UnitBallSampler {
+ /** The source of randomness. */
+ private final UniformRandomProvider rng;
+
+ /**
+ * @param rng Source of randomness.
+ */
+ UnitBallSampler1D(UniformRandomProvider rng) {
+ this.rng = rng;
+ }
+
+ @Override
+ public double[] sample() {
+ return new double[] {makeSignedDouble(rng.nextLong())};
+ }
+
+ @Override
+ public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
+ return new UnitBallSampler1D(rng);
+ }
+ }
+
+ /**
+ * Sample uniformly from a 2D unit disk.
+ */
+ private static class UnitBallSampler2D extends UnitBallSampler {
+ /** The source of randomness. */
+ private final UniformRandomProvider rng;
+
+ /**
+ * @param rng Source of randomness.
+ */
+ UnitBallSampler2D(UniformRandomProvider rng) {
+ this.rng = rng;
+ }
+
+ @Override
+ public double[] sample() {
+ // Generate via rejection method of a circle inside a square of edge length 2.
+ // This should compute approximately 2^2 / pi = 1.27 square positions per sample.
+ double x;
+ double y;
+ do {
+ x = makeSignedDouble(rng.nextLong());
+ y = makeSignedDouble(rng.nextLong());
+ } while (x * x + y * y > 1.0);
+ return new double[] {x, y};
+ }
+
+ @Override
+ public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
+ return new UnitBallSampler2D(rng);
+ }
+ }
+
+ /**
+ * Sample uniformly from a 3D unit ball. This is an non-array based specialisation of
+ * {@link UnitBallSamplerND} for performance.
+ */
+ private static class UnitBallSampler3D extends UnitBallSampler {
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+
+ /**
+ * @param rng Source of randomness.
+ */
+ UnitBallSampler3D(UniformRandomProvider rng) {
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ // Discard 2 samples from the coordinate but include in the sum
+ final double x0 = normal.sample();
+ final double x1 = normal.sample();
+ final double x = normal.sample();
+ final double y = normal.sample();
+ final double z = normal.sample();
+ final double sum = x0 * x0 + x1 * x1 + x * x + y * y + z * z;
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ final double f = 1.0 / Math.sqrt(sum);
+ return new double[] {x * f, y * f, z * f};
+ }
+
+ @Override
+ public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
+ return new UnitBallSampler3D(rng);
+ }
+ }
+
+ /**
+ * Sample uniformly from a unit n-ball.
+ * Take a random point on the (n+1)-dimensional hypersphere and drop two coordinates.
+ * Remember that the (n+1)-hypersphere is the unit sphere of R^(n+2), i.e. the surface
+ * of the (n+2)-dimensional ball.
+ * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
+ * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
+ */
+ private static class UnitBallSamplerND extends UnitBallSampler {
+ /** The dimension. */
+ private final int dimension;
+ /** The normal distribution. */
+ private final NormalizedGaussianSampler normal;
+
+ /**
+ * @param dimension Space dimension.
+ * @param rng Source of randomness.
+ */
+ UnitBallSamplerND(int dimension, UniformRandomProvider rng) {
+ this.dimension = dimension;
+ normal = new ZigguratNormalizedGaussianSampler(rng);
+ }
+
+ @Override
+ public double[] sample() {
+ final double[] sample = new double[dimension];
+ // Discard 2 samples from the coordinate but include in the sum
+ final double x0 = normal.sample();
+ final double x1 = normal.sample();
+ double sum = x0 * x0 + x1 * x1;
+ for (int i = 0; i < dimension; i++) {
+ final double x = normal.sample();
+ sum += x * x;
+ sample[i] = x;
+ }
+ // Note: Handle the possibility of a zero sum and invalid inverse
+ if (sum == 0) {
+ return sample();
+ }
+ final double f = 1.0 / Math.sqrt(sum);
+ for (int i = 0; i < dimension; i++) {
+ sample[i] *= f;
+ }
+ return sample;
+ }
+
+ @Override
+ public UnitBallSampler withUniformRandomProvider(UniformRandomProvider rng) {
+ return new UnitBallSamplerND(dimension, rng);
+ }
+ }
+
+ /**
+ * @return a random Cartesian coordinate within the unit n-ball.
+ */
+ public abstract double[] sample();
+
+ /**
+ * Create a unit n-ball sampler for the given dimension.
+ *
+ * @param dimension Space dimension.
+ * @param rng Generator for the individual components of the coordinates. A shallow
+ * copy will be stored in this instance.
+ * @return the sampler
+ * @throws IllegalArgumentException If {@code dimension <= 0}
+ */
+ public static UnitBallSampler of(int dimension,
+ UniformRandomProvider rng) {
+ if (dimension <= 0) {
+ throw new IllegalArgumentException("Dimension must be strictly positive");
+ } else if (dimension == ONE_D) {
+ return new UnitBallSampler1D(rng);
+ } else if (dimension == TWO_D) {
+ return new UnitBallSampler2D(rng);
+ } else if (dimension == THREE_D) {
+ return new UnitBallSampler3D(rng);
+ }
+ return new UnitBallSamplerND(dimension, rng);
+ }
+
+ /**
+ * Creates a signed double in the range {@code [-1, 1)}. The magnitude is sampled evenly
+ * from the 2<sup>54</sup> dyadic rationals in the range.
+ *
+ * <p>Note: This method will not return samples for both -0.0 and 0.0.
+ *
+ * @param bits the bits
+ * @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);
+ }
+}
diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/UnitBallSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/UnitBallSamplerTest.java
new file mode 100644
index 0000000..6f2b9b0
--- /dev/null
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/UnitBallSamplerTest.java
@@ -0,0 +1,423 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.rng.sampling;
+
+import org.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;
+
+/**
+ * Test for {@link UnitBallSampler}.
+ */
+public class UnitBallSamplerTest {
+ /**
+ * Test a non-positive dimension.
+ */
+ @Test(expected = IllegalArgumentException.class)
+ public void testInvalidDimensionThrows() {
+ UnitBallSampler.of(0, null);
+ }
+
+ /**
+ * Test the distribution of points in one dimension.
+ */
+ @Test
+ public void testDistribution1D() {
+ testDistributionND(1);
+ }
+
+ /**
+ * Test the distribution of points in two dimensions.
+ */
+ @Test
+ public void testDistribution2D() {
+ testDistributionND(2);
+ }
+
+ /**
+ * Test the distribution of points in three dimensions.
+ */
+ @Test
+ public void testDistribution3D() {
+ testDistributionND(3);
+ }
+
+ /**
+ * Test the distribution of points in four dimensions.
+ */
+ @Test
+ public void testDistribution4D() {
+ testDistributionND(4);
+ }
+
+ /**
+ * Test the distribution of points in five dimensions.
+ */
+ @Test
+ public void testDistribution5D() {
+ testDistributionND(5);
+ }
+
+ /**
+ * Test the distribution of points in six dimensions.
+ */
+ @Test
+ public void testDistribution6D() {
+ testDistributionND(6);
+ }
+
+ /**
+ * Test the distribution of points in n dimensions. The output coordinates
+ * should be uniform in the unit n-ball. The unit n-ball is divided into inner
+ * n-balls. The radii of the internal n-balls are varied to ensure that successive layers
+ * have the same volume. This assigns each coordinate to an inner n-ball layer and an
+ * orthant using the sign bits of the coordinates. The number of samples in each bin
+ * should be the same.
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Volume_of_an_n-ball">Volume of an n-ball</a>
+ * @see <a href="https://en.wikipedia.org/wiki/Orthant">Orthant</a>
+ */
+ private static void testDistributionND(int dimension) {
+ // The number of inner layers and samples has been selected by trial and error using
+ // random seeds and multiple runs to ensure correctness of the test (i.e. it fails with
+ // approximately the fraction expected for the test p-value).
+ // A fixed seed is used to make the test suite robust.
+ final int layers = 10;
+ final int samplesPerBin = 20;
+ final int orthants = 1 << dimension;
+
+ // Compute the radius for each layer to have the same volume.
+ final double volume = createVolumeFunction(dimension).apply(1);
+ final DoubleUnaryOperator radius = createRadiusFunction(dimension);
+ final double[] r = new double[layers];
+ for (int i = 1; i < layers; i++) {
+ r[i - 1] = radius.apply(volume * ((double) i / layers));
+ }
+ // The final radius should be 1.0. Any coordinates with a radius above 1
+ // should fail so explicitly set the value as 1.
+ r[layers - 1] = 1.0;
+
+ // Expect a uniform distribution
+ final double[] expected = new double[layers * orthants];
+ final int samples = samplesPerBin * expected.length;
+ Arrays.fill(expected, (double) samples / layers);
+
+ // Increase the loops and use a null seed (i.e. randomly generated) to verify robustness
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_SHI_RO_512_PP, 0xa1b2c3d4L);
+ final UnitBallSampler sampler = UnitBallSampler.of(dimension, rng);
+ for (int loop = 0; loop < 1; loop++) {
+ // Assign each coordinate to a layer inside the ball and an orthant using the sign
+ final long[] observed = new long[layers * orthants];
+ NEXT:
+ for (int i = 0; i < samples; i++) {
+ final double[] v = sampler.sample();
+ final double length = length(v);
+ for (int layer = 0; layer < layers; layer++) {
+ if (length <= r[layer]) {
+ final int orthant = orthant(v);
+ observed[layer * orthants + orthant]++;
+ continue NEXT;
+ }
+ }
+ // Radius above 1
+ Assert.fail("Invalid sample length: " + length);
+ }
+ final double p = new ChiSquareTest().chiSquareTest(expected, observed);
+ Assert.assertFalse("p-value too small: " + p, p < 0.001);
+ }
+ }
+
+ /**
+ * 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.
+ */
+ 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(0x1a2b3cL) {
+ private int count = -2 * dimension;
+
+ @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 = UnitBallSampler.of(dimension, bad).sample();
+ Assert.assertEquals(dimension, vector.length);
+ // A non-zero coordinate should occur with a SplitMix which returns 0 only once.
+ Assert.assertNotEquals(0.0, length(vector));
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for 1D.
+ */
+ @Test
+ public void testSharedStateSampler1D() {
+ testSharedStateSampler(1);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for 2D.
+ */
+ @Test
+ public void testSharedStateSampler2D() {
+ testSharedStateSampler(2);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for 3D.
+ */
+ @Test
+ public void testSharedStateSampler3D() {
+ testSharedStateSampler(3);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for 4D.
+ */
+ @Test
+ public void testSharedStateSampler4D() {
+ testSharedStateSampler(4);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for the given dimension.
+ */
+ private static void testSharedStateSampler(int dimension) {
+ final UniformRandomProvider rng1 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
+ final UniformRandomProvider rng2 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
+ final UnitBallSampler sampler1 = UnitBallSampler.of(dimension, rng1);
+ final UnitBallSampler sampler2 = sampler1.withUniformRandomProvider(rng2);
+ RandomAssert.assertProduceSameSequence(
+ new RandomAssert.Sampler<double[]>() {
+ @Override
+ public double[] sample() {
+ return sampler1.sample();
+ }
+ },
+ new RandomAssert.Sampler<double[]>() {
+ @Override
+ public double[] sample() {
+ return sampler2.sample();
+ }
+ });
+ }
+
+ /**
+ * @return the length (L2-norm) of given vector.
+ */
+ private static double length(double[] vector) {
+ double total = 0;
+ for (double d : vector) {
+ total += d * d;
+ }
+ return Math.sqrt(total);
+ }
+
+ /**
+ * Assign an orthant to the vector using the sign of each component.
+ * The i<sup>th</sup> bit is set in the orthant for the i<sup>th</sup> component
+ * if the component is negative.
+ *
+ * @return the orthant in the range [0, vector.length)
+ * @see <a href="https://en.wikipedia.org/wiki/Orthant">Orthant</a>
+ */
+ private static int orthant(double[] vector) {
+ int orthant = 0;
+ for (int i = 0; i < vector.length; i++) {
+ if (vector[i] < 0) {
+ orthant |= 1 << i;
+ }
+ }
+ return orthant;
+ }
+
+ /**
+ * Check the n-ball volume functions can map the radius to the volume and back.
+ * These functions are used to divide the n-ball into uniform volume bins to test sampling
+ * within the n-ball.
+ */
+ @Test
+ public void checkVolumeFunctions() {
+ final double[] radii = {0, 0.1, 0.25, 0.5, 0.75, 1.0};
+ for (int n = 1; n <= 6; n++) {
+ final DoubleUnaryOperator volume = createVolumeFunction(n);
+ final DoubleUnaryOperator radius = createRadiusFunction(n);
+ for (final double r : radii) {
+ Assert.assertEquals(r, radius.apply(volume.apply(r)), 1e-10);
+ }
+ }
+ }
+
+ /**
+ * Specify computation of a double-valued result for a given double.
+ *
+ * <p>This can be replaced with java.util.function.DoubleUnaryOperator when the source
+ * requires Java 1.8.
+ */
+ private interface DoubleUnaryOperator {
+ /**
+ * Compute the result.
+ *
+ * @param value the input value
+ * @return the result
+ */
+ double apply(double value);
+ }
+
+ /**
+ * Creates a function to compute the volume of a ball of the given dimension
+ * from the radius.
+ *
+ * @param dimension the dimension
+ * @return the volume function
+ * @see <a href="https://en.wikipedia.org/wiki/Volume_of_an_n-ball">Volume of an n-ball</a>
+ */
+ private static DoubleUnaryOperator createVolumeFunction(final int dimension) {
+ // This can be simplified using lambda functions when the source requires Java 1.8
+ if (dimension == 1) {
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double r) {
+ return r * 2;
+ }
+ };
+ } else if (dimension == 2) {
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double r) {
+ return Math.PI * r * r;
+ }
+ };
+ } else if (dimension == 3) {
+ final double factor = 4 * Math.PI / 3;
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double r) {
+ return factor * Math.pow(r, 3);
+ }
+ };
+ } else if (dimension == 4) {
+ final double factor = Math.PI * Math.PI / 2;
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double r) {
+ return factor * Math.pow(r, 4);
+ }
+ };
+ } else if (dimension == 5) {
+ final double factor = 8 * Math.PI * Math.PI / 15;
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double r) {
+ return factor * Math.pow(r, 5);
+ }
+ };
+ } else if (dimension == 6) {
+ final double factor = Math.pow(Math.PI, 3) / 6;
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double r) {
+ return factor * Math.pow(r, 6);
+ }
+ };
+ }
+ throw new IllegalStateException("Unsupported dimension: " + dimension);
+ }
+
+ /**
+ * Creates a function to compute the radius of a ball of the given dimension
+ * from the volume.
+ *
+ * @param dimension the dimension
+ * @return the radius function
+ * @see <a href="https://en.wikipedia.org/wiki/Volume_of_an_n-ball">Volume of an n-ball</a>
+ */
+ private static DoubleUnaryOperator createRadiusFunction(final int dimension) {
+ // This can be simplified using lambda functions when the source requires Java 1.8
+ if (dimension == 1) {
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double v) {
+ return v * 0.5;
+ }
+ };
+ } else if (dimension == 2) {
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double v) {
+ return Math.sqrt(v / Math.PI);
+ }
+ };
+ } else if (dimension == 3) {
+ final double factor = 3.0 / (4 * Math.PI);
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double v) {
+ return Math.cbrt(v * factor);
+ }
+ };
+ } else if (dimension == 4) {
+ final double factor = 2.0 / (Math.PI * Math.PI);
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double v) {
+ return Math.pow(v * factor, 0.25);
+ }
+ };
+ } else if (dimension == 5) {
+ final double factor = 15.0 / (8 * Math.PI * Math.PI);
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double v) {
+ return Math.pow(v * factor, 0.2);
+ }
+ };
+ } else if (dimension == 6) {
+ final double factor = 6.0 / Math.pow(Math.PI, 3);
+ return new DoubleUnaryOperator() {
+ @Override
+ public double apply(double v) {
+ return Math.pow(v * factor, 1.0 / 6);
+ }
+ };
+ }
+ throw new IllegalStateException("Unsupported dimension: " + dimension);
+ }
+}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index fe0fe49..a8bbd88 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -75,6 +75,9 @@ 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="add" issue="128">
+ New "UnitBallSampler" to generate coordinates uniformly within an n-unit ball.
+ </action>
<action dev="aherbert" type="add" issue="126">
"PoissonSamplerCache": Method to return a SharedStateDiscreteSampler.
</action>