All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.apache.commons.math3.util.MathArrays Maven / Gradle / Ivy

There is a newer version: 5.0.70
Show newest version
/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math3.util;

import java.lang.reflect.Array;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.Comparator;
import java.util.Iterator;
import java.util.List;
import java.util.TreeSet;

import org.apache.commons.math3.Field;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;
import org.apache.commons.math3.distribution.UniformIntegerDistribution;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MathArithmeticException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.MathInternalError;
import org.apache.commons.math3.exception.NoDataException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NotPositiveException;
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.exception.NotANumberException;
import org.apache.commons.math3.exception.util.LocalizedFormats;

/**
 * Arrays utilities.
 *
 * @since 3.0
 */
public class MathArrays {

    /**
     * Private constructor.
     */
    private MathArrays() {}

    /**
     * Real-valued function that operate on an array or a part of it.
     * @since 3.1
     */
    public interface Function {
        /**
         * Operates on an entire array.
         *
         * @param array Array to operate on.
         * @return the result of the operation.
         */
        double evaluate(double[] array);
        /**
         * @param array Array to operate on.
         * @param startIndex Index of the first element to take into account.
         * @param numElements Number of elements to take into account.
         * @return the result of the operation.
         */
        double evaluate(double[] array,
                        int startIndex,
                        int numElements);
    }

    /**
     * Create a copy of an array scaled by a value.
     *
     * @param arr Array to scale.
     * @param val Scalar.
     * @return scaled copy of array with each entry multiplied by val.
     * @since 3.2
     */
    public static double[] scale(double val, final double[] arr) {
        double[] newArr = new double[arr.length];
        for (int i = 0; i < arr.length; i++) {
            newArr[i] = arr[i] * val;
        }
        return newArr;
    }

