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/08/25 22:35:14 UTC
[commons-rng] branch master updated (2286528 -> c3d7e60)
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 2286528 Update observed frequencies with expected where appropriate.
new e9a67d3 Update ziggurat performance benchmark with more methods
new c3d7e60 PMD: Allow long class for ZigguratSampler
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:
.../distribution/ZigguratSamplerPerformance.java | 271 ++++++++++++++++++++-
src/main/resources/pmd/pmd-ruleset.xml | 2 +-
2 files changed, 269 insertions(+), 4 deletions(-)
[commons-rng] 02/02: PMD: Allow long class for ZigguratSampler
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 c3d7e6044b6585b696285422c3c6a6aca33c9174
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Wed Aug 25 23:35:11 2021 +0100
PMD: Allow long class for ZigguratSampler
---
src/main/resources/pmd/pmd-ruleset.xml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/main/resources/pmd/pmd-ruleset.xml b/src/main/resources/pmd/pmd-ruleset.xml
index ddce4af..ae7f871 100644
--- a/src/main/resources/pmd/pmd-ruleset.xml
+++ b/src/main/resources/pmd/pmd-ruleset.xml
@@ -161,7 +161,7 @@
<properties>
<!-- The length is due to multiple implementations as inner classes -->
<property name="violationSuppressXPath" value="//ClassOrInterfaceDeclaration[@SimpleName='MarsagliaTsangWangDiscreteSampler'
- or @SimpleName='CompositeSamplers' or @SimpleName='StableSampler']"/>
+ or @SimpleName='CompositeSamplers' or @SimpleName='StableSampler' or @SimpleName='ZigguratSampler']"/>
</properties>
</rule>
<rule ref="category/java/design.xml/LogicInversion">
[commons-rng] 01/02: Update ziggurat performance benchmark with
more methods
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 e9a67d339fe4032101305086f89085bc1c18486d
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Wed Aug 25 21:02:41 2021 +0100
Update ziggurat performance benchmark with more methods
Use only bit shifts for random positive longs.
Use a different method to sample y.
Benchmark methods to interpolate X or Y
---
.../distribution/ZigguratSamplerPerformance.java | 271 ++++++++++++++++++++-
1 file changed, 268 insertions(+), 3 deletions(-)
diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java
index f681ecc..58e1f41 100644
--- a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java
+++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java
@@ -35,6 +35,7 @@ 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 java.util.Arrays;
import java.util.concurrent.TimeUnit;
/**
@@ -79,6 +80,8 @@ public class ZigguratSamplerPerformance {
private static final String MOD_GAUSSIAN_SIMPLE_OVERHANGS = "ModGaussianSimpleOverhangs";
/** The name for the {@link ModifiedZigguratNormalizedGaussianSamplerInlining}. */
private static final String MOD_GAUSSIAN_INLINING = "ModGaussianInlining";
+ /** The name for the {@link ModifiedZigguratNormalizedGaussianSamplerInliningShift}. */
+ private static final String MOD_GAUSSIAN_INLINING_SHIFT = "ModGaussianInliningShift";
/** The name for the {@link ModifiedZigguratNormalizedGaussianSamplerInliningSimpleOverhangs}. */
private static final String MOD_GAUSSIAN_INLINING_SIMPLE_OVERHANGS = "ModGaussianInliningSimpleOverhangs";
/** The name for the {@link ModifiedZigguratNormalizedGaussianSamplerIntMap}. */
@@ -248,6 +251,87 @@ public class ZigguratSamplerPerformance {
}
/**
+ * Defines method to use for interpolating the X or Y tables from unsigned {@code long} values.
+ */
+ @State(Scope.Benchmark)
+ public static class InterpolationSources {
+ /** The method to obtain the long. */
+ @Param({"U1", "1minusU2", "U_1minusU"})
+ private String method;
+
+ /** The sampler. */
+ private ContinuousSampler sampler;
+
+ /**
+ * @return the sampler.
+ */
+ public ContinuousSampler getSampler() {
+ return sampler;
+ }
+
+ /** Instantiates sampler. */
+ @Setup
+ public void setup() {
+ // Use a fast generator
+ final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create();
+ // Get an x table. This is length 254.
+ // We will sample from this internally to avoid index out-of-bounds issues.
+ final int start = 42;
+ final double[] x = Arrays.copyOfRange(
+ ModifiedZigguratNormalizedGaussianSampler.getX(), start, start + 129);
+ final int mask = 127;
+ if ("U1".equals(method)) {
+ sampler = () -> {
+ final long u = rng.nextLong() >>> 1;
+ final int j = 1 + (((int) u) & mask);
+ // double multiply
+ // double add
+ // double subtract
+ // long convert to double
+ // double multiply
+ // int subtract
+ // three index loads
+ return x[j] * TWO_POW_63 + (x[j - 1] - x[j]) * u;
+ };
+ } else if ("1minusU2".equals(method)) {
+ sampler = () -> {
+ final long u = rng.nextLong() >>> 1;
+ final int j = 1 + (((int) u) & mask);
+ // Since u is in [0, 2^63) to create (1 - u) using Long.MIN_VALUE
+ // as an unsigned integer of 2^63.
+ // double multiply
+ // double add
+ // double subtract
+ // long subtract
+ // long convert to double
+ // double multiply
+ // int subtract * 2
+ // three index loads
+ return x[j - 1] * TWO_POW_63 + (x[j] - x[j - 1]) * (Long.MIN_VALUE - u);
+ };
+ } else if ("U_1minusU".equals(method)) {
+ sampler = () -> {
+ final long u = rng.nextLong() >>> 1;
+ final int j = 1 + (((int) u) & mask);
+ // Interpolation between bounds a and b using:
+ // a * u + b * (1 - u) == b + u * (a - b)
+ // long convert to double
+ // double multiply
+ // double add
+ // long subtract
+ // long convert to double
+ // double multiply
+ // int subtract
+ // two index loads
+ return x[j - 1] * u + x[j] * (Long.MIN_VALUE - u);
+ };
+ } else {
+ throwIllegalStateException(method);
+ }
+ }
+ }
+
+ /**
* The samplers to use for testing the ziggurat method.
*/
@State(Scope.Benchmark)
@@ -260,7 +344,8 @@ public class ZigguratSamplerPerformance {
// Experimental Marsaglia exponential ziggurat sampler
EXPONENTIAL,
// Experimental McFarland Gaussian ziggurat samplers
- MOD_GAUSSIAN2, MOD_GAUSSIAN_SIMPLE_OVERHANGS, MOD_GAUSSIAN_INLINING,
+ MOD_GAUSSIAN2, MOD_GAUSSIAN_SIMPLE_OVERHANGS,
+ MOD_GAUSSIAN_INLINING, MOD_GAUSSIAN_INLINING_SHIFT,
MOD_GAUSSIAN_INLINING_SIMPLE_OVERHANGS, MOD_GAUSSIAN_INT_MAP, MOD_GAUSSIAN_512,
// Experimental McFarland Gaussian ziggurat samplers
MOD_EXPONENTIAL2, MOD_EXPONENTIAL_SIMPLE_OVERHANGS, MOD_EXPONENTIAL_INLINING,
@@ -292,6 +377,8 @@ public class ZigguratSamplerPerformance {
return new ModifiedZigguratNormalizedGaussianSamplerSimpleOverhangs(rng);
} else if (MOD_GAUSSIAN_INLINING.equals(type)) {
return new ModifiedZigguratNormalizedGaussianSamplerInlining(rng);
+ } else if (MOD_GAUSSIAN_INLINING_SHIFT.equals(type)) {
+ return new ModifiedZigguratNormalizedGaussianSamplerInliningShift(rng);
} else if (MOD_GAUSSIAN_INLINING_SIMPLE_OVERHANGS.equals(type)) {
return new ModifiedZigguratNormalizedGaussianSamplerInliningSimpleOverhangs(rng);
} else if (MOD_GAUSSIAN_INT_MAP.equals(type)) {
@@ -1101,6 +1188,17 @@ public class ZigguratSamplerPerformance {
exponential = new ModifiedZigguratExponentialSampler(rng);
}
+ /**
+ * Provide access to the precomputed ziggurat lengths.
+ *
+ * <p>This is package-private to allow usage in the interpolation test.
+ *
+ * @return x
+ */
+ static double[] getX() {
+ return X;
+ }
+
/** {@inheritDoc} */
@Override
public double sample() {
@@ -1478,6 +1576,160 @@ public class ZigguratSamplerPerformance {
*
* <p>Uses the algorithm from McFarland, C.D. (2016).
*
+ * <p>This implementation separates sampling of the main ziggurat and sampling from the edge
+ * into different methods. This allows inlining of the main sample method.
+ *
+ * <p>Positive longs are created using bit shifts (not masking). The y coordinate is
+ * generated with u2 not (1-u2) which avoids a subtraction.
+ */
+ static class ModifiedZigguratNormalizedGaussianSamplerInliningShift
+ extends ModifiedZigguratNormalizedGaussianSampler {
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ */
+ ModifiedZigguratNormalizedGaussianSamplerInliningShift(UniformRandomProvider rng) {
+ super(rng);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double sample() {
+ // Ideally this method byte code size should be below -XX:MaxInlineSize
+ // (which defaults to 35 bytes). This compiles to 33 bytes.
+
+ final long xx = nextLong();
+ // Float multiplication squashes these last 8 bits, so they can be used to sample i
+ final int i = ((int) xx) & 0xff;
+
+ if (i < I_MAX) {
+ // Early exit.
+ // Branch frequency: 0.988283
+ return X[i] * xx;
+ }
+
+ return edgeSample(xx);
+ }
+
+ /**
+ * Create the sample from the edge of the ziggurat.
+ *
+ * <p>This method has been extracted to fit the main sample method within 35 bytes (the
+ * default size for a JVM to inline a method).
+ *
+ * @param xx Initial random deviate
+ * @return a sample
+ */
+ private double edgeSample(long xx) {
+ // Recycle bits.
+ // Remove sign bit to create u1.
+ long u1 = (xx << 1) >>> 1;
+ // Extract the sign bit:
+ // Use 2 - 1 or 0 - 1
+ final double signBit = ((xx >>> 62) & 0x2) - 1.0;
+ final int j = normSampleA();
+ // Four kinds of overhangs:
+ // j = 0 : Sample from tail
+ // 0 < j < J_INFLECTION : Overhang is concave; only sample from Lower-Left triangle
+ // j = J_INFLECTION : Must sample from entire overhang rectangle
+ // j > J_INFLECTION : Overhangs are convex; implicitly accept point in Lower-Left triangle
+ //
+ // Conditional statements are arranged such that the more likely outcomes are first.
+ double x;
+ if (j > J_INFLECTION) {
+ // Convex overhang
+ // Branch frequency: 0.00892897
+ // Loop repeat frequency: 0.389804
+ for (;;) {
+ x = interpolateSample(X, j, u1);
+ final long uDiff = (nextLong() >>> 1) - u1;
+ if (uDiff >= 0) {
+ // Lower-left triangle
+ break;
+ }
+ if (uDiff >= MAX_IE &&
+ // Within maximum distance of f(x) from the triangle hypotenuse.
+ // Frequency (per upper-right triangle): 0.431497
+ // Reject frequency: 0.489630
+ // u2 = (u1 + uDiff)
+ interpolateSample(Y, j, u1 + uDiff) < Math.exp(-0.5 * x * x)) {
+ break;
+ }
+ // uDiff < MAX_IE (upper-right triangle) or rejected as above the curve
+ u1 = randomInt63();
+ }
+ } else if (j < J_INFLECTION) {
+ if (j == 0) {
+ // Tail
+ // Branch frequency: 0.000276358
+ // Note: Although less frequent than the next branch, j == 0 is a subset of
+ // j < J_INFLECTION and must be first.
+ // Loop repeat frequency: 0.0634786
+ do {
+ x = ONE_OVER_X_0 * exponential.sample();
+ } while (exponential.sample() < 0.5 * x * x);
+ x += X_0;
+ } else {
+ // Concave overhang
+ // Branch frequency: 0.00249563
+ // Loop repeat frequency: 0.0123784
+ for (;;) {
+ // U_x <- min(U_1, U_2)
+ // distance <- | U_1 - U_2 |
+ // U_y <- 1 - (U_x + distance)
+ long uDiff = (nextLong() >>> 1) - u1;
+ if (uDiff < 0) {
+ uDiff = -uDiff;
+ u1 -= uDiff;
+ }
+ x = interpolateSample(X, j, u1);
+ if (uDiff > MIN_IE ||
+ interpolateSample(Y, j, u1 + uDiff) < Math.exp(-0.5 * x * x)) {
+ break;
+ }
+ u1 = nextLong() >>> 1;
+ }
+ }
+ } else {
+ // Inflection point
+ // Branch frequency: 0.0000159359
+ // Loop repeat frequency: 0.500213
+ for (;;) {
+ x = interpolateSample(X, j, u1);
+ if (interpolateSample(Y, j, nextLong() >>> 1) < Math.exp(-0.5 * x * x)) {
+ break;
+ }
+ u1 = nextLong() >>> 1;
+ }
+ }
+ return signBit * x;
+ }
+
+ /**
+ * Interpolate x from the uniform deviate.
+ * <pre>
+ * value = x[j] + u * (x[j-1] - x[j])
+ * </pre>
+ * <p>If x is the precomputed lengths of the ziggurat (X) then x[j-1] is larger and this adds
+ * a delta to x[j].
+ * <p>If x is the precomputed pdf for the lengths of the ziggurat (Y) then x[j-1] is smaller
+ * larger and this subtracts a delta from x[j].
+ *
+ * @param x x
+ * @param j j
+ * @param u uniform deviate
+ * @return the sample
+ */
+ private static double interpolateSample(double[] x, int j, long u) {
+ return x[j] * TWO_POW_63 + (x[j - 1] - x[j]) * u;
+ }
+ }
+
+ /**
+ * Modified Ziggurat method for sampling from a Gaussian distribution with mean 0 and standard deviation 1.
+ *
+ * <p>Uses the algorithm from McFarland, C.D. (2016).
+ *
* <p>This implementation extracts the separates the sample method for the main ziggurat
* from the edge sampling. This allows inlining of the main sample method.
*
@@ -2899,12 +3151,12 @@ public class ZigguratSamplerPerformance {
protected double expOverhang(int j, long xx) {
// Recycle the initial random deviate.
// Shift right to make an unsigned long.
- for (long ux = xx >>> 1;; ux = randomInt63()) {
+ for (long ux = xx >>> 1;; ux = nextLong() >>> 1) {
// To sample a unit right-triangle:
// U_x <- min(U_1, U_2)
// distance <- | U_1 - U_2 |
// U_y <- 1 - (U_x + distance)
- long uDistance = randomInt63() - ux;
+ long uDistance = (nextLong() >>> 1) - ux;
if (uDistance < 0) {
uDistance = -uDistance;
ux -= uDistance;
@@ -3668,6 +3920,19 @@ public class ZigguratSamplerPerformance {
}
/**
+ * Benchmark methods for obtaining an interpolated value from an unsigned long.
+ *
+ * <p>Note: To be disabled. The speed is: U1, 1minusU2, U_1minusU.
+ *
+ * @param sources Source of randomness.
+ * @return the sample value
+ */
+ @Benchmark
+ public double interpolate(InterpolationSources sources) {
+ return sources.getSampler().sample();
+ }
+
+ /**
* Run the sampler.
*
* @param sources Source of randomness.