| /* |
| * 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; |
| } |
| } |
| } |