    /**
     * 

Multiply each element of an array by a value.

* *

The array is modified in place (no copy is created).

* * @param arr Array to scale * @param val Scalar * @since 3.2 */ public static void scaleInPlace(double val, final double[] arr) { for (int i = 0; i < arr.length; i++) { arr[i] *= val; } } /** * Creates an array whose contents will be the element-by-element * addition of the arguments. * * @param a First term of the addition. * @param b Second term of the addition. * @return a new array {@code r} where {@code r[i] = a[i] + b[i]}. * @throws DimensionMismatchException if the array lengths differ. * @since 3.1 */ public static double[] ebeAdd(double[] a, double[] b) throws DimensionMismatchException { checkEqualLength(a, b); final double[] result = a.clone(); for (int i = 0; i < a.length; i++) { result[i] += b[i]; } return result; } /** * Creates an array whose contents will be the element-by-element * subtraction of the second argument from the first. * * @param a First term. * @param b Element to be subtracted. * @return a new array {@code r} where {@code r[i] = a[i] - b[i]}. * @throws DimensionMismatchException if the array lengths differ. * @since 3.1 */ public static double[] ebeSubtract(double[] a, double[] b) throws DimensionMismatchException { checkEqualLength(a, b); final double[] result = a.clone(); for (int i = 0; i < a.length; i++) { result[i] -= b[i]; } return result; } /** * Creates an array whose contents will be the element-by-element * multiplication of the arguments. * * @param a First factor of the multiplication. * @param b Second factor of the multiplication. * @return a new array {@code r} where {@code r[i] = a[i] * b[i]}. * @throws DimensionMismatchException if the array lengths differ. * @since 3.1 */ public static double[] ebeMultiply(double[] a, double[] b) throws DimensionMismatchException { checkEqualLength(a, b); final double[] result = a.clone(); for (int i = 0; i < a.length; i++) { result[i] *= b[i]; } return result; } /** * Creates an array whose contents will be the element-by-element * division of the first argument by the second. * * @param a Numerator of the division. * @param b Denominator of the division. * @return a new array {@code r} where {@code r[i] = a[i] / b[i]}. * @throws DimensionMismatchException if the array lengths differ. * @since 3.1 */ public static double[] ebeDivide(double[] a, double[] b) throws DimensionMismatchException { checkEqualLength(a, b); final double[] result = a.clone(); for (int i = 0; i < a.length; i++) { result[i] /= b[i]; } return result; } /** * Calculates the L1 (sum of abs) distance between two points. * * @param p1 the first point * @param p2 the second point * @return the L1 distance between the two points * @throws DimensionMismatchException if the array lengths differ. */ public static double distance1(double[] p1, double[] p2) throws DimensionMismatchException { checkEqualLength(p1, p2); double sum = 0; for (int i = 0; i < p1.length; i++) { sum += FastMath.abs(p1[i] - p2[i]); } return sum; } /** * Calculates the L1 (sum of abs) distance between two points. * * @param p1 the first point * @param p2 the second point * @return the L1 distance between the two points * @throws DimensionMismatchException if the array lengths differ. */ public static int distance1(int[] p1, int[] p2) throws DimensionMismatchException { checkEqualLength(p1, p2); int sum = 0; for (int i = 0; i < p1.length; i++) { sum += FastMath.abs(p1[i] - p2[i]); } return sum; } /** * Calculates the L2 (Euclidean) distance between two points. * * @param p1 the first point * @param p2 the second point * @return the L2 distance between the two points * @throws DimensionMismatchException if the array lengths differ. */ public static double distance(double[] p1, double[] p2) throws DimensionMismatchException { checkEqualLength(p1, p2); double sum = 0; for (int i = 0; i < p1.length; i++) { final double dp = p1[i] - p2[i]; sum += dp * dp; } return FastMath.sqrt(sum); } /** * Calculates the cosine of the angle between two vectors. * * @param v1 Cartesian coordinates of the first vector. * @param v2 Cartesian coordinates of the second vector. * @return the cosine of the angle between the vectors. * @since 3.6 */ public static double cosAngle(double[] v1, double[] v2) { return linearCombination(v1, v2) / (safeNorm(v1) * safeNorm(v2)); } /** * Calculates the L2 (Euclidean) distance between two points. * * @param p1 the first point * @param p2 the second point * @return the L2 distance between the two points * @throws DimensionMismatchException if the array lengths differ. */ public static double distance(int[] p1, int[] p2) throws DimensionMismatchException { checkEqualLength(p1, p2); double sum = 0; for (int i = 0; i < p1.length; i++) { final double dp = p1[i] - p2[i]; sum += dp * dp; } return FastMath.sqrt(sum); } /** * Calculates the L (max of abs) distance between two points. * * @param p1 the first point * @param p2 the second point * @return the L distance between the two points * @throws DimensionMismatchException if the array lengths differ. */ public static double distanceInf(double[] p1, double[] p2) throws DimensionMismatchException { checkEqualLength(p1, p2); double max = 0; for (int i = 0; i < p1.length; i++) { max = FastMath.max(max, FastMath.abs(p1[i] - p2[i])); } return max; } /** * Calculates the L (max of abs) distance between two points. * * @param p1 the first point * @param p2 the second point * @return the L distance between the two points * @throws DimensionMismatchException if the array lengths differ. */ public static int distanceInf(int[] p1, int[] p2) throws DimensionMismatchException { checkEqualLength(p1, p2); int max = 0; for (int i = 0; i < p1.length; i++) { max = FastMath.max(max, FastMath.abs(p1[i] - p2[i])); } return max; } /** * Specification of ordering direction. */ public enum OrderDirection { /** Constant for increasing direction. */ INCREASING, /** Constant for decreasing direction. */ DECREASING } /** * Check that an array is monotonically increasing or decreasing. * * @param the type of the elements in the specified array * @param val Values. * @param dir Ordering direction. * @param strict Whether the order should be strict. * @return {@code true} if sorted, {@code false} otherwise. */ public static > boolean isMonotonic(T[] val, OrderDirection dir, boolean strict) { T previous = val[0]; final int max = val.length; for (int i = 1; i < max; i++) { final int comp; switch (dir) { case INCREASING: comp = previous.compareTo(val[i]); if (strict) { if (comp >= 0) { return false; } } else { if (comp > 0) { return false; } } break; case DECREASING: comp = val[i].compareTo(previous); if (strict) { if (comp >= 0) { return false; } } else { if (comp > 0) { return false; } } break; default: // Should never happen. throw new MathInternalError(); } previous = val[i]; } return true; } /** * Check that an array is monotonically increasing or decreasing. * * @param val Values. * @param dir Ordering direction. * @param strict Whether the order should be strict. * @return {@code true} if sorted, {@code false} otherwise. */ public static boolean isMonotonic(double[] val, OrderDirection dir, boolean strict) { return checkOrder(val, dir, strict, false); } /** * Check that both arrays have the same length. * * @param a Array. * @param b Array. * @param abort Whether to throw an exception if the check fails. * @return {@code true} if the arrays have the same length. * @throws DimensionMismatchException if the lengths differ and * {@code abort} is {@code true}. * @since 3.6 */ public static boolean checkEqualLength(double[] a, double[] b, boolean abort) { if (a.length == b.length) { return true; } else { if (abort) { throw new DimensionMismatchException(a.length, b.length); } return false; } } /** * Check that both arrays have the same length. * * @param a Array. * @param b Array. * @throws DimensionMismatchException if the lengths differ. * @since 3.6 */ public static void checkEqualLength(double[] a, double[] b) { checkEqualLength(a, b, true); } /** * Check that both arrays have the same length. * * @param a Array. * @param b Array. * @param abort Whether to throw an exception if the check fails. * @return {@code true} if the arrays have the same length. * @throws DimensionMismatchException if the lengths differ and * {@code abort} is {@code true}. * @since 3.6 */ public static boolean checkEqualLength(int[] a, int[] b, boolean abort) { if (a.length == b.length) { return true; } else { if (abort) { throw new DimensionMismatchException(a.length, b.length); } return false; } } /** * Check that both arrays have the same length. * * @param a Array. * @param b Array. * @throws DimensionMismatchException if the lengths differ. * @since 3.6 */ public static void checkEqualLength(int[] a, int[] b) { checkEqualLength(a, b, true); } /** * Check that the given array is sorted. * * @param val Values. * @param dir Ordering direction. * @param strict Whether the order should be strict. * @param abort Whether to throw an exception if the check fails. * @return {@code true} if the array is sorted. * @throws NonMonotonicSequenceException if the array is not sorted * and {@code abort} is {@code true}. */ public static boolean checkOrder(double[] val, OrderDirection dir, boolean strict, boolean abort) throws NonMonotonicSequenceException { double previous = val[0]; final int max = val.length; int index; ITEM: for (index = 1; index < max; index++) { switch (dir) { case INCREASING: if (strict) { if (val[index] <= previous) { break ITEM; } } else { if (val[index] < previous) { break ITEM; } } break; case DECREASING: if (strict) { if (val[index] >= previous) { break ITEM; } } else { if (val[index] > previous) { break ITEM; } } break; default: // Should never happen. throw new MathInternalError(); } previous = val[index]; } if (index == max) { // Loop completed. return true; } // Loop early exit means wrong ordering. if (abort) { throw new NonMonotonicSequenceException(val[index], previous, index, dir, strict); } else { return false; } } /** * Check that the given array is sorted. * * @param val Values. * @param dir Ordering direction. * @param strict Whether the order should be strict. * @throws NonMonotonicSequenceException if the array is not sorted. * @since 2.2 */ public static void checkOrder(double[] val, OrderDirection dir, boolean strict) throws NonMonotonicSequenceException { checkOrder(val, dir, strict, true); } /** * Check that the given array is sorted in strictly increasing order. * * @param val Values. * @throws NonMonotonicSequenceException if the array is not sorted. * @since 2.2 */ public static void checkOrder(double[] val) throws NonMonotonicSequenceException { checkOrder(val, OrderDirection.INCREASING, true); } /** * Throws DimensionMismatchException if the input array is not rectangular. * * @param in array to be tested * @throws NullArgumentException if input array is null * @throws DimensionMismatchException if input array is not rectangular * @since 3.1 */ public static void checkRectangular(final long[][] in) throws NullArgumentException, DimensionMismatchException { MathUtils.checkNotNull(in); for (int i = 1; i < in.length; i++) { if (in[i].length != in[0].length) { throw new DimensionMismatchException( LocalizedFormats.DIFFERENT_ROWS_LENGTHS, in[i].length, in[0].length); } } } /** * Check that all entries of the input array are strictly positive. * * @param in Array to be tested * @throws NotStrictlyPositiveException if any entries of the array are not * strictly positive. * @since 3.1 */ public static void checkPositive(final double[] in) throws NotStrictlyPositiveException { for (int i = 0; i < in.length; i++) { if (in[i] <= 0) { throw new NotStrictlyPositiveException(in[i]); } } } /** * Check that no entry of the input array is {@code NaN}. * * @param in Array to be tested. * @throws NotANumberException if an entry is {@code NaN}. * @since 3.4 */ public static void checkNotNaN(final double[] in) throws NotANumberException { for(int i = 0; i < in.length; i++) { if (Double.isNaN(in[i])) { throw new NotANumberException(); } } } /** * Check that all entries of the input array are >= 0. * * @param in Array to be tested * @throws NotPositiveException if any array entries are less than 0. * @since 3.1 */ public static void checkNonNegative(final long[] in) throws NotPositiveException { for (int i = 0; i < in.length; i++) { if (in[i] < 0) { throw new NotPositiveException(in[i]); } } } /** * Check all entries of the input array are >= 0. * * @param in Array to be tested * @throws NotPositiveException if any array entries are less than 0. * @since 3.1 */ public static void checkNonNegative(final long[][] in) throws NotPositiveException { for (int i = 0; i < in.length; i ++) { for (int j = 0; j < in[i].length; j++) { if (in[i][j] < 0) { throw new NotPositiveException(in[i][j]); } } } } /** * Returns the Cartesian norm (2-norm), handling both overflow and underflow. * Translation of the minpack enorm subroutine. * * The redistribution policy for MINPACK is available * here, for * convenience, it is reproduced below.

* * * * *
* Minpack Copyright Notice (1999) University of Chicago. * All rights reserved *
* Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: *
    *
  1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer.
  2. *
  3. Redistributions in binary form must reproduce the above * copyright notice, this list of conditions and the following * disclaimer in the documentation and/or other materials provided * with the distribution.
  4. *
  5. The end-user documentation included with the redistribution, if any, * must include the following acknowledgment: * {@code This product includes software developed by the University of * Chicago, as Operator of Argonne National Laboratory.} * Alternately, this acknowledgment may appear in the software itself, * if and wherever such third-party acknowledgments normally appear.
  6. *
  7. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL * BE CORRECTED.
  8. *
  9. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE * POSSIBILITY OF SUCH LOSS OR DAMAGES.
  10. *
    * * @param v Vector of doubles. * @return the 2-norm of the vector. * @since 2.2 */ public static double safeNorm(double[] v) { double rdwarf = 3.834e-20; double rgiant = 1.304e+19; double s1 = 0; double s2 = 0; double s3 = 0; double x1max = 0; double x3max = 0; double floatn = v.length; double agiant = rgiant / floatn; for (int i = 0; i < v.length; i++) { double xabs = FastMath.abs(v[i]); if (xabs < rdwarf || xabs > agiant) { if (xabs > rdwarf) { if (xabs > x1max) { double r = x1max / xabs; s1= 1 + s1 * r * r; x1max = xabs; } else { double r = xabs / x1max; s1 += r * r; } } else { if (xabs > x3max) { double r = x3max / xabs; s3= 1 + s3 * r * r; x3max = xabs; } else { if (xabs != 0) { double r = xabs / x3max; s3 += r * r; } } } } else { s2 += xabs * xabs; } } double norm; if (s1 != 0) { norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max); } else { if (s2 == 0) { norm = x3max * Math.sqrt(s3); } else { if (s2 >= x3max) { norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3))); } else { norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3))); } } } return norm; } /** * A helper data structure holding a double and an integer value. */ private static class PairDoubleInteger { /** Key */ private final double key; /** Value */ private final int value; /** * @param key Key. * @param value Value. */ PairDoubleInteger(double key, int value) { this.key = key; this.value = value; } /** @return the key. */ public double getKey() { return key; } /** @return the value. */ public int getValue() { return value; } } /** * Sort an array in ascending order in place and perform the same reordering * of entries on other arrays. For example, if * {@code x = [3, 1, 2], y = [1, 2, 3]} and {@code z = [0, 5, 7]}, then * {@code sortInPlace(x, y, z)} will update {@code x} to {@code [1, 2, 3]}, * {@code y} to {@code [2, 3, 1]} and {@code z} to {@code [5, 7, 0]}. * * @param x Array to be sorted and used as a pattern for permutation * of the other arrays. * @param yList Set of arrays whose permutations of entries will follow * those performed on {@code x}. * @throws DimensionMismatchException if any {@code y} is not the same * size as {@code x}. * @throws NullArgumentException if {@code x} or any {@code y} is null. * @since 3.0 */ public static void sortInPlace(double[] x, double[] ... yList) throws DimensionMismatchException, NullArgumentException { sortInPlace(x, OrderDirection.INCREASING, yList); } /** * Sort an array in place and perform the same reordering of entries on * other arrays. This method works the same as the other * {@link #sortInPlace(double[], double[][]) sortInPlace} method, but * allows the order of the sort to be provided in the {@code dir} * parameter. * * @param x Array to be sorted and used as a pattern for permutation * of the other arrays. * @param dir Order direction. * @param yList Set of arrays whose permutations of entries will follow * those performed on {@code x}. * @throws DimensionMismatchException if any {@code y} is not the same * size as {@code x}. * @throws NullArgumentException if {@code x} or any {@code y} is null * @since 3.0 */ public static void sortInPlace(double[] x, final OrderDirection dir, double[] ... yList) throws NullArgumentException, DimensionMismatchException { // Consistency checks. if (x == null) { throw new NullArgumentException(); } final int yListLen = yList.length; final int len = x.length; for (int j = 0; j < yListLen; j++) { final double[] y = yList[j]; if (y == null) { throw new NullArgumentException(); } if (y.length != len) { throw new DimensionMismatchException(y.length, len); } } // Associate each abscissa "x[i]" with its index "i". final List list = new ArrayList(len); for (int i = 0; i < len; i++) { list.add(new PairDoubleInteger(x[i], i)); } // Create comparators for increasing and decreasing orders. final Comparator comp = dir == MathArrays.OrderDirection.INCREASING ? new Comparator() { /** {@inheritDoc} */ public int compare(PairDoubleInteger o1, PairDoubleInteger o2) { return Double.compare(o1.getKey(), o2.getKey()); } } : new Comparator() { /** {@inheritDoc} */ public int compare(PairDoubleInteger o1, PairDoubleInteger o2) { return Double.compare(o2.getKey(), o1.getKey()); } }; // Sort. Collections.sort(list, comp); // Modify the original array so that its elements are in // the prescribed order. // Retrieve indices of original locations. final int[] indices = new int[len]; for (int i = 0; i < len; i++) { final PairDoubleInteger e = list.get(i); x[i] = e.getKey(); indices[i] = e.getValue(); } // In each of the associated arrays, move the // elements to their new location. for (int j = 0; j < yListLen; j++) { // Input array will be modified in place. final double[] yInPlace = yList[j]; final double[] yOrig = yInPlace.clone(); for (int i = 0; i < len; i++) { yInPlace[i] = yOrig[indices[i]]; } } } /** * Creates a copy of the {@code source} array. * * @param source Array to be copied. * @return the copied array. */ public static int[] copyOf(int[] source) { return copyOf(source, source.length); } /** * Creates a copy of the {@code source} array. * * @param source Array to be copied. * @return the copied array. */ public static double[] copyOf(double[] source) { return copyOf(source, source.length); } /** * Creates a copy of the {@code source} array. * * @param source Array to be copied. * @param len Number of entries to copy. If smaller then the source * length, the copy will be truncated, if larger it will padded with * zeroes. * @return the copied array. */ public static int[] copyOf(int[] source, int len) { final int[] output = new int[len]; System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length)); return output; } /** * Creates a copy of the {@code source} array. * * @param source Array to be copied. * @param len Number of entries to copy. If smaller then the source * length, the copy will be truncated, if larger it will padded with * zeroes. * @return the copied array. */ public static double[] copyOf(double[] source, int len) { final double[] output = new double[len]; System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length)); return output; } /** * Creates a copy of the {@code source} array. * * @param source Array to be copied. * @param from Initial index of the range to be copied, inclusive. * @param to Final index of the range to be copied, exclusive. (This index may lie outside the array.) * @return the copied array. */ public static double[] copyOfRange(double[] source, int from, int to) { final int len = to - from; final double[] output = new double[len]; System.arraycopy(source, from, output, 0, FastMath.min(len, source.length - from)); return output; } /** * Compute a linear combination accurately. * This method computes the sum of the products * ai bi to high accuracy. * It does so by using specific multiplication and addition algorithms to * preserve accuracy and reduce cancellation effects. *
    * It is based on the 2005 paper * * Accurate Sum and Dot Product by Takeshi Ogita, Siegfried M. Rump, * and Shin'ichi Oishi published in SIAM J. Sci. Comput. * * @param a Factors. * @param b Factors. * @return Σi ai bi. * @throws DimensionMismatchException if arrays dimensions don't match */ public static double linearCombination(final double[] a, final double[] b) throws DimensionMismatchException { checkEqualLength(a, b); final int len = a.length; if (len == 1) { // Revert to scalar multiplication. return a[0] * b[0]; } final double[] prodHigh = new double[len]; double prodLowSum = 0; for (int i = 0; i < len; i++) { final double ai = a[i]; final double aHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(ai) & ((-1L) << 27)); final double aLow = ai - aHigh; final double bi = b[i]; final double bHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(bi) & ((-1L) << 27)); final double bLow = bi - bHigh; prodHigh[i] = ai * bi; final double prodLow = aLow * bLow - (((prodHigh[i] - aHigh * bHigh) - aLow * bHigh) - aHigh * bLow); prodLowSum += prodLow; } final double prodHighCur = prodHigh[0]; double prodHighNext = prodHigh[1]; double sHighPrev = prodHighCur + prodHighNext; double sPrime = sHighPrev - prodHighNext; double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime); final int lenMinusOne = len - 1; for (int i = 1; i < lenMinusOne; i++) { prodHighNext = prodHigh[i + 1]; final double sHighCur = sHighPrev + prodHighNext; sPrime = sHighCur - prodHighNext; sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime); sHighPrev = sHighCur; } double result = sHighPrev + (prodLowSum + sLowSum); if (Double.isNaN(result)) { // either we have split infinite numbers or some coefficients were NaNs, // just rely on the naive implementation and let IEEE754 handle this result = 0; for (int i = 0; i < len; ++i) { result += a[i] * b[i]; } } return result; } /** * Compute a linear combination accurately. *

    * This method computes a1×b1 + * a2×b2 to high accuracy. It does * so by using specific multiplication and addition algorithms to * preserve accuracy and reduce cancellation effects. It is based * on the 2005 paper * Accurate Sum and Dot Product by Takeshi Ogita, * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput. *

    * @param a1 first factor of the first term * @param b1 second factor of the first term * @param a2 first factor of the second term * @param b2 second factor of the second term * @return a1×b1 + * a2×b2 * @see #linearCombination(double, double, double, double, double, double) * @see #linearCombination(double, double, double, double, double, double, double, double) */ public static double linearCombination(final double a1, final double b1, final double a2, final double b2) { // the code below is split in many additions/subtractions that may // appear redundant. However, they should NOT be simplified, as they // use IEEE754 floating point arithmetic rounding properties. // The variable naming conventions are that xyzHigh contains the most significant // bits of xyz and xyzLow contains its least significant bits. So theoretically // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot // be represented in only one double precision number so we preserve two numbers // to hold it as long as we can, combining the high and low order bits together // only at the end, after cancellation may have occurred on high order bits // split a1 and b1 as one 26 bits number and one 27 bits number final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27)); final double a1Low = a1 - a1High; final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27)); final double b1Low = b1 - b1High; // accurate multiplication a1 * b1 final double prod1High = a1 * b1; final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low); // split a2 and b2 as one 26 bits number and one 27 bits number final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27)); final double a2Low = a2 - a2High; final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27)); final double b2Low = b2 - b2High; // accurate multiplication a2 * b2 final double prod2High = a2 * b2; final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low); // accurate addition a1 * b1 + a2 * b2 final double s12High = prod1High + prod2High; final double s12Prime = s12High - prod2High; final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime); // final rounding, s12 may have suffered many cancellations, we try // to recover some bits from the extra words we have saved up to now double result = s12High + (prod1Low + prod2Low + s12Low); if (Double.isNaN(result)) { // either we have split infinite numbers or some coefficients were NaNs, // just rely on the naive implementation and let IEEE754 handle this result = a1 * b1 + a2 * b2; } return result; } /** * Compute a linear combination accurately. *

    * This method computes a1×b1 + * a2×b2 + a3×b3 * to high accuracy. It does so by using specific multiplication and * addition algorithms to preserve accuracy and reduce cancellation effects. * It is based on the 2005 paper * Accurate Sum and Dot Product by Takeshi Ogita, * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput. *

    * @param a1 first factor of the first term * @param b1 second factor of the first term * @param a2 first factor of the second term * @param b2 second factor of the second term * @param a3 first factor of the third term * @param b3 second factor of the third term * @return a1×b1 + * a2×b2 + a3×b3 * @see #linearCombination(double, double, double, double) * @see #linearCombination(double, double, double, double, double, double, double, double) */ public static double linearCombination(final double a1, final double b1, final double a2, final double b2, final double a3, final double b3) { // the code below is split in many additions/subtractions that may // appear redundant. However, they should NOT be simplified, as they // do use IEEE754 floating point arithmetic rounding properties. // The variables naming conventions are that xyzHigh contains the most significant // bits of xyz and xyzLow contains its least significant bits. So theoretically // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot // be represented in only one double precision number so we preserve two numbers // to hold it as long as we can, combining the high and low order bits together // only at the end, after cancellation may have occurred on high order bits // split a1 and b1 as one 26 bits number and one 27 bits number final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27)); final double a1Low = a1 - a1High; final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27)); final double b1Low = b1 - b1High; // accurate multiplication a1 * b1 final double prod1High = a1 * b1; final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low); // split a2 and b2 as one 26 bits number and one 27 bits number final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27)); final double a2Low = a2 - a2High; final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27)); final double b2Low = b2 - b2High; // accurate multiplication a2 * b2 final double prod2High = a2 * b2; final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low); // split a3 and b3 as one 26 bits number and one 27 bits number final double a3High = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27)); final double a3Low = a3 - a3High; final double b3High = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27)); final double b3Low = b3 - b3High; // accurate multiplication a3 * b3 final double prod3High = a3 * b3; final double prod3Low = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low); // accurate addition a1 * b1 + a2 * b2 final double s12High = prod1High + prod2High; final double s12Prime = s12High - prod2High; final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime); // accurate addition a1 * b1 + a2 * b2 + a3 * b3 final double s123High = s12High + prod3High; final double s123Prime = s123High - prod3High; final double s123Low = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime); // final rounding, s123 may have suffered many cancellations, we try // to recover some bits from the extra words we have saved up to now double result = s123High + (prod1Low + prod2Low + prod3Low + s12Low + s123Low); if (Double.isNaN(result)) { // either we have split infinite numbers or some coefficients were NaNs, // just rely on the naive implementation and let IEEE754 handle this result = a1 * b1 + a2 * b2 + a3 * b3; } return result; } /** * Compute a linear combination accurately. *

    * This method computes a1×b1 + * a2×b2 + a3×b3 + * a4×b4 * to high accuracy. It does so by using specific multiplication and * addition algorithms to preserve accuracy and reduce cancellation effects. * It is based on the 2005 paper * Accurate Sum and Dot Product by Takeshi Ogita, * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput. *

    * @param a1 first factor of the first term * @param b1 second factor of the first term * @param a2 first factor of the second term * @param b2 second factor of the second term * @param a3 first factor of the third term * @param b3 second factor of the third term * @param a4 first factor of the third term * @param b4 second factor of the third term * @return a1×b1 + * a2×b2 + a3×b3 + * a4×b4 * @see #linearCombination(double, double, double, double) * @see #linearCombination(double, double, double, double, double, double) */ public static double linearCombination(final double a1, final double b1, final double a2, final double b2, final double a3, final double b3, final double a4, final double b4) { // the code below is split in many additions/subtractions that may // appear redundant. However, they should NOT be simplified, as they // do use IEEE754 floating point arithmetic rounding properties. // The variables naming conventions are that xyzHigh contains the most significant // bits of xyz and xyzLow contains its least significant bits. So theoretically // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot // be represented in only one double precision number so we preserve two numbers // to hold it as long as we can, combining the high and low order bits together // only at the end, after cancellation may have occurred on high order bits // split a1 and b1 as one 26 bits number and one 27 bits number final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27)); final double a1Low = a1 - a1High; final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27)); final double b1Low = b1 - b1High; // accurate multiplication a1 * b1 final double prod1High = a1 * b1; final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low); // split a2 and b2 as one 26 bits number and one 27 bits number final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27)); final double a2Low = a2 - a2High; final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27)); final double b2Low = b2 - b2High; // accurate multiplication a2 * b2 final double prod2High = a2 * b2; final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low); // split a3 and b3 as one 26 bits number and one 27 bits number final double a3High = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27)); final double a3Low = a3 - a3High; final double b3High = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27)); final double b3Low = b3 - b3High; // accurate multiplication a3 * b3 final double prod3High = a3 * b3; final double prod3Low = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low); // split a4 and b4 as one 26 bits number and one 27 bits number final double a4High = Double.longBitsToDouble(Double.doubleToRawLongBits(a4) & ((-1L) << 27)); final double a4Low = a4 - a4High; final double b4High = Double.longBitsToDouble(Double.doubleToRawLongBits(b4) & ((-1L) << 27)); final double b4Low = b4 - b4High; // accurate multiplication a4 * b4 final double prod4High = a4 * b4; final double prod4Low = a4Low * b4Low - (((prod4High - a4High * b4High) - a4Low * b4High) - a4High * b4Low); // accurate addition a1 * b1 + a2 * b2 final double s12High = prod1High + prod2High; final double s12Prime = s12High - prod2High; final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime); // accurate addition a1 * b1 + a2 * b2 + a3 * b3 final double s123High = s12High + prod3High; final double s123Prime = s123High - prod3High; final double s123Low = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime); // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4 final double s1234High = s123High + prod4High; final double s1234Prime = s1234High - prod4High; final double s1234Low = (prod4High - (s1234High - s1234Prime)) + (s123High - s1234Prime); // final rounding, s1234 may have suffered many cancellations, we try // to recover some bits from the extra words we have saved up to now double result = s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + s12Low + s123Low + s1234Low); if (Double.isNaN(result)) { // either we have split infinite numbers or some coefficients were NaNs, // just rely on the naive implementation and let IEEE754 handle this result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4; } return result; } /** * Returns true iff both arguments are null or have same dimensions and all * their elements are equal as defined by * {@link Precision#equals(float,float)}. * * @param x first array * @param y second array * @return true if the values are both null or have same dimension * and equal elements. */ public static boolean equals(float[] x, float[] y) { if ((x == null) || (y == null)) { return !((x == null) ^ (y == null)); } if (x.length != y.length) { return false; } for (int i = 0; i < x.length; ++i) { if (!Precision.equals(x[i], y[i])) { return false; } } return true; } /** * Returns true iff both arguments are null or have same dimensions and all * their elements are equal as defined by * {@link Precision#equalsIncludingNaN(double,double) this method}. * * @param x first array * @param y second array * @return true if the values are both null or have same dimension and * equal elements * @since 2.2 */ public static boolean equalsIncludingNaN(float[] x, float[] y) { if ((x == null) || (y == null)) { return !((x == null) ^ (y == null)); } if (x.length != y.length) { return false; } for (int i = 0; i < x.length; ++i) { if (!Precision.equalsIncludingNaN(x[i], y[i])) { return false; } } return true; } /** * Returns {@code true} iff both arguments are {@code null} or have same * dimensions and all their elements are equal as defined by * {@link Precision#equals(double,double)}. * * @param x First array. * @param y Second array. * @return {@code true} if the values are both {@code null} or have same * dimension and equal elements. */ public static boolean equals(double[] x, double[] y) { if ((x == null) || (y == null)) { return !((x == null) ^ (y == null)); } if (x.length != y.length) { return false; } for (int i = 0; i < x.length; ++i) { if (!Precision.equals(x[i], y[i])) { return false; } } return true; } /** * Returns {@code true} iff both arguments are {@code null} or have same * dimensions and all their elements are equal as defined by * {@link Precision#equalsIncludingNaN(double,double) this method}. * * @param x First array. * @param y Second array. * @return {@code true} if the values are both {@code null} or have same * dimension and equal elements. * @since 2.2 */ public static boolean equalsIncludingNaN(double[] x, double[] y) { if ((x == null) || (y == null)) { return !((x == null) ^ (y == null)); } if (x.length != y.length) { return false; } for (int i = 0; i < x.length; ++i) { if (!Precision.equalsIncludingNaN(x[i], y[i])) { return false; } } return true; } /** * Normalizes an array to make it sum to a specified value. * Returns the result of the transformation *
         *    x |-> x * normalizedSum / sum
         * 
    * applied to each non-NaN element x of the input array, where sum is the * sum of the non-NaN entries in the input array. *

    * Throws IllegalArgumentException if {@code normalizedSum} is infinite * or NaN and ArithmeticException if the input array contains any infinite elements * or sums to 0. *

    * Ignores (i.e., copies unchanged to the output array) NaNs in the input array. * * @param values Input array to be normalized * @param normalizedSum Target sum for the normalized array * @return the normalized array. * @throws MathArithmeticException if the input array contains infinite * elements or sums to zero. * @throws MathIllegalArgumentException if the target sum is infinite or {@code NaN}. * @since 2.1 */ public static double[] normalizeArray(double[] values, double normalizedSum) throws MathIllegalArgumentException, MathArithmeticException { if (Double.isInfinite(normalizedSum)) { throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_INFINITE); } if (Double.isNaN(normalizedSum)) { throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_NAN); } double sum = 0d; final int len = values.length; double[] out = new double[len]; for (int i = 0; i < len; i++) { if (Double.isInfinite(values[i])) { throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i); } if (!Double.isNaN(values[i])) { sum += values[i]; } } if (sum == 0) { throw new MathArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO); } for (int i = 0; i < len; i++) { if (Double.isNaN(values[i])) { out[i] = Double.NaN; } else { out[i] = values[i] * normalizedSum / sum; } } return out; } /** Build an array of elements. *

    * Arrays are filled with field.getZero() * * @param the type of the field elements * @param field field to which array elements belong * @param length of the array * @return a new array * @since 3.2 */ public static T[] buildArray(final Field field, final int length) { @SuppressWarnings("unchecked") // OK because field must be correct class T[] array = (T[]) Array.newInstance(field.getRuntimeClass(), length); Arrays.fill(array, field.getZero()); return array; } /** Build a double dimension array of elements. *

    * Arrays are filled with field.getZero() * * @param the type of the field elements * @param field field to which array elements belong * @param rows number of rows in the array * @param columns number of columns (may be negative to build partial * arrays in the same way new Field[rows][] works) * @return a new array * @since 3.2 */ @SuppressWarnings("unchecked") public static T[][] buildArray(final Field field, final int rows, final int columns) { final T[][] array; if (columns < 0) { T[] dummyRow = buildArray(field, 0); array = (T[][]) Array.newInstance(dummyRow.getClass(), rows); } else { array = (T[][]) Array.newInstance(field.getRuntimeClass(), new int[] { rows, columns }); for (int i = 0; i < rows; ++i) { Arrays.fill(array[i], field.getZero()); } } return array; } /** * Calculates the * convolution between two sequences. *

    * The solution is obtained via straightforward computation of the * convolution sum (and not via FFT). Whenever the computation needs * an element that would be located at an index outside the input arrays, * the value is assumed to be zero. * * @param x First sequence. * Typically, this sequence will represent an input signal to a system. * @param h Second sequence. * Typically, this sequence will represent the impulse response of the system. * @return the convolution of {@code x} and {@code h}. * This array's length will be {@code x.length + h.length - 1}. * @throws NullArgumentException if either {@code x} or {@code h} is {@code null}. * @throws NoDataException if either {@code x} or {@code h} is empty. * * @since 3.3 */ public static double[] convolve(double[] x, double[] h) throws NullArgumentException, NoDataException { MathUtils.checkNotNull(x); MathUtils.checkNotNull(h); final int xLen = x.length; final int hLen = h.length; if (xLen == 0 || hLen == 0) { throw new NoDataException(); } // initialize the output array final int totalLength = xLen + hLen - 1; final double[] y = new double[totalLength]; // straightforward implementation of the convolution sum for (int n = 0; n < totalLength; n++) { double yn = 0; int k = FastMath.max(0, n + 1 - xLen); int j = n - k; while (k < hLen && j >= 0) { yn += x[j--] * h[k++]; } y[n] = yn; } return y; } /** * Specification for indicating that some operation applies * before or after a given index. */ public enum Position { /** Designates the beginning of the array (near index 0). */ HEAD, /** Designates the end of the array. */ TAIL } /** * Shuffle the entries of the given array. * The {@code start} and {@code pos} parameters select which portion * of the array is randomized and which is left untouched. * * @see #shuffle(int[],int,Position,RandomGenerator) * * @param list Array whose entries will be shuffled (in-place). * @param start Index at which shuffling begins. * @param pos Shuffling is performed for index positions between * {@code start} and either the end (if {@link Position#TAIL}) * or the beginning (if {@link Position#HEAD}) of the array. */ public static void shuffle(int[] list, int start, Position pos) { shuffle(list, start, pos, new Well19937c()); } /** * Shuffle the entries of the given array, using the * * Fisher–Yates algorithm. * The {@code start} and {@code pos} parameters select which portion * of the array is randomized and which is left untouched. * * @param list Array whose entries will be shuffled (in-place). * @param start Index at which shuffling begins. * @param pos Shuffling is performed for index positions between * {@code start} and either the end (if {@link Position#TAIL}) * or the beginning (if {@link Position#HEAD}) of the array. * @param rng Random number generator. */ public static void shuffle(int[] list, int start, Position pos, RandomGenerator rng) { switch (pos) { case TAIL: { for (int i = list.length - 1; i >= start; i--) { final int target; if (i == start) { target = start; } else { // NumberIsTooLargeException cannot occur. target = new UniformIntegerDistribution(rng, start, i).sample(); } final int temp = list[target]; list[target] = list[i]; list[i] = temp; } } break; case HEAD: { for (int i = 0; i <= start; i++) { final int target; if (i == start) { target = start; } else { // NumberIsTooLargeException cannot occur. target = new UniformIntegerDistribution(rng, i, start).sample(); } final int temp = list[target]; list[target] = list[i]; list[i] = temp; } } break; default: throw new MathInternalError(); // Should never happen. } } /** * Shuffle the entries of the given array. * * @see #shuffle(int[],int,Position,RandomGenerator) * * @param list Array whose entries will be shuffled (in-place). * @param rng Random number generator. */ public static void shuffle(int[] list, RandomGenerator rng) { shuffle(list, 0, Position.TAIL, rng); } /** * Shuffle the entries of the given array. * * @see #shuffle(int[],int,Position,RandomGenerator) * * @param list Array whose entries will be shuffled (in-place). */ public static void shuffle(int[] list) { shuffle(list, new Well19937c()); } /** * Returns an array representing the natural number {@code n}. * * @param n Natural number. * @return an array whose entries are the numbers 0, 1, ..., {@code n}-1. * If {@code n == 0}, the returned array is empty. */ public static int[] natural(int n) { return sequence(n, 0, 1); } /** * Returns an array of {@code size} integers starting at {@code start}, * skipping {@code stride} numbers. * * @param size Natural number. * @param start Natural number. * @param stride Natural number. * @return an array whose entries are the numbers * {@code start, start + stride, ..., start + (size - 1) * stride}. * If {@code size == 0}, the returned array is empty. * * @since 3.4 */ public static int[] sequence(int size, int start, int stride) { final int[] a = new int[size]; for (int i = 0; i < size; i++) { a[i] = start + i * stride; } return a; } /** * This method is used * to verify that the input parameters designate a subarray of positive length. *

    *

      *
    • returns true iff the parameters designate a subarray of * positive length
    • *
    • throws MathIllegalArgumentException if the array is null or * or the indices are invalid
    • *
    • returns false
    • if the array is non-null, but * length is 0. *

    * * @param values the input array * @param begin index of the first array element to include * @param length the number of elements to include * @return true if the parameters are valid and designate a subarray of positive length * @throws MathIllegalArgumentException if the indices are invalid or the array is null * @since 3.3 */ public static boolean verifyValues(final double[] values, final int begin, final int length) throws MathIllegalArgumentException { return verifyValues(values, begin, length, false); } /** * This method is used * to verify that the input parameters designate a subarray of positive length. *

    *

      *
    • returns true iff the parameters designate a subarray of * non-negative length
    • *
    • throws IllegalArgumentException if the array is null or * or the indices are invalid
    • *
    • returns false
    • if the array is non-null, but * length is 0 unless allowEmpty is true *

    * * @param values the input array * @param begin index of the first array element to include * @param length the number of elements to include * @param allowEmpty if true then zero length arrays are allowed * @return true if the parameters are valid * @throws MathIllegalArgumentException if the indices are invalid or the array is null * @since 3.3 */ public static boolean verifyValues(final double[] values, final int begin, final int length, final boolean allowEmpty) throws MathIllegalArgumentException { if (values == null) { throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY); } if (begin < 0) { throw new NotPositiveException(LocalizedFormats.START_POSITION, Integer.valueOf(begin)); } if (length < 0) { throw new NotPositiveException(LocalizedFormats.LENGTH, Integer.valueOf(length)); } if (begin + length > values.length) { throw new NumberIsTooLargeException(LocalizedFormats.SUBARRAY_ENDS_AFTER_ARRAY_END, Integer.valueOf(begin + length), Integer.valueOf(values.length), true); } if (length == 0 && !allowEmpty) { return false; } return true; } /** * This method is used * to verify that the begin and length parameters designate a subarray of positive length * and the weights are all non-negative, non-NaN, finite, and not all zero. *

    *

      *
    • returns true iff the parameters designate a subarray of * positive length and the weights array contains legitimate values.
    • *
    • throws IllegalArgumentException if any of the following are true: *
      • the values array is null
      • *
      • the weights array is null
      • *
      • the weights array does not have the same length as the values array
      • *
      • the weights array contains one or more infinite values
      • *
      • the weights array contains one or more NaN values
      • *
      • the weights array contains negative values
      • *
      • the start and length arguments do not determine a valid array
      *
    • *
    • returns false
    • if the array is non-null, but * length is 0. *

    * * @param values the input array * @param weights the weights array * @param begin index of the first array element to include * @param length the number of elements to include * @return true if the parameters are valid and designate a subarray of positive length * @throws MathIllegalArgumentException if the indices are invalid or the array is null * @since 3.3 */ public static boolean verifyValues( final double[] values, final double[] weights, final int begin, final int length) throws MathIllegalArgumentException { return verifyValues(values, weights, begin, length, false); } /** * This method is used * to verify that the begin and length parameters designate a subarray of positive length * and the weights are all non-negative, non-NaN, finite, and not all zero. *

    *

      *
    • returns true iff the parameters designate a subarray of * non-negative length and the weights array contains legitimate values.
    • *
    • throws MathIllegalArgumentException if any of the following are true: *
      • the values array is null
      • *
      • the weights array is null
      • *
      • the weights array does not have the same length as the values array
      • *
      • the weights array contains one or more infinite values
      • *
      • the weights array contains one or more NaN values
      • *
      • the weights array contains negative values
      • *
      • the start and length arguments do not determine a valid array
      *
    • *
    • returns false
    • if the array is non-null, but * length is 0 unless allowEmpty is true. *

    * * @param values the input array. * @param weights the weights array. * @param begin index of the first array element to include. * @param length the number of elements to include. * @param allowEmpty if {@code true} than allow zero length arrays to pass. * @return {@code true} if the parameters are valid. * @throws NullArgumentException if either of the arrays are null * @throws MathIllegalArgumentException if the array indices are not valid, * the weights array contains NaN, infinite or negative elements, or there * are no positive weights. * @since 3.3 */ public static boolean verifyValues(final double[] values, final double[] weights, final int begin, final int length, final boolean allowEmpty) throws MathIllegalArgumentException { if (weights == null || values == null) { throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY); } checkEqualLength(weights, values); boolean containsPositiveWeight = false; for (int i = begin; i < begin + length; i++) { final double weight = weights[i]; if (Double.isNaN(weight)) { throw new MathIllegalArgumentException(LocalizedFormats.NAN_ELEMENT_AT_INDEX, Integer.valueOf(i)); } if (Double.isInfinite(weight)) { throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT, Double.valueOf(weight), Integer.valueOf(i)); } if (weight < 0) { throw new MathIllegalArgumentException(LocalizedFormats.NEGATIVE_ELEMENT_AT_INDEX, Integer.valueOf(i), Double.valueOf(weight)); } if (!containsPositiveWeight && weight > 0.0) { containsPositiveWeight = true; } } if (!containsPositiveWeight) { throw new MathIllegalArgumentException(LocalizedFormats.WEIGHT_AT_LEAST_ONE_NON_ZERO); } return verifyValues(values, begin, length, allowEmpty); } /** * Concatenates a sequence of arrays. The return array consists of the * entries of the input arrays concatenated in the order they appear in * the argument list. Null arrays cause NullPointerExceptions; zero * length arrays are allowed (contributing nothing to the output array). * * @param x list of double[] arrays to concatenate * @return a new array consisting of the entries of the argument arrays * @throws NullPointerException if any of the arrays are null * @since 3.6 */ public static double[] concatenate(double[] ...x) { int combinedLength = 0; for (double[] a : x) { combinedLength += a.length; } int offset = 0; int curLength = 0; final double[] combined = new double[combinedLength]; for (int i = 0; i < x.length; i++) { curLength = x[i].length; System.arraycopy(x[i], 0, combined, offset, curLength); offset += curLength; } return combined; } /** * Returns an array consisting of the unique values in {@code data}. * The return array is sorted in descending order. Empty arrays * are allowed, but null arrays result in NullPointerException. * Infinities are allowed. NaN values are allowed with maximum * sort order - i.e., if there are NaN values in {@code data}, * {@code Double.NaN} will be the first element of the output array, * even if the array also contains {@code Double.POSITIVE_INFINITY}. * * @param data array to scan * @return descending list of values included in the input array * @throws NullPointerException if data is null * @since 3.6 */ public static double[] unique(double[] data) { TreeSet values = new TreeSet(); for (int i = 0; i < data.length; i++) { values.add(data[i]); } final int count = values.size(); final double[] out = new double[count]; Iterator iterator = values.iterator(); int i = 0; while (iterator.hasNext()) { out[count - ++i] = iterator.next(); } return out; } }