/* * 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 org.apache.commons.math3.Field; 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.NotANumberException; 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.util.LocalizedFormats; import org.apache.commons.math3.random.RandomGenerator; import org.apache.commons.math3.random.Well19937c; 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; /** * 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 redistribution policy for MINPACK is available here, for convenience, it is reproduced
* below.
*
* 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
*
* 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 Arrays are filled with field.getZero()
*
* @param 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.
*
*
*
*
*
*
* @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
*
* 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:
*
*
* 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.
*
*
* 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.
*
* new Field[rows][]
works)
* @return a new array
* @since 3.2
*/
@SuppressWarnings("unchecked")
public static
*
true
iff the parameters designate a subarray of
* positive lengthMathIllegalArgumentException
if the array is null or
* or the indices are invalidfalse
length
is 0.
*
*
true
iff the parameters designate a subarray of
* non-negative lengthIllegalArgumentException
if the array is null or
* or the indices are invalidfalse
length
is 0 unless allowEmpty
is true
* 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.
*
*
true
iff the parameters designate a subarray of
* positive length and the weights array contains legitimate values.IllegalArgumentException
if any of the following are true:
*
* false
length
is 0.
*
*
true
iff the parameters designate a subarray of
* non-negative length and the weights array contains legitimate values.MathIllegalArgumentException
if any of the following are true:
* false
length
is 0 unless allowEmpty
is true
.
*