View Javadoc
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.numbers.examples.jmh.complex;
19  
20  import org.apache.commons.numbers.complex.Complex;
21  import org.apache.commons.rng.UniformRandomProvider;
22  import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
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.Measurement;
28  import org.openjdk.jmh.annotations.Mode;
29  import org.openjdk.jmh.annotations.OutputTimeUnit;
30  import org.openjdk.jmh.annotations.Param;
31  import org.openjdk.jmh.annotations.Scope;
32  import org.openjdk.jmh.annotations.Setup;
33  import org.openjdk.jmh.annotations.State;
34  import org.openjdk.jmh.annotations.Warmup;
35  import org.openjdk.jmh.infra.Blackhole;
36  import java.util.Arrays;
37  import java.util.concurrent.TimeUnit;
38  import java.util.function.BiFunction;
39  import java.util.function.Predicate;
40  import java.util.function.Supplier;
41  import java.util.function.ToDoubleFunction;
42  import java.util.function.UnaryOperator;
43  import java.util.stream.Stream;
44  
45  /**
46   * Executes a benchmark to measure the speed of operations in the {@link Complex} class.
47   */
48  @BenchmarkMode(Mode.AverageTime)
49  @OutputTimeUnit(TimeUnit.NANOSECONDS)
50  @Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
51  @Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
52  @State(Scope.Benchmark)
53  @Fork(value = 1, jvmArgs = {"-server", "-Xms512M", "-Xmx512M"})
54  public class ComplexPerformance {
55      /**
56       * An array of edge numbers that will produce edge case results from functions:
57       * {@code +/-inf, +/-max, +/-min, +/-0, nan}.
58       */
59      private static final double[] EDGE_NUMBERS = {
60          Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.MAX_VALUE,
61          -Double.MAX_VALUE, Double.MIN_VALUE, -Double.MIN_VALUE, 0.0, -0.0, Double.NaN};
62  
63      /** The range to use for uniform random numbers. */
64      private static final double RANGE = 3.456789;
65  
66      /**
67       * Contains the size of numbers.
68       */
69      @State(Scope.Benchmark)
70      public static class ComplexNumberSize {
71          /**
72           * The size of the data.
73           */
74          @Param({"10000"})
75          private int size;
76  
77          /**
78           * Gets the size.
79           *
80           * @return the size
81           */
82          public int getSize() {
83              return size;
84          }
85      }
86  
87      /**
88       * Contains an array of complex numbers.
89       */
90      @State(Scope.Benchmark)
91      public static class ComplexNumbers extends ComplexNumberSize {
92          /** The numbers. */
93          protected Complex[] numbers;
94  
95          /**
96           * The type of the data.
97           */
98          @Param({"cis", "vector", "log-uniform", "uniform", "edge"})
99          private String type;
100 
101         /**
102          * Gets the numbers.
103          *
104          * @return the numbers
105          */
106         public Complex[] getNumbers() {
107             return numbers;
108         }
109 
110         /**
111          * Create the complex numbers.
112          */
113         @Setup
114         public void setup() {
115             numbers = createNumbers(RandomSource.XO_RO_SHI_RO_128_PP.create());
116         }
117 
118         /**
119          * Creates the numbers.
120          *
121          * @param rng Random number generator.
122          * @return the random complex number
123          */
124         Complex[] createNumbers(UniformRandomProvider rng) {
125             Supplier<Complex> generator;
126             if ("cis".equals(type)) {
127                 generator = () -> Complex.ofCis(rng.nextDouble() * 2 * Math.PI);
128             } else if ("vector".equals(type)) {
129                 // An unnormalised random vector is created using a Gaussian sample
130                 // for each dimension. Normalisation would create a cis number.
131                 // This is effectively a polar complex number with random modulus
132                 // in [-pi, pi] and random magnitude in a range defined by a Chi-squared
133                 // distribution with 2 degrees of freedom.
134                 final ZigguratNormalizedGaussianSampler s = ZigguratNormalizedGaussianSampler.of(rng);
135                 generator = () -> Complex.ofCartesian(s.sample(), s.sample());
136             } else if ("log-uniform".equals(type)) {
137                 generator = () -> Complex.ofCartesian(createLogUniformNumber(rng), createLogUniformNumber(rng));
138             } else if ("uniform".equals(type)) {
139                 generator = () -> Complex.ofCartesian(createUniformNumber(rng), createUniformNumber(rng));
140             } else if ("edge".equals(type)) {
141                 generator = () -> Complex.ofCartesian(createEdgeNumber(rng), createEdgeNumber(rng));
142             } else {
143                 throw new IllegalStateException("Unknown number type: " + type);
144             }
145             return Stream.generate(generator).limit(getSize()).toArray(Complex[]::new);
146         }
147     }
148 
149     /**
150      * Contains two arrays of complex numbers.
151      */
152     @State(Scope.Benchmark)
153     public static class TwoComplexNumbers extends ComplexNumbers {
154         /** The numbers. */
155         private Complex[] numbers2;
156 
157         /**
158          * Gets the second set of numbers.
159          *
160          * @return the numbers
161          */
162         public Complex[] getNumbers2() {
163             return numbers2;
164         }
165 
166         /**
167          * Create the complex numbers.
168          */
169         @Override
170         @Setup
171         public void setup() {
172             // Do not call super.setup() so we recycle the RNG and avoid duplicates
173             final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create();
174             numbers = createNumbers(rng);
175             numbers2 = createNumbers(rng);
176         }
177     }
178 
179     /**
180      * Contains an array of complex numbers and an array of real numbers.
181      */
182     @State(Scope.Benchmark)
183     public static class ComplexAndRealNumbers extends ComplexNumbers {
184         /** The numbers. */
185         private double[] numbers2;
186 
187         /**
188          * Gets the second set of numbers.
189          *
190          * @return the numbers
191          */
192         public double[] getNumbers2() {
193             return numbers2;
194         }
195 
196         /**
197          * Create the complex numbers.
198          */
199         @Override
200         @Setup
201         public void setup() {
202             // Do not call super.setup() so we recycle the RNG and avoid duplicates
203             final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create();
204             numbers = createNumbers(rng);
205             numbers2 = Arrays.stream(createNumbers(rng)).mapToDouble(Complex::real).toArray();
206         }
207     }
208 
209     /**
210      * Define a function between a complex and real number.
211      */
212     private interface ComplexRealFunction {
213         /**
214          * Applies this function to the given arguments.
215          *
216          * @param z the complex argument
217          * @param x the real argument
218          * @return the function result
219          */
220         Complex apply(Complex z, double x);
221     }
222 
223     /**
224      * Creates a random double number with a random sign and mantissa and a large range for
225      * the exponent. The numbers will not be uniform over the range. This samples randomly
226      * using the components of a double. The limiting distribution is the log-uniform distribution.
227      *
228      * @param rng Random number generator.
229      * @return the random number
230      * @see <a href="https://en.wikipedia.org/wiki/Reciprocal_distribution">Reciprocal (log-uniform) distribution</a>
231      */
232     private static double createLogUniformNumber(UniformRandomProvider rng) {
233         // Create random doubles using random bits in the sign bit and the mantissa.
234         // Then create an exponent in the range -64 to 64. Thus the sum product
235         // of 4 max or min values will not over or underflow.
236         final long mask = ((1L << 52) - 1) | 1L << 63;
237         final long bits = rng.nextLong() & mask;
238         // The exponent must be unsigned so + 1023 to the signed exponent
239         final long exp = rng.nextInt(129) - 64 + 1023;
240         return Double.longBitsToDouble(bits | (exp << 52));
241     }
242 
243     /**
244      * Creates a random double number with a random sign and uniform range.
245      *
246      * @param rng Random number generator.
247      * @return the random number
248      */
249     private static double createUniformNumber(UniformRandomProvider rng) {
250         // Note: [0, 1) - 1 is [-1, 0).
251         // Since the 1 is a 50/50 sample the result is the interval [-1, 1)
252         // using the 2^54 dyadic rationals in the interval.
253         // The range is not critical. The numbers will have approximately 50%
254         // with the same exponent, max, matching that of RANGE and the rest smaller
255         // exponents down to (max - 53) since the uniform deviate is limited to 2^-53.
256         return (rng.nextDouble() - rng.nextInt(1)) * RANGE;
257     }
258 
259     /**
260      * Creates a random double number that will be an edge case:
261      * {@code +/-inf, +/-max, +/-min, +/-0, nan}.
262      *
263      * @param rng Random number generator.
264      * @return the random number
265      */
266     private static double createEdgeNumber(UniformRandomProvider rng) {
267         return EDGE_NUMBERS[rng.nextInt(EDGE_NUMBERS.length)];
268     }
269 
270     /**
271      * Apply the function to all the numbers.
272      *
273      * @param numbers Numbers.
274      * @param fun Function.
275      * @param bh Data sink.
276      */
277     private static void apply(Complex[] numbers, Predicate<Complex> fun, Blackhole bh) {
278         for (int i = 0; i < numbers.length; i++) {
279             bh.consume(fun.test(numbers[i]));
280         }
281     }
282 
283     /**
284      * Apply the function to all the numbers.
285      *
286      * @param numbers Numbers.
287      * @param fun Function.
288      * @param bh Data sink.
289      */
290     private static void apply(Complex[] numbers, ToDoubleFunction<Complex> fun, Blackhole bh) {
291         for (int i = 0; i < numbers.length; i++) {
292             bh.consume(fun.applyAsDouble(numbers[i]));
293         }
294     }
295 
296     /**
297      * Apply the function to all the numbers.
298      *
299      * @param numbers Numbers.
300      * @param fun Function.
301      * @param bh Data sink.
302      */
303     private static void apply(Complex[] numbers, UnaryOperator<Complex> fun, Blackhole bh) {
304         for (int i = 0; i < numbers.length; i++) {
305             bh.consume(fun.apply(numbers[i]));
306         }
307     }
308 
309     /**
310      * Apply the function to the paired numbers.
311      *
312      * @param numbers First numbers of the pairs.
313      * @param numbers2 Second numbers of the pairs.
314      * @param fun Function.
315      * @param bh Data sink.
316      */
317     private static void apply(Complex[] numbers, Complex[] numbers2,
318             BiFunction<Complex, Complex, Complex> fun, Blackhole bh) {
319         for (int i = 0; i < numbers.length; i++) {
320             bh.consume(fun.apply(numbers[i], numbers2[i]));
321         }
322     }
323 
324     /**
325      * Apply the function to the paired numbers.
326      *
327      * @param numbers First numbers of the pairs.
328      * @param numbers2 Second numbers of the pairs.
329      * @param fun Function.
330      * @param bh Data sink.
331      */
332     private static void apply(Complex[] numbers, double[] numbers2,
333             ComplexRealFunction fun, Blackhole bh) {
334         for (int i = 0; i < numbers.length; i++) {
335             bh.consume(fun.apply(numbers[i], numbers2[i]));
336         }
337     }
338 
339     /**
340      * Identity function. This can be used to measure overhead of object array creation.
341      *
342      * @param z Complex number.
343      * @return the complex number
344      */
345     private static Complex identity(Complex z) {
346         return z;
347     }
348 
349     /**
350      * Copy function. This can be used to measure overhead of object array creation plus
351      * new Complex creation.
352      *
353      * @param z Complex number.
354      * @return a copy of the complex number
355      */
356     private static Complex copy(Complex z) {
357         return Complex.ofCartesian(z.real(), z.imag());
358     }
359 
360     // Benchmark methods.
361     //
362     // The methods are partially documented as the names are self-documenting.
363     // CHECKSTYLE: stop JavadocMethod
364     // CHECKSTYLE: stop DesignForExtension
365     //
366     // Benchmarks use function references to perform different operations on the complex numbers.
367     // Tests show that explicit programming of the same benchmarks run in the same time.
368     // For reference examples are provided for the fastest operations: real() and conj().
369 
370     /**
371      * Explicit benchmark without using a method reference.
372      * This should run in the same time as {@link #real(ComplexNumbers, Blackhole)}.
373      * This is commented out as it exists for reference purposes.
374      *
375      * @param numbers Numbers.
376      * @param bh Data sink.
377      */
378     //@Benchmark
379     public void real2(ComplexNumbers numbers, Blackhole bh) {
380         final Complex[] z = numbers.getNumbers();
381         for (int i = 0; i < z.length; i++) {
382             bh.consume(z[i].real());
383         }
384     }
385 
386     /**
387      * Explicit benchmark without using a method reference.
388      * This should run in the same time as {@link #conj(ComplexNumbers, Blackhole)}.
389      * This is commented out as it exists for reference purposes.
390      *
391      * @param numbers Numbers.
392      * @param bh Data sink.
393      */
394     //@Benchmark
395     public void conj2(ComplexNumbers numbers, Blackhole bh) {
396         final Complex[] z = numbers.getNumbers();
397         for (int i = 0; i < z.length; i++) {
398             bh.consume(z[i].conj());
399         }
400     }
401 
402     /**
403      * Baseline the JMH overhead for the loop execute to consume Complex objects.
404      *
405      * @param numbers Numbers.
406      * @param bh Data sink.
407      */
408     @Benchmark
409     public void baselineIdentity(ComplexNumbers numbers, Blackhole bh) {
410         apply(numbers.getNumbers(), ComplexPerformance::identity, bh);
411     }
412 
413     /**
414      * Baseline the creation of a copy complex number. This
415      * measures the overhead of creation of new complex numbers including field access
416      * to the real and imaginary parts.
417      *
418      * @param numbers Numbers.
419      * @param bh Data sink.
420      */
421     @Benchmark
422     public void baselineCopy(ComplexNumbers numbers, Blackhole bh) {
423         apply(numbers.getNumbers(), ComplexPerformance::copy, bh);
424     }
425 
426     // Unary operations that a boolean
427 
428     @Benchmark
429     public void isNaN(ComplexNumbers numbers, Blackhole bh) {
430         apply(numbers.getNumbers(), Complex::isNaN, bh);
431     }
432 
433     @Benchmark
434     public void isInfinite(ComplexNumbers numbers, Blackhole bh) {
435         apply(numbers.getNumbers(), Complex::isInfinite, bh);
436     }
437 
438     @Benchmark
439     public void isFinite(ComplexNumbers numbers, Blackhole bh) {
440         apply(numbers.getNumbers(), Complex::isFinite, bh);
441     }
442 
443     // Unary operations that a double
444 
445     @Benchmark
446     public void real(ComplexNumbers numbers, Blackhole bh) {
447         apply(numbers.getNumbers(), Complex::real, bh);
448     }
449 
450     @Benchmark
451     public void imag(ComplexNumbers numbers, Blackhole bh) {
452         apply(numbers.getNumbers(), Complex::imag, bh);
453     }
454 
455     @Benchmark
456     public void abs(ComplexNumbers numbers, Blackhole bh) {
457         apply(numbers.getNumbers(), Complex::abs, bh);
458     }
459 
460     @Benchmark
461     public void arg(ComplexNumbers numbers, Blackhole bh) {
462         apply(numbers.getNumbers(), Complex::arg, bh);
463     }
464 
465     @Benchmark
466     public void norm(ComplexNumbers numbers, Blackhole bh) {
467         apply(numbers.getNumbers(), Complex::norm, bh);
468     }
469 
470     /**
471      * This test demonstrates that the method used in abs() is not as fast as using square
472      * root of the norm. The C99 standard for the abs() function requires over/underflow
473      * protection in the intermediate computation and infinity edge case handling. This
474      * has a performance overhead.
475      *
476      * @param numbers Numbers.
477      * @param bh Data sink.
478      */
479     @Benchmark
480     public void sqrtNorm(ComplexNumbers numbers, Blackhole bh) {
481         apply(numbers.getNumbers(), (ToDoubleFunction<Complex>) z -> Math.sqrt(z.norm()), bh);
482     }
483 
484     /**
485      * This test demonstrates that the {@link Math#hypot(double, double)} method
486      * is not as fast as the custom implementation in abs().
487      *
488      * @param numbers Numbers.
489      * @param bh Data sink.
490      */
491     @Benchmark
492     public void absMathHypot(ComplexNumbers numbers, Blackhole bh) {
493         apply(numbers.getNumbers(), (ToDoubleFunction<Complex>) z -> Math.hypot(z.real(), z.imag()), bh);
494     }
495 
496     // Unary operations that a complex number
497 
498     @Benchmark
499     public void conj(ComplexNumbers numbers, Blackhole bh) {
500         apply(numbers.getNumbers(), Complex::conj, bh);
501     }
502 
503     @Benchmark
504     public void negate(ComplexNumbers numbers, Blackhole bh) {
505         apply(numbers.getNumbers(), Complex::negate, bh);
506     }
507 
508     @Benchmark
509     public void proj(ComplexNumbers numbers, Blackhole bh) {
510         apply(numbers.getNumbers(), Complex::proj, bh);
511     }
512 
513     @Benchmark
514     public void cos(ComplexNumbers numbers, Blackhole bh) {
515         apply(numbers.getNumbers(), Complex::cos, bh);
516     }
517 
518     @Benchmark
519     public void cosh(ComplexNumbers numbers, Blackhole bh) {
520         apply(numbers.getNumbers(), Complex::cosh, bh);
521     }
522 
523     @Benchmark
524     public void exp(ComplexNumbers numbers, Blackhole bh) {
525         apply(numbers.getNumbers(), Complex::exp, bh);
526     }
527 
528     @Benchmark
529     public void log(ComplexNumbers numbers, Blackhole bh) {
530         apply(numbers.getNumbers(), Complex::log, bh);
531     }
532 
533     @Benchmark
534     public void log10(ComplexNumbers numbers, Blackhole bh) {
535         apply(numbers.getNumbers(), Complex::log10, bh);
536     }
537 
538     @Benchmark
539     public void sin(ComplexNumbers numbers, Blackhole bh) {
540         apply(numbers.getNumbers(), Complex::sin, bh);
541     }
542 
543     @Benchmark
544     public void sinh(ComplexNumbers numbers, Blackhole bh) {
545         apply(numbers.getNumbers(), Complex::sinh, bh);
546     }
547 
548     @Benchmark
549     public void sqrt(ComplexNumbers numbers, Blackhole bh) {
550         apply(numbers.getNumbers(), Complex::sqrt, bh);
551     }
552 
553     @Benchmark
554     public void tan(ComplexNumbers numbers, Blackhole bh) {
555         apply(numbers.getNumbers(), Complex::tan, bh);
556     }
557 
558     @Benchmark
559     public void tanh(ComplexNumbers numbers, Blackhole bh) {
560         apply(numbers.getNumbers(), Complex::tanh, bh);
561     }
562 
563     @Benchmark
564     public void acos(ComplexNumbers numbers, Blackhole bh) {
565         apply(numbers.getNumbers(), Complex::acos, bh);
566     }
567 
568     @Benchmark
569     public void acosh(ComplexNumbers numbers, Blackhole bh) {
570         apply(numbers.getNumbers(), Complex::acosh, bh);
571     }
572 
573     @Benchmark
574     public void asin(ComplexNumbers numbers, Blackhole bh) {
575         apply(numbers.getNumbers(), Complex::asin, bh);
576     }
577 
578     @Benchmark
579     public void asinh(ComplexNumbers numbers, Blackhole bh) {
580         apply(numbers.getNumbers(), Complex::asinh, bh);
581     }
582 
583     @Benchmark
584     public void atan(ComplexNumbers numbers, Blackhole bh) {
585         apply(numbers.getNumbers(), Complex::atan, bh);
586     }
587 
588     @Benchmark
589     public void atanh(ComplexNumbers numbers, Blackhole bh) {
590         apply(numbers.getNumbers(), Complex::atanh, bh);
591     }
592 
593     // Binary operations on two complex numbers.
594 
595     @Benchmark
596     public void pow(TwoComplexNumbers numbers, Blackhole bh) {
597         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::pow, bh);
598     }
599 
600     @Benchmark
601     public void multiply(TwoComplexNumbers numbers, Blackhole bh) {
602         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::multiply, bh);
603     }
604 
605     @Benchmark
606     public void divide(TwoComplexNumbers numbers, Blackhole bh) {
607         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::divide, bh);
608     }
609 
610     @Benchmark
611     public void add(TwoComplexNumbers numbers, Blackhole bh) {
612         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::add, bh);
613     }
614 
615     @Benchmark
616     public void subtract(TwoComplexNumbers numbers, Blackhole bh) {
617         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::subtract, bh);
618     }
619 
620     // Binary operations on a complex and a real number.
621     // These only benchmark methods on the real component as the
622     // following are expected to be the same speed as the real-only operations
623     // given the equivalent primitive operations:
624     // - multiplyImaginary
625     // - divideImaginary
626     // - addImaginary
627     // - subtractImaginary
628     // - subtractFrom
629     // - subtractFromImaginary
630 
631     @Benchmark
632     public void powReal(ComplexAndRealNumbers numbers, Blackhole bh) {
633         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::pow, bh);
634     }
635 
636     @Benchmark
637     public void multiplyReal(ComplexAndRealNumbers numbers, Blackhole bh) {
638         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::multiply, bh);
639     }
640 
641     @Benchmark
642     public void divideReal(ComplexAndRealNumbers numbers, Blackhole bh) {
643         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::divide, bh);
644     }
645 
646     @Benchmark
647     public void addReal(ComplexAndRealNumbers numbers, Blackhole bh) {
648         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::add, bh);
649     }
650 
651     @Benchmark
652     public void subtractReal(ComplexAndRealNumbers numbers, Blackhole bh) {
653         apply(numbers.getNumbers(), numbers.getNumbers2(), Complex::subtract, bh);
654     }
655 }