1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.rng.examples.jmh.sampling.shape;
19
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.sampling.ObjectSampler;
22 import org.apache.commons.rng.sampling.UnitSphereSampler;
23 import org.apache.commons.rng.simple.RandomSource;
24 import org.openjdk.jmh.annotations.Benchmark;
25 import org.openjdk.jmh.annotations.BenchmarkMode;
26 import org.openjdk.jmh.annotations.Fork;
27 import org.openjdk.jmh.annotations.Level;
28 import org.openjdk.jmh.annotations.Measurement;
29 import org.openjdk.jmh.annotations.Mode;
30 import org.openjdk.jmh.annotations.OutputTimeUnit;
31 import org.openjdk.jmh.annotations.Param;
32 import org.openjdk.jmh.annotations.Scope;
33 import org.openjdk.jmh.annotations.Setup;
34 import org.openjdk.jmh.annotations.State;
35 import org.openjdk.jmh.annotations.Warmup;
36 import org.openjdk.jmh.infra.Blackhole;
37 import java.util.concurrent.TimeUnit;
38
39 /**
40 * Executes benchmark to compare the speed of generating samples within a tetrahedron.
41 */
42 @BenchmarkMode(Mode.AverageTime)
43 @OutputTimeUnit(TimeUnit.NANOSECONDS)
44 @Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
45 @Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
46 @State(Scope.Benchmark)
47 @Fork(value = 1, jvmArgs = { "-server", "-Xms512M", "-Xmx512M" })
48 public class TetrahedronSamplerBenchmark {
49 /** Name for the baseline method. */
50 private static final String BASELINE = "Baseline";
51 /** Name for the method using array coordinates. */
52 private static final String ARRAY = "Array";
53 /** Name for the method using non-array (primitive) coordinates. */
54 private static final String NON_ARRAY = "NonArray";
55 /** Name for the method using array coordinates and inline sample algorithm. */
56 private static final String ARRAY_INLINE = "ArrayInline";
57 /** Name for the method using non-array (primitive) coordinates and inline sample algorithm. */
58 private static final String NON_ARRAY_INLINE = "NonArrayInline";
59 /** Error message for an unknown sampler type. */
60 private static final String UNKNOWN_SAMPLER = "Unknown sampler type: ";
61
62 /**
63 * The base class for sampling from a tetrahedron.
64 *
65 * <ul>
66 * <li>
67 * Uses the algorithm described in:
68 * <blockquote>
69 * Rocchini, C. and Cignoni, P. (2001)<br>
70 * <i>Generating Random Points in a Tetrahedron</i>.<br>
71 * <strong>Journal of Graphics Tools</strong> 5(4), pp. 9-12.
72 * </blockquote>
73 * </li>
74 * </ul>
75 *
76 * @see <a href="https://doi.org/10.1080/10867651.2000.10487528">
77 * Rocchini, C. & Cignoni, P. (2001) Journal of Graphics Tools 5, pp. 9-12</a>
78 */
79 private abstract static class TetrahedronSampler implements ObjectSampler<double[]> {
80 /** The source of randomness. */
81 private final UniformRandomProvider rng;
82
83 /**
84 * @param rng Source of randomness.
85 */
86 TetrahedronSampler(UniformRandomProvider rng) {
87 this.rng = rng;
88 }
89
90 @Override
91 public double[] sample() {
92 double s = rng.nextDouble();
93 double t = rng.nextDouble();
94 final double u = rng.nextDouble();
95 // Care is taken to ensure the 3 deviates remain in the 2^53 dyadic rationals in [0, 1).
96 // The following are exact for all the 2^53 dyadic rationals:
97 // 1 - u; u in [0, 1]
98 // u - 1; u in [0, 1]
99 // u + 1; u in [-1, 0]
100 // u + v; u in [-1, 0], v in [0, 1]
101 // u + v; u, v in [0, 1], u + v <= 1
102
103 // Cut and fold with the plane s + t = 1
104 if (s + t > 1) {
105 // (s, t, u) = (1 - s, 1 - t, u) if s + t > 1
106 s = 1 - s;
107 t = 1 - t;
108 }
109 // Now s + t <= 1.
110 // Cut and fold with the planes t + u = 1 and s + t + u = 1.
111 final double tpu = t + u;
112 final double sptpu = s + tpu;
113 if (sptpu > 1) {
114 if (tpu > 1) {
115 // (s, t, u) = (s, 1 - u, 1 - s - t) if t + u > 1
116 // 1 - s - (1-u) - (1-s-t) == u - 1 + t
117 return createSample(u - 1 + t, s, 1 - u, 1 - s - t);
118 }
119 // (s, t, u) = (1 - t - u, t, s + t + u - 1) if t + u <= 1
120 // 1 - (1-t-u) - t - (s+t+u-1) == 1 - s - t
121 return createSample(1 - s - t, 1 - tpu, t, s - 1 + tpu);
122 }
123 return createSample(1 - sptpu, s, t, u);
124 }
125
126 /**
127 * Creates the sample given the random variates {@code s}, {@code t} and {@code u} in the
128 * interval {@code [0, 1]} and {@code s + t + u <= 1}. The sum {@code 1 - s - t - u} is
129 * provided. The sample can be obtained from the tetrahedron {@code abcd} using:
130 *
131 * <pre>
132 * p = (1 - s - t - u)a + sb + tc + ud
133 * </pre>
134 *
135 * @param p1msmtmu plus 1 minus s minus t minus u (1 - s - t - u)
136 * @param s the first variate s
137 * @param t the second variate t
138 * @param u the third variate u
139 * @return the sample
140 */
141 protected abstract double[] createSample(double p1msmtmu, double s, double t, double u);
142 }
143
144 /**
145 * Sample from a tetrahedron using array coordinates.
146 */
147 private static class ArrayTetrahedronSampler extends TetrahedronSampler {
148 // CHECKSTYLE: stop JavadocVariableCheck
149 private final double[] a;
150 private final double[] b;
151 private final double[] c;
152 private final double[] d;
153 // CHECKSTYLE: resume JavadocVariableCheck
154
155 /**
156 * @param rng the source of randomness
157 * @param a The first vertex.
158 * @param b The second vertex.
159 * @param c The third vertex.
160 * @param d The fourth vertex.
161 */
162 ArrayTetrahedronSampler(UniformRandomProvider rng,
163 double[] a, double[] b, double[] c, double[] d) {
164 super(rng);
165 this.a = a.clone();
166 this.b = b.clone();
167 this.c = c.clone();
168 this.d = d.clone();
169 }
170
171 @Override
172 protected double[] createSample(double p1msmtmu, double s, double t, double u) {
173 return new double[] {p1msmtmu * a[0] + s * b[0] + t * c[0] + u * d[0],
174 p1msmtmu * a[1] + s * b[1] + t * c[1] + u * d[1],
175 p1msmtmu * a[2] + s * b[2] + t * c[2] + u * d[2]};
176 }
177 }
178
179 /**
180 * Sample from a tetrahedron using non-array coordinates.
181 */
182 private static class NonArrayTetrahedronSampler extends TetrahedronSampler {
183 // CHECKSTYLE: stop JavadocVariableCheck
184 private final double ax;
185 private final double bx;
186 private final double cx;
187 private final double dx;
188 private final double ay;
189 private final double by;
190 private final double cy;
191 private final double dy;
192 private final double az;
193 private final double bz;
194 private final double cz;
195 private final double dz;
196 // CHECKSTYLE: resume JavadocVariableCheck
197
198 /**
199 * @param rng the source of randomness
200 * @param a The first vertex.
201 * @param b The second vertex.
202 * @param c The third vertex.
203 * @param d The fourth vertex.
204 */
205 NonArrayTetrahedronSampler(UniformRandomProvider rng,
206 double[] a, double[] b, double[] c, double[] d) {
207 super(rng);
208 ax = a[0];
209 ay = a[1];
210 az = a[2];
211 bx = b[0];
212 by = b[1];
213 bz = b[2];
214 cx = c[0];
215 cy = c[1];
216 cz = c[2];
217 dx = d[0];
218 dy = d[1];
219 dz = d[2];
220 }
221
222 @Override
223 protected double[] createSample(double p1msmtmu, double s, double t, double u) {
224 return new double[] {p1msmtmu * ax + s * bx + t * cx + u * dx,
225 p1msmtmu * ay + s * by + t * cy + u * dy,
226 p1msmtmu * az + s * bz + t * cz + u * dz};
227 }
228 }
229
230 /**
231 * Sample from a tetrahedron using array coordinates with an inline sample algorithm
232 * in-place of a method call.
233 */
234 private static class ArrayInlineTetrahedronSampler implements ObjectSampler<double[]> {
235 // CHECKSTYLE: stop JavadocVariableCheck
236 private final double[] a;
237 private final double[] b;
238 private final double[] c;
239 private final double[] d;
240 private final UniformRandomProvider rng;
241 // CHECKSTYLE: resume JavadocVariableCheck
242
243 /**
244 * @param rng the source of randomness
245 * @param a The first vertex.
246 * @param b The second vertex.
247 * @param c The third vertex.
248 * @param d The fourth vertex.
249 */
250 ArrayInlineTetrahedronSampler(UniformRandomProvider rng,
251 double[] a, double[] b, double[] c, double[] d) {
252 this.a = a.clone();
253 this.b = b.clone();
254 this.c = c.clone();
255 this.d = d.clone();
256 this.rng = rng;
257 }
258
259 @Override
260 public double[] sample() {
261 double s = rng.nextDouble();
262 double t = rng.nextDouble();
263 double u = rng.nextDouble();
264
265 if (s + t > 1) {
266 // (s, t, u) = (1 - s, 1 - t, u) if s + t > 1
267 s = 1 - s;
268 t = 1 - t;
269 }
270
271 double p1msmtmu;
272 final double tpu = t + u;
273 final double sptpu = s + tpu;
274 if (sptpu > 1) {
275 if (tpu > 1) {
276 // (s, t, u) = (s, 1 - u, 1 - s - t) if t + u > 1
277 final double tt = t;
278 p1msmtmu = u - 1 + t;
279 t = 1 - u;
280 u = 1 - s - tt;
281 } else {
282 // (s, t, u) = (1 - t - u, t, s + t + u - 1) if t + u <= 1
283 final double ss = s;
284 p1msmtmu = 1 - s - t;
285 s = 1 - tpu;
286 u = ss - 1 + tpu;
287 }
288 } else {
289 p1msmtmu = 1 - sptpu;
290 }
291
292 return new double[] {p1msmtmu * a[0] + s * b[0] + t * c[0] + u * d[0],
293 p1msmtmu * a[1] + s * b[1] + t * c[1] + u * d[1],
294 p1msmtmu * a[2] + s * b[2] + t * c[2] + u * d[2]};
295 }
296 }
297
298 /**
299 * Sample from a tetrahedron using non-array coordinates with an inline sample algorithm
300 * in-place of a method call.
301 */
302 private static class NonArrayInlineTetrahedronSampler implements ObjectSampler<double[]> {
303 // CHECKSTYLE: stop JavadocVariableCheck
304 private final double ax;
305 private final double bx;
306 private final double cx;
307 private final double dx;
308 private final double ay;
309 private final double by;
310 private final double cy;
311 private final double dy;
312 private final double az;
313 private final double bz;
314 private final double cz;
315 private final double dz;
316 private final UniformRandomProvider rng;
317 // CHECKSTYLE: resume JavadocVariableCheck
318
319 /**
320 * @param rng the source of randomness
321 * @param a The first vertex.
322 * @param b The second vertex.
323 * @param c The third vertex.
324 * @param d The fourth vertex.
325 */
326 NonArrayInlineTetrahedronSampler(UniformRandomProvider rng,
327 double[] a, double[] b, double[] c, double[] d) {
328 ax = a[0];
329 ay = a[1];
330 az = a[2];
331 bx = b[0];
332 by = b[1];
333 bz = b[2];
334 cx = c[0];
335 cy = c[1];
336 cz = c[2];
337 dx = d[0];
338 dy = d[1];
339 dz = d[2];
340 this.rng = rng;
341 }
342
343 @Override
344 public double[] sample() {
345 double s = rng.nextDouble();
346 double t = rng.nextDouble();
347 double u = rng.nextDouble();
348
349 if (s + t > 1) {
350 // (s, t, u) = (1 - s, 1 - t, u) if s + t > 1
351 s = 1 - s;
352 t = 1 - t;
353 }
354
355 double p1msmtmu;
356 final double tpu = t + u;
357 final double sptpu = s + tpu;
358 if (sptpu > 1) {
359 if (tpu > 1) {
360 // (s, t, u) = (s, 1 - u, 1 - s - t) if t + u > 1
361 final double tt = t;
362 p1msmtmu = u - 1 + t;
363 t = 1 - u;
364 u = 1 - s - tt;
365 } else {
366 // (s, t, u) = (1 - t - u, t, s + t + u - 1) if t + u <= 1
367 final double ss = s;
368 p1msmtmu = 1 - s - t;
369 s = 1 - tpu;
370 u = ss - 1 + tpu;
371 }
372 } else {
373 p1msmtmu = 1 - sptpu;
374 }
375
376 return new double[] {p1msmtmu * ax + s * bx + t * cx + u * dx,
377 p1msmtmu * ay + s * by + t * cy + u * dy,
378 p1msmtmu * az + s * bz + t * cz + u * dz};
379 }
380 }
381
382 /**
383 * Contains the sampler and the number of samples.
384 */
385 @State(Scope.Benchmark)
386 public static class SamplerData {
387 /** The sampler. */
388 private ObjectSampler<double[]> sampler;
389
390 /** The number of samples. */
391 @Param({"1", "10", "100", "1000", "10000"})
392 private int size;
393
394 /** The sampler type. */
395 @Param({BASELINE, ARRAY, NON_ARRAY, ARRAY_INLINE, NON_ARRAY_INLINE})
396 private String type;
397
398 /**
399 * Gets the size.
400 *
401 * @return the size
402 */
403 public int getSize() {
404 return size;
405 }
406
407 /**
408 * Gets the sampler.
409 *
410 * @return the sampler
411 */
412 public ObjectSampler<double[]> getSampler() {
413 return sampler;
414 }
415
416 /**
417 * Create the source of randomness and the sampler.
418 */
419 @Setup(Level.Iteration)
420 public void setup() {
421 // This could be configured using @Param
422 final UniformRandomProvider rng = RandomSource.XO_SHI_RO_256_PP.create();
423 final UnitSphereSampler s = UnitSphereSampler.of(rng, 3);
424 final double[] a = s.sample();
425 final double[] b = s.sample();
426 final double[] c = s.sample();
427 final double[] d = s.sample();
428 sampler = createSampler(rng, a, b, c, d);
429 }
430
431 /**
432 * Creates the tetrahedron sampler.
433 *
434 * @param a The first vertex.
435 * @param b The second vertex.
436 * @param c The third vertex.
437 * @param d The fourth vertex.
438 * @param rng the source of randomness
439 * @return the sampler
440 */
441 private ObjectSampler<double[]> createSampler(final UniformRandomProvider rng,
442 double[] a, double[] b, double[] c, double[] d) {
443 if (BASELINE.equals(type)) {
444 return () -> {
445 final double s = rng.nextDouble();
446 final double t = rng.nextDouble();
447 final double u = rng.nextDouble();
448 return new double[] {s, t, u};
449 };
450 } else if (ARRAY.equals(type)) {
451 return new ArrayTetrahedronSampler(rng, a, b, c, d);
452 } else if (NON_ARRAY.equals(type)) {
453 return new NonArrayTetrahedronSampler(rng, a, b, c, d);
454 } else if (ARRAY_INLINE.equals(type)) {
455 return new ArrayInlineTetrahedronSampler(rng, a, b, c, d);
456 } else if (NON_ARRAY_INLINE.equals(type)) {
457 return new NonArrayInlineTetrahedronSampler(rng, a, b, c, d);
458 }
459 throw new IllegalStateException(UNKNOWN_SAMPLER + type);
460 }
461 }
462
463 /**
464 * Run the sampler for the configured number of samples.
465 *
466 * @param bh Data sink
467 * @param data Input data.
468 */
469 @Benchmark
470 public void sample(Blackhole bh, SamplerData data) {
471 final ObjectSampler<double[]> sampler = data.getSampler();
472 for (int i = data.getSize() - 1; i >= 0; i--) {
473 bh.consume(sampler.sample());
474 }
475 }
476 }