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.math4.legacy.core;
19  
20  import java.lang.reflect.Array;
21  import java.util.Arrays;
22  import java.util.Iterator;
23  import java.util.TreeSet;
24  
25  import org.apache.commons.numbers.core.Precision;
26  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
27  import org.apache.commons.math4.legacy.exception.MathArithmeticException;
28  import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
29  import org.apache.commons.math4.legacy.exception.MathInternalError;
30  import org.apache.commons.math4.legacy.exception.NoDataException;
31  import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
32  import org.apache.commons.math4.legacy.exception.NotANumberException;
33  import org.apache.commons.math4.legacy.exception.NotPositiveException;
34  import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
35  import org.apache.commons.math4.legacy.exception.NullArgumentException;
36  import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
37  import org.apache.commons.math4.legacy.exception.NotFiniteNumberException;
38  import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
39  import org.apache.commons.math4.core.jdkmath.JdkMath;
40  
41  /**
42   * Arrays utilities.
43   *
44   * @since 3.0
45   */
46  public final class MathArrays {
47  
48      /**
49       * Private constructor.
50       */
51      private MathArrays() {}
52  
53      /**
54       * Real-valued function that operate on an array or a part of it.
55       * @since 3.1
56       */
57      public interface Function {
58          /**
59           * Operates on an entire array.
60           *
61           * @param array Array to operate on.
62           * @return the result of the operation.
63           */
64          double evaluate(double[] array);
65          /**
66           * @param array Array to operate on.
67           * @param startIndex Index of the first element to take into account.
68           * @param numElements Number of elements to take into account.
69           * @return the result of the operation.
70           */
71          double evaluate(double[] array,
72                          int startIndex,
73                          int numElements);
74      }
75  
76      /**
77       * Create a copy of an array scaled by a value.
78       *
79       * @param arr Array to scale.
80       * @param val Scalar.
81       * @return scaled copy of array with each entry multiplied by val.
82       * @since 3.2
83       */
84      public static double[] scale(double val, final double[] arr) {
85          double[] newArr = new double[arr.length];
86          for (int i = 0; i < arr.length; i++) {
87              newArr[i] = arr[i] * val;
88          }
89          return newArr;
90      }
91  
92      /**
93       * <p>Multiply each element of an array by a value.</p>
94       *
95       * <p>The array is modified in place (no copy is created).</p>
96       *
97       * @param arr Array to scale
98       * @param val Scalar
99       * @since 3.2
100      */
101     public static void scaleInPlace(double val, final double[] arr) {
102         for (int i = 0; i < arr.length; i++) {
103             arr[i] *= val;
104         }
105     }
106 
107     /**
108      * Creates an array whose contents will be the element-by-element
109      * addition of the arguments.
110      *
111      * @param a First term of the addition.
112      * @param b Second term of the addition.
113      * @return a new array {@code r} where {@code r[i] = a[i] + b[i]}.
114      * @throws DimensionMismatchException if the array lengths differ.
115      * @since 3.1
116      */
117     public static double[] ebeAdd(double[] a, double[] b) {
118         checkEqualLength(a, b);
119 
120         final double[] result = a.clone();
121         for (int i = 0; i < a.length; i++) {
122             result[i] += b[i];
123         }
124         return result;
125     }
126     /**
127      * Creates an array whose contents will be the element-by-element
128      * subtraction of the second argument from the first.
129      *
130      * @param a First term.
131      * @param b Element to be subtracted.
132      * @return a new array {@code r} where {@code r[i] = a[i] - b[i]}.
133      * @throws DimensionMismatchException if the array lengths differ.
134      * @since 3.1
135      */
136     public static double[] ebeSubtract(double[] a, double[] b) {
137         checkEqualLength(a, b);
138 
139         final double[] result = a.clone();
140         for (int i = 0; i < a.length; i++) {
141             result[i] -= b[i];
142         }
143         return result;
144     }
145     /**
146      * Creates an array whose contents will be the element-by-element
147      * multiplication of the arguments.
148      *
149      * @param a First factor of the multiplication.
150      * @param b Second factor of the multiplication.
151      * @return a new array {@code r} where {@code r[i] = a[i] * b[i]}.
152      * @throws DimensionMismatchException if the array lengths differ.
153      * @since 3.1
154      */
155     public static double[] ebeMultiply(double[] a, double[] b) {
156         checkEqualLength(a, b);
157 
158         final double[] result = a.clone();
159         for (int i = 0; i < a.length; i++) {
160             result[i] *= b[i];
161         }
162         return result;
163     }
164     /**
165      * Creates an array whose contents will be the element-by-element
166      * division of the first argument by the second.
167      *
168      * @param a Numerator of the division.
169      * @param b Denominator of the division.
170      * @return a new array {@code r} where {@code r[i] = a[i] / b[i]}.
171      * @throws DimensionMismatchException if the array lengths differ.
172      * @since 3.1
173      */
174     public static double[] ebeDivide(double[] a, double[] b) {
175         checkEqualLength(a, b);
176 
177         final double[] result = a.clone();
178         for (int i = 0; i < a.length; i++) {
179             result[i] /= b[i];
180         }
181         return result;
182     }
183 
184     /**
185      * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
186      *
187      * @param p1 the first point
188      * @param p2 the second point
189      * @return the L<sub>1</sub> distance between the two points
190      * @throws DimensionMismatchException if the array lengths differ.
191      */
192     public static double distance1(double[] p1, double[] p2) {
193         checkEqualLength(p1, p2);
194         double sum = 0;
195         for (int i = 0; i < p1.length; i++) {
196             sum += JdkMath.abs(p1[i] - p2[i]);
197         }
198         return sum;
199     }
200 
201     /**
202      * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
203      *
204      * @param p1 the first point
205      * @param p2 the second point
206      * @return the L<sub>1</sub> distance between the two points
207      * @throws DimensionMismatchException if the array lengths differ.
208      */
209     public static int distance1(int[] p1, int[] p2) {
210         checkEqualLength(p1, p2);
211         int sum = 0;
212         for (int i = 0; i < p1.length; i++) {
213             sum += JdkMath.abs(p1[i] - p2[i]);
214         }
215         return sum;
216     }
217 
218     /**
219      * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
220      *
221      * @param p1 the first point
222      * @param p2 the second point
223      * @return the L<sub>2</sub> distance between the two points
224      * @throws DimensionMismatchException if the array lengths differ.
225      */
226     public static double distance(double[] p1, double[] p2) {
227         checkEqualLength(p1, p2);
228         double sum = 0;
229         for (int i = 0; i < p1.length; i++) {
230             final double dp = p1[i] - p2[i];
231             sum += dp * dp;
232         }
233         return JdkMath.sqrt(sum);
234     }
235 
236     /**
237      * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
238      *
239      * @param p1 the first point
240      * @param p2 the second point
241      * @return the L<sub>2</sub> distance between the two points
242      * @throws DimensionMismatchException if the array lengths differ.
243      */
244     public static double distance(int[] p1, int[] p2) {
245         checkEqualLength(p1, p2);
246         double sum = 0;
247         for (int i = 0; i < p1.length; i++) {
248             final double dp = (double) p1[i] - p2[i];
249             sum += dp * dp;
250         }
251         return JdkMath.sqrt(sum);
252     }
253 
254     /**
255      * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
256      *
257      * @param p1 the first point
258      * @param p2 the second point
259      * @return the L<sub>&infin;</sub> distance between the two points
260      * @throws DimensionMismatchException if the array lengths differ.
261      */
262     public static double distanceInf(double[] p1, double[] p2) {
263         checkEqualLength(p1, p2);
264         double max = 0;
265         for (int i = 0; i < p1.length; i++) {
266             max = JdkMath.max(max, JdkMath.abs(p1[i] - p2[i]));
267         }
268         return max;
269     }
270 
271     /**
272      * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
273      *
274      * @param p1 the first point
275      * @param p2 the second point
276      * @return the L<sub>&infin;</sub> distance between the two points
277      * @throws DimensionMismatchException if the array lengths differ.
278      */
279     public static int distanceInf(int[] p1, int[] p2) {
280         checkEqualLength(p1, p2);
281         int max = 0;
282         for (int i = 0; i < p1.length; i++) {
283             max = JdkMath.max(max, JdkMath.abs(p1[i] - p2[i]));
284         }
285         return max;
286     }
287 
288     /**
289      * Specification of ordering direction.
290      */
291     public enum OrderDirection {
292         /** Constant for increasing direction. */
293         INCREASING,
294         /** Constant for decreasing direction. */
295         DECREASING
296     }
297 
298     /**
299      * Check that an array is monotonically increasing or decreasing.
300      *
301      * @param <T> the type of the elements in the specified array
302      * @param val Values.
303      * @param dir Ordering direction.
304      * @param strict Whether the order should be strict.
305      * @return {@code true} if sorted, {@code false} otherwise.
306      */
307     public static  <T extends Comparable<? super T>> boolean isMonotonic(T[] val,
308                                       OrderDirection dir,
309                                       boolean strict) {
310         T previous = val[0];
311         final int max = val.length;
312         for (int i = 1; i < max; i++) {
313             final int comp;
314             switch (dir) {
315             case INCREASING:
316                 comp = previous.compareTo(val[i]);
317                 if (strict) {
318                     if (comp >= 0) {
319                         return false;
320                     }
321                 } else {
322                     if (comp > 0) {
323                         return false;
324                     }
325                 }
326                 break;
327             case DECREASING:
328                 comp = val[i].compareTo(previous);
329                 if (strict) {
330                     if (comp >= 0) {
331                         return false;
332                     }
333                 } else {
334                     if (comp > 0) {
335                         return false;
336                     }
337                 }
338                 break;
339             default:
340                 // Should never happen.
341                 throw new MathInternalError();
342             }
343 
344             previous = val[i];
345         }
346         return true;
347     }
348 
349     /**
350      * Check that an array is monotonically increasing or decreasing.
351      *
352      * @param val Values.
353      * @param dir Ordering direction.
354      * @param strict Whether the order should be strict.
355      * @return {@code true} if sorted, {@code false} otherwise.
356      */
357     public static boolean isMonotonic(double[] val, OrderDirection dir, boolean strict) {
358         return checkOrder(val, dir, strict, false);
359     }
360 
361     /**
362      * Check that both arrays have the same length.
363      *
364      * @param a Array.
365      * @param b Array.
366      * @param abort Whether to throw an exception if the check fails.
367      * @return {@code true} if the arrays have the same length.
368      * @throws DimensionMismatchException if the lengths differ and
369      * {@code abort} is {@code true}.
370      * @since 3.6
371      */
372     public static boolean checkEqualLength(double[] a,
373                                            double[] b,
374                                            boolean abort) {
375         if (a.length == b.length) {
376             return true;
377         } else {
378             if (abort) {
379                 throw new DimensionMismatchException(a.length, b.length);
380             }
381             return false;
382         }
383     }
384 
385     /**
386      * Check that both arrays have the same length.
387      *
388      * @param a Array.
389      * @param b Array.
390      * @throws DimensionMismatchException if the lengths differ.
391      * @since 3.6
392      */
393     public static void checkEqualLength(double[] a,
394                                         double[] b) {
395         checkEqualLength(a, b, true);
396     }
397 
398 
399     /**
400      * Check that both arrays have the same length.
401      *
402      * @param a Array.
403      * @param b Array.
404      * @param abort Whether to throw an exception if the check fails.
405      * @return {@code true} if the arrays have the same length.
406      * @throws DimensionMismatchException if the lengths differ and
407      * {@code abort} is {@code true}.
408      * @since 3.6
409      */
410     public static boolean checkEqualLength(int[] a,
411                                            int[] b,
412                                            boolean abort) {
413         if (a.length == b.length) {
414             return true;
415         } else {
416             if (abort) {
417                 throw new DimensionMismatchException(a.length, b.length);
418             }
419             return false;
420         }
421     }
422 
423     /**
424      * Check that both arrays have the same length.
425      *
426      * @param a Array.
427      * @param b Array.
428      * @throws DimensionMismatchException if the lengths differ.
429      * @since 3.6
430      */
431     public static void checkEqualLength(int[] a,
432                                         int[] b) {
433         checkEqualLength(a, b, true);
434     }
435 
436     /**
437      * Check that the given array is sorted.
438      *
439      * @param val Values.
440      * @param dir Ordering direction.
441      * @param strict Whether the order should be strict.
442      * @param abort Whether to throw an exception if the check fails.
443      * @return {@code true} if the array is sorted.
444      * @throws NonMonotonicSequenceException if the array is not sorted
445      * and {@code abort} is {@code true}.
446      */
447     public static boolean checkOrder(double[] val, OrderDirection dir,
448                                      boolean strict, boolean abort) {
449         double previous = val[0];
450         final int max = val.length;
451 
452         int index;
453         ITEM:
454         for (index = 1; index < max; index++) {
455             switch (dir) {
456             case INCREASING:
457                 if (strict) {
458                     if (val[index] <= previous) {
459                         break ITEM;
460                     }
461                 } else {
462                     if (val[index] < previous) {
463                         break ITEM;
464                     }
465                 }
466                 break;
467             case DECREASING:
468                 if (strict) {
469                     if (val[index] >= previous) {
470                         break ITEM;
471                     }
472                 } else {
473                     if (val[index] > previous) {
474                         break ITEM;
475                     }
476                 }
477                 break;
478             default:
479                 // Should never happen.
480                 throw new MathInternalError();
481             }
482 
483             previous = val[index];
484         }
485 
486         if (index == max) {
487             // Loop completed.
488             return true;
489         }
490 
491         // Loop early exit means wrong ordering.
492         if (abort) {
493             throw new NonMonotonicSequenceException(val[index],
494                                                     previous,
495                                                     index,
496                                                     dir == OrderDirection.INCREASING,
497                                                     strict);
498         } else {
499             return false;
500         }
501     }
502 
503     /**
504      * Check that the given array is sorted.
505      *
506      * @param val Values.
507      * @param dir Ordering direction.
508      * @param strict Whether the order should be strict.
509      * @throws NonMonotonicSequenceException if the array is not sorted.
510      * @since 2.2
511      */
512     public static void checkOrder(double[] val, OrderDirection dir, boolean strict) {
513         checkOrder(val, dir, strict, true);
514     }
515 
516     /**
517      * Check that the given array is sorted in strictly increasing order.
518      *
519      * @param val Values.
520      * @throws NonMonotonicSequenceException if the array is not sorted.
521      * @since 2.2
522      */
523     public static void checkOrder(double[] val) {
524         checkOrder(val, OrderDirection.INCREASING, true);
525     }
526 
527     /**
528      * Throws DimensionMismatchException if the input array is not rectangular.
529      *
530      * @param in array to be tested
531      * @throws NullArgumentException if input array is null
532      * @throws DimensionMismatchException if input array is not rectangular
533      * @since 3.1
534      */
535     public static void checkRectangular(final long[][] in) {
536         NullArgumentException.check(in);
537         for (int i = 1; i < in.length; i++) {
538             if (in[i].length != in[0].length) {
539                 throw new DimensionMismatchException(
540                         LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
541                         in[i].length, in[0].length);
542             }
543         }
544     }
545 
546     /**
547      * Check that all entries of the input array are strictly positive.
548      *
549      * @param in Array to be tested
550      * @throws NotStrictlyPositiveException if any entries of the array are not
551      * strictly positive.
552      * @since 3.1
553      */
554     public static void checkPositive(final double[] in) {
555         for (double x : in) {
556             if (x <= 0) {
557                 throw new NotStrictlyPositiveException(x);
558             }
559         }
560     }
561 
562     /**
563      * Check that no entry of the input array is {@code NaN}.
564      *
565      * @param in Array to be tested.
566      * @throws NotANumberException if an entry is {@code NaN}.
567      * @since 3.4
568      */
569     public static void checkNotNaN(final double[] in) {
570         for (double x : in) {
571             if (Double.isNaN(x)) {
572                 throw new NotANumberException();
573             }
574         }
575     }
576 
577     /**
578      * Check that all the elements are real numbers.
579      *
580      * @param val Arguments.
581      * @throws NotFiniteNumberException if any values of the array is not a
582      * finite real number.
583      */
584     public static void checkFinite(final double[] val) {
585         for (double x : val) {
586             if (!Double.isFinite(x)) {
587                 throw new NotFiniteNumberException(x);
588             }
589         }
590     }
591 
592     /**
593      * Check that all entries of the input array are &gt;= 0.
594      *
595      * @param in Array to be tested
596      * @throws NotPositiveException if any array entries are less than 0.
597      * @since 3.1
598      */
599     public static void checkNonNegative(final long[] in) {
600         for (long i : in) {
601             if (i < 0) {
602                 throw new NotPositiveException(i);
603             }
604         }
605     }
606 
607     /**
608      * Check all entries of the input array are &gt;= 0.
609      *
610      * @param in Array to be tested
611      * @throws NotPositiveException if any array entries are less than 0.
612      * @since 3.1
613      */
614     public static void checkNonNegative(final long[][] in) {
615         for (int i = 0; i < in.length; i++) {
616             for (int j = 0; j < in[i].length; j++) {
617                 if (in[i][j] < 0) {
618                     throw new NotPositiveException(in[i][j]);
619                 }
620             }
621         }
622     }
623 
624     /**
625      * Returns true iff both arguments are null or have same dimensions and all
626      * their elements are equal as defined by
627      * {@link Precision#equals(float,float)}.
628      *
629      * @param x first array
630      * @param y second array
631      * @return true if the values are both null or have same dimension
632      * and equal elements.
633      */
634     public static boolean equals(float[] x, float[] y) {
635         if (x == null || y == null) {
636             return (x == null) == (y == null);
637         }
638         if (x.length != y.length) {
639             return false;
640         }
641         for (int i = 0; i < x.length; ++i) {
642             if (!Precision.equals(x[i], y[i])) {
643                 return false;
644             }
645         }
646         return true;
647     }
648 
649     /**
650      * Returns true iff both arguments are null or have same dimensions and all
651      * their elements are equal as defined by
652      * {@link Precision#equalsIncludingNaN(double,double) this method}.
653      *
654      * @param x first array
655      * @param y second array
656      * @return true if the values are both null or have same dimension and
657      * equal elements
658      * @since 2.2
659      */
660     public static boolean equalsIncludingNaN(float[] x, float[] y) {
661         if (x == null || y == null) {
662             return (x == null) == (y == null);
663         }
664         if (x.length != y.length) {
665             return false;
666         }
667         for (int i = 0; i < x.length; ++i) {
668             if (!Precision.equalsIncludingNaN(x[i], y[i])) {
669                 return false;
670             }
671         }
672         return true;
673     }
674 
675     /**
676      * Returns {@code true} iff both arguments are {@code null} or have same
677      * dimensions and all their elements are equal as defined by
678      * {@link Precision#equals(double,double)}.
679      *
680      * @param x First array.
681      * @param y Second array.
682      * @return {@code true} if the values are both {@code null} or have same
683      * dimension and equal elements.
684      */
685     public static boolean equals(double[] x, double[] y) {
686         if (x == null || y == null) {
687             return (x == null) == (y == null);
688         }
689         if (x.length != y.length) {
690             return false;
691         }
692         for (int i = 0; i < x.length; ++i) {
693             if (!Precision.equals(x[i], y[i])) {
694                 return false;
695             }
696         }
697         return true;
698     }
699 
700     /**
701      * Returns {@code true} iff both arguments are {@code null} or have same
702      * dimensions and all their elements are equal as defined by
703      * {@link Precision#equalsIncludingNaN(double,double) this method}.
704      *
705      * @param x First array.
706      * @param y Second array.
707      * @return {@code true} if the values are both {@code null} or have same
708      * dimension and equal elements.
709      * @since 2.2
710      */
711     public static boolean equalsIncludingNaN(double[] x, double[] y) {
712         if (x == null || y == null) {
713             return (x == null) == (y == null);
714         }
715         if (x.length != y.length) {
716             return false;
717         }
718         for (int i = 0; i < x.length; ++i) {
719             if (!Precision.equalsIncludingNaN(x[i], y[i])) {
720                 return false;
721             }
722         }
723         return true;
724     }
725 
726     /**
727      * Normalizes an array to make it sum to a specified value.
728      * Returns the result of the transformation
729      * <pre>
730      *    x |-&gt; x * normalizedSum / sum
731      * </pre>
732      * applied to each non-NaN element x of the input array, where sum is the
733      * sum of the non-NaN entries in the input array.
734      * <p>
735      * Throws IllegalArgumentException if {@code normalizedSum} is infinite
736      * or NaN and ArithmeticException if the input array contains any infinite elements
737      * or sums to 0.
738      * <p>
739      * Ignores (i.e., copies unchanged to the output array) NaNs in the input array.
740      *
741      * @param values Input array to be normalized
742      * @param normalizedSum Target sum for the normalized array
743      * @return the normalized array.
744      * @throws MathArithmeticException if the input array contains infinite
745      * elements or sums to zero.
746      * @throws MathIllegalArgumentException if the target sum is infinite or {@code NaN}.
747      * @since 2.1
748      */
749     public static double[] normalizeArray(double[] values, double normalizedSum) {
750         if (Double.isInfinite(normalizedSum)) {
751             throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_INFINITE);
752         }
753         if (Double.isNaN(normalizedSum)) {
754             throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_NAN);
755         }
756         double sum = 0d;
757         final int len = values.length;
758         double[] out = new double[len];
759         for (int i = 0; i < len; i++) {
760             if (Double.isInfinite(values[i])) {
761                 throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i);
762             }
763             if (!Double.isNaN(values[i])) {
764                 sum += values[i];
765             }
766         }
767         if (sum == 0) {
768             throw new MathArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO);
769         }
770         final double scale = normalizedSum / sum;
771         for (int i = 0; i < len; i++) {
772             if (Double.isNaN(values[i])) {
773                 out[i] = Double.NaN;
774             } else {
775                 out[i] = values[i] * scale;
776             }
777         }
778         return out;
779     }
780 
781     /** Build an array of elements.
782      * <p>
783      * Arrays are filled with field.getZero()
784      *
785      * @param <T> the type of the field elements
786      * @param field field to which array elements belong
787      * @param length of the array
788      * @return a new array
789      * @since 3.2
790      */
791     public static <T> T[] buildArray(final Field<T> field, final int length) {
792         @SuppressWarnings("unchecked") // OK because field must be correct class
793         T[] array = (T[]) Array.newInstance(field.getRuntimeClass(), length);
794         Arrays.fill(array, field.getZero());
795         return array;
796     }
797 
798     /** Build a double dimension  array of elements.
799      * <p>
800      * Arrays are filled with field.getZero()
801      *
802      * @param <T> the type of the field elements
803      * @param field field to which array elements belong
804      * @param rows number of rows in the array
805      * @param columns number of columns (may be negative to build partial
806      * arrays in the same way <code>new Field[rows][]</code> works)
807      * @return a new array
808      * @since 3.2
809      */
810     @SuppressWarnings("unchecked")
811     public static <T> T[][] buildArray(final Field<T> field, final int rows, final int columns) {
812         final T[][] array;
813         if (columns < 0) {
814             T[] dummyRow = buildArray(field, 0);
815             array = (T[][]) Array.newInstance(dummyRow.getClass(), rows);
816         } else {
817             array = (T[][]) Array.newInstance(field.getRuntimeClass(),
818                                               rows, columns);
819             for (int i = 0; i < rows; ++i) {
820                 Arrays.fill(array[i], field.getZero());
821             }
822         }
823         return array;
824     }
825 
826     /**
827      * Calculates the <a href="http://en.wikipedia.org/wiki/Convolution">
828      * convolution</a> between two sequences.
829      * <p>
830      * The solution is obtained via straightforward computation of the
831      * convolution sum (and not via FFT). Whenever the computation needs
832      * an element that would be located at an index outside the input arrays,
833      * the value is assumed to be zero.
834      *
835      * @param x First sequence.
836      * Typically, this sequence will represent an input signal to a system.
837      * @param h Second sequence.
838      * Typically, this sequence will represent the impulse response of the system.
839      * @return the convolution of {@code x} and {@code h}.
840      * This array's length will be {@code x.length + h.length - 1}.
841      * @throws NullArgumentException if either {@code x} or {@code h} is {@code null}.
842      * @throws NoDataException if either {@code x} or {@code h} is empty.
843      *
844      * @since 3.3
845      */
846     public static double[] convolve(double[] x, double[] h) {
847         NullArgumentException.check(x);
848         NullArgumentException.check(h);
849 
850         final int xLen = x.length;
851         final int hLen = h.length;
852 
853         if (xLen == 0 || hLen == 0) {
854             throw new NoDataException();
855         }
856 
857         // initialize the output array
858         final int totalLength = xLen + hLen - 1;
859         final double[] y = new double[totalLength];
860 
861         // straightforward implementation of the convolution sum
862         for (int n = 0; n < totalLength; n++) {
863             double yn = 0;
864             int k = JdkMath.max(0, n + 1 - xLen);
865             int j = n - k;
866             while (k < hLen && j >= 0) {
867                 yn += x[j--] * h[k++];
868             }
869             y[n] = yn;
870         }
871 
872         return y;
873     }
874 
875     /**
876      * Returns an array representing the natural number {@code n}.
877      *
878      * @param n Natural number.
879      * @return an array whose entries are the numbers 0, 1, ..., {@code n}-1.
880      * If {@code n == 0}, the returned array is empty.
881      */
882     public static int[] natural(int n) {
883         return sequence(n, 0, 1);
884     }
885     /**
886      * Returns an array of {@code size} integers starting at {@code start},
887      * skipping {@code stride} numbers.
888      *
889      * @param size Natural number.
890      * @param start Natural number.
891      * @param stride Natural number.
892      * @return an array whose entries are the numbers
893      * {@code start, start + stride, ..., start + (size - 1) * stride}.
894      * If {@code size == 0}, the returned array is empty.
895      *
896      * @since 3.4
897      */
898     public static int[] sequence(int size,
899                                  int start,
900                                  int stride) {
901         final int[] a = new int[size];
902         for (int i = 0; i < size; i++) {
903             a[i] = start + i * stride;
904         }
905         return a;
906     }
907     /**
908      * This method is used
909      * to verify that the input parameters designate a subarray of positive length.
910      * <ul>
911      * <li>returns <code>true</code> iff the parameters designate a subarray of
912      * positive length</li>
913      * <li>throws <code>MathIllegalArgumentException</code> if the array is null or
914      * or the indices are invalid</li>
915      * <li>returns <code>false</code> if the array is non-null, but
916      * <code>length</code> is 0.</li>
917      * </ul>
918      *
919      * @param values the input array
920      * @param begin index of the first array element to include
921      * @param length the number of elements to include
922      * @return true if the parameters are valid and designate a subarray of positive length
923      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
924      * @since 3.3
925      */
926     public static boolean verifyValues(final double[] values, final int begin, final int length) {
927         return verifyValues(values, begin, length, false);
928     }
929 
930     /**
931      * This method is used
932      * to verify that the input parameters designate a subarray of positive length.
933      * <ul>
934      * <li>returns <code>true</code> iff the parameters designate a subarray of
935      * non-negative length</li>
936      * <li>throws <code>IllegalArgumentException</code> if the array is null or
937      * or the indices are invalid</li>
938      * <li>returns <code>false</code> if the array is non-null, but
939      * <code>length</code> is 0 unless <code>allowEmpty</code> is <code>true</code></li>
940      * </ul>
941      *
942      * @param values the input array
943      * @param begin index of the first array element to include
944      * @param length the number of elements to include
945      * @param allowEmpty if <code>true</code> then zero length arrays are allowed
946      * @return true if the parameters are valid
947      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
948      * @since 3.3
949      */
950     public static boolean verifyValues(final double[] values, final int begin,
951                                        final int length, final boolean allowEmpty) {
952 
953         if (values == null) {
954             throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
955         }
956 
957         if (begin < 0) {
958             throw new NotPositiveException(LocalizedFormats.START_POSITION, Integer.valueOf(begin));
959         }
960 
961         if (length < 0) {
962             throw new NotPositiveException(LocalizedFormats.LENGTH, Integer.valueOf(length));
963         }
964 
965         if (begin + length > values.length) {
966             throw new NumberIsTooLargeException(LocalizedFormats.SUBARRAY_ENDS_AFTER_ARRAY_END,
967                     Integer.valueOf(begin + length), Integer.valueOf(values.length), true);
968         }
969 
970         return !(length == 0 && !allowEmpty);
971     }
972 
973     /**
974      * This method is used
975      * to verify that the begin and length parameters designate a subarray of positive length
976      * and the weights are all non-negative, non-NaN, finite, and not all zero.
977      * <ul>
978      * <li>returns <code>true</code> iff the parameters designate a subarray of
979      * positive length and the weights array contains legitimate values.</li>
980      * <li>throws <code>IllegalArgumentException</code> if any of the following are true:
981      * <ul><li>the values array is null</li>
982      *     <li>the weights array is null</li>
983      *     <li>the weights array does not have the same length as the values array</li>
984      *     <li>the weights array contains one or more infinite values</li>
985      *     <li>the weights array contains one or more NaN values</li>
986      *     <li>the weights array contains negative values</li>
987      *     <li>the weights array does not contain at least one non-zero value (applies when length is non zero)</li>
988      *     <li>the start and length arguments do not determine a valid array</li></ul>
989      * </li>
990      * <li>returns <code>false</code> if the array is non-null, but
991      * <code>length</code> is 0.</li>
992      * </ul>
993      *
994      * @param values the input array
995      * @param weights the weights array
996      * @param begin index of the first array element to include
997      * @param length the number of elements to include
998      * @return true if the parameters are valid and designate a subarray of positive length
999      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
1000      * @since 3.3
1001      */
1002     public static boolean verifyValues(
1003         final double[] values,
1004         final double[] weights,
1005         final int begin,
1006         final int length) {
1007         return verifyValues(values, weights, begin, length, false);
1008     }
1009 
1010     /**
1011      * This method is used
1012      * to verify that the begin and length parameters designate a subarray of positive length
1013      * and the weights are all non-negative, non-NaN, finite, and not all zero.
1014      * <ul>
1015      * <li>returns <code>true</code> iff the parameters designate a subarray of
1016      * non-negative length and the weights array contains legitimate values.</li>
1017      * <li>throws <code>MathIllegalArgumentException</code> if any of the following are true:
1018      * <ul><li>the values array is null</li>
1019      *     <li>the weights array is null</li>
1020      *     <li>the weights array does not have the same length as the values array</li>
1021      *     <li>the weights array contains one or more infinite values</li>
1022      *     <li>the weights array contains one or more NaN values</li>
1023      *     <li>the weights array contains negative values</li>
1024      *     <li>the weights array does not contain at least one non-zero value (applies when length is non zero)</li>
1025      *     <li>the start and length arguments do not determine a valid array</li></ul>
1026      * </li>
1027      * <li>returns <code>false</code> if the array is non-null, but
1028      * <code>length</code> is 0 unless <code>allowEmpty</code> is <code>true</code>.</li>
1029      * </ul>
1030      *
1031      * @param values the input array.
1032      * @param weights the weights array.
1033      * @param begin index of the first array element to include.
1034      * @param length the number of elements to include.
1035      * @param allowEmpty if {@code true} than allow zero length arrays to pass.
1036      * @return {@code true} if the parameters are valid.
1037      * @throws NullArgumentException if either of the arrays are null
1038      * @throws MathIllegalArgumentException if the array indices are not valid,
1039      * the weights array contains NaN, infinite or negative elements, or there
1040      * are no positive weights.
1041      * @since 3.3
1042      */
1043     public static boolean verifyValues(final double[] values, final double[] weights,
1044                                        final int begin, final int length, final boolean allowEmpty) {
1045 
1046         if (weights == null || values == null) {
1047             throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
1048         }
1049 
1050         checkEqualLength(weights, values);
1051 
1052         if (length != 0) {
1053             boolean containsPositiveWeight = false;
1054             for (int i = begin; i < begin + length; i++) {
1055                 final double weight = weights[i];
1056                 if (Double.isNaN(weight)) {
1057                     throw new MathIllegalArgumentException(LocalizedFormats.NAN_ELEMENT_AT_INDEX, Integer.valueOf(i));
1058                 }
1059                 if (Double.isInfinite(weight)) {
1060                     throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT,
1061                         Double.valueOf(weight), Integer.valueOf(i));
1062                 }
1063                 if (weight < 0) {
1064                     throw new MathIllegalArgumentException(LocalizedFormats.NEGATIVE_ELEMENT_AT_INDEX,
1065                         Integer.valueOf(i), Double.valueOf(weight));
1066                 }
1067                 if (!containsPositiveWeight && weight > 0.0) {
1068                     containsPositiveWeight = true;
1069                 }
1070             }
1071 
1072             if (!containsPositiveWeight) {
1073                 throw new MathIllegalArgumentException(LocalizedFormats.WEIGHT_AT_LEAST_ONE_NON_ZERO);
1074             }
1075         }
1076 
1077         return verifyValues(values, begin, length, allowEmpty);
1078     }
1079 
1080     /**
1081      * Concatenates a sequence of arrays. The return array consists of the
1082      * entries of the input arrays concatenated in the order they appear in
1083      * the argument list.  Null arrays cause NullPointerExceptions; zero
1084      * length arrays are allowed (contributing nothing to the output array).
1085      *
1086      * @param x list of double[] arrays to concatenate
1087      * @return a new array consisting of the entries of the argument arrays
1088      * @throws NullPointerException if any of the arrays are null
1089      * @since 3.6
1090      */
1091     public static double[] concatenate(double[]... x) {
1092         int combinedLength = 0;
1093         for (double[] a : x) {
1094             combinedLength += a.length;
1095         }
1096         int offset = 0;
1097         int curLength = 0;
1098         final double[] combined = new double[combinedLength];
1099         for (int i = 0; i < x.length; i++) {
1100             curLength = x[i].length;
1101             System.arraycopy(x[i], 0, combined, offset, curLength);
1102             offset += curLength;
1103         }
1104         return combined;
1105     }
1106 
1107     /**
1108      * Returns an array consisting of the unique values in {@code data}.
1109      * The return array is sorted in descending order.  Empty arrays
1110      * are allowed, but null arrays result in NullPointerException.
1111      * Infinities are allowed.  NaN values are allowed with maximum
1112      * sort order - i.e., if there are NaN values in {@code data},
1113      * {@code Double.NaN} will be the first element of the output array,
1114      * even if the array also contains {@code Double.POSITIVE_INFINITY}.
1115      *
1116      * @param data array to scan
1117      * @return descending list of values included in the input array
1118      * @throws NullPointerException if data is null
1119      * @since 3.6
1120      */
1121     public static double[] unique(double[] data) {
1122         TreeSet<Double> values = new TreeSet<>();
1123         for (int i = 0; i < data.length; i++) {
1124             values.add(data[i]);
1125         }
1126         final int count = values.size();
1127         final double[] out = new double[count];
1128         Iterator<Double> iterator = values.descendingIterator();
1129         int i = 0;
1130         while (iterator.hasNext()) {
1131             out[i++] = iterator.next();
1132         }
1133         return out;
1134     }
1135 }