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.distribution.ContinuousSampler;
23 import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
24 import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
25 import org.apache.commons.rng.simple.RandomSource;
26 import org.openjdk.jmh.annotations.Benchmark;
27 import org.openjdk.jmh.annotations.BenchmarkMode;
28 import org.openjdk.jmh.annotations.Fork;
29 import org.openjdk.jmh.annotations.Measurement;
30 import org.openjdk.jmh.annotations.Mode;
31 import org.openjdk.jmh.annotations.OutputTimeUnit;
32 import org.openjdk.jmh.annotations.Param;
33 import org.openjdk.jmh.annotations.Scope;
34 import org.openjdk.jmh.annotations.Setup;
35 import org.openjdk.jmh.annotations.State;
36 import org.openjdk.jmh.annotations.Warmup;
37 import org.openjdk.jmh.infra.Blackhole;
38 import java.util.concurrent.TimeUnit;
39
40 /**
41 * Executes benchmark to compare the speed of generating samples within an N-dimension unit ball.
42 */
43 @BenchmarkMode(Mode.AverageTime)
44 @OutputTimeUnit(TimeUnit.NANOSECONDS)
45 @Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
46 @Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
47 @State(Scope.Benchmark)
48 @Fork(value = 1, jvmArgs = { "-server", "-Xms128M", "-Xmx128M" })
49 public class UnitBallSamplerBenchmark {
50 /** Name for the baseline method. */
51 private static final String BASELINE = "Baseline";
52 /** Name for the rejection method. */
53 private static final String REJECTION = "Rejection";
54 /** Name for the disk-point picking method. */
55 private static final String DISK_POINT = "DiskPoint";
56 /** Name for the ball-point picking method. */
57 private static final String BALL_POINT = "BallPoint";
58 /** Name for the method moving from the surface of the hypersphere to an internal point. */
59 private static final String HYPERSPHERE_INTERNAL = "HypersphereInternal";
60 /** Name for the picking from the surface of a greater dimension hypersphere and discarding 2 points. */
61 private static final String HYPERSPHERE_DISCARD = "HypersphereDiscard";
62 /** Error message for an unknown sampler type. */
63 private static final String UNKNOWN_SAMPLER = "Unknown sampler type: ";
64 /** The mask to extract the lower 53-bits from a long. */
65 private static final long LOWER_53_BITS = -1L >>> 11;
66
67 /**
68 * Base class for a sampler using a provided source of randomness.
69 */
70 private abstract static class BaseSampler implements ObjectSampler<double[]> {
71 /** The source of randomness. */
72 protected UniformRandomProvider rng;
73
74 /**
75 * Create an instance.
76 *
77 * @param rng the source of randomness
78 */
79 BaseSampler(UniformRandomProvider rng) {
80 this.rng = rng;
81 }
82 }
83
84 /**
85 * Base class for the sampler data.
86 * Contains the source of randomness and the number of samples.
87 * The sampler should be created by a sub-class of the data.
88 */
89 @State(Scope.Benchmark)
90 public abstract static class SamplerData {
91 /** The sampler. */
92 private ObjectSampler<double[]> sampler;
93
94 /** The number of samples. */
95 @Param({//"1",
96 "100"})
97 private int size;
98
99 /**
100 * Gets the size.
101 *
102 * @return the size
103 */
104 public int getSize() {
105 return size;
106 }
107
108 /**
109 * Gets the sampler.
110 *
111 * @return the sampler
112 */
113 public ObjectSampler<double[]> getSampler() {
114 return sampler;
115 }
116
117 /**
118 * Create the source of randomness.
119 */
120 @Setup
121 public void setup() {
122 // This could be configured using @Param
123 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create();
124 sampler = createSampler(rng);
125 }
126
127 /**
128 * Creates the sampler.
129 *
130 * @param rng the source of randomness
131 * @return the sampler
132 */
133 protected abstract ObjectSampler<double[]> createSampler(UniformRandomProvider rng);
134 }
135
136 /**
137 * The 1D unit line sampler.
138 */
139 @State(Scope.Benchmark)
140 public static class Sampler1D extends SamplerData {
141 /** Name for the signed double method. */
142 private static final String SIGNED_DOUBLE = "signedDouble";
143 /** Name for the signed double method, version 2. */
144 private static final String SIGNED_DOUBLE2 = "signedDouble2";
145 /** Name for the two doubles method. */
146 private static final String TWO_DOUBLES = "twoDoubles";
147 /** Name for the boolean double method. */
148 private static final String BOOLEAN_DOUBLE = "booleanDouble";
149
150 /** The sampler type. */
151 @Param({BASELINE, SIGNED_DOUBLE, SIGNED_DOUBLE2, TWO_DOUBLES, BOOLEAN_DOUBLE})
152 private String type;
153
154 /** {@inheritDoc} */
155 @Override
156 protected ObjectSampler<double[]> createSampler(final UniformRandomProvider rng) {
157 if (BASELINE.equals(type)) {
158 return () -> new double[] {0.5};
159 } else if (SIGNED_DOUBLE.equals(type)) {
160 // Sample [-1, 1) uniformly
161 return () -> new double[] {makeSignedDouble(rng.nextLong())};
162 } else if (SIGNED_DOUBLE2.equals(type)) {
163 // Sample [-1, 1) uniformly
164 return () -> new double[] {makeSignedDouble2(rng.nextLong())};
165 } else if (TWO_DOUBLES.equals(type)) {
166 // Sample [-1, 1) excluding -0.0 but also missing the final 1.0 - 2^-53.
167 // The 1.0 could be adjusted to 1.0 - 2^-53 to create the interval (-1, 1).
168 return () -> new double[] {rng.nextDouble() + rng.nextDouble() - 1.0};
169 } else if (BOOLEAN_DOUBLE.equals(type)) {
170 // This will sample (-1, 1) including -0.0 and 0.0
171 return () -> new double[] {rng.nextBoolean() ? -rng.nextDouble() : rng.nextDouble()};
172 }
173 throw new IllegalStateException(UNKNOWN_SAMPLER + type);
174 }
175 }
176
177 /**
178 * The 2D unit disk sampler.
179 */
180 @State(Scope.Benchmark)
181 public static class Sampler2D extends SamplerData {
182 /** The sampler type. */
183 @Param({BASELINE, REJECTION, DISK_POINT, HYPERSPHERE_INTERNAL, HYPERSPHERE_DISCARD})
184 private String type;
185
186 /** {@inheritDoc} */
187 @Override
188 protected ObjectSampler<double[]> createSampler(final UniformRandomProvider rng) {
189 if (BASELINE.equals(type)) {
190 return () -> new double[] {0.5, 0};
191 } else if (REJECTION.equals(type)) {
192 return new RejectionSampler(rng);
193 } else if (DISK_POINT.equals(type)) {
194 return new DiskPointSampler(rng);
195 } else if (HYPERSPHERE_INTERNAL.equals(type)) {
196 return new HypersphereInternalSampler(rng);
197 } else if (HYPERSPHERE_DISCARD.equals(type)) {
198 return new HypersphereDiscardSampler(rng);
199 }
200 throw new IllegalStateException(UNKNOWN_SAMPLER + type);
201 }
202
203 /**
204 * Sample using a simple rejection method.
205 */
206 private static class RejectionSampler extends BaseSampler {
207 /**
208 * @param rng the source of randomness
209 */
210 RejectionSampler(UniformRandomProvider rng) {
211 super(rng);
212 }
213
214 @Override
215 public double[] sample() {
216 // Generate via rejection method of a circle inside a square of edge length 2.
217 // This should compute approximately 2^2 / pi = 1.27 square positions per sample.
218 double x;
219 double y;
220 do {
221 x = makeSignedDouble(rng.nextLong());
222 y = makeSignedDouble(rng.nextLong());
223 } while (x * x + y * y > 1.0);
224 return new double[] {x, y};
225 }
226 }
227
228 /**
229 * Sample using disk point picking.
230 * @see <a href="https://mathworld.wolfram.com/DiskPointPicking.html">Disk point picking</a>
231 */
232 private static class DiskPointSampler extends BaseSampler {
233 /** 2 pi. */
234 private static final double TWO_PI = 2 * Math.PI;
235
236 /**
237 * @param rng the source of randomness
238 */
239 DiskPointSampler(UniformRandomProvider rng) {
240 super(rng);
241 }
242
243 @Override
244 public double[] sample() {
245 final double t = TWO_PI * rng.nextDouble();
246 final double r = Math.sqrt(rng.nextDouble());
247 final double x = r * Math.cos(t);
248 final double y = r * Math.sin(t);
249 return new double[] {x, y};
250 }
251 }
252
253 /**
254 * Choose a uniform point X on the unit hypersphere, then multiply it by U<sup>1/n</sup>
255 * where U in [0, 1].
256 * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
257 * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
258 */
259 private static class HypersphereInternalSampler extends BaseSampler {
260 /** The normal distribution. */
261 private final NormalizedGaussianSampler normal;
262
263 /**
264 * @param rng the source of randomness
265 */
266 HypersphereInternalSampler(UniformRandomProvider rng) {
267 super(rng);
268 normal = ZigguratSampler.NormalizedGaussian.of(rng);
269 }
270
271 @Override
272 public double[] sample() {
273 final double x = normal.sample();
274 final double y = normal.sample();
275 final double sum = x * x + y * y;
276 // Note: Handle the possibility of a zero sum and invalid inverse
277 if (sum == 0) {
278 return sample();
279 }
280 // Take a point on the unit hypersphere then multiply it by U^1/n
281 final double f = Math.sqrt(rng.nextDouble()) / Math.sqrt(sum);
282 return new double[] {x * f, y * f};
283 }
284 }
285
286 /**
287 * Take a random point on the (n+1)-dimensional hypersphere and drop two coordinates.
288 * Remember that the (n+1)-hypersphere is the unit sphere of R^(n+2).
289 * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
290 * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
291 */
292 private static class HypersphereDiscardSampler extends BaseSampler {
293 /** The normal distribution. */
294 private final NormalizedGaussianSampler normal;
295
296 /**
297 * @param rng the source of randomness
298 */
299 HypersphereDiscardSampler(UniformRandomProvider rng) {
300 super(rng);
301 normal = ZigguratSampler.NormalizedGaussian.of(rng);
302 }
303
304 @Override
305 public double[] sample() {
306 // Discard 2 samples from the coordinate but include in the sum
307 final double x0 = normal.sample();
308 final double x1 = normal.sample();
309 final double x = normal.sample();
310 final double y = normal.sample();
311 final double sum = x0 * x0 + x1 * x1 + x * x + y * y;
312 // Note: Handle the possibility of a zero sum and invalid inverse
313 if (sum == 0) {
314 return sample();
315 }
316 final double f = 1.0 / Math.sqrt(sum);
317 return new double[] {x * f, y * f};
318 }
319 }
320 }
321
322 /**
323 * The 3D unit ball sampler.
324 */
325 @State(Scope.Benchmark)
326 public static class Sampler3D extends SamplerData {
327 /** The sampler type. */
328 @Param({BASELINE, REJECTION, BALL_POINT, HYPERSPHERE_INTERNAL, HYPERSPHERE_DISCARD})
329 private String type;
330
331 /** {@inheritDoc} */
332 @Override
333 protected ObjectSampler<double[]> createSampler(final UniformRandomProvider rng) {
334 if (BASELINE.equals(type)) {
335 return () -> new double[] {0.5, 0, 0};
336 } else if (REJECTION.equals(type)) {
337 return new RejectionSampler(rng);
338 } else if (BALL_POINT.equals(type)) {
339 return new BallPointSampler(rng);
340 } else if (HYPERSPHERE_INTERNAL.equals(type)) {
341 return new HypersphereInternalSampler(rng);
342 } else if (HYPERSPHERE_DISCARD.equals(type)) {
343 return new HypersphereDiscardSampler(rng);
344 }
345 throw new IllegalStateException(UNKNOWN_SAMPLER + type);
346 }
347
348 /**
349 * Sample using a simple rejection method.
350 */
351 private static class RejectionSampler extends BaseSampler {
352 /**
353 * @param rng the source of randomness
354 */
355 RejectionSampler(UniformRandomProvider rng) {
356 super(rng);
357 }
358
359 @Override
360 public double[] sample() {
361 // Generate via rejection method of a ball inside a cube of edge length 2.
362 // This should compute approximately 2^3 / (4pi/3) = 1.91 cube positions per sample.
363 double x;
364 double y;
365 double z;
366 do {
367 x = makeSignedDouble(rng.nextLong());
368 y = makeSignedDouble(rng.nextLong());
369 z = makeSignedDouble(rng.nextLong());
370 } while (x * x + y * y + z * z > 1.0);
371 return new double[] {x, y, z};
372 }
373 }
374
375 /**
376 * Sample using ball point picking.
377 * @see <a href="https://mathworld.wolfram.com/BallPointPicking.html">Ball point picking</a>
378 */
379 private static class BallPointSampler extends BaseSampler {
380 /** The normal distribution. */
381 private final NormalizedGaussianSampler normal;
382 /** The exponential distribution. */
383 private final ContinuousSampler exp;
384
385 /**
386 * @param rng the source of randomness
387 */
388 BallPointSampler(UniformRandomProvider rng) {
389 super(rng);
390 normal = ZigguratSampler.NormalizedGaussian.of(rng);
391 // Exponential(mean=2) == Chi-squared distribution(degrees freedom=2)
392 // thus is the equivalent of the HypersphereDiscardSampler.
393 // Here we use mean = 1 and scale the output later.
394 exp = ZigguratSampler.Exponential.of(rng);
395 }
396
397 @Override
398 public double[] sample() {
399 final double x = normal.sample();
400 final double y = normal.sample();
401 final double z = normal.sample();
402 // Include the exponential sample. It has mean 1 so multiply by 2.
403 final double sum = exp.sample() * 2 + x * x + y * y + z * z;
404 // Note: Handle the possibility of a zero sum and invalid inverse
405 if (sum == 0) {
406 return sample();
407 }
408 final double f = 1.0 / Math.sqrt(sum);
409 return new double[] {x * f, y * f, z * f};
410 }
411 }
412
413 /**
414 * Choose a uniform point X on the unit hypersphere, then multiply it by U<sup>1/n</sup>
415 * where U in [0, 1].
416 * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
417 * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
418 */
419 private static class HypersphereInternalSampler extends BaseSampler {
420 /** The normal distribution. */
421 private final NormalizedGaussianSampler normal;
422
423 /**
424 * @param rng the source of randomness
425 */
426 HypersphereInternalSampler(UniformRandomProvider rng) {
427 super(rng);
428 normal = ZigguratSampler.NormalizedGaussian.of(rng);
429 }
430
431 @Override
432 public double[] sample() {
433 final double x = normal.sample();
434 final double y = normal.sample();
435 final double z = normal.sample();
436 final double sum = x * x + y * y + z * z;
437 // Note: Handle the possibility of a zero sum and invalid inverse
438 if (sum == 0) {
439 return sample();
440 }
441 // Take a point on the unit hypersphere then multiply it by U^1/n
442 final double f = Math.cbrt(rng.nextDouble()) / Math.sqrt(sum);
443 return new double[] {x * f, y * f, z * f};
444 }
445 }
446
447 /**
448 * Take a random point on the (n+1)-dimensional hypersphere and drop two coordinates.
449 * Remember that the (n+1)-hypersphere is the unit sphere of R^(n+2).
450 * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
451 * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
452 */
453 private static class HypersphereDiscardSampler extends BaseSampler {
454 /** The normal distribution. */
455 private final NormalizedGaussianSampler normal;
456
457 /**
458 * @param rng the source of randomness
459 */
460 HypersphereDiscardSampler(UniformRandomProvider rng) {
461 super(rng);
462 normal = ZigguratSampler.NormalizedGaussian.of(rng);
463 }
464
465 @Override
466 public double[] sample() {
467 // Discard 2 samples from the coordinate but include in the sum
468 final double x0 = normal.sample();
469 final double x1 = normal.sample();
470 final double x = normal.sample();
471 final double y = normal.sample();
472 final double z = normal.sample();
473 final double sum = x0 * x0 + x1 * x1 + x * x + y * y + z * z;
474 // Note: Handle the possibility of a zero sum and invalid inverse
475 if (sum == 0) {
476 return sample();
477 }
478 final double f = 1.0 / Math.sqrt(sum);
479 return new double[] {x * f, y * f, z * f};
480 }
481 }
482 }
483
484 /**
485 * The ND unit ball sampler.
486 */
487 @State(Scope.Benchmark)
488 public static class SamplerND extends SamplerData {
489 /** The sampler type. */
490 @Param({BASELINE, REJECTION, BALL_POINT, HYPERSPHERE_INTERNAL, HYPERSPHERE_DISCARD})
491 private String type;
492 /** The number of dimensions. */
493 @Param({"3", "4", "5"})
494 private int dimension;
495
496 /** {@inheritDoc} */
497 @Override
498 protected ObjectSampler<double[]> createSampler(final UniformRandomProvider rng) {
499 if (BASELINE.equals(type)) {
500 return new ObjectSampler<double[]>() {
501 private final int dim = dimension;
502 @Override
503 public double[] sample() {
504 final double[] sample = new double[dim];
505 for (int i = 0; i < dim; i++) {
506 sample[i] = 0.01;
507 }
508 return sample;
509 }
510 };
511 } else if (REJECTION.equals(type)) {
512 return new RejectionSampler(rng, dimension);
513 } else if (BALL_POINT.equals(type)) {
514 return new BallPointSampler(rng, dimension);
515 } else if (HYPERSPHERE_INTERNAL.equals(type)) {
516 return new HypersphereInternalSampler(rng, dimension);
517 } else if (HYPERSPHERE_DISCARD.equals(type)) {
518 return new HypersphereDiscardSampler(rng, dimension);
519 }
520 throw new IllegalStateException(UNKNOWN_SAMPLER + type);
521 }
522
523 /**
524 * Sample using a simple rejection method.
525 */
526 private static class RejectionSampler extends BaseSampler {
527 /** The dimension. */
528 private final int dimension;
529
530 /**
531 * @param rng the source of randomness
532 * @param dimension the dimension
533 */
534 RejectionSampler(UniformRandomProvider rng, int dimension) {
535 super(rng);
536 this.dimension = dimension;
537 }
538
539 @Override
540 public double[] sample() {
541 // Generate via rejection method of a ball inside a hypercube of edge length 2.
542 final double[] sample = new double[dimension];
543 double sum;
544 do {
545 sum = 0;
546 for (int i = 0; i < dimension; i++) {
547 final double x = makeSignedDouble(rng.nextLong());
548 sum += x * x;
549 sample[i] = x;
550 }
551 } while (sum > 1.0);
552 return sample;
553 }
554 }
555
556 /**
557 * Sample using ball point picking.
558 * @see <a href="https://mathworld.wolfram.com/BallPointPicking.html">Ball point picking</a>
559 */
560 private static class BallPointSampler extends BaseSampler {
561 /** The dimension. */
562 private final int dimension;
563 /** The normal distribution. */
564 private final NormalizedGaussianSampler normal;
565 /** The exponential distribution. */
566 private final ContinuousSampler exp;
567
568 /**
569 * @param rng the source of randomness
570 * @param dimension the dimension
571 */
572 BallPointSampler(UniformRandomProvider rng, int dimension) {
573 super(rng);
574 this.dimension = dimension;
575 normal = ZigguratSampler.NormalizedGaussian.of(rng);
576 // Exponential(mean=2) == Chi-squared distribution(degrees freedom=2)
577 // thus is the equivalent of the HypersphereDiscardSampler.
578 // Here we use mean = 1 and scale the output later.
579 exp = ZigguratSampler.Exponential.of(rng);
580 }
581
582 @Override
583 public double[] sample() {
584 final double[] sample = new double[dimension];
585 // Include the exponential sample. It has mean 1 so multiply by 2.
586 double sum = exp.sample() * 2;
587 for (int i = 0; i < dimension; i++) {
588 final double x = normal.sample();
589 sum += x * x;
590 sample[i] = x;
591 }
592 // Note: Handle the possibility of a zero sum and invalid inverse
593 if (sum == 0) {
594 return sample();
595 }
596 final double f = 1.0 / Math.sqrt(sum);
597 for (int i = 0; i < dimension; i++) {
598 sample[i] *= f;
599 }
600 return sample;
601 }
602 }
603
604 /**
605 * Choose a uniform point X on the unit hypersphere, then multiply it by U<sup>1/n</sup>
606 * where U in [0, 1].
607 * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
608 * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
609 */
610 private static class HypersphereInternalSampler extends BaseSampler {
611 /** The dimension. */
612 private final int dimension;
613 /** The normal distribution. */
614 private final NormalizedGaussianSampler normal;
615 /** Reciprocal of the dimension. */
616 private final double power;
617
618 /**
619 * @param rng the source of randomness
620 * @param dimension the dimension
621 */
622 HypersphereInternalSampler(UniformRandomProvider rng, int dimension) {
623 super(rng);
624 this.dimension = dimension;
625 power = 1.0 / dimension;
626 normal = ZigguratSampler.NormalizedGaussian.of(rng);
627 }
628
629 @Override
630 public double[] sample() {
631 final double[] sample = new double[dimension];
632 double sum = 0;
633 for (int i = 0; i < dimension; i++) {
634 final double x = normal.sample();
635 sum += x * x;
636 sample[i] = x;
637 }
638 // Note: Handle the possibility of a zero sum and invalid inverse
639 if (sum == 0) {
640 return sample();
641 }
642 // Take a point on the unit hypersphere then multiply it by U^1/n
643 final double f = Math.pow(rng.nextDouble(), power) / Math.sqrt(sum);
644 for (int i = 0; i < dimension; i++) {
645 sample[i] *= f;
646 }
647 return sample;
648 }
649 }
650
651 /**
652 * Take a random point on the (n+1)-dimensional hypersphere and drop two coordinates.
653 * Remember that the (n+1)-hypersphere is the unit sphere of R^(n+2).
654 * @see <a href="https://mathoverflow.net/questions/309567/sampling-a-uniformly-distributed-point-inside-a-hypersphere">
655 * Sampling a uniformly distributed point INSIDE a hypersphere?</a>
656 */
657 private static class HypersphereDiscardSampler extends BaseSampler {
658 /** The dimension. */
659 private final int dimension;
660 /** The normal distribution. */
661 private final NormalizedGaussianSampler normal;
662
663 /**
664 * @param rng the source of randomness
665 * @param dimension the dimension
666 */
667 HypersphereDiscardSampler(UniformRandomProvider rng, int dimension) {
668 super(rng);
669 this.dimension = dimension;
670 normal = ZigguratSampler.NormalizedGaussian.of(rng);
671 }
672
673 @Override
674 public double[] sample() {
675 final double[] sample = new double[dimension];
676 // Discard 2 samples from the coordinate but include in the sum
677 final double x0 = normal.sample();
678 final double x1 = normal.sample();
679 double sum = x0 * x0 + x1 * x1;
680 for (int i = 0; i < dimension; i++) {
681 final double x = normal.sample();
682 sum += x * x;
683 sample[i] = x;
684 }
685 // Note: Handle the possibility of a zero sum and invalid inverse
686 if (sum == 0) {
687 return sample();
688 }
689 final double f = 1.0 / Math.sqrt(sum);
690 for (int i = 0; i < dimension; i++) {
691 sample[i] *= f;
692 }
693 return sample;
694 }
695 }
696 }
697
698 /**
699 * Creates a signed double in the range {@code [-1, 1)}. The magnitude is sampled evenly
700 * from the 2<sup>54</sup> dyadic rationals in the range.
701 *
702 * <p>Note: This method will not return samples for both -0.0 and 0.0.
703 *
704 * @param bits the bits
705 * @return the double
706 */
707 private static double makeSignedDouble(long bits) {
708 // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a signed
709 // shift of 10 in place of an unsigned shift of 11.
710 return (bits >> 10) * 0x1.0p-53d;
711 }
712
713 /**
714 * Creates a signed double in the range {@code [-1, 1)}. The magnitude is sampled evenly
715 * from the 2<sup>54</sup> dyadic rationals in the range.
716 *
717 * <p>Note: This method will not return samples for both -0.0 and 0.0.
718 *
719 * @param bits the bits
720 * @return the double
721 */
722 private static double makeSignedDouble2(long bits) {
723 // Use the upper 54 bits on the assumption they are more random.
724 // The sign bit generates a value of 0 or 1 for subtraction.
725 // The next 53 bits generates a positive number in the range [0, 1).
726 // [0, 1) - (0 or 1) => [-1, 1)
727 return (((bits >>> 10) & LOWER_53_BITS) * 0x1.0p-53d) - (bits >>> 63);
728 }
729
730 /**
731 * Run the sampler for the configured number of samples.
732 *
733 * @param bh Data sink
734 * @param data Input data.
735 */
736 private static void runSampler(Blackhole bh, SamplerData data) {
737 final ObjectSampler<double[]> sampler = data.getSampler();
738 for (int i = data.getSize() - 1; i >= 0; i--) {
739 bh.consume(sampler.sample());
740 }
741 }
742
743 /**
744 * Generation of uniform samples on a 1D unit line.
745 *
746 * @param bh Data sink
747 * @param data Input data.
748 */
749 @Benchmark
750 public void create1D(Blackhole bh, Sampler1D data) {
751 runSampler(bh, data);
752 }
753
754 /**
755 * Generation of uniform samples from a 2D unit disk.
756 *
757 * @param bh Data sink
758 * @param data Input data.
759 */
760 @Benchmark
761 public void create2D(Blackhole bh, Sampler2D data) {
762 runSampler(bh, data);
763 }
764
765 /**
766 * Generation of uniform samples from a 3D unit ball.
767 *
768 * @param bh Data sink
769 * @param data Input data.
770 */
771 @Benchmark
772 public void create3D(Blackhole bh, Sampler3D data) {
773 runSampler(bh, data);
774 }
775
776 /**
777 * Generation of uniform samples from an ND unit ball.
778 *
779 * @param bh Data sink
780 * @param data Input data.
781 */
782 @Benchmark
783 public void createND(Blackhole bh, SamplerND data) {
784 runSampler(bh, data);
785 }
786 }