[prev in list] [next in list] [prev in thread] [next in thread]
List: jakarta-commons-dev
Subject: [commons-rng] branch master updated: RNG-131: TriangleSampler to sample uniformly from a triangle.
From: aherbert () apache ! org
Date: 2021-04-29 12:30:47
Message-ID: 161969944720.32213.15811547232796477103 () gitbox ! apache ! org
[Download RAW message or body]
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
The following commit(s) were added to refs/heads/master by this push:
new 6e346fe RNG-131: TriangleSampler to sample uniformly from a triangle.
6e346fe is described below
commit 6e346fefff881e9fcc0509a5a0af93283089c53a
Author: Alex Herbert <aherbert@apache.org>
AuthorDate: Wed Apr 21 22:32:04 2021 +0100
RNG-131: TriangleSampler to sample uniformly from a triangle.
Use a method that samples from any triangle with finite valued
coordinates.
A TriangleSampler JMH benchmark is included to evaluate different
sampling methods.
---
.../jmh/sampling/TriangleSamplerBenchmark.java | 672 ++++++++++++++++++
.../commons/rng/sampling/shape/Coordinates.java | 65 ++
.../rng/sampling/shape/TriangleSampler.java | 354 ++++++++++
.../rng/sampling/shape/CoordinatesTest.java | 49 ++
.../rng/sampling/shape/TriangleSamplerTest.java | 786 +++++++++++++++++++++
src/changes/changes.xml | 3 +
6 files changed, 1929 insertions(+)
diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/TriangleSamplerBenchmark.java \
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/TriangleSamplerBenchmark.java
new file mode 100644
index 0000000..68086a9
--- /dev/null
+++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/TriangleSamplerBenchmark.java
@@ -0,0 +1,672 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.rng.examples.jmh.sampling;
+
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.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 an \
N-dimension triangle. + */
+@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 TriangleSamplerBenchmark {
+ /** Name for the baseline method. */
+ private static final String BASELINE = "Baseline";
+ /** Name for the baseline method including a 50:50 if statement. */
+ private static final String BASELINE_IF = "BaselineIf";
+ /** Name for the method using vectors. */
+ private static final String VECTORS = "Vectors";
+ /** Name for the method using the coordinates. */
+ private static final String COORDINATES = "Coordinates";
+ /** Error message for an unknown sampler type. */
+ private static final String UNKNOWN_SAMPLER = "Unknown sampler type: ";
+
+ /**
+ * The sampler.
+ */
+ private interface Sampler {
+ /**
+ * Gets the next sample.
+ *
+ * @return the sample
+ */
+ double[] sample();
+ }
+
+ /**
+ * Base class for a sampler using vectors: {@code p = a + sv + tw}.
+ */
+ private abstract static class VectorTriangleSampler implements Sampler {
+ /** The source of randomness. */
+ private final UniformRandomProvider rng;
+
+ /**
+ * Create an instance.
+ *
+ * @param rng the source of randomness
+ */
+ VectorTriangleSampler(UniformRandomProvider rng) {
+ this.rng = rng;
+ }
+
+ /**
+ * @return a random Cartesian point inside the triangle.
+ */
+ @Override
+ public double[] sample() {
+ final double s = rng.nextDouble();
+ final double t = rng.nextDouble();
+ if (s + t > 1) {
+ return createSample(1.0 - s, 1.0 - t);
+ }
+ return createSample(s, t);
+ }
+
+ /**
+ * Creates the sample given the random variates {@code s} and {@code t} in \
the + * interval {@code [0, 1]} and {@code s + t <= 1}.
+ * The sample can be obtained from the triangle abc with v = b-a and w = c-a \
using: + * <pre>
+ * p = a + sv + tw
+ * </pre>
+ *
+ * @param s the first variate s
+ * @param t the second variate t
+ * @return the sample
+ */
+ protected abstract double[] createSample(double s, double t);
+ }
+
+ /**
+ * Base class for a sampler using coordinates: {@code a(1 - s - t) + sb + tc}.
+ */
+ private abstract static class CoordinateTriangleSampler implements Sampler {
+ /** The source of randomness. */
+ private final UniformRandomProvider rng;
+
+ /**
+ * Create an instance.
+ *
+ * @param rng the source of randomness
+ */
+ CoordinateTriangleSampler(UniformRandomProvider rng) {
+ this.rng = rng;
+ }
+
+ /**
+ * @return a random Cartesian point inside the triangle.
+ */
+ @Override
+ public double[] sample() {
+ final double s = rng.nextDouble();
+ final double t = rng.nextDouble();
+ final double spt = s + t;
+ if (spt > 1) {
+ // Transform: s1 = 1 - s; t1 = 1 - t.
+ // Compute: 1 - s1 - t1
+ // Do not assume (1 - (1-s) - (1-t)) is (s + t - 1), i.e. (spt - \
1.0), + // to avoid loss of a random bit due to rounding when s + t > \
1. + // An exact sum is (s - 1 + t).
+ return createSample(s - 1.0 + t, 1.0 - s, 1.0 - t);
+ }
+ // Here s + t is exact so can be subtracted to make 1 - s - t
+ return createSample(1.0 - spt, s, t);
+ }
+
+ /**
+ * Creates the sample given the random variates {@code s} and {@code t} in \
the + * interval {@code [0, 1]} and {@code s + t <= 1}. The sum {@code 1 - s \
- t} is provided. + * The sample can be obtained from the triangle abc using:
+ * <pre>
+ * p = a(1 - s - t) + sb + tc
+ * </pre>
+ *
+ * @param p1msmt plus 1 minus s minus t (1 - s - t)
+ * @param s the first variate s
+ * @param t the second variate t
+ * @return the sample
+ */
+ protected abstract double[] createSample(double p1msmt, double s, double t);
+ }
+
+ /**
+ * Base class for the sampler data.
+ * Contains the source of randomness and the number of samples.
+ * The sampler should be created by a sub-class of the data.
+ */
+ @State(Scope.Benchmark)
+ public abstract static class SamplerData {
+ /** The sampler. */
+ private Sampler sampler;
+
+ /** The number of samples. */
+ @Param({"1000"})
+ private int size;
+
+ /**
+ * Gets the size.
+ *
+ * @return the size
+ */
+ public int getSize() {
+ return size;
+ }
+
+ /**
+ * Gets the sampler.
+ *
+ * @return the sampler
+ */
+ public Sampler getSampler() {
+ return sampler;
+ }
+
+ /**
+ * Create the source of randomness and the sampler.
+ */
+ @Setup(Level.Iteration)
+ public void setup() {
+ // This could be configured using @Param
+ final UniformRandomProvider rng = \
RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP); + final int \
dimension = getDimension(); + final UnitSphereSampler s = \
UnitSphereSampler.of(dimension, rng); + final double[] a = s.nextVector();
+ final double[] b = s.nextVector();
+ final double[] c = s.nextVector();
+ sampler = createSampler(a, b, c, rng);
+ }
+
+ /**
+ * Gets the dimension of the triangle vertices.
+ *
+ * @return the dimension
+ */
+ protected abstract int getDimension();
+
+ /**
+ * Creates the triangle sampler.
+ *
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng the source of randomness
+ * @return the sampler
+ */
+ protected abstract Sampler createSampler(double[] a, double[] b, double[] c, \
UniformRandomProvider rng); + }
+
+ /**
+ * The 2D triangle sampler.
+ */
+ @State(Scope.Benchmark)
+ public static class Sampler2D extends SamplerData {
+ /** The sampler type. */
+ @Param({BASELINE, BASELINE_IF, VECTORS, COORDINATES})
+ private String type;
+
+ /** {@inheritDoc} */
+ @Override
+ protected int getDimension() {
+ return 2;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected Sampler createSampler(final double[] a, final double[] b, final \
double[] c, + final UniformRandomProvider rng) \
{ + if (BASELINE.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ final double s = rng.nextDouble();
+ final double t = rng.nextDouble();
+ return new double[] {s, t};
+ }
+ };
+ } else if (BASELINE_IF.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ final double s = rng.nextDouble();
+ final double t = rng.nextDouble();
+ if (s + t > 1) {
+ return new double[] {s, t};
+ }
+ return new double[] {t, s};
+ }
+ };
+ } else if (VECTORS.equals(type)) {
+ return new VectorTriangleSampler2D(a, b, c, rng);
+ } else if (COORDINATES.equals(type)) {
+ return new CoordinateTriangleSampler2D(a, b, c, rng);
+ }
+ throw new IllegalStateException(UNKNOWN_SAMPLER + type);
+ }
+
+ /**
+ * Sample vectors in 2D.
+ */
+ private static class VectorTriangleSampler2D extends VectorTriangleSampler {
+ // CHECKSTYLE: stop JavadocVariableCheck
+ private final double ax;
+ private final double ay;
+ private final double vx;
+ private final double vy;
+ private final double wx;
+ private final double wy;
+ // CHECKSTYLE: resume JavadocVariableCheck
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng the source of randomness
+ */
+ VectorTriangleSampler2D(double[] a, double[] b, double[] c, \
UniformRandomProvider rng) { + super(rng);
+ ax = a[0];
+ ay = a[1];
+ vx = b[0] - ax;
+ vy = b[1] - ay;
+ wx = c[0] - ax;
+ wy = c[1] - ay;
+ }
+
+ @Override
+ protected double[] createSample(double s, double t) {
+ return new double[] {ax + s * vx + t * wx,
+ ay + s * vy + t * wy};
+ }
+ }
+
+ /**
+ * Sample using coordinates in 2D.
+ */
+ private static class CoordinateTriangleSampler2D extends \
CoordinateTriangleSampler { + // CHECKSTYLE: stop JavadocVariableCheck
+ private final double ax;
+ private final double ay;
+ private final double bx;
+ private final double by;
+ private final double cx;
+ private final double cy;
+ // CHECKSTYLE: resume JavadocVariableCheck
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng the source of randomness
+ */
+ CoordinateTriangleSampler2D(double[] a, double[] b, double[] c, \
UniformRandomProvider rng) { + super(rng);
+ ax = a[0];
+ ay = a[1];
+ bx = b[0];
+ by = b[1];
+ cx = c[0];
+ cy = c[1];
+ }
+
+ @Override
+ protected double[] createSample(double p1msmt, double s, double t) {
+ return new double[] {p1msmt * ax + s * bx + t * cx,
+ p1msmt * ay + s * by + t * cy};
+ }
+ } }
+
+ /**
+ * The 3D triangle sampler.
+ */
+ @State(Scope.Benchmark)
+ public static class Sampler3D extends SamplerData {
+ /** The sampler type. */
+ @Param({BASELINE, BASELINE_IF, VECTORS, COORDINATES})
+ private String type;
+
+ /** {@inheritDoc} */
+ @Override
+ protected int getDimension() {
+ return 3;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected Sampler createSampler(final double[] a, final double[] b, final \
double[] c, + final UniformRandomProvider rng) \
{ + if (BASELINE.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ final double s = rng.nextDouble();
+ final double t = rng.nextDouble();
+ return new double[] {s, t, s};
+ }
+ };
+ } else if (BASELINE_IF.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ final double s = rng.nextDouble();
+ final double t = rng.nextDouble();
+ if (s + t > 1) {
+ return new double[] {s, t, s};
+ }
+ return new double[] {t, s, t};
+ }
+ };
+ } else if (VECTORS.equals(type)) {
+ return new VectorTriangleSampler3D(a, b, c, rng);
+ } else if (COORDINATES.equals(type)) {
+ return new CoordinateTriangleSampler3D(a, b, c, rng);
+ }
+ throw new IllegalStateException(UNKNOWN_SAMPLER + type);
+ }
+
+ /**
+ * Sample vectors in 3D.
+ */
+ private static class VectorTriangleSampler3D extends VectorTriangleSampler {
+ // CHECKSTYLE: stop JavadocVariableCheck
+ private final double ax;
+ private final double ay;
+ private final double az;
+ private final double vx;
+ private final double vy;
+ private final double vz;
+ private final double wx;
+ private final double wy;
+ private final double wz;
+ // CHECKSTYLE: resume JavadocVariableCheck
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng the source of randomness
+ */
+ VectorTriangleSampler3D(double[] a, double[] b, double[] c, \
UniformRandomProvider rng) { + super(rng);
+ ax = a[0];
+ ay = a[1];
+ az = a[2];
+ vx = b[0] - ax;
+ vy = b[1] - ay;
+ vz = b[2] - az;
+ wx = c[0] - ax;
+ wy = c[1] - ay;
+ wz = c[2] - az;
+ }
+
+ @Override
+ protected double[] createSample(double s, double t) {
+ return new double[] {ax + s * vx + t * wx,
+ ay + s * vy + t * wy,
+ az + s * vz + t * wz};
+ }
+ }
+
+ /**
+ * Sample using coordinates in 3D.
+ */
+ private static class CoordinateTriangleSampler3D extends \
CoordinateTriangleSampler { + // CHECKSTYLE: stop JavadocVariableCheck
+ private final double ax;
+ private final double ay;
+ private final double az;
+ private final double bx;
+ private final double by;
+ private final double bz;
+ private final double cx;
+ private final double cy;
+ private final double cz;
+ // CHECKSTYLE: resume JavadocVariableCheck
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng the source of randomness
+ */
+ CoordinateTriangleSampler3D(double[] a, double[] b, double[] c, \
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];
+ }
+
+ @Override
+ protected double[] createSample(double p1msmt, double s, double t) {
+ return new double[] {p1msmt * ax + s * bx + t * cx,
+ p1msmt * ay + s * by + t * cy,
+ p1msmt * az + s * bz + t * cz};
+ }
+ }
+ }
+
+ /**
+ * The ND triangle sampler.
+ */
+ @State(Scope.Benchmark)
+ public static class SamplerND extends SamplerData {
+ /** The number of dimensions. */
+ @Param({"2", "3", "4", "8", "16", "32"})
+ private int dimension;
+ /** The sampler type. */
+ @Param({BASELINE, BASELINE_IF, VECTORS, COORDINATES})
+ private String type;
+
+ /** {@inheritDoc} */
+ @Override
+ protected int getDimension() {
+ return dimension;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected Sampler createSampler(final double[] a, final double[] b, final \
double[] c, + final UniformRandomProvider rng) \
{ + if (BASELINE.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ double s = rng.nextDouble();
+ double t = rng.nextDouble();
+ final double[] x = new double[a.length];
+ for (int i = 0; i < x.length; i++) {
+ x[i] = s;
+ s = t;
+ t = x[i];
+ }
+ return x;
+ }
+ };
+ } else if (BASELINE_IF.equals(type)) {
+ return new Sampler() {
+ @Override
+ public double[] sample() {
+ double s = rng.nextDouble();
+ double t = rng.nextDouble();
+ final double[] x = new double[a.length];
+ if (s + t > 1) {
+ for (int i = 0; i < x.length; i++) {
+ x[i] = t;
+ t = s;
+ s = x[i];
+ }
+ return x;
+ }
+ for (int i = 0; i < x.length; i++) {
+ x[i] = s;
+ s = t;
+ t = x[i];
+ }
+ return x;
+ }
+ };
+ } else if (VECTORS.equals(type)) {
+ return new VectorTriangleSamplerND(a, b, c, rng);
+ } else if (COORDINATES.equals(type)) {
+ return new CoordinateTriangleSamplerND(a, b, c, rng);
+ }
+ throw new IllegalStateException(UNKNOWN_SAMPLER + type);
+ }
+
+ /**
+ * Sample vectors in ND.
+ */
+ private static class VectorTriangleSamplerND extends VectorTriangleSampler {
+ // CHECKSTYLE: stop JavadocVariableCheck
+ private final double[] a;
+ private final double[] v;
+ private final double[] w;
+ // CHECKSTYLE: resume JavadocVariableCheck
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng the source of randomness
+ */
+ VectorTriangleSamplerND(double[] a, double[] b, double[] c, \
UniformRandomProvider rng) { + super(rng);
+ this.a = a.clone();
+ v = new double[a.length];
+ w = new double[a.length];
+ for (int i = 0; i < a.length; i++) {
+ v[i] = b[i] - a[i];
+ w[i] = c[i] - a[i];
+ }
+ }
+
+ @Override
+ protected double[] createSample(double s, double t) {
+ final double[] x = new double[a.length];
+ for (int i = 0; i < x.length; i++) {
+ x[i] = a[i] + s * v[i] + t * w[i];
+ }
+ return x;
+ }
+ }
+
+ /**
+ * Sample using coordinates in ND.
+ */
+ private static class CoordinateTriangleSamplerND extends \
CoordinateTriangleSampler { + // CHECKSTYLE: stop JavadocVariableCheck
+ private final double[] a;
+ private final double[] b;
+ private final double[] c;
+ // CHECKSTYLE: resume JavadocVariableCheck
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng the source of randomness
+ */
+ CoordinateTriangleSamplerND(double[] a, double[] b, double[] c, \
UniformRandomProvider rng) { + super(rng);
+ this.a = a.clone();
+ this.b = b.clone();
+ this.c = c.clone();
+ }
+
+ @Override
+ protected double[] createSample(double p1msmt, double s, double t) {
+ final double[] x = new double[a.length];
+ for (int i = 0; i < x.length; i++) {
+ x[i] = p1msmt * a[i] + s * b[i] + t * c[i];
+ }
+ return x;
+ }
+ }
+ }
+
+ /**
+ * Run the sampler for the configured number of samples.
+ *
+ * @param bh Data sink
+ * @param data Input data.
+ */
+ private static void runSampler(Blackhole bh, SamplerData data) {
+ final Sampler sampler = data.getSampler();
+ for (int i = data.getSize() - 1; i >= 0; i--) {
+ bh.consume(sampler.sample());
+ }
+ }
+
+ /**
+ * Generation of uniform samples from a 2D triangle.
+ *
+ * @param bh Data sink
+ * @param data Input data.
+ */
+ @Benchmark
+ public void create2D(Blackhole bh, Sampler2D data) {
+ runSampler(bh, data);
+ }
+
+ /**
+ * Generation of uniform samples from a 3D triangle.
+ *
+ * @param bh Data sink
+ * @param data Input data.
+ */
+ @Benchmark
+ public void create3D(Blackhole bh, Sampler3D data) {
+ runSampler(bh, data);
+ }
+
+ /**
+ * Generation of uniform samples from an ND triangle.
+ *
+ * @param bh Data sink
+ * @param data Input data.
+ */
+ @Benchmark
+ public void createND(Blackhole bh, SamplerND data) {
+ runSampler(bh, data);
+ }
+}
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
new file mode 100644
index 0000000..d729e66
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/Coordinates.java
@@ -0,0 +1,65 @@
+/*
+ * 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;
+
+/**
+ * Utility class for common coordinate operations for shape samplers.
+ *
+ * @since 1.4
+ */
+final class Coordinates {
+
+ /** No public construction. */
+ private Coordinates() {}
+
+ /**
+ * Check that the values are finite. This method is primarily for parameter
+ * validation in methods and constructors, for example:
+ *
+ * <pre>
+ * public Line(double[] start, double[] end) {
+ * this.start = Coordinates.requireFinite(start, "start");
+ * this.end = Coordinates.requireFinite(end, "end");
+ * }
+ * </pre>
+ *
+ * @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
+ */
+ static double[] requireFinite(double[] values, String message) {
+ for (final double value : values) {
+ if (!isFinite(value)) {
+ throw new IllegalArgumentException(message + " contains non-finite \
value: " + value); + }
+ }
+ return values;
+ }
+
+ /**
+ * Checks if the value is finite.
+ * To be replaced by {@code Double.isFinite(double)} when source requires Java \
8. + *
+ * @param value the value
+ * @return true if finite
+ */
+ private static boolean isFinite(double value) {
+ return Math.abs(value) <= Double.MAX_VALUE;
+ }
+}
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TriangleSampler.java \
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TriangleSampler.java
new file mode 100644
index 0000000..77c25f1
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/TriangleSampler.java
@@ -0,0 +1,354 @@
+/*
+ * 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 <a \
href="https://mathworld.wolfram.com/TrianglePointPicking.html"> + * inside a \
triangle</a>. + *
+ * <ul>
+ * <li>
+ * Uses the algorithm described in:
+ * <blockquote>
+ * Turk, G. <i>Generating random points in triangles</i>. Glassner, A. S. (ed) \
(1990).<br> + * <strong>Graphic Gems</strong> Academic Press, pp. 24-28.
+ * </blockquote>
+ * </li>
+ * </ul>
+ *
+ * <p>Sampling uses:</p>
+ *
+ * <ul>
+ * <li>{@link UniformRandomProvider#nextDouble()}
+ * </ul>
+ *
+ * @since 1.4
+ */
+public abstract class TriangleSampler implements SharedStateSampler<TriangleSampler> \
{ + /** The dimension for 2D sampling. */
+ private static final int TWO_D = 2;
+ /** The dimension for 3D sampling. */
+ private static final int THREE_D = 3;
+ /** The source of randomness. */
+ private final UniformRandomProvider rng;
+
+ // The following code defines a plane as the set of all points r:
+ // r = r0 + sv + tw
+ // where s and t range over all real numbers, v and w are given linearly \
independent + // vectors defining the plane, and r0 is an arbitrary (but fixed) \
point in the plane. + //
+ // Sampling from a triangle (a,b,c) is performed when:
+ // s and t are in [0, 1] and s+t <= 1;
+ // r0 is one triangle vertex (a);
+ // and v (b-a) and w (c-a) are vectors from the other two vertices to r0.
+ //
+ // For robustness with large value coordinates the point r is computed without
+ // the requirement to compute v and w which can overflow:
+ //
+ // a + s(b-a) + t(c-a) == a + sb - sa + tc - ta
+ // == a(1 - s - t) + sb + tc
+ //
+ // Assuming the uniform deviates are from the 2^53 dyadic rationals in [0, 1) if \
s+t <= 1 + // then 1 - (s+t) is exact. Sampling is then done using:
+ //
+ // if (s + t <= 1):
+ // p = a(1 - (s + t)) + sb + tc
+ // else:
+ // p = a(1 - (1 - s) - (1 - t)) + (1 - s)b + (1 - t)c
+ // p = a(s - 1 + t) + (1 - s)b + (1 - t)c
+ //
+ // Note do not simplify (1 - (1 - s) - (1 - t)) to (s + t - 1) as s+t > 1 and \
has potential + // loss of a single bit of randomness due to rounding. An exact \
sum is s - 1 + t. +
+ /**
+ * Sample uniformly from a triangle in 2D. This is an non-array based \
specialisation of + * {@link TriangleSamplerND} for performance.
+ */
+ private static class TriangleSampler2D extends TriangleSampler {
+ /** The x component of vertex a. */
+ private final double ax;
+ /** The y component of vertex a. */
+ private final double ay;
+ /** The x component of vertex b. */
+ private final double bx;
+ /** The y component of vertex b. */
+ private final double by;
+ /** The x component of vertex c. */
+ private final double cx;
+ /** The y component of vertex c. */
+ private final double cy;
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng Source of randomness.
+ */
+ TriangleSampler2D(double[] a, double[] b, double[] c, UniformRandomProvider \
rng) { + super(rng);
+ ax = a[0];
+ ay = a[1];
+ bx = b[0];
+ by = b[1];
+ cx = c[0];
+ cy = c[1];
+ }
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers
+ * @param source Source to copy.
+ */
+ TriangleSampler2D(UniformRandomProvider rng, TriangleSampler2D source) {
+ super(rng);
+ ax = source.ax;
+ ay = source.ay;
+ bx = source.bx;
+ by = source.by;
+ cx = source.cx;
+ cy = source.cy;
+ }
+
+ @Override
+ public double[] createSample(double p1msmt, double s, double t) {
+ return new double[] {p1msmt * ax + s * bx + t * cx,
+ p1msmt * ay + s * by + t * cy};
+ }
+
+ @Override
+ public TriangleSampler withUniformRandomProvider(UniformRandomProvider rng) \
{ + return new TriangleSampler2D(rng, this);
+ }
+ }
+
+ /**
+ * Sample uniformly from a triangle in 3D. This is an non-array based \
specialisation of + * {@link TriangleSamplerND} for performance.
+ */
+ private static class TriangleSampler3D extends TriangleSampler {
+ /** The x component of vertex a. */
+ private final double ax;
+ /** The y component of vertex a. */
+ private final double ay;
+ /** The z component of vertex a. */
+ private final double az;
+ /** The x component of vertex b. */
+ private final double bx;
+ /** The y component of vertex b. */
+ private final double by;
+ /** The z component of vertex b. */
+ private final double bz;
+ /** The x component of vertex c. */
+ private final double cx;
+ /** The y component of vertex c. */
+ private final double cy;
+ /** The z component of vertex c. */
+ private final double cz;
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng Source of randomness.
+ */
+ TriangleSampler3D(double[] a, double[] b, double[] c, 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];
+ }
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers
+ * @param source Source to copy.
+ */
+ TriangleSampler3D(UniformRandomProvider rng, TriangleSampler3D source) {
+ super(rng);
+ ax = source.ax;
+ ay = source.ay;
+ az = source.az;
+ bx = source.bx;
+ by = source.by;
+ bz = source.bz;
+ cx = source.cx;
+ cy = source.cy;
+ cz = source.cz;
+ }
+
+ @Override
+ public double[] createSample(double p1msmt, double s, double t) {
+ return new double[] {p1msmt * ax + s * bx + t * cx,
+ p1msmt * ay + s * by + t * cy,
+ p1msmt * az + s * bz + t * cz};
+ }
+
+ @Override
+ public TriangleSampler withUniformRandomProvider(UniformRandomProvider rng) \
{ + return new TriangleSampler3D(rng, this);
+ }
+ }
+
+ /**
+ * Sample uniformly from a triangle in ND.
+ */
+ private static class TriangleSamplerND extends TriangleSampler {
+ /** The first vertex. */
+ private final double[] a;
+ /** The second vertex. */
+ private final double[] b;
+ /** The third vertex. */
+ private final double[] c;
+
+ /**
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng Source of randomness.
+ */
+ TriangleSamplerND(double[] a, double[] b, double[] c, UniformRandomProvider \
rng) { + super(rng);
+ // Defensive copy
+ this.a = a.clone();
+ this.b = b.clone();
+ this.c = c.clone();
+ }
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers
+ * @param source Source to copy.
+ */
+ TriangleSamplerND(UniformRandomProvider rng, TriangleSamplerND source) {
+ super(rng);
+ // Shared state is immutable
+ a = source.a;
+ b = source.b;
+ c = source.c;
+ }
+
+ @Override
+ public double[] createSample(double p1msmt, double s, double t) {
+ final double[] x = new double[a.length];
+ for (int i = 0; i < x.length; i++) {
+ x[i] = p1msmt * a[i] + s * b[i] + t * c[i];
+ }
+ return x;
+ }
+
+ @Override
+ public TriangleSampler withUniformRandomProvider(UniformRandomProvider rng) \
{ + return new TriangleSamplerND(rng, this);
+ }
+ }
+
+ /**
+ * @param rng Source of randomness.
+ */
+ TriangleSampler(UniformRandomProvider rng) {
+ this.rng = rng;
+ }
+
+ /**
+ * @return a random Cartesian point inside the triangle.
+ */
+ public double[] sample() {
+ final double s = rng.nextDouble();
+ final double t = rng.nextDouble();
+ final double spt = s + t;
+ if (spt > 1) {
+ // Transform: s1 = 1 - s; t1 = 1 - t.
+ // Compute: 1 - s1 - t1
+ // Do not assume (1 - (1-s) - (1-t)) is (s + t - 1), i.e. (spt - 1.0),
+ // to avoid loss of a random bit due to rounding when s + t > 1.
+ // An exact sum is (s - 1 + t).
+ return createSample(s - 1.0 + t, 1.0 - s, 1.0 - t);
+ }
+ // Here s + t is exact so can be subtracted to make 1 - s - t
+ return createSample(1.0 - spt, s, t);
+ }
+
+ /**
+ * Creates the sample given the random variates {@code s} and {@code t} in the
+ * interval {@code [0, 1]} and {@code s + t <= 1}. The sum {@code 1 - s - t} is \
provided. + * The sample can be obtained from the triangle abc using:
+ * <pre>
+ * p = a(1 - s - t) + sb + tc
+ * </pre>
+ *
+ * @param p1msmt plus 1 minus s minus t (1 - s - t)
+ * @param s the first variate s
+ * @param t the second variate t
+ * @return the sample
+ */
+ protected abstract double[] createSample(double p1msmt, double s, double t);
+
+ /**
+ * Create a triangle sampler with vertices {@code a}, {@code b} and {@code c}.
+ * Points are returned within the triangle.
+ *
+ * <p>Sampling is supported in dimensions of 2 or above. Samples will lie in the
+ * plane (2D Euclidean space) defined by using the three triangle vertices to
+ * create two vectors starting at a point in the plane and orientated in
+ * different directions along the plane.
+ *
+ * <p>No test for collinear points is performed. If the points are collinear the \
sampling + * distribution is undefined.
+ *
+ * @param a The first vertex.
+ * @param b The second vertex.
+ * @param c The third vertex.
+ * @param rng Source of randomness.
+ * @return the sampler
+ * @throws IllegalArgumentException If the vertices do not have the same
+ * dimension; the dimension is less than 2; or vertices have non-finite \
coordinates + */
+ public static TriangleSampler of(double[] a,
+ double[] b,
+ double[] c,
+ UniformRandomProvider rng) {
+ final int dimension = a.length;
+ if (dimension != b.length || dimension != c.length) {
+ throw new IllegalArgumentException(
+ new StringBuilder("Mismatch of vertex dimensions: \
").append(dimension).append(',') + \
.append(b.length).append(',') + \
.append(c.length).toString()); + }
+ // Detect non-finite vertices
+ Coordinates.requireFinite(a, "Vertex a");
+ Coordinates.requireFinite(b, "Vertex b");
+ Coordinates.requireFinite(c, "Vertex c");
+ // Low dimension specialisations
+ if (dimension == TWO_D) {
+ return new TriangleSampler2D(a, b, c, rng);
+ } else if (dimension == THREE_D) {
+ return new TriangleSampler3D(a, b, c, rng);
+ } else if (dimension > THREE_D) {
+ return new TriangleSamplerND(a, b, c, rng);
+ }
+ // Less than 2D
+ throw new IllegalArgumentException(
+ "Unsupported dimension: " + dimension);
+ }
+}
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
new file mode 100644
index 0000000..27c8ea0
--- /dev/null
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/CoordinatesTest.java
@@ -0,0 +1,49 @@
+/*
+ * 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.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test for {@link Coordinates} utility class.
+ */
+public class CoordinatesTest {
+ /**
+ * Test {@link Coordinates#requireFinite(double[], String)} detects infinite and \
NaN. + */
+ @Test
+ public void testRequireFiniteWithMessageThrows() {
+ final double[] c = {0, 1, 2};
+ final String message = "This should be prepended";
+ Assert.assertSame(c, Coordinates.requireFinite(c, message));
+ final double[] bad = {Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, \
Double.NaN}; + for (int i = 0; i < c.length; i++) {
+ for (final double d : bad) {
+ final double value = c[i];
+ c[i] = d;
+ try {
+ Coordinates.requireFinite(c, message);
+ Assert.fail(String.format("Did not detect non-finite coordinate: \
%d = %s", i, d)); + } catch (IllegalArgumentException ex) {
+ Assert.assertTrue("Missing message prefix", \
ex.getMessage().startsWith(message)); + }
+ c[i] = value;
+ }
+ }
+ }
+}
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
new file mode 100644
index 0000000..ea87a35
--- /dev/null
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/shape/TriangleSamplerTest.java
@@ -0,0 +1,786 @@
+/*
+ * 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.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;
+
+/**
+ * Test for {@link TriangleSampler}.
+ */
+public class TriangleSamplerTest {
+ // Precomputed 3D and 4D rotation matrices for forward and backward transforms
+ // created using the matlab function by jodag:
+ // https://stackoverflow.com/questions/50337642/how-to-calculate-a-rotation-matrix-in-n-dimensions-given-the-point-to-rotate-an
+
+ /** 3D rotation around the vector [1; 2; 3] by pi/4. */
+ private static final double[][] F3 = {
+ {0.728027725387508, -0.525104821111919, 0.440727305612110},
+ {0.608788597915763, 0.790790557990391, -0.063456571298848},
+ {-0.315201640406345, 0.314507901710379, 0.895395278995195}};
+ /** 3D rotation around the vector [1; 2; 3] by -pi/4. */
+ private static final double[][] R3 = {
+ {0.728027725387508, 0.608788597915762, -0.315201640406344},
+ {-0.525104821111919, 0.790790557990391, 0.314507901710379},
+ {0.440727305612110, -0.063456571298848, 0.895395278995195}};
+ /** 4D rotation around the orthogonal vectors [1 0; 0 1; 1 0; 0 1] by pi/4. */
+ private static final double[][] F4 = {
+ {0.853553390593274, -0.353553390593274, 0.146446609406726, \
0.353553390593274}, + {0.353553390593274, 0.853553390593274, \
-0.353553390593274, 0.146446609406726}, + {0.146446609406726, \
0.353553390593274, 0.853553390593274, -0.353553390593274}, + \
{-0.353553390593274, 0.146446609406726, 0.353553390593274, 0.853553390593274}}; + \
/** 4D rotation around the orthogonal vectors [1 0; 0 1; 1 0; 0 1] by -pi/4. */ + \
private static final double[][] R4 = { + {0.853553390593274, \
0.353553390593274, 0.146446609406726, -0.353553390593274}, + \
{-0.353553390593274, 0.853553390593274, 0.353553390593274, 0.146446609406726}, + \
{0.146446609406726, -0.353553390593274, 0.853553390593274, 0.353553390593274}, + \
{0.353553390593274, 0.146446609406726, -0.353553390593274, 0.853553390593274}}; +
+ /**
+ * Test the sampling assumptions used to transform coordinates outside the \
triangle + * back inside the triangle.
+ */
+ @Test
+ public void testSamplingAssumptions() {
+ // The separation between the 2^53 dyadic rationals in the interval [0, 1)
+ final double delta = 0x1.0p-53;
+ double s = 0.5;
+ double t = 0.5 + delta;
+ // This value cannot be exactly represented and is rounded
+ final double spt = s + t;
+ // Test that (1 - (1-s) - (1-t)) is not equal to (s + t - 1).
+ // This is due to the rounding to store s + t as a double.
+ final double expected = 1 - (1 - s) - (1 - t);
+ Assert.assertNotEquals(expected, spt - 1, 0.0);
+ Assert.assertNotEquals(expected, s + t - 1, 0.0);
+ // For any uniform deviate u in [0, 1], u - 1 is exact, thus s - 1 is exact
+ // and s - 1 + t is exact.
+ Assert.assertEquals(expected, s - 1 + t, 0.0);
+
+ // Test that a(1 - s - t) + sb + tc does not overflow is s+t = 1
+ final double max = Double.MAX_VALUE;
+ s -= delta;
+ final UniformRandomProvider rng = \
RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP); + for (int n = 0; n < \
100; n++) { + Assert.assertNotEquals(Double.POSITIVE_INFINITY, (1 - s - t) \
* max + s * max + t * max, 0.0); + s = rng.nextDouble();
+ t = 1.0 - s;
+ }
+ }
+
+ /**
+ * Test an unsupported dimension.
+ */
+ @Test(expected = IllegalArgumentException.class)
+ public void testInvalidDimensionThrows() {
+ TriangleSampler.of(new double[1], new double[1], new double[1], null);
+ }
+
+ /**
+ * Test a dimension mismatch between vertices.
+ */
+ @Test
+ public void testDimensionMismatchThrows() {
+ final double[] c2 = new double[2];
+ final double[] c3 = new double[3];
+ for (double[][] c : new double[][][] {
+ {c2, c2, c3},
+ {c2, c3, c2},
+ {c3, c2, c2},
+ {c2, c3, c3},
+ {c3, c3, c2},
+ {c3, c2, c3},
+ }) {
+ try {
+ TriangleSampler.of(c[0], c[1], c[2], null);
+ Assert.fail(String.format("Did not detect dimension mismatch: \
%d,%d,%d", + c[0].length, c[1].length, c[2].length));
+ } catch (IllegalArgumentException ex) {
+ // Expected
+ }
+ }
+ }
+
+ /**
+ * Test non-finite vertices.
+ */
+ @Test
+ public void testNonFiniteVertexCoordinates() {
+ // A valid triangle
+ final double[][] c = new double[][] {
+ {0, 0, 1}, {2, 1, 0}, {-1, 2, 3}
+ };
+ Assert.assertNotNull(TriangleSampler.of(c[0], c[1], c[2], null));
+ 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 {
+ TriangleSampler.of(c[0], c[1], c[2], null);
+ 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 triangle with coordinates that are separated by more than
+ * {@link Double#MAX_VALUE} in two dimensions.
+ */
+ @Test
+ public void testExtremeValueCoordinates2D() {
+ testExtremeValueCoordinates(2);
+ }
+
+ /**
+ * Test a triangle with coordinates that are separated by more than
+ * {@link Double#MAX_VALUE} in three dimensions.
+ */
+ @Test
+ public void testExtremeValueCoordinates3D() {
+ testExtremeValueCoordinates(3);
+ }
+
+ /**
+ * Test a triangle with coordinates that are separated by more than
+ * {@link Double#MAX_VALUE} in four dimensions.
+ */
+ @Test
+ public void testExtremeValueCoordinates4D() {
+ testExtremeValueCoordinates(4);
+ }
+
+ /**
+ * Test a triangle with coordinates that are separated by more than
+ * {@link Double#MAX_VALUE}.
+ *
+ * @param dimension the dimension
+ */
+ private static void testExtremeValueCoordinates(int dimension) {
+ // Object seed so use Long not long
+ final Long seed = 999666333L;
+ final double[][] c1 = new double[3][dimension];
+ final double[][] c2 = new double[3][dimension];
+ // Create a valid triangle that can be scaled
+ c1[0][0] = -1;
+ c1[0][1] = -1;
+ c1[1][0] = 1;
+ c1[1][1] = 1;
+ c1[2][0] = 1;
+ c1[2][1] = -1;
+ // Extremely large value for scaling. Use a power of 2 for exact scaling.
+ final double scale = 0x1.0p1023;
+ for (int i = 0; i < 3; i++) {
+ // Fill the remaining dimensions with 1 or -1
+ for (int j = 2; j < dimension; j++) {
+ c1[i][j] = 1 - (j & 0x2);
+ }
+ // Scale the second triangle
+ for (int j = 0; j < dimension; j++) {
+ c2[i][j] = c1[i][j] * scale;
+ }
+ }
+ // Show the triangle is too big to compute vectors between points.
+ Assert.assertEquals("Expect vector c - a to be infinite in the x dimension",
+ Double.POSITIVE_INFINITY, c2[2][0] - c2[0][0], 0.0);
+ Assert.assertEquals("Expect vector c - b to be infinite in the y dimension",
+ Double.NEGATIVE_INFINITY, c2[2][1] - c2[1][1], 0.0);
+
+ final TriangleSampler sampler1 = TriangleSampler.of(c1[0], c1[1], c1[2],
+ RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP, seed));
+ final TriangleSampler sampler2 = TriangleSampler.of(c2[0], c2[1], c2[2],
+ 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 two dimensions.
+ */
+ @Test
+ public void testDistribution2D() {
+ testDistributionND(2);
+ }
+
+ /**
+ * Test the distribution of points in three dimensions.
+ */
+ @Test
+ public void testDistribution3D() {
+ testDistributionND(3);
+ }
+
+ /**
+ * Test the distribution of points in four dimensions.
+ */
+ @Test
+ public void testDistribution4D() {
+ testDistributionND(4);
+ }
+
+ /**
+ * Test the distribution of points in N dimensions. The output coordinates
+ * should be uniform in the triangle.
+ *
+ * <p>ND is supported by transforming 2D triangles to ND, sampling in ND then
+ * transforming back to 2D. The transformed results should have a zero
+ * coordinate for indices above 1. Supports 2 to 4 dimensions.
+ *
+ * @param dimension the dimension in the range [2, 4]
+ */
+ private static void testDistributionND(int dimension) {
+ // Create 3 triangles that form a 3x4 rectangle:
+ // b-----c
+ // |\ /|
+ // | \ / |
+ // | e |
+ // | \ |
+ // | \|
+ // a-----d
+ //
+ // Sample from them proportional to the area:
+ // adb = 4 * 3 / 2 = 6
+ // bce = 3 * 2 / 2 = 3
+ // cde = 4 * 1.5 / 2 = 3
+ // The samples in the ratio 2:1:1 should be uniform within the rectangle.
+ final double[] a = {0, 0};
+ final double[] b = {0, 4};
+ final double[] c = {3, 4};
+ final double[] d = {3, 0};
+ final double[] e = {1.5, 2};
+
+ // Assign bins
+ final int bins = 20;
+ final int samplesPerBin = 10;
+ // Scale factors to assign x,y to a bin
+ final double sx = bins / 3.0;
+ final double sy = bins / 4.0;
+
+ // Expect a uniform distribution (this is rescaled by the ChiSquareTest)
+ final double[] expected = new double[bins * bins];
+ Arrays.fill(expected, 1);
+
+ // Support ND by applying a rotation transform to the 2D triangles to \
generate + // n-coordinates. The samples are transformed back and should be \
uniform in 2D. + Transform forward;
+ Transform reverse;
+ if (dimension == 4) {
+ forward = new ForwardTransform(F4);
+ reverse = new ReverseTransform(R4);
+ } else if (dimension == 3) {
+ forward = new ForwardTransform(F3);
+ reverse = new ReverseTransform(R3);
+ } else if (dimension == 2) {
+ forward = reverse = new Transform2Dto2D();
+ } else {
+ throw new AssertionError("Unsupported dimension: " + dimension);
+ }
+
+ // 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, 0xfabcab); + final \
TriangleSampler sampler1 = TriangleSampler.of(forward.apply(a), forward.apply(d), \
forward.apply(b), rng); + final TriangleSampler sampler2 = \
TriangleSampler.of(forward.apply(b), forward.apply(c), forward.apply(e), rng); + \
final TriangleSampler sampler3 = TriangleSampler.of(forward.apply(c), \
forward.apply(d), forward.apply(e), rng); + final InsideTriangle test1 = new \
InsideTriangle(a, d, b); + final InsideTriangle test2 = new InsideTriangle(b, \
c, e); + final InsideTriangle test3 = new InsideTriangle(c, d, e);
+ final int samples = expected.length * samplesPerBin;
+ for (int n = 0; n < 1; n++) {
+ // Assign each coordinate to a region inside the combined rectangle
+ final long[] observed = new long[expected.length];
+ // Sample according to the area of each triangle (ratio 2:1:1)
+ for (int i = 0; i < samples; i += 4) {
+ addObservation(reverse.apply(sampler1.sample()), observed, bins, sx, \
sy, test1); + addObservation(reverse.apply(sampler1.sample()), \
observed, bins, sx, sy, test1); + \
addObservation(reverse.apply(sampler2.sample()), observed, bins, sx, sy, test2); + \
addObservation(reverse.apply(sampler3.sample()), observed, bins, sx, sy, test3); + \
} + 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 * sx)
+ * </pre>
+ *
+ * @param v the sample (2D coordinate xy)
+ * @param observed the observations
+ * @param bins the numbers of bins in the x dimension
+ * @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 inside the inside triangle test
+ */
+ private static void addObservation(double[] v, long[] observed,
+ int bins, double sx, double sy, InsideTriangle inside) {
+ final double x = v[0];
+ final double y = v[1];
+ // Test the point is inside the triangle
+ Assert.assertTrue(inside.test(x, y));
+ // Add to the correct bin after using the offset
+ final int binx = (int) (x * sx);
+ final int biny = (int) (y * sy);
+ observed[biny * bins + binx]++;
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for 2D.
+ */
+ @Test
+ public void testSharedStateSampler2D() {
+ testSharedStateSampler(2);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for 3D.
+ */
+ @Test
+ public void testSharedStateSampler3D() {
+ testSharedStateSampler(3);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for 4D.
+ */
+ @Test
+ public void testSharedStateSampler4D() {
+ testSharedStateSampler(4);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation for the given dimension.
+ */
+ private static void testSharedStateSampler(int dimension) {
+ 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, dimension); + final double[] c2 = \
createCoordinate(2, dimension); + final double[] c3 = createCoordinate(-3, \
dimension); + final TriangleSampler sampler1 = TriangleSampler.of(c1, c2, c3, \
rng1); + final TriangleSampler 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 for 2D.
+ */
+ @Test
+ public void testChangedInputCoordinates2D() {
+ testChangedInputCoordinates(2);
+ }
+
+ /**
+ * Test the input vectors are copied and not used by reference for 3D.
+ */
+ @Test
+ public void testChangedInputCoordinates3D() {
+ testChangedInputCoordinates(3);
+ }
+
+ /**
+ * Test the input vectors are copied and not used by reference for 4D.
+ */
+ @Test
+ public void testChangedInputCoordinates4D() {
+ testChangedInputCoordinates(4);
+ }
+
+ /**
+ * Test the input vectors are copied and not used by reference for the given
+ * dimension.
+ *
+ * @param dimension the dimension
+ */
+ private static void testChangedInputCoordinates(int dimension) {
+ 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, dimension); + final double[] c2 = \
createCoordinate(2, dimension); + final double[] c3 = createCoordinate(-3, \
dimension); + final TriangleSampler sampler1 = TriangleSampler.of(c1, c2, c3, \
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 < dimension; i++) {
+ c1[i] += offset;
+ c2[i] += offset;
+ c3[i] += offset;
+ }
+ final TriangleSampler sampler2 = TriangleSampler.of(c1, c2, c3, rng2);
+ for (int n = 0; n < 3; n++) {
+ final double[] s1 = sampler1.sample();
+ final double[] s2 = sampler2.sample();
+ Assert.assertEquals(s1.length, s2.length);
+ Assert.assertFalse("First sampler has used the vertices by reference",
+ Arrays.equals(s1, s2));
+ for (int i = 0; i < dimension; i++) {
+ Assert.assertEquals(s1[i] + offset, s2[i], 1e-10);
+ }
+ }
+ }
+
+ /**
+ * Creates the coordinate of length specified by the dimension filled with
+ * the given value and the dimension index: x + i.
+ *
+ * @param x the value for index 0
+ * @param dimension the dimension
+ * @return the coordinate
+ */
+ private static double[] createCoordinate(double x, int dimension) {
+ final double[] coord = new double[dimension];
+ for (int i = 0; i < dimension; i++) {
+ coord[0] = x + i;
+ }
+ return coord;
+ }
+
+ /**
+ * Test the transform generates 3D coordinates and can reverse them.
+ */
+ @Test
+ public void testTransform3D() {
+ testTransformND(3);
+ }
+
+ /**
+ * Test the transform generates 4D coordinates and can reverse them.
+ */
+ @Test
+ public void testTransform4D() {
+ testTransformND(4);
+ }
+
+ /**
+ * Test the transform generates ND coordinates and can reverse them.
+ *
+ * @param dimension the dimension
+ */
+ private static void testTransformND(int dimension) {
+ Transform forward;
+ Transform reverse;
+ if (dimension == 4) {
+ forward = new ForwardTransform(F4);
+ reverse = new ReverseTransform(R4);
+ } else if (dimension == 3) {
+ forward = new ForwardTransform(F3);
+ reverse = new ReverseTransform(R3);
+ } else if (dimension == 2) {
+ forward = reverse = new Transform2Dto2D();
+ } else {
+ throw new AssertionError("Unsupported dimension: " + dimension);
+ }
+
+ final UniformRandomProvider rng = \
RandomSource.create(RandomSource.SPLIT_MIX_64, 789L); + double sum = 0;
+ for (int n = 0; n < 10; n++) {
+ final double[] a = new double[] {rng.nextDouble(), rng.nextDouble()};
+ final double[] b = forward.apply(a);
+ Assert.assertEquals(dimension, b.length);
+ for (int i = 2; i < dimension; i++) {
+ sum += Math.abs(b[i]);
+ }
+ final double[] c = reverse.apply(b);
+ Assert.assertArrayEquals(a, c, 1e-10);
+ }
+ // Check that higher dimension coordinates are generated
+ Assert.assertTrue(sum > 0.5);
+ }
+
+ /**
+ * Test 3D rotations forward and reverse.
+ */
+ @Test
+ public void testRotations3D() {
+ final double[] x = {1, 0.5, 0};
+ final double[] y = multiply(F3, x);
+ Assert.assertArrayEquals(new double[] {0.465475314831549, 1.004183876910958, \
-0.157947689551155}, y, 1e-10); + Assert.assertEquals(length(x), length(y), \
1e-10); + final double[] x2 = multiply(R3, y);
+ Assert.assertArrayEquals(x, x2, 1e-10);
+ }
+
+ /**
+ * Test 4D rotations forward and reverse.
+ */
+ @Test
+ public void testRotations4D() {
+ final double[] x = {1, 0.5, 0, 0};
+ final double[] y = multiply(F4, x);
+ Assert.assertArrayEquals(
+ new double[] {0.676776695296637, 0.780330085889911, \
0.323223304703363, -0.280330085889911}, y, 1e-10); + \
Assert.assertEquals(length(x), length(y), 1e-10); + final double[] x2 = \
multiply(R4, y); + Assert.assertArrayEquals(x, x2, 1e-10);
+ }
+
+ /**
+ * Matrix multiplication. It is assumed the matrix is square and matches (or \
exceeds) + * the vector length.
+ *
+ * @param m the matrix
+ * @param v the vector
+ * @return the result
+ */
+ private static double[] multiply(double[][] matrix, double[] v) {
+ final int n = matrix.length;
+ final int m = v.length;
+ final double[] r = new double[n];
+ for (int i = 0; i < n; i++) {
+ double sum = 0;
+ for (int j = 0; j < m; j++) {
+ sum += matrix[i][j] * v[j];
+ }
+ r[i] = sum;
+ }
+ return r;
+ }
+
+ /**
+ * @return the length (L2-norm) of given vector.
+ */
+ private static double length(double[] vector) {
+ double total = 0;
+ for (final double d : vector) {
+ total += d * d;
+ }
+ return Math.sqrt(total);
+ }
+
+ /**
+ * Test the inside triangle predicate.
+ */
+ @Test
+ public void testInsideTriangle() {
+ final InsideTriangle inside = new InsideTriangle(1, 2, 3, 1, 0.5, 6);
+ // Vertices
+ Assert.assertTrue(inside.test(1, 2));
+ Assert.assertTrue(inside.test(3, 1));
+ Assert.assertTrue(inside.test(0.5, 6));
+ // Edge
+ Assert.assertTrue(inside.test(0.75, 4));
+ // Inside
+ Assert.assertTrue(inside.test(1.5, 3));
+ // Outside
+ Assert.assertFalse(inside.test(0, 20));
+ Assert.assertFalse(inside.test(-20, 0));
+ Assert.assertFalse(inside.test(6, 6));
+ // Just outside
+ Assert.assertFalse(inside.test(0.75, 4 - 1e-10));
+
+ // Note:
+ // Touching triangles can both have the point inside.
+ // This predicate is not suitable for assigning points uniquely to
+ // non-overlapping triangles that share an edge.
+ final InsideTriangle inside2 = new InsideTriangle(1, 2, 3, 1, 0, -2);
+ Assert.assertTrue(inside.test(2, 1.5));
+ Assert.assertTrue(inside2.test(2, 1.5));
+ }
+
+ /**
+ * Define a transform on coordinates.
+ */
+ private interface Transform {
+ /**
+ * Apply the transform.
+ *
+ * @param coord the coordinates
+ * @return the new coordinates
+ */
+ double[] apply(double[] coord);
+ }
+
+ /**
+ * Transform coordinates from 2D to a higher dimension using the rotation \
matrix. + */
+ private static class ForwardTransform implements Transform {
+ private final double[][] r;
+
+ /**
+ * @param r the rotation matrix
+ */
+ ForwardTransform(double[][] r) {
+ this.r = r;
+ }
+
+ @Override
+ public double[] apply(double[] coord) {
+ return multiply(r, coord);
+ }
+ }
+
+ /**
+ * Transform coordinates from a higher dimension to 2D using the rotation \
matrix. + * The result should be in the 2D plane (i.e. higher dimensions of the \
transformed vector + * are asserted to be zero).
+ */
+ private static class ReverseTransform implements Transform {
+ private final double[][] r;
+ private final int n;
+
+ /**
+ * @param r the rotation matrix
+ */
+ ReverseTransform(double[][] r) {
+ this.r = r;
+ n = r.length;
+ }
+
+ @Override
+ public double[] apply(double[] coord) {
+ Assert.assertEquals(n, coord.length);
+ final double[] x = multiply(r, coord);
+ // This should reverse the 2D transform and return to the XY plane.
+ for (int i = 2; i < x.length; i++) {
+ Assert.assertEquals(0.0, x[i], 1e-14);
+ }
+ return new double[] {x[0], x[1]};
+ }
+ }
+
+ /**
+ * No-operation transform on 2D input. Asserts the input coordinates are length \
2. + */
+ private static class Transform2Dto2D implements Transform {
+ @Override
+ public double[] apply(double[] coord) {
+ Assert.assertEquals(2, coord.length);
+ return coord;
+ }
+ }
+
+ /**
+ * Class to test if a point is inside the triangle.
+ *
+ * <p>This function has been adapted from a StackOverflow post by Cédric \
Dufour. It converts the + * point to unscaled barycentric coordinates (s,t) and \
tests they are within the triangle. + * (Scaling would be done by dividing by \
twice the area of the triangle.) + *
+ * <h2>Warning</h2>
+ *
+ * <p>This assigns points on the edges as inside irrespective of edge \
orientation. Thus + * back-to-back touching triangles can both have the point \
inside them. A predicate for geometry + * applications where the point must be \
within a unique non-overlapping triangle should use + * a different solution, \
e.g. assigning new points to the result of a triangulation. + * For testing \
sampling within the triangle this predicate is acceptable. + *
+ * @see <a
+ * href="https://stackoverflow.com/questions/2049582/how-to-determine-if-a-point-is-in-a-2d-triangle?lq=1">
+ * Point in a triangle</a>
+ * @see <a href="https://stackoverflow.com/a/34093754">Point inside triangle by \
Cédric Dufour</a> + */
+ private static class InsideTriangle {
+ private final double p2x;
+ private final double p2y;
+ private final double dX21;
+ private final double dY12;
+ private final double dY20;
+ private final double dX02;
+ private final double d;
+
+ /**
+ * Create an instance.
+ *
+ * @param p0 triangle vertex 0
+ * @param p1 triangle vertex 1
+ * @param p2 triangle vertex 2
+ */
+ InsideTriangle(double[] p0, double[] p1, double[] p2) {
+ this(p0[0], p0[1], p1[0], p1[1], p2[0], p2[1]);
+ }
+ /**
+ * Create an instance.
+ *
+ * @param p0x triangle vertex 0 x coordinate
+ * @param p0y triangle vertex 0 y coordinate
+ * @param p1x triangle vertex 1 x coordinate
+ * @param p1y triangle vertex 1 y coordinate
+ * @param p2x triangle vertex 2 x coordinate
+ * @param p2y triangle vertex 2 y coordinate
+ */
+ InsideTriangle(double p0x, double p0y, double p1x, double p1y, double p2x, \
double p2y) { + this.p2x = p2x;
+ this.p2y = p2y;
+ // Precompute factors
+ dX21 = p2x - p1x;
+ dY12 = p1y - p2y;
+ dY20 = p2y - p0y;
+ dX02 = p0x - p2x;
+ // d = twice the signed area of the triangle
+ d = dY12 * (p0x - p2x) + dX21 * (p0y - p2y);
+ }
+
+ /**
+ * Check the point is inside the triangle.
+ *
+ * @param px the point x coordinate
+ * @param py the point y coordinate
+ * @return true if inside the triangle
+ */
+ boolean test(double px, double py) {
+ // Barycentric coordinates:
+ // p = p0 + (p1 - p0) * s + (p2 - p0) * t
+ // The point p is inside the triangle if 0 <= s <= 1 and 0 <= t <= 1 and \
s + t <= 1 + //
+ // The following solves s and t.
+ // Some factors are precomputed.
+ final double dX = px - p2x;
+ final double dY = py - p2y;
+ final double s = dY12 * dX + dX21 * dY;
+ final double t = dY20 * dX + dX02 * dY;
+ if (d < 0) {
+ return s <= 0 && t <= 0 && s + t >= d;
+ }
+ return s >= 0 && t >= 0 && s + t <= d;
+ }
+ }
+}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 1ab4aa9..9026445 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="131">
+ New "TriangleSampler" to sample uniformly from a triangle.
+ </action>
<action dev="aherbert" type="add" issue="132">
New "o.a.c.rng.sampling.shape" package for sampling coordinates from shapes.
</action>
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic