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/05/11 21:45:39 UTC
[commons-rng] 01/02: RNG-135: TetrahedronSampler to sample
uniformly from a tetrahedron
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 <ah...@apache.org>
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>