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 b8f4d82fc0eb5f00aafac90e4a040b9032be458b Author: Alex Herbert <[email protected]> AuthorDate: Sun May 2 12:29:28 2021 +0100 RNG-135: TetrahedronSampler to sample uniformly from a tetrahedron --- .../shape/TetrahedronSamplerBenchmark.java | 490 ++++++++++++++++++ .../commons/rng/sampling/shape/Coordinates.java | 28 +- .../rng/sampling/shape/TetrahedronSampler.java | 203 ++++++++ .../rng/sampling/shape/CoordinatesTest.java | 24 + .../rng/sampling/shape/TetrahedronSamplerTest.java | 551 +++++++++++++++++++++ src/changes/changes.xml | 3 + src/main/resources/pmd/pmd-ruleset.xml | 3 +- 7 files changed, 1300 insertions(+), 2 deletions(-) diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/TetrahedronSamplerBenchmark.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/TetrahedronSamplerBenchmark.java new file mode 100644 index 0000000..490a850 --- /dev/null +++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/TetrahedronSamplerBenchmark.java @@ -0,0 +1,490 @@ +/* + * 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.shape; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.UnitSphereSampler; +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.Level; +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 a tetrahedron. + */ +@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", "-Xms512M", "-Xmx512M" }) +public class TetrahedronSamplerBenchmark { + /** Name for the baseline method. */ + private static final String BASELINE = "Baseline"; + /** Name for the method using array coordinates. */ + private static final String ARRAY = "Array"; + /** Name for the method using non-array (primitive) coordinates. */ + private static final String NON_ARRAY = "NonArray"; + /** Name for the method using array coordinates and inline sample algorithm. */ + private static final String ARRAY_INLINE = "ArrayInline"; + /** Name for the method using non-array (primitive) coordinates and inline sample algorithm. */ + private static final String NON_ARRAY_INLINE = "NonArrayInline"; + /** 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 a random Cartesian point within the tetrahedron. + */ + double[] sample(); + } + + /** + * The base class for sampling from a tetrahedron. + * + * <ul> + * <li> + * Uses the algorithm described in: + * <blockquote> + * Rocchini, C. and Cignoni, P. (2001)<br> + * <i>Generating Random Points in a Tetrahedron</i>.<br> + * <strong>Journal of Graphics Tools</strong> 5(4), pp. 9-12. + * </blockquote> + * </li> + * </ul> + * + * @see <a href="https://doi.org/10.1080/10867651.2000.10487528";> + * Rocchini, C. & Cignoni, P. (2001) Journal of Graphics Tools 5, pp. 9-12</a> + */ + private abstract static class TetrahedronSampler implements Sampler { + /** The source of randomness. */ + private final UniformRandomProvider rng; + + /** + * @param rng Source of randomness. + */ + TetrahedronSampler(UniformRandomProvider rng) { + this.rng = rng; + } + + @Override + public double[] sample() { + double s = rng.nextDouble(); + double t = rng.nextDouble(); + final double u = rng.nextDouble(); + // Care is taken to ensure the 3 deviates remain in the 2^53 dyadic rationals in [0, 1). + // The following are exact for all the 2^53 dyadic rationals: + // 1 - u; u in [0, 1] + // u - 1; u in [0, 1] + // u + 1; u in [-1, 0] + // u + v; u in [-1, 0], v in [0, 1] + // u + v; u, v in [0, 1], u + v <= 1 + + // Cut and fold with the plane s + t = 1 + if (s + t > 1) { + // (s, t, u) = (1 - s, 1 - t, u) if s + t > 1 + s = 1 - s; + t = 1 - t; + } + // Now s + t <= 1. + // Cut and fold with the planes t + u = 1 and s + t + u = 1. + final double tpu = t + u; + final double sptpu = s + tpu; + if (sptpu > 1) { + if (tpu > 1) { + // (s, t, u) = (s, 1 - u, 1 - s - t) if t + u > 1 + // 1 - s - (1-u) - (1-s-t) == u - 1 + t + return createSample(u - 1 + t, s, 1 - u, 1 - s - t); + } + // (s, t, u) = (1 - t - u, t, s + t + u - 1) if t + u <= 1 + // 1 - (1-t-u) - t - (s+t+u-1) == 1 - s - t + return createSample(1 - s - t, 1 - tpu, t, s - 1 + tpu); + } + return createSample(1 - sptpu, s, t, u); + } + + /** + * Creates the sample given the random variates {@code s}, {@code t} and {@code u} in the + * interval {@code [0, 1]} and {@code s + t + u <= 1}. The sum {@code 1 - s - t - u} is + * provided. The sample can be obtained from the tetrahedron {@code abcd} using: + * + * <pre> + * p = (1 - s - t - u)a + sb + tc + ud + * </pre> + * + * @param p1msmtmu plus 1 minus s minus t minus u (1 - s - t - u) + * @param s the first variate s + * @param t the second variate t + * @param u the third variate u + * @return the sample + */ + protected abstract double[] createSample(double p1msmtmu, double s, double t, double u); + } + + /** + * Sample from a tetrahedron using array coordinates. + */ + private static class ArrayTetrahedronSampler extends TetrahedronSampler { + // CHECKSTYLE: stop JavadocVariableCheck + private final double[] a; + private final double[] b; + private final double[] c; + private final double[] d; + // CHECKSTYLE: resume JavadocVariableCheck + + /** + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @param d The fourth vertex. + * @param rng Source of randomness. + */ + ArrayTetrahedronSampler(double[] a, double[] b, double[] c, double[] d, + UniformRandomProvider rng) { + super(rng); + this.a = a.clone(); + this.b = b.clone(); + this.c = c.clone(); + this.d = d.clone(); + } + + @Override + protected double[] createSample(double p1msmtmu, double s, double t, double u) { + return new double[] {p1msmtmu * a[0] + s * b[0] + t * c[0] + u * d[0], + p1msmtmu * a[1] + s * b[1] + t * c[1] + u * d[1], + p1msmtmu * a[2] + s * b[2] + t * c[2] + u * d[2]}; + } + } + + /** + * Sample from a tetrahedron using non-array coordinates. + */ + private static class NonArrayTetrahedronSampler extends TetrahedronSampler { + // CHECKSTYLE: stop JavadocVariableCheck + private final double ax; + private final double bx; + private final double cx; + private final double dx; + private final double ay; + private final double by; + private final double cy; + private final double dy; + private final double az; + private final double bz; + private final double cz; + private final double dz; + // CHECKSTYLE: resume JavadocVariableCheck + + /** + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @param d The fourth vertex. + * @param rng Source of randomness. + */ + NonArrayTetrahedronSampler(double[] a, double[] b, double[] c, double[] d, + UniformRandomProvider rng) { + super(rng); + ax = a[0]; + ay = a[1]; + az = a[2]; + bx = b[0]; + by = b[1]; + bz = b[2]; + cx = c[0]; + cy = c[1]; + cz = c[2]; + dx = d[0]; + dy = d[1]; + dz = d[2]; + } + + @Override + protected double[] createSample(double p1msmtmu, double s, double t, double u) { + return new double[] {p1msmtmu * ax + s * bx + t * cx + u * dx, + p1msmtmu * ay + s * by + t * cy + u * dy, + p1msmtmu * az + s * bz + t * cz + u * dz}; + } + } + + /** + * Sample from a tetrahedron using array coordinates with an inline sample algorithm + * in-place of a method call. + */ + private static class ArrayInlineTetrahedronSampler implements Sampler { + // CHECKSTYLE: stop JavadocVariableCheck + private final double[] a; + private final double[] b; + private final double[] c; + private final double[] d; + private final UniformRandomProvider rng; + // CHECKSTYLE: resume JavadocVariableCheck + + /** + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @param d The fourth vertex. + * @param rng Source of randomness. + */ + ArrayInlineTetrahedronSampler(double[] a, double[] b, double[] c, double[] d, + UniformRandomProvider rng) { + this.a = a.clone(); + this.b = b.clone(); + this.c = c.clone(); + this.d = d.clone(); + this.rng = rng; + } + + @Override + public double[] sample() { + double s = rng.nextDouble(); + double t = rng.nextDouble(); + double u = rng.nextDouble(); + + if (s + t > 1) { + // (s, t, u) = (1 - s, 1 - t, u) if s + t > 1 + s = 1 - s; + t = 1 - t; + } + + double p1msmtmu; + final double tpu = t + u; + final double sptpu = s + tpu; + if (sptpu > 1) { + if (tpu > 1) { + // (s, t, u) = (s, 1 - u, 1 - s - t) if t + u > 1 + final double tt = t; + p1msmtmu = u - 1 + t; + t = 1 - u; + u = 1 - s - tt; + } else { + // (s, t, u) = (1 - t - u, t, s + t + u - 1) if t + u <= 1 + final double ss = s; + p1msmtmu = 1 - s - t; + s = 1 - tpu; + u = ss - 1 + tpu; + } + } else { + p1msmtmu = 1 - sptpu; + } + + return new double[] {p1msmtmu * a[0] + s * b[0] + t * c[0] + u * d[0], + p1msmtmu * a[1] + s * b[1] + t * c[1] + u * d[1], + p1msmtmu * a[2] + s * b[2] + t * c[2] + u * d[2]}; + } + } + + /** + * Sample from a tetrahedron using non-array coordinates with an inline sample algorithm + * in-place of a method call. + */ + private static class NonArrayInlineTetrahedronSampler implements Sampler { + // CHECKSTYLE: stop JavadocVariableCheck + private final double ax; + private final double bx; + private final double cx; + private final double dx; + private final double ay; + private final double by; + private final double cy; + private final double dy; + private final double az; + private final double bz; + private final double cz; + private final double dz; + private final UniformRandomProvider rng; + // CHECKSTYLE: resume JavadocVariableCheck + + /** + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @param d The fourth vertex. + * @param rng Source of randomness. + */ + NonArrayInlineTetrahedronSampler(double[] a, double[] b, double[] c, double[] d, + UniformRandomProvider rng) { + ax = a[0]; + ay = a[1]; + az = a[2]; + bx = b[0]; + by = b[1]; + bz = b[2]; + cx = c[0]; + cy = c[1]; + cz = c[2]; + dx = d[0]; + dy = d[1]; + dz = d[2]; + this.rng = rng; + } + + @Override + public double[] sample() { + double s = rng.nextDouble(); + double t = rng.nextDouble(); + double u = rng.nextDouble(); + + if (s + t > 1) { + // (s, t, u) = (1 - s, 1 - t, u) if s + t > 1 + s = 1 - s; + t = 1 - t; + } + + double p1msmtmu; + final double tpu = t + u; + final double sptpu = s + tpu; + if (sptpu > 1) { + if (tpu > 1) { + // (s, t, u) = (s, 1 - u, 1 - s - t) if t + u > 1 + final double tt = t; + p1msmtmu = u - 1 + t; + t = 1 - u; + u = 1 - s - tt; + } else { + // (s, t, u) = (1 - t - u, t, s + t + u - 1) if t + u <= 1 + final double ss = s; + p1msmtmu = 1 - s - t; + s = 1 - tpu; + u = ss - 1 + tpu; + } + } else { + p1msmtmu = 1 - sptpu; + } + + return new double[] {p1msmtmu * ax + s * bx + t * cx + u * dx, + p1msmtmu * ay + s * by + t * cy + u * dy, + p1msmtmu * az + s * bz + t * cz + u * dz}; + } + } + + /** + * Contains the sampler and the number of samples. + */ + @State(Scope.Benchmark) + public static class SamplerData { + /** The sampler. */ + private Sampler sampler; + + /** The number of samples. */ + @Param({"1", "10", "100", "1000", "10000"}) + private int size; + + /** The sampler type. */ + @Param({BASELINE, ARRAY, NON_ARRAY, ARRAY_INLINE, NON_ARRAY_INLINE}) + private String type; + + /** + * 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 and the sampler. + */ + @Setup(Level.Iteration) + public void setup() { + // This could be configured using @Param + final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_SHI_RO_256_PP); + final UnitSphereSampler s = UnitSphereSampler.of(3, rng); + final double[] a = s.nextVector(); + final double[] b = s.nextVector(); + final double[] c = s.nextVector(); + final double[] d = s.nextVector(); + sampler = createSampler(a, b, c, d, rng); + } + + /** + * Creates the tetrahedron sampler. + * + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @param d The fourth vertex. + * @param rng the source of randomness + * @return the sampler + */ + private Sampler createSampler(double[] a, double[] b, double[] c, double[] d, + final UniformRandomProvider rng) { + if (BASELINE.equals(type)) { + return new Sampler() { + @Override + public double[] sample() { + final double s = rng.nextDouble(); + final double t = rng.nextDouble(); + final double u = rng.nextDouble(); + return new double[] {s, t, u}; + } + }; + } else if (ARRAY.equals(type)) { + return new ArrayTetrahedronSampler(a, b, c, d, rng); + } else if (NON_ARRAY.equals(type)) { + return new NonArrayTetrahedronSampler(a, b, c, d, rng); + } else if (ARRAY_INLINE.equals(type)) { + return new ArrayInlineTetrahedronSampler(a, b, c, d, rng); + } else if (NON_ARRAY_INLINE.equals(type)) { + return new NonArrayInlineTetrahedronSampler(a, b, c, d, rng); + } + throw new IllegalStateException(UNKNOWN_SAMPLER + type); + } + } + + /** + * Run the sampler for the configured number of samples. + * + * @param bh Data sink + * @param data Input data. + */ + @Benchmark + public void sample(Blackhole bh, SamplerData data) { + final Sampler sampler = data.getSampler(); + for (int i = data.getSize() - 1; i >= 0; i--) { + bh.consume(sampler.sample()); + } + } +} diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/Coordinates.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/Coordinates.java index d729e66..a1f1fc2 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/Coordinates.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/Coordinates.java @@ -41,7 +41,7 @@ final class Coordinates { * @param values the values * @param message the message detail to prepend to the message in the event an exception is thrown * @return the values - * @throws IllegalArgumentException If a non-finite value is found + * @throws IllegalArgumentException if a non-finite value is found */ static double[] requireFinite(double[] values, String message) { for (final double value : values) { @@ -62,4 +62,30 @@ final class Coordinates { private static boolean isFinite(double value) { return Math.abs(value) <= Double.MAX_VALUE; } + + /** + * Check that the values is the specified length. This method is primarily for + * parameter validation in methods and constructors, for example: + * + * <pre> + * public Square(double[] topLeft, double[] bottomRight) { + * this.topLeft = Coordinates.requireLength(topLeft, 2, "topLeft"); + * this.bottomRight = Coordinates.requireLength(bottomRight, 2, "bottomRight"); + * } + * </pre> + * + * @param values the values + * @param length the length + * @param message the message detail to prepend to the message in the event an + * exception is thrown + * @return the values + * @throws IllegalArgumentException if the array length is not the specified length + */ + static double[] requireLength(double[] values, int length, String message) { + if (values.length != length) { + throw new IllegalArgumentException(String.format("%s length mismatch: %d != %d", + message, values.length, length)); + } + return values; + } } diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TetrahedronSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TetrahedronSampler.java new file mode 100644 index 0000000..00128c3 --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TetrahedronSampler.java @@ -0,0 +1,203 @@ +/* + * 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.shape; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.SharedStateSampler; + +/** + * Generate points uniformly distributed within a + * <a href="https://en.wikipedia.org/wiki/Tetrahedron";>tetrahedron</a>. + * + * <ul> + * <li> + * Uses the algorithm described in: + * <blockquote> + * Rocchini, C. and Cignoni, P. (2001)<br> + * <i>Generating Random Points in a Tetrahedron</i>.<br> + * <strong>Journal of Graphics Tools</strong> 5(4), pp. 9-12. + * </blockquote> + * </li> + * </ul> + * + * <p>Sampling uses:</p> + * + * <ul> + * <li>{@link UniformRandomProvider#nextDouble()} + * </ul> + * + * @see <a href="https://doi.org/10.1080/10867651.2000.10487528";> + * Rocchini, C. & Cignoni, P. (2001) Journal of Graphics Tools 5, pp. 9-12</a> + * @since 1.4 + */ +public class TetrahedronSampler implements SharedStateSampler<TetrahedronSampler> { + /** The dimension for 3D sampling. */ + private static final int THREE_D = 3; + /** The name of vertex a. */ + private static final String VERTEX_A = "Vertex a"; + /** The name of vertex b. */ + private static final String VERTEX_B = "Vertex b"; + /** The name of vertex c. */ + private static final String VERTEX_C = "Vertex c"; + /** The name of vertex d. */ + private static final String VERTEX_D = "Vertex d"; + + /** The first vertex. */ + private final double[] a; + /** The second vertex. */ + private final double[] b; + /** The third vertex. */ + private final double[] c; + /** The fourth vertex. */ + private final double[] d; + /** The source of randomness. */ + private final UniformRandomProvider rng; + + /** + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @param d The fourth vertex. + * @param rng Source of randomness. + */ + TetrahedronSampler(double[] a, double[] b, double[] c, double[] d, UniformRandomProvider rng) { + // Defensive copy + this.a = a.clone(); + this.b = b.clone(); + this.c = c.clone(); + this.d = d.clone(); + this.rng = rng; + } + + /** + * @param rng Generator of uniformly distributed random numbers + * @param source Source to copy. + */ + TetrahedronSampler(UniformRandomProvider rng, TetrahedronSampler source) { + // Shared state is immutable + a = source.a; + b = source.b; + c = source.c; + d = source.d; + this.rng = rng; + } + + /** + * @return a random Cartesian point within the tetrahedron. + */ + public double[] sample() { + double s = rng.nextDouble(); + double t = rng.nextDouble(); + final double u = rng.nextDouble(); + // Care is taken to ensure the 3 deviates remain in the 2^53 dyadic rationals in [0, 1). + // The following are exact for all the 2^53 dyadic rationals: + // 1 - u; u in [0, 1] + // u - 1; u in [0, 1] + // u + 1; u in [-1, 0] + // u + v; u in [-1, 0], v in [0, 1] + // u + v; u, v in [0, 1], u + v <= 1 + + // Cut and fold with the plane s + t = 1 + if (s + t > 1) { + // (s, t, u) = (1 - s, 1 - t, u) if s + t > 1 + s = 1 - s; + t = 1 - t; + } + // Now s + t <= 1. + // Cut and fold with the planes t + u = 1 and s + t + u = 1. + final double tpu = t + u; + final double sptpu = s + tpu; + if (sptpu > 1) { + if (tpu > 1) { + // (s, t, u) = (s, 1 - u, 1 - s - t) if t + u > 1 + // 1 - s - (1-u) - (1-s-t) == u - 1 + t + return createSample(u - 1 + t, s, 1 - u, 1 - s - t); + } + // (s, t, u) = (1 - t - u, t, s + t + u - 1) if t + u <= 1 + // 1 - (1-t-u) - t - (s+t+u-1) == 1 - s - t + return createSample(1 - s - t, 1 - tpu, t, s - 1 + tpu); + } + return createSample(1 - sptpu, s, t, u); + } + + /** + * Creates the sample given the random variates {@code s}, {@code t} and {@code u} in the + * interval {@code [0, 1]} and {@code s + t + u <= 1}. The sum {@code 1 - s - t - u} is + * provided. The sample can be obtained from the tetrahedron {@code abcd} using: + * + * <pre> + * p = (1 - s - t - u)a + sb + tc + ud + * </pre> + * + * @param p1msmtmu plus 1 minus s minus t minus u ({@code 1 - s - t - u}) + * @param s the first variate s + * @param t the second variate t + * @param u the third variate u + * @return the sample + */ + private double[] createSample(double p1msmtmu, double s, double t, double u) { + // From the barycentric coordinates s,t,u create the point by moving along + // vectors ab, ac and ad. + // Here we do not compute the vectors and use the original vertices: + // p = a + s(b-a) + t(c-a) + u(d-a) + // = (1-s-t-u)a + sb + tc + ud + return new double[] {p1msmtmu * a[0] + s * b[0] + t * c[0] + u * d[0], + p1msmtmu * a[1] + s * b[1] + t * c[1] + u * d[1], + p1msmtmu * a[2] + s * b[2] + t * c[2] + u * d[2]}; + } + + /** {@inheritDoc} */ + @Override + public TetrahedronSampler withUniformRandomProvider(UniformRandomProvider rng) { + return new TetrahedronSampler(rng, this); + } + + /** + * Create a tetrahedron sampler with vertices {@code a}, {@code b}, {@code c} and {@code d}. + * Sampled points are uniformly distributed within the tetrahedron. + * + * <p>No test for a volume is performed. If the vertices are coplanar the sampling + * distribution is undefined. + * + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @param d The fourth vertex. + * @param rng Source of randomness. + * @return the sampler + * @throws IllegalArgumentException If the vertices do not have length 3; + * or vertices have non-finite coordinates + */ + public static TetrahedronSampler of(double[] a, + double[] b, + double[] c, + double[] d, + UniformRandomProvider rng) { + // Must be 3D + Coordinates.requireLength(a, THREE_D, VERTEX_A); + Coordinates.requireLength(b, THREE_D, VERTEX_B); + Coordinates.requireLength(c, THREE_D, VERTEX_C); + Coordinates.requireLength(d, THREE_D, VERTEX_D); + // Detect non-finite vertices + Coordinates.requireFinite(a, VERTEX_A); + Coordinates.requireFinite(b, VERTEX_B); + Coordinates.requireFinite(c, VERTEX_C); + Coordinates.requireFinite(d, VERTEX_D); + return new TetrahedronSampler(a, b, c, d, rng); + } +} diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/CoordinatesTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/CoordinatesTest.java index 27c8ea0..cbfc336 100644 --- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/CoordinatesTest.java +++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/CoordinatesTest.java @@ -46,4 +46,28 @@ public class CoordinatesTest { } } } + + /** + * Test {@link Coordinates#requireLength(double[], int, String)} detects invalid lengths. + */ + @Test + public void testRequireLengthWithMessageThrows() { + final String message = "This should be prepended"; + for (final double[] c : new double[][] {{0, 1}, {0, 1, 2}}) { + final int length = c.length; + Assert.assertSame(c, Coordinates.requireLength(c, length, message)); + try { + Coordinates.requireLength(c, length - 1, message); + Assert.fail("Did not detect length was too long: " + (length - 1)); + } catch (IllegalArgumentException ex) { + Assert.assertTrue("Missing message prefix", ex.getMessage().startsWith(message)); + } + try { + Coordinates.requireLength(c, length + 1, message); + Assert.fail("Did not detect length was too short: " + (length + 1)); + } catch (IllegalArgumentException ex) { + Assert.assertTrue("Missing message prefix", ex.getMessage().startsWith(message)); + } + } + } } diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TetrahedronSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TetrahedronSamplerTest.java new file mode 100644 index 0000000..e879b2b --- /dev/null +++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TetrahedronSamplerTest.java @@ -0,0 +1,551 @@ +/* + * 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.shape; + +import java.util.Arrays; + +import org.junit.Assert; +import org.junit.Test; + +import org.apache.commons.math3.stat.inference.ChiSquareTest; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.RandomAssert; +import org.apache.commons.rng.simple.RandomSource; + +/** + * Test for {@link TetrahedronSampler}. + */ +public class TetrahedronSamplerTest { + /** + * Test invalid vertex dimensions (i.e. not 3D coordinates). + */ + @Test + public void testInvalidDimensionThrows() { + final UniformRandomProvider rng = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); + final double[] ok = new double[3]; + final double[] bad = new double[2]; + final double[][] c = {ok, ok, ok, ok}; + for (int i = 0; i < c.length; i++) { + c[i] = bad; + try { + TetrahedronSampler.of(c[0], c[1], c[2], c[3], rng); + Assert.fail(String.format("Did not detect invalid dimension for vertex: %d", i)); + } catch (IllegalArgumentException ex) { + // Expected + } + c[i] = ok; + } + } + + /** + * Test non-finite vertices. + */ + @Test + public void testNonFiniteVertexCoordinates() { + final UniformRandomProvider rng = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); + // A valid tetrahedron + final double[][] c = new double[][] { + {1, 1, 1}, {1, -1, 1}, {-1, 1, 1}, {1, 1, -1} + }; + Assert.assertNotNull(TetrahedronSampler.of(c[0], c[1], c[2], c[3], rng)); + final double[] bad = {Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NaN}; + for (int i = 0; i < c.length; i++) { + for (int j = 0; j < c[0].length; j++) { + for (final double d : bad) { + final double value = c[i][j]; + c[i][j] = d; + try { + TetrahedronSampler.of(c[0], c[1], c[2], c[3], rng); + Assert.fail(String.format("Did not detect non-finite coordinate: %d,%d = %s", i, j, d)); + } catch (IllegalArgumentException ex) { + // Expected + } + c[i][j] = value; + } + } + } + } + + /** + * Test a tetrahedron with coordinates that are separated by more than + * {@link Double#MAX_VALUE}. + */ + @Test + public void testExtremeValueCoordinates() { + // Object seed so use Long not long + final Long seed = 876543L; + // Create a valid tetrahedron that can be scaled + final double[][] c1 = new double[][] { + {1, 1, 1}, {1, -1, -1}, {-1, -1, 1}, {-1, 1, -1} + }; + final double[][] c2 = new double[4][3]; + // Extremely large value for scaling. Use a power of 2 for exact scaling. + final double scale = 0x1.0p1023; + for (int i = 0; i < c1.length; i++) { + // Scale the second tetrahedron + for (int j = 0; j < 3; j++) { + c2[i][j] = c1[i][j] * scale; + } + } + // Show the tetrahedron is too big to compute vectors between points. + Assert.assertEquals("Expect vector c - a to be infinite in the x dimension", + Double.NEGATIVE_INFINITY, c2[2][0] - c2[0][0], 0.0); + Assert.assertEquals("Expect vector c - d to be infinite in the y dimension", + Double.NEGATIVE_INFINITY, c2[2][1] - c2[3][1], 0.0); + Assert.assertEquals("Expect vector c - b to be infinite in the z dimension", + Double.POSITIVE_INFINITY, c2[2][2] - c2[1][2], 0.0); + + final TetrahedronSampler sampler1 = TetrahedronSampler.of(c1[0], c1[1], c1[2], c1[3], + RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP, seed)); + final TetrahedronSampler sampler2 = TetrahedronSampler.of(c2[0], c2[1], c2[2], c2[3], + RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP, seed)); + + for (int n = 0; n < 10; n++) { + final double[] a = sampler1.sample(); + final double[] b = sampler2.sample(); + for (int i = 0; i < a.length; i++) { + a[i] *= scale; + } + Assert.assertArrayEquals(a, b, 0.0); + } + } + + /** + * Test the distribution of points in three dimensions. 6 tetrahedra are used to create + * a box. The distribution should be uniform inside the box. + */ + @Test + public void testDistribution() { + // Create the lower and upper limits of the box + final double lx = -1; + final double ly = -2; + final double lz = 1; + final double ux = 3; + final double uy = 4; + final double uz = 5; + // Create vertices abcd and efgh for the lower and upper rectangles + // (in the XY plane) of the box + final double[] a = {lx, ly, lz}; + final double[] b = {ux, ly, lz}; + final double[] c = {ux, uy, lz}; + final double[] d = {lx, uy, lz}; + final double[] e = {lx, ly, uz}; + final double[] f = {ux, ly, uz}; + final double[] g = {ux, uy, uz}; + final double[] h = {lx, uy, uz}; + + // Assign bins + final int bins = 10; + // Samples should be a multiple of 6 (due to combining 6 equal volume tetrahedra) + final int samplesPerBin = 12; + // Scale factors to assign x,y,z to a bin + final double sx = bins / (ux - lx); + final double sy = bins / (uy - ly); + final double sz = bins / (uz - lz); + + // Compute factor to allocate bin index: + // index = x + y * binsX + z * binsX * binsY + final int binsXy = bins * bins; + + // Expect a uniform distribution (this is rescaled by the ChiSquareTest) + final double[] expected = new double[binsXy * bins]; + Arrays.fill(expected, 1); + + // 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, 0xaabbccddeeffL); + + // Cut the box into 6 equal volume tetrahedra by cutting the box in half three times, + // cutting diagonally through each of the three pairs of opposing faces. In this way, + // the tetrahedra all share one of the main diagonals of the box (d-f). + // See the cuts used for the marching tetrahedra algorithm: + // https://en.wikipedia.org/wiki/Marching_tetrahedra + final TetrahedronSampler[] samplers = {TetrahedronSampler.of(d, f, b, c, rng), + TetrahedronSampler.of(d, f, c, g, rng), + TetrahedronSampler.of(d, f, g, h, rng), + TetrahedronSampler.of(d, f, h, e, rng), + TetrahedronSampler.of(d, f, e, a, rng), + TetrahedronSampler.of(d, f, a, b, rng)}; + // To determine the sample is inside the correct tetrahedron it is projected to the + // 4 faces of the tetrahedron along the face normals. The distance should be negative + // when the face normals are orientated outwards. + final InsideTetrahedron[] insides = {new InsideTetrahedron(d, f, b, c), + new InsideTetrahedron(d, f, c, g), + new InsideTetrahedron(d, f, g, h), + new InsideTetrahedron(d, f, h, e), + new InsideTetrahedron(d, f, e, a), + new InsideTetrahedron(d, f, a, b)}; + + final int samples = expected.length * samplesPerBin; + for (int n = 0; n < 1; n++) { + // Assign each coordinate to a region inside the combined box + final long[] observed = new long[expected.length]; + // Equal volume tetrahedra so sample from each one + for (int i = 0; i < samples; i += 6) { + for (int j = 0; j < 6; j++) { + addObservation(samplers[j].sample(), observed, bins, binsXy, + lx, ly, lz, sx, sy, sz, insides[j]); + } + } + final double p = new ChiSquareTest().chiSquareTest(expected, observed); + Assert.assertFalse("p-value too small: " + p, p < 0.001); + } + } + + /** + * Adds the observation. Coordinates are mapped using the offsets, scaled and + * then cast to an integer bin. + * + * <pre> + * binx = (int) ((x - lx) * sx) + * </pre> + * + * @param v the sample (3D coordinate xyz) + * @param observed the observations + * @param binsX the numbers of bins in the x dimension + * @param binsXy the numbers of bins in the combined x and y dimensions + * @param lx the lower limit to convert the x coordinate to the x bin + * @param ly the lower limit to convert the y coordinate to the y bin + * @param lz the lower limit to convert the z coordinate to the z bin + * @param sx the scale to convert the x coordinate to the x bin + * @param sy the scale to convert the y coordinate to the y bin + * @param sz the scale to convert the z coordinate to the z bin + * @param inside the inside tetrahedron test + */ + // CHECKSTYLE: stop ParameterNumberCheck + private static void addObservation(double[] v, long[] observed, + int binsX, int binsXy, + double lx, double ly, double lz, + double sx, double sy, double sz, + InsideTetrahedron inside) { + Assert.assertEquals(3, v.length); + // Test the point is inside the correct tetrahedron + Assert.assertTrue("Not inside the tetrahedron", inside.test(v)); + final double x = v[0]; + final double y = v[1]; + final double z = v[2]; + // Add to the correct bin after using the offset + final int binx = (int) ((x - lx) * sx); + final int biny = (int) ((y - ly) * sy); + final int binz = (int) ((z - lz) * sz); + observed[binz * binsXy + biny * binsX + binx]++; + } + // CHECKSTYLE: resume ParameterNumberCheck + + /** + * Test the SharedStateSampler implementation. + */ + @Test + public void testSharedStateSampler() { + final UniformRandomProvider rng1 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); + final UniformRandomProvider rng2 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); + final double[] c1 = createCoordinate(-1); + final double[] c2 = createCoordinate(2); + final double[] c3 = createCoordinate(-3); + final double[] c4 = createCoordinate(4); + final TetrahedronSampler sampler1 = TetrahedronSampler.of(c1, c2, c3, c4, rng1); + final TetrahedronSampler 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(); + } + }); + } + + /** + * Test the input vectors are copied and not used by reference. + */ + @Test + public void testChangedInputCoordinates() { + final UniformRandomProvider rng1 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); + final UniformRandomProvider rng2 = RandomSource.create(RandomSource.SPLIT_MIX_64, 0L); + final double[] c1 = createCoordinate(1); + final double[] c2 = createCoordinate(2); + final double[] c3 = createCoordinate(-3); + final double[] c4 = createCoordinate(-4); + final TetrahedronSampler sampler1 = TetrahedronSampler.of(c1, c2, c3, c4, rng1); + // Check the input vectors are copied and not used by reference. + // Change them in place and create a new sampler. It should have different output + // translated by the offset. + final double offset = 10; + for (int i = 0; i < 3; i++) { + c1[i] += offset; + c2[i] += offset; + c3[i] += offset; + c4[i] += offset; + } + final TetrahedronSampler sampler2 = TetrahedronSampler.of(c1, c2, c3, c4, rng2); + for (int n = 0; n < 5; n++) { + final double[] s1 = sampler1.sample(); + final double[] s2 = sampler2.sample(); + Assert.assertEquals(3, s1.length); + Assert.assertEquals(3, s2.length); + for (int i = 0; i < 3; i++) { + Assert.assertEquals(s1[i] + offset, s2[i], 1e-10); + } + } + } + + /** + * Test the inside tetrahedron predicate. + */ + @Test + public void testInsideTetrahedron() { + final double[][] c1 = new double[][] { + {1, 1, 1}, {1, -1, -1}, {-1, -1, 1}, {-1, 1, -1} + }; + final InsideTetrahedron inside = new InsideTetrahedron(c1[0], c1[1], c1[2], c1[3]); + // Testing points on the vertices, edges or faces are subject to floating point error + final double epsilon = 1e-14; + // Vertices + for (int i = 0; i < 4; i++) { + Assert.assertTrue(inside.test(c1[i], epsilon)); + } + // Edge + Assert.assertTrue(inside.test(new double[] {1, 0, 0}, epsilon)); + Assert.assertTrue(inside.test(new double[] {0.5, 0.5, 1}, epsilon)); + // Just inside the edge + Assert.assertTrue(inside.test(new double[] {1 - 1e-10, 0, 0})); + Assert.assertTrue(inside.test(new double[] {0.5, 0.5, 1 - 1e-10})); + // Just outside the edge + Assert.assertFalse(inside.test(new double[] {1, 0, 1e-10}, epsilon)); + Assert.assertFalse(inside.test(new double[] {0.5, 0.5 + 1e-10, 1}, epsilon)); + // Face + double x = 1.0 / 3; + Assert.assertTrue(inside.test(new double[] {x, -x, x}, epsilon)); + Assert.assertTrue(inside.test(new double[] {-x, -x, -x}, epsilon)); + Assert.assertTrue(inside.test(new double[] {x, x, -x}, epsilon)); + Assert.assertTrue(inside.test(new double[] {-x, x, x}, epsilon)); + // Just outside the face + x += 1e-10; + Assert.assertFalse(inside.test(new double[] {x, -x, x}, epsilon)); + Assert.assertFalse(inside.test(new double[] {-x, -x, -x}, epsilon)); + Assert.assertFalse(inside.test(new double[] {x, x, -x}, epsilon)); + Assert.assertFalse(inside.test(new double[] {-x, x, x}, epsilon)); + // Inside + Assert.assertTrue(inside.test(new double[] {0, 0, 0})); + Assert.assertTrue(inside.test(new double[] {0.5, 0.25, -0.1})); + // Outside + Assert.assertFalse(inside.test(new double[] {0, 20, 0})); + Assert.assertFalse(inside.test(new double[] {-20, 0, 0})); + Assert.assertFalse(inside.test(new double[] {6, 6, 4})); + } + + /** + * Creates the coordinate of length 3 filled with + * the given value and the dimension index: x + i. + * + * @param x the value for index 0 + * @return the coordinate + */ + private static double[] createCoordinate(double x) { + final double[] coord = new double[3]; + for (int i = 0; i < 3; i++) { + coord[0] = x + i; + } + return coord; + } + + /** + * Class to test if a point is inside the tetrahedron. + * + * <p>Computes the outer pointing face normals for the tetrahedron. A point is inside + * if the point lies below each of the face planes of the shape. + * + * @see <a href="https://mathworld.wolfram.com/Point-PlaneDistance.html";>Point-Plane distance</a> + */ + private static class InsideTetrahedron { + /** The face normals. */ + private final double[][] n; + /** The distance of each face from the origin. */ + private final double[] d; + + /** + * Create an instance. + * + * @param v1 The first vertex. + * @param v2 The second vertex. + * @param v3 The third vertex. + * @param v4 The fourth vertex. + */ + InsideTetrahedron(double[] v1, double[] v2, double[] v3, double[] v4) { + // Compute the centre of each face + final double[][] x = new double[][] { + centre(v1, v2, v3), + centre(v2, v3, v4), + centre(v3, v4, v1), + centre(v4, v1, v2) + }; + + // Compute the normal for each face + n = new double[][] { + normal(v1, v2, v3), + normal(v2, v3, v4), + normal(v3, v4, v1), + normal(v4, v1, v2) + }; + + // Given the plane: + // 0 = ax + by + cz + d + // Where abc is the face normal and d is the distance of the plane from the origin. + // Compute d: + // d = -(ax + by + cz) + d = new double[] { + -dot(n[0], x[0]), + -dot(n[1], x[1]), + -dot(n[2], x[2]), + -dot(n[3], x[3]), + }; + + // Compute the distance of the other vertex from each face plane. + // When below the distance should be negative. Orient each normal so this is true. + // + // This distance D of a point xyz to the plane is: + // D = ax + by + cz + d + // Above plane: + // ax + by + cz + d > 0 + // ax + by + cz > -d + final double[][] other = {v4, v1, v2, v3}; + for (int i = 0; i < 4; i++) { + if (dot(n[i], other[i]) > -d[i]) { + // Swap orientation + n[i][0] = -n[i][0]; + n[i][1] = -n[i][1]; + n[i][2] = -n[i][2]; + d[i] = -d[i]; + } + } + } + + /** + * Compute the centre of the triangle face. + * + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @return the centre + */ + private static double[] centre(double[] a, double[] b, double[] c) { + return new double[] { + (a[0] + b[0] + c[0]) / 3, + (a[1] + b[1] + c[1]) / 3, + (a[2] + b[2] + c[2]) / 3 + }; + } + + /** + * Compute the normal of the triangle face. + * + * @param a The first vertex. + * @param b The second vertex. + * @param c The third vertex. + * @return the normal + */ + private static double[] normal(double[] a, double[] b, double[] c) { + final double[] v1 = subtract(b, a); + final double[] v2 = subtract(c, a); + // Cross product + final double[] normal = { + v1[1] * v2[2] - v1[2] * v2[1], + v1[2] * v2[0] - v1[0] * v2[2], + v1[0] * v2[1] - v1[1] * v2[0] + }; + // Normalise + final double scale = 1.0 / Math.sqrt(dot(normal, normal)); + normal[0] *= scale; + normal[1] *= scale; + normal[2] *= scale; + return normal; + } + + /** + * Compute the dot product of vector {@code a} and {@code b}. + * + * <pre> + * a.b = a.x * b.x + a.y * b.y + a.z * b.z + * </pre> + * + * @param a the first vector + * @param b the second vector + * @return the dot product + */ + private static double dot(double[] a, double[] b) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; + } + + /** + * Subtract the second term from the first: {@code a - b}. + * + * @param a The first term. + * @param b The second term. + * @return the vector {@code a - b} + */ + private static double[] subtract(double[] a, double[] b) { + return new double[] { + a[0] - b[0], + a[1] - b[1], + a[2] - b[2] + }; + } + + /** + * Check the point is inside the tetrahedron. + * + * @param x the coordinate + * @return true if inside the tetrahedron + */ + boolean test(double[] x) { + // Must be below all the face planes + for (int i = 0; i < 4; i++) { + // This distance D of a point xyz to the plane is: + // D = ax + by + cz + d + // Above plane: + // ax + by + cz + d > 0 + // ax + by + cz > -d + if (dot(n[i], x) > -d[i]) { + return false; + } + } + return true; + } + + /** + * Check the point is inside the tetrahedron within the given absolute epsilon. + * + * @param x the coordinate + * @param epsilon the epsilon + * @return true if inside the tetrahedron + */ + boolean test(double[] x, double epsilon) { + for (int i = 0; i < 4; i++) { + // As above but with an epsilon above zero + if (dot(n[i], x) > epsilon - d[i]) { + return false; + } + } + return true; + } + } +} diff --git a/src/changes/changes.xml b/src/changes/changes.xml index d323cfb..3c97cd0 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="135"> + New "TetrahedronSampler" to sample uniformly from a tetrahedron. + </action> <action dev="aherbert" type="add" issue="134"> New "BoxSampler" to sample uniformly from a box (or hyperrectangle). </action> diff --git a/src/main/resources/pmd/pmd-ruleset.xml b/src/main/resources/pmd/pmd-ruleset.xml index 9853424..2ef758b 100644 --- a/src/main/resources/pmd/pmd-ruleset.xml +++ b/src/main/resources/pmd/pmd-ruleset.xml @@ -84,7 +84,8 @@ <!-- Do not require Utils/Helper suffix --> <property name="violationSuppressXPath" value="//ClassOrInterfaceDeclaration[@SimpleName='ListSampler' or @SimpleName='ProviderBuilder' - or @SimpleName='ThreadLocalRandomSource' or @SimpleName='SeedFactory']"/> + or @SimpleName='ThreadLocalRandomSource' or @SimpleName='SeedFactory' + or @SimpleName='Coordinates']"/> <!-- Allow samplers to have only factory constructors --> <property name="utilityClassPattern" value="[A-Z][a-zA-Z0-9]+(Utils?|Helper|Sampler)" /> </properties>
