package jsci.maths;
import jsci.GlobalSettings;
/**
* Arrays are faster than object, so this class is here to take full
* advantage of arrays without encapsulation.
* All methods are safe, that is, they create copies of arrays whenever
* necessary
* This makes for slower methods, but the programmer can always
* define his own methods optimized for performance.
* @author Daniel Lemire
* @author Mark Hale
*/
public final class ArrayMath extends AbstractMath {
private ArrayMath() {}
public static boolean isSquare(double[][] a) {
for(int i=0; i1)) {
throw new IllegalArgumentException("Percentile must be between 0 and 1 : "+p);
}
final double[] ans = (isSorted ? v : sortMinToMax(v));
final int pos = (int) Math.floor(p*(ans.length-1));
final double dif = p*(ans.length-1) - Math.floor(p*(ans.length-1));
if(pos == (ans.length-1))
return(ans[ans.length-1]);
else
return(ans[pos]*(1.0-dif) + ans[pos+1]*dif);
}
/**
* Computes the median of an array.
* If the array contains an even number of elements,
* then the mean of the lower and upper medians is returned.
*/
public static double median(double[] v) {
if(v.length == 3) {
return median3(copy(v));
} else if(v.length == 5) {
return median5(copy(v));
} else if(v.length == 7) {
return median7(copy(v));
} else if(v.length % 2 == 1) { // odd
return quickSelect(copy(v), false);
} else { // even
double[] tmp = copy(v);
double lowerMedian = quickSelect(tmp, false);
double upperMedian = quickSelect(tmp, true);
return (lowerMedian+upperMedian)/2.0;
}
}
/*
* http://ndevilla.free.fr/median/median/index.html
* optmed.c
* @author N. Devillard
*/
private static void sortInPlace(final double[] v, final int i, final int j) {
if(v[i] > v[j]) {
final double tmp = v[i];
v[i] = v[j];
v[j] = tmp;
}
}
private static double median3(double[] v) {
sortInPlace(v,0,1);
sortInPlace(v,1,2);
sortInPlace(v,0,1);
return v[1];
}
private static double median5(double[] v) {
sortInPlace(v,0,1);
sortInPlace(v,3,4);
sortInPlace(v,0,3);
sortInPlace(v,1,4);
sortInPlace(v,1,2);
sortInPlace(v,2,3);
sortInPlace(v,1,2);
return v[2];
}
private static double median7(double[] v) {
sortInPlace(v,0,5);
sortInPlace(v,0,3);
sortInPlace(v,1,6);
sortInPlace(v,2,4);
sortInPlace(v,0,1);
sortInPlace(v,3,5);
sortInPlace(v,2,6);
sortInPlace(v,2,3);
sortInPlace(v,3,6);
sortInPlace(v,4,5);
sortInPlace(v,1,4);
sortInPlace(v,1,3);
sortInPlace(v,3,4);
return v[3];
}
/*
* http://ndevilla.free.fr/median/median/index.html
* quickselect.c
* @author N. Devillard
*/
private static double quickSelect(double[] v, boolean upperMedian) {
int low=0, high=v.length-1;
final int median = upperMedian ? high/2+1 : high/2;
int middle, ll, hh;
for (;;) {
if (high <= low) { /* One element only */
return v[median];
}
if (high == low + 1) { /* Two elements only */
if (v[low] > v[high])
swap(v, low, high);
return v[median];
}
/* Find median of low, middle and high items; swap into position low */
middle = (low + high) / 2;
if (v[middle] > v[high]) swap(v, middle, high);
if (v[low] > v[high]) swap(v, low, high);
if (v[middle] > v[low]) swap(v, middle, low);
/* Swap low item (now in position middle) into position (low+1) */
swap(v, middle, low+1);
/* Nibble from each end towards middle, swapping items when stuck */
ll = low + 1;
hh = high;
for (;;) {
do ll++; while (v[low] > v[ll]);
do hh--; while (v[hh] > v[low]);
if (hh < ll)
break;
swap(v, ll, hh);
}
/* Swap middle item (in position low) back into correct position */
swap(v, low, hh);
/* Re-set active partition */
if (hh <= median)
low = ll;
if (hh >= median)
high = hh - 1;
}
}
/**
* Checks if two arrays are equal within a tolerance.
*/
public static boolean equals(double[] a, double[] b) {
if(a.length!=b.length)
return(false);
for(int i=0;iGlobalSettings.ZERO_TOL)
return(false);
}
return(true);
}
/**
* Takes the absolute value of each component of an array.
*/
public static double[] abs(double[] v) {
double[] ans=new double[v.length];
for(int i=0;iv[i]) {
min=v[i];
}
}
return(min);
}
/**
* Returns the minimum of an array.
*/
public static double min(double[][] v) {
double min=min(v[0]);
for(int i=1;imin(v[i])) {
min=min(v[i]);
}
}
return(min);
}
/**
* Returns the componentwise modulus of
* an array of Complex numbers.
*/
public static double[] mod(Complex[] v) {
double[] ans=new double[v.length];
for(int k=0;k k1).
* @param k0 location of the first component
* @param k1 location of the last component
*/
public static double[] extract(int k0, int k1, double[] invect) {
if ((k0<0)||(k1<0)||(k0>invect.length-1)||(k1>invect.length-1)) {
throw new IllegalArgumentException("The parameters are incorrect : "+k0+", "+k1+", "+invect.length);
}
if(k1>k0) {
double[] ans=new double[k1-k0+1];
System.arraycopy(invect,k0,ans,0,k1-k0+1);
return(ans);
}
double[] ans=new double[-k1+k0+1];
for(int k=k1;k<=k0;k++) {
ans[k0-k]=invect[k];
}
return(ans);
}
/**
* Inverts an array from left to right.
*/
public static double [] invert (double[] v) {
double[] w=new double[v.length];
for(int k=0;kn0)||(pos<0)) {
throw new IllegalArgumentException("Array is to long for this : "+n0+", "+pos+", "+v.length);
}
double[] w=new double[n0];
System.arraycopy(v,0,w,pos,v.length);
return(w);
}
/**
* Adds to an array w, a times v where a is a scalar.
* Since v can be smaller then w, we must specified the position at which v will be added.
* @param a scalar.
* @param p position.
*/
public static double[] add(double[]w, double a, double[] v,int p) {
if(v.length>w.length) {
throw new IllegalArgumentException("Second array must be shorter or equal to the first one : "+w.length+", "+v.length);
}
double[] ans=copy(w);
for(int k=p;kb) {
step=-step;
}
for(int k=1;k1)) {
throw new IllegalArgumentException("Percentile must be between 0 and 1 : "+p);
}
final int[] ans = (isSorted ? v : sortMinToMax(v));
final int pos = (int) Math.floor(p*(ans.length-1));
final double dif = p*(ans.length-1)-Math.floor(p*(ans.length-1));
if(pos == (ans.length-1))
return(ans[ans.length-1]);
else
return(ans[pos]*(1.0-dif)+ans[pos+1]*dif);
}
/**
* Computes the median of an array.
* If the array contains an even number of elements,
* the lower median is returned.
*/
public static int median(int[] v) {
if(v.length == 3) {
return median3(copy(v));
} else if(v.length == 5) {
return median5(copy(v));
} else if(v.length == 7) {
return median7(copy(v));
} else {
return quickSelect(copy(v));
}
}
/*
* http://ndevilla.free.fr/median/median/index.html
* optmed.c
* @author N. Devillard
*/
private static void sortInPlace(final int[] v, final int i, final int j) {
if(v[i] > v[j]) {
final int tmp = v[i];
v[i] = v[j];
v[j] = tmp;
}
}
private static int median3(int[] v) {
sortInPlace(v,0,1);
sortInPlace(v,1,2);
sortInPlace(v,0,1);
return v[1];
}
private static int median5(int[] v) {
sortInPlace(v,0,1);
sortInPlace(v,3,4);
sortInPlace(v,0,3);
sortInPlace(v,1,4);
sortInPlace(v,1,2);
sortInPlace(v,2,3);
sortInPlace(v,1,2);
return v[2];
}
private static int median7(int[] v) {
sortInPlace(v,0,5);
sortInPlace(v,0,3);
sortInPlace(v,1,6);
sortInPlace(v,2,4);
sortInPlace(v,0,1);
sortInPlace(v,3,5);
sortInPlace(v,2,6);
sortInPlace(v,2,3);
sortInPlace(v,3,6);
sortInPlace(v,4,5);
sortInPlace(v,1,4);
sortInPlace(v,1,3);
sortInPlace(v,3,4);
return v[3];
}
/*
* http://ndevilla.free.fr/median/median/index.html
* quickselect.c
* @author N. Devillard
*/
private static int quickSelect(int[] v) {
int low=0, high=v.length-1;
final int median=high/2;
int middle, ll, hh;
for (;;) {
if (high <= low) { /* One element only */
return v[median];
}
if (high == low + 1) { /* Two elements only */
if (v[low] > v[high])
swap(v, low, high);
return v[median];
}
/* Find median of low, middle and high items; swap into position low */
middle = (low + high) / 2;
if (v[middle] > v[high]) swap(v, middle, high);
if (v[low] > v[high]) swap(v, low, high);
if (v[middle] > v[low]) swap(v, middle, low);
/* Swap low item (now in position middle) into position (low+1) */
swap(v, middle, low+1);
/* Nibble from each end towards middle, swapping items when stuck */
ll = low + 1;
hh = high;
for (;;) {
do ll++; while (v[low] > v[ll]);
do hh--; while (v[hh] > v[low]);
if (hh < ll)
break;
swap(v, ll, hh);
}
/* Swap middle item (in position low) back into correct position */
swap(v, low, hh);
/* Re-set active partition */
if (hh <= median)
low = ll;
if (hh >= median)
high = hh - 1;
}
}
/**
* Use java.util.Arrays instead.
* Check if two arrays are equal.
* @deprecated Use java.util.Arrays instead.
*/
public static boolean equals(int[] a, int[] b) {
if(a.length!=b.length) {
return(false);
}
for(int i=0;iv[i]) {
min=v[i];
}
}
return(min);
}
/**
* Returns the minimum of an array.
*/
public static int min(int[][] v) {
int min=min(v[0]);
for(int i=1;imin(v[i])) {
min=min(v[i]);
}
}
return(min);
}
/**
* Computes the L2 norm of an array (Euclidean norm or "length").
*/
public static double norm (int[] data) {
return(Math.sqrt(sumSquares(data)));
}
/**
* Sums the squares of all components;
* also called the energy of the array.
*/
public static int sumSquares (int[] data) {
int ans=0;
for(int k=0;k k1).
* @param k0 location of the first component.
* @param k1 location of the last component.
*/
public static int[] extract(int k0, int k1, int[] invect) {
if ((k0<0)||(k1<0)||(k0>invect.length-1)||(k1>invect.length-1)) {
throw new IllegalArgumentException("The parameters are incorrect : "+k0+", "+k1+", "+invect.length);
}
if(k1>k0) {
int[] ans=new int[k1-k0+1];
System.arraycopy(invect,k0,ans,0,k1-k0+1);
return(ans);
}
int[] ans=new int[-k1+k0+1];
for(int k=k1;k<=k0;k++) {
ans[k0-k]=invect[k];
}
return(ans);
}
/**
* Inverts an array from left to right.
*/
public static int [] invert (int[] v) {
int[] w=new int[v.length];
for(int k=0;kn0)||(pos<0)) {
throw new IllegalArgumentException("The array is too long for this : "+n0+", "+pos+", "+v.length);
}
int[] w=new int[n0];
System.arraycopy(v,0,w,pos,v.length);
return(w);
}
/**
* Adds to an array w, a times v where a is a scalar.
* Since v can be smaller then w, we must specified the position at which v will be added.
* @param a scalar.
* @param p position.
* @param w longer array.
* @param v shorter array.
* @exception IllegalArgumentException if the second array
* is not shorter than the first one.
*/
public static int[] add(int[]w, double a, int[] v,int p) {
if(v.length>w.length) {
throw new IllegalArgumentException("Second array must be shorter or equal to the first one : "+w.length+", "+v.length);
}
int[] ans=copy(w);
for(int k=p;kb) {
step=-step;
}
for(int k=1;k
Ads help maintain this website.