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:38 UTC
[commons-rng] branch master updated (838295a -> d365385)
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 838295a Javadoc: Use the same description for all shape samplers.
new b8f4d82 RNG-135: TetrahedronSampler to sample uniformly from a tetrahedron
new d365385 Organise imports
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:
.../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 +++++++++++++++++++++
.../rng/sampling/shape/TriangleSamplerTest.java | 7 +-
.../rng/sampling/shape/UnitBallSamplerTest.java | 8 +-
src/changes/changes.xml | 3 +
src/main/resources/pmd/pmd-ruleset.xml | 3 +-
9 files changed, 1310 insertions(+), 7 deletions(-)
create mode 100644 commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/TetrahedronSamplerBenchmark.java
create mode 100644 commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TetrahedronSampler.java
create mode 100644 commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TetrahedronSamplerTest.java
[commons-rng] 01/02: RNG-135: TetrahedronSampler to sample
uniformly from a tetrahedron
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 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>
[commons-rng] 02/02: Organise imports
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 d36538538e0dd425e365c8b9b2e39d93a4a10f41
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Tue May 11 22:43:37 2021 +0100
Organise imports
---
.../apache/commons/rng/sampling/shape/TriangleSamplerTest.java | 7 +++++--
.../apache/commons/rng/sampling/shape/UnitBallSamplerTest.java | 8 +++++---
2 files changed, 10 insertions(+), 5 deletions(-)
diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TriangleSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TriangleSamplerTest.java
index ea87a35..dfb1a2d 100644
--- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TriangleSamplerTest.java
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TriangleSamplerTest.java
@@ -16,13 +16,16 @@
*/
package org.apache.commons.rng.sampling.shape;
+import java.util.Arrays;
+
import org.junit.Assert;
import org.junit.Test;
-import org.apache.commons.rng.simple.RandomSource;
-import java.util.Arrays;
+
import org.apache.commons.math3.stat.inference.ChiSquareTest;
+
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.RandomAssert;
+import org.apache.commons.rng.simple.RandomSource;
/**
* Test for {@link TriangleSampler}.
diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/UnitBallSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/UnitBallSamplerTest.java
index 67c6bba..6e4f00c 100644
--- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/UnitBallSamplerTest.java
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/UnitBallSamplerTest.java
@@ -16,15 +16,17 @@
*/
package org.apache.commons.rng.sampling.shape;
+import java.util.Arrays;
+
import org.junit.Assert;
import org.junit.Test;
-import org.apache.commons.rng.simple.RandomSource;
-import java.util.Arrays;
+
import org.apache.commons.math3.stat.inference.ChiSquareTest;
+
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.core.source64.SplitMix64;
import org.apache.commons.rng.sampling.RandomAssert;
-import org.apache.commons.rng.sampling.RandomAssert.Sampler;
+import org.apache.commons.rng.simple.RandomSource;
/**
* Test for {@link UnitBallSampler}.