| /* |
| * 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); |
| } |
| } |