/*
Copyright (C) 1999 CERN - European Organization for Nuclear Research.
Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
is hereby granted without fee, provided that the above copyright notice appear in all copies and
that both that copyright notice and this permission notice appear in supporting documentation.
CERN makes no representations about the suitability of this software for any purpose.
It is provided "as is" without expressed or implied warranty.
*/
package cern.colt.matrix.tfloat.algo;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.Future;
import cern.colt.matrix.AbstractFormatter;
import cern.colt.matrix.tfloat.FloatMatrix1D;
import cern.colt.matrix.tfloat.FloatMatrix2D;
import cern.colt.matrix.tfloat.FloatMatrix3D;
import cern.colt.matrix.tfloat.impl.DenseColumnFloatMatrix2D;
import cern.colt.matrix.tfloat.impl.DenseFloatMatrix1D;
import cern.colt.matrix.tfloat.impl.DenseFloatMatrix2D;
import cern.colt.matrix.tfloat.impl.SparseCCFloatMatrix2D;
import cern.colt.matrix.tfloat.impl.SparseFloatMatrix1D;
import cern.colt.matrix.tfloat.impl.SparseRCFloatMatrix2D;
import cern.jet.math.tfloat.FloatFunctions;
import edu.emory.utils.ConcurrencyUtils;
/**
* Tests matrices for linear algebraic properties (equality, tridiagonality,
* symmetry, singularity, etc).
*
* Except where explicitly indicated, all methods involving equality tests (
* ==) allow for numerical instability, to a degree specified upon
* instance construction and returned by method {@link #tolerance()}. The public
* static final variable DEFAULT represents a default Property object
* with a tolerance of 1.0E-5. The public static final variable
* ZERO represents a Property object with a tolerance of 0.0.
* The public static final variable SEVEN represents a Property object
* with a tolerance of 1.0E-7. As long as you are happy with these
* tolerances, there is no need to construct Property objects. Simply use idioms
* like Property.DEFAULT.equals(A,B),
* Property.ZERO.equals(A,B), Property.SEVEN.equals(A,B).
*
* To work with a different tolerance (e.g. 1.0E-2) use the constructor
* and/or method {@link #setTolerance(float)}. Note that the public static final
* Property objects are immutable: Is is not possible to alter their tolerance.
* Any attempt to do so will throw an Exception.
*
* Note that this implementation is not synchronized.
*
* Example: equals(FloatMatrix2D A, FloatMatrix2D B) is defined as
* follows
*
*
*
*
* { some other tests not related to tolerance go here }
* float epsilon = tolerance();
* for (int row=rows; --row >= 0;) {
* for (int column=columns; --column >= 0;) {
* //if (!(A.getQuick(row,column) == B.getQuick(row,column))) return false;
* if (Math.abs(A.getQuick(row,column) - B.getQuick(row,column)) > epsilon) return false;
* }
* }
* return true;
*
*
*
*
* Here are some example properties
*
*
* matrix
* 4 x 4
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
* 4 x 4
1 0 0 0
0 0 0 0
0 0 0 0
0 0 0 1
* 4 x 4
1 1 0 0
1 1 1 0
0 1 1 1
0 0 1 1
* 4 x 4
0 1 1 1
0 1 1 1
0 0 0 1
0 0 0 1
* 4 x 4
0 0 0 0
1 1 0 0
1 1 0 0
1 1 1 1
* 4 x 4
1 1 0 0
0 1 1 0
0 1 0 1
1 0 1 1
* 4 x 4
1 1 1 0
0 1 0 0
1 1 0 1
0 0 1 1
*
*
* upperBandwidth
* 0
* 0
* 1
* 3
* 0
* 1
* 2
*
*
* lowerBandwidth
* 0
* 0
* 1
* 0
* 3
* 3
* 2
*
*
* semiBandwidth
* 1
* 1
* 2
* 4
* 4
* 4
* 3
*
*
* description
* zero
* diagonal
* tridiagonal
* upper triangular
* lower triangular
* unstructured
*
* unstructured
*
*
*
*
* @author wolfgang.hoschek@cern.ch
* @version 1.1, 28/May/2000 (fixed strange bugs involving NaN, -inf, inf)
*
* @author Piotr Wendykier (piotr.wendykier@gmail.com)
*/
public class FloatProperty extends cern.colt.PersistentObject {
private static final long serialVersionUID = 1L;
/**
* The default Property object; currently has tolerance()==1.0E-5.
*/
public static final FloatProperty DEFAULT = new FloatProperty(1.0E-5f);
/**
* A Property object with tolerance()==0.0.
*/
public static final FloatProperty ZERO = new FloatProperty(0.0f);
/**
* A Property object with tolerance()==1.0E-7.
*/
public static final FloatProperty SEVEN = new FloatProperty(1.0E-7f);
protected float tolerance;
/**
* Not instantiable by no-arg constructor.
*/
private FloatProperty() {
this(1.0E-5f); // just to be on the safe side
}
/**
* Constructs an instance with a tolerance of
* Math.abs(newTolerance).
*/
public FloatProperty(float newTolerance) {
tolerance = Math.abs(newTolerance);
}
/**
* Returns a String with length blanks.
*/
protected static String blanks(int length) {
if (length < 0)
length = 0;
StringBuffer buf = new StringBuffer(length);
for (int k = 0; k < length; k++) {
buf.append(' ');
}
return buf.toString();
}
/**
* Checks whether the given matrix A is rectangular.
*
* @throws IllegalArgumentException
* if A.rows() < A.columns().
*/
public void checkRectangular(FloatMatrix2D A) {
if (A.rows() < A.columns()) {
throw new IllegalArgumentException("Matrix must be rectangular: " + AbstractFormatter.shape(A));
}
}
/**
* Checks whether the given matrix A is square.
*
* @throws IllegalArgumentException
* if A.rows() != A.columns().
*/
public void checkSquare(FloatMatrix2D A) {
if (A.rows() != A.columns())
throw new IllegalArgumentException("Matrix must be square: " + AbstractFormatter.shape(A));
}
public void checkDense(FloatMatrix2D A) {
if (!(A instanceof DenseFloatMatrix2D) && !(A instanceof DenseColumnFloatMatrix2D))
throw new IllegalArgumentException("Matrix must be dense");
}
public void checkDense(FloatMatrix1D A) {
if (!(A instanceof DenseFloatMatrix1D))
throw new IllegalArgumentException("Matrix must be dense");
}
public void checkSparse(FloatMatrix1D A) {
if (!(A instanceof SparseFloatMatrix1D))
throw new IllegalArgumentException("Matrix must be sparse");
}
public void checkSparse(FloatMatrix2D A) {
// if (!(A instanceof SparseFloatMatrix2D) && !(A instanceof RCFloatMatrix2D) && !(A instanceof RCMFloatMatrix2D)
// && !(A instanceof CCFloatMatrix2D) && !(A instanceof CCMFloatMatrix2D))
if (!(A instanceof SparseCCFloatMatrix2D) && !(A instanceof SparseRCFloatMatrix2D))
throw new IllegalArgumentException("Matrix must be sparse");
}
/**
* Returns the matrix's fraction of non-zero cells;
* A.cardinality() / A.size().
*/
public float density(FloatMatrix2D A) {
return A.cardinality() / (float) A.size();
}
/**
* Returns whether all cells of the given matrix A are equal to the
* given value. The result is true if and only if
* A != null and ! (Math.abs(value - A[i]) > tolerance())
* holds for all coordinates.
*
* @param A
* the first matrix to compare.
* @param value
* the value to compare against.
* @return true if the matrix is equal to the value; false
* otherwise.
*/
public boolean equals(final FloatMatrix1D A, final float value) {
if (A == null)
return false;
int size = (int) A.size();
final float epsilon = tolerance();
boolean result = false;
int nthreads = ConcurrencyUtils.getNumberOfThreads();
if ((nthreads > 1) && (size >= ConcurrencyUtils.getThreadsBeginN_1D())) {
nthreads = Math.min(nthreads, size);
Future[] futures = new Future[nthreads];
Boolean[] results = new Boolean[nthreads];
int k = size / nthreads;
for (int j = 0; j < nthreads; j++) {
final int firstIdx = j * k;
final int lastIdx = (j == nthreads - 1) ? size : firstIdx + k;
futures[j] = ConcurrencyUtils.submit(new Callable() {
public Boolean call() throws Exception {
for (int i = firstIdx; i < lastIdx; i++) {
float x = A.getQuick(i);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
return true;
}
});
}
try {
for (int j = 0; j < nthreads; j++) {
results[j] = (Boolean) futures[j].get();
}
result = results[0].booleanValue();
for (int j = 1; j < nthreads; j++) {
result = result && results[j].booleanValue();
}
} catch (ExecutionException ex) {
ex.printStackTrace();
} catch (InterruptedException e) {
e.printStackTrace();
}
return result;
} else {
for (int i = 0; i < size; i++) {
float x = A.getQuick(i);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
return true;
}
}
/**
* Returns whether both given matrices A and B are equal.
* The result is true if A==B. Otherwise, the result is
* true if and only if both arguments are != null, have
* the same size and ! (Math.abs(A[i] - B[i]) > tolerance()) holds
* for all indexes.
*
* @param A
* the first matrix to compare.
* @param B
* the second matrix to compare.
* @return true if both matrices are equal; false
* otherwise.
*/
public boolean equals(final FloatMatrix1D A, final FloatMatrix1D B) {
if (A == B)
return true;
if (!(A != null && B != null))
return false;
int size = (int) A.size();
if (size != B.size())
return false;
final float epsilon = tolerance();
boolean result = false;
int nthreads = ConcurrencyUtils.getNumberOfThreads();
if ((nthreads > 1) && (size >= ConcurrencyUtils.getThreadsBeginN_1D())) {
nthreads = Math.min(nthreads, size);
Future[] futures = new Future[nthreads];
Boolean[] results = new Boolean[nthreads];
int k = size / nthreads;
for (int j = 0; j < nthreads; j++) {
final int firstIdx = j * k;
final int lastIdx = (j == nthreads - 1) ? size : firstIdx + k;
futures[j] = ConcurrencyUtils.submit(new Callable() {
public Boolean call() throws Exception {
for (int i = firstIdx; i < lastIdx; i++) {
float x = A.getQuick(i);
float value = B.getQuick(i);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
return true;
}
});
}
try {
for (int j = 0; j < nthreads; j++) {
results[j] = (Boolean) futures[j].get();
}
result = results[0].booleanValue();
for (int j = 1; j < nthreads; j++) {
result = result && results[j].booleanValue();
}
} catch (ExecutionException ex) {
ex.printStackTrace();
} catch (InterruptedException e) {
e.printStackTrace();
}
return result;
} else {
for (int i = 0; i < size; i++) {
float x = A.getQuick(i);
float value = B.getQuick(i);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
return true;
}
}
/**
* Returns whether all cells of the given matrix A are equal to the
* given value. The result is true if and only if
* A != null and
* ! (Math.abs(value - A[row,col]) > tolerance()) holds for all
* coordinates.
*
* @param A
* the first matrix to compare.
* @param value
* the value to compare against.
* @return true if the matrix is equal to the value; false
* otherwise.
*/
public boolean equals(final FloatMatrix2D A, final float value) {
if (A == null)
return false;
final int rows = A.rows();
final int columns = A.columns();
boolean result = false;
final float epsilon = tolerance();
int nthreads = ConcurrencyUtils.getNumberOfThreads();
if ((nthreads > 1) && (A.size() >= ConcurrencyUtils.getThreadsBeginN_2D())) {
nthreads = Math.min(nthreads, A.rows());
Future[] futures = new Future[nthreads];
Boolean[] results = new Boolean[nthreads];
int k = A.rows() / nthreads;
for (int j = 0; j < nthreads; j++) {
final int firstRow = j * k;
final int lastRow = (j == nthreads - 1) ? A.rows() : firstRow + k;
futures[j] = ConcurrencyUtils.submit(new Callable() {
public Boolean call() throws Exception {
for (int r = firstRow; r < lastRow; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
return true;
}
});
}
try {
for (int j = 0; j < nthreads; j++) {
results[j] = (Boolean) futures[j].get();
}
result = results[0].booleanValue();
for (int j = 1; j < nthreads; j++) {
result = result && results[j].booleanValue();
}
} catch (ExecutionException ex) {
ex.printStackTrace();
} catch (InterruptedException e) {
e.printStackTrace();
}
return result;
} else {
for (int r = 0; r < rows; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
return true;
}
}
/**
* Returns whether both given matrices A and B are equal.
* The result is true if A==B. Otherwise, the result is
* true if and only if both arguments are != null, have
* the same number of columns and rows and
* ! (Math.abs(A[row,col] - B[row,col]) > tolerance()) holds for
* all coordinates.
*
* @param A
* the first matrix to compare.
* @param B
* the second matrix to compare.
* @return true if both matrices are equal; false
* otherwise.
*/
public boolean equals(final FloatMatrix2D A, final FloatMatrix2D B) {
if (A == B)
return true;
if (!(A != null && B != null))
return false;
final int rows = A.rows();
final int columns = A.columns();
if (columns != B.columns() || rows != B.rows())
return false;
boolean result = false;
final float epsilon = tolerance();
int nthreads = ConcurrencyUtils.getNumberOfThreads();
if ((nthreads > 1) && (A.size() >= ConcurrencyUtils.getThreadsBeginN_2D())) {
nthreads = Math.min(nthreads, A.rows());
Future[] futures = new Future[nthreads];
Boolean[] results = new Boolean[nthreads];
int k = A.rows() / nthreads;
for (int j = 0; j < nthreads; j++) {
final int firstRow = j * k;
final int lastRow = (j == nthreads - 1) ? A.rows() : firstRow + k;
futures[j] = ConcurrencyUtils.submit(new Callable() {
public Boolean call() throws Exception {
for (int r = firstRow; r < lastRow; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(r, c);
float value = B.getQuick(r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
return true;
}
});
}
try {
for (int j = 0; j < nthreads; j++) {
results[j] = (Boolean) futures[j].get();
}
result = results[0].booleanValue();
for (int j = 1; j < nthreads; j++) {
result = result && results[j].booleanValue();
}
} catch (ExecutionException ex) {
ex.printStackTrace();
} catch (InterruptedException e) {
e.printStackTrace();
}
return result;
} else {
for (int r = 0; r < rows; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(r, c);
float value = B.getQuick(r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
return true;
}
}
/**
* Returns whether all cells of the given matrix A are equal to the
* given value. The result is true if and only if
* A != null and
* ! (Math.abs(value - A[slice,row,col]) > tolerance()) holds for
* all coordinates.
*
* @param A
* the first matrix to compare.
* @param value
* the value to compare against.
* @return true if the matrix is equal to the value; false
* otherwise.
*/
public boolean equals(final FloatMatrix3D A, final float value) {
if (A == null)
return false;
final int slices = A.slices();
final int rows = A.rows();
final int columns = A.columns();
boolean result = false;
final float epsilon = tolerance();
int nthreads = ConcurrencyUtils.getNumberOfThreads();
if ((nthreads > 1) && (A.size() >= ConcurrencyUtils.getThreadsBeginN_3D())) {
nthreads = Math.min(nthreads, slices);
Future[] futures = new Future[nthreads];
Boolean[] results = new Boolean[nthreads];
int k = slices / nthreads;
for (int j = 0; j < nthreads; j++) {
final int firstSlice = j * k;
final int lastSlice = (j == nthreads - 1) ? slices : firstSlice + k;
futures[j] = ConcurrencyUtils.submit(new Callable() {
public Boolean call() throws Exception {
for (int s = firstSlice; s < lastSlice; s++) {
for (int r = 0; r < rows; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(s, r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
}
return true;
}
});
}
try {
for (int j = 0; j < nthreads; j++) {
results[j] = (Boolean) futures[j].get();
}
result = results[0].booleanValue();
for (int j = 1; j < nthreads; j++) {
result = result && results[j].booleanValue();
}
} catch (ExecutionException ex) {
ex.printStackTrace();
} catch (InterruptedException e) {
e.printStackTrace();
}
return result;
} else {
for (int s = 0; s < slices; s++) {
for (int r = 0; r < rows; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(s, r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
}
return true;
}
}
/**
* Returns whether both given matrices A and B are equal.
* The result is true if A==B. Otherwise, the result is
* true if and only if both arguments are != null, have
* the same number of columns, rows and slices, and
* ! (Math.abs(A[slice,row,col] - B[slice,row,col]) > tolerance())
* holds for all coordinates.
*
* @param A
* the first matrix to compare.
* @param B
* the second matrix to compare.
* @return true if both matrices are equal; false
* otherwise.
*/
public boolean equals(final FloatMatrix3D A, final FloatMatrix3D B) {
if (A == B)
return true;
if (!(A != null && B != null))
return false;
final int slices = A.slices();
final int rows = A.rows();
final int columns = A.columns();
if (columns != B.columns() || rows != B.rows() || slices != B.slices())
return false;
boolean result = false;
final float epsilon = tolerance();
int nthreads = ConcurrencyUtils.getNumberOfThreads();
if ((nthreads > 1) && (A.size() >= ConcurrencyUtils.getThreadsBeginN_3D())) {
nthreads = Math.min(nthreads, slices);
Future[] futures = new Future[nthreads];
Boolean[] results = new Boolean[nthreads];
int k = slices / nthreads;
for (int j = 0; j < nthreads; j++) {
final int firstSlice = j * k;
final int lastSlice = (j == nthreads - 1) ? slices : firstSlice + k;
futures[j] = ConcurrencyUtils.submit(new Callable() {
public Boolean call() throws Exception {
for (int s = firstSlice; s < lastSlice; s++) {
for (int r = 0; r < rows; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(s, r, c);
float value = B.getQuick(s, r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
}
return true;
}
});
}
try {
for (int j = 0; j < nthreads; j++) {
results[j] = (Boolean) futures[j].get();
}
result = results[0].booleanValue();
for (int j = 1; j < nthreads; j++) {
result = result && results[j].booleanValue();
}
} catch (ExecutionException ex) {
ex.printStackTrace();
} catch (InterruptedException e) {
e.printStackTrace();
}
return result;
} else {
for (int s = 0; s < slices; s++) {
for (int r = 0; r < rows; r++) {
for (int c = 0; c < columns; c++) {
float x = A.getQuick(s, r, c);
float value = B.getQuick(s, r, c);
float diff = Math.abs(value - x);
if ((diff != diff) && ((value != value && x != x) || value == x))
diff = 0;
if (!(diff <= epsilon)) {
return false;
}
}
}
}
return true;
}
}
/**
* Modifies the given matrix square matrix A such that it is
* diagonally dominant by row and column, hence non-singular, hence
* invertible. For testing purposes only.
*
* @param A
* the square matrix to modify.
* @throws IllegalArgumentException
* if !isSquare(A).
*/
public void generateNonSingular(FloatMatrix2D A) {
checkSquare(A);
cern.jet.math.tfloat.FloatFunctions F = cern.jet.math.tfloat.FloatFunctions.functions;
int min = Math.min(A.rows(), A.columns());
for (int i = min; --i >= 0;) {
A.setQuick(i, i, 0);
}
for (int i = min; --i >= 0;) {
float rowSum = A.viewRow(i).aggregate(FloatFunctions.plus, FloatFunctions.abs);
float colSum = A.viewColumn(i).aggregate(FloatFunctions.plus, FloatFunctions.abs);
A.setQuick(i, i, Math.max(rowSum, colSum) + i + 1);
}
}
/**
*/
protected static String get(cern.colt.list.tobject.ObjectArrayList list, int index) {
return ((String) list.get(index));
}
/**
* A matrix A is diagonal if A[i,j] == 0 whenever
* i != j. Matrix may but need not be square.
*/
public boolean isDiagonal(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int row = rows; --row >= 0;) {
for (int column = columns; --column >= 0;) {
if (row != column && !(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
return true;
}
/**
* A matrix A is diagonally dominant by column if the
* absolute value of each diagonal element is larger than the sum of the
* absolute values of the off-diagonal elements in the corresponding column.
*
* returns true if for all i: abs(A[i,i]) > Sum(abs(A[j,i])); j != i.
* Matrix may but need not be square.
*
* Note: Ignores tolerance.
*/
public boolean isDiagonallyDominantByColumn(FloatMatrix2D A) {
cern.jet.math.tfloat.FloatFunctions F = cern.jet.math.tfloat.FloatFunctions.functions;
int min = Math.min(A.rows(), A.columns());
for (int i = min; --i >= 0;) {
float diag = Math.abs(A.getQuick(i, i));
diag += diag;
if (diag <= A.viewColumn(i).aggregate(FloatFunctions.plus, FloatFunctions.abs))
return false;
}
return true;
}
/**
* A matrix A is diagonally dominant by row if the absolute
* value of each diagonal element is larger than the sum of the absolute
* values of the off-diagonal elements in the corresponding row.
* returns true if for all i: abs(A[i,i]) > Sum(abs(A[i,j])); j != i.
* Matrix may but need not be square.
*
* Note: Ignores tolerance.
*/
public boolean isDiagonallyDominantByRow(FloatMatrix2D A) {
cern.jet.math.tfloat.FloatFunctions F = cern.jet.math.tfloat.FloatFunctions.functions;
int min = Math.min(A.rows(), A.columns());
for (int i = min; --i >= 0;) {
float diag = Math.abs(A.getQuick(i, i));
diag += diag;
if (diag <= A.viewRow(i).aggregate(FloatFunctions.plus, FloatFunctions.abs))
return false;
}
return true;
}
/**
* A matrix A is an identity matrix if A[i,i] == 1
* and all other cells are zero. Matrix may but need not be square.
*/
public boolean isIdentity(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int row = rows; --row >= 0;) {
for (int column = columns; --column >= 0;) {
float v = A.getQuick(row, column);
if (row == column) {
if (!(Math.abs(1 - v) < epsilon))
return false;
} else if (!(Math.abs(v) <= epsilon))
return false;
}
}
return true;
}
/**
* A matrix A is lower bidiagonal if A[i,j]==0
* unless i==j || i==j+1. Matrix may but need not be square.
*/
public boolean isLowerBidiagonal(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int row = rows; --row >= 0;) {
for (int column = columns; --column >= 0;) {
if (!(row == column || row == column + 1)) {
if (!(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
}
return true;
}
/**
* A matrix A is lower triangular if A[i,j]==0
* whenever i < j. Matrix may but need not be square.
*/
public boolean isLowerTriangular(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int column = columns; --column >= 0;) {
for (int row = Math.min(column, rows); --row >= 0;) {
if (!(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
return true;
}
/**
* A matrix A is non-negative if A[i,j] >= 0
* holds for all cells.
*
* Note: Ignores tolerance.
*/
public boolean isNonNegative(FloatMatrix2D A) {
int rows = A.rows();
int columns = A.columns();
for (int row = rows; --row >= 0;) {
for (int column = columns; --column >= 0;) {
if (!(A.getQuick(row, column) >= 0))
return false;
}
}
return true;
}
/**
* A square matrix A is orthogonal if
* A*transpose(A) = I.
*
* @throws IllegalArgumentException
* if !isSquare(A).
*/
public boolean isOrthogonal(FloatMatrix2D A) {
checkSquare(A);
return equals(A.zMult(A, null, 1, 0, false, true), cern.colt.matrix.tfloat.FloatFactory2D.dense.identity(A
.rows()));
}
/**
* A matrix A is positive if A[i,j] > 0 holds
* for all cells.
*
* Note: Ignores tolerance.
*/
public boolean isPositive(FloatMatrix2D A) {
int rows = A.rows();
int columns = A.columns();
for (int row = rows; --row >= 0;) {
for (int column = columns; --column >= 0;) {
if (!(A.getQuick(row, column) > 0))
return false;
}
}
return true;
}
/**
* A matrix A is singular if it has no inverse, that is, iff
* det(A)==0.
*/
public boolean isSingular(FloatMatrix2D A) {
return !(Math.abs(DenseFloatAlgebra.DEFAULT.det(A)) >= tolerance());
}
/**
* A square matrix A is skew-symmetric if
* A = -transpose(A), that is A[i,j] == -A[j,i].
*
* @throws IllegalArgumentException
* if !isSquare(A).
*/
public boolean isSkewSymmetric(FloatMatrix2D A) {
checkSquare(A);
float epsilon = tolerance();
int rows = A.rows();
for (int row = rows; --row >= 0;) {
for (int column = rows; --column >= 0;) {
if (!(Math.abs(A.getQuick(row, column) + A.getQuick(column, row)) <= epsilon))
return false;
}
}
return true;
}
/**
* A matrix A is square if it has the same number of rows
* and columns.
*/
public boolean isSquare(FloatMatrix2D A) {
return A.rows() == A.columns();
}
/**
* A matrix A is strictly lower triangular if
* A[i,j]==0 whenever i <= j. Matrix may but need not
* be square.
*/
public boolean isStrictlyLowerTriangular(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int column = columns; --column >= 0;) {
for (int row = Math.min(rows, column + 1); --row >= 0;) {
if (!(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
return true;
}
/**
* A matrix A is strictly triangular if it is triangular and
* its diagonal elements all equal 0. Matrix may but need not be square.
*/
public boolean isStrictlyTriangular(FloatMatrix2D A) {
if (!isTriangular(A))
return false;
float epsilon = tolerance();
for (int i = Math.min(A.rows(), A.columns()); --i >= 0;) {
if (!(Math.abs(A.getQuick(i, i)) <= epsilon))
return false;
}
return true;
}
/**
* A matrix A is strictly upper triangular if
* A[i,j]==0 whenever i >= j. Matrix may but need not
* be square.
*/
public boolean isStrictlyUpperTriangular(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int column = columns; --column >= 0;) {
for (int row = rows; --row >= column;) {
if (!(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
return true;
}
/**
* A matrix A is symmetric if A = tranpose(A), that
* is A[i,j] == A[j,i].
*
* @throws IllegalArgumentException
* if !isSquare(A).
*/
public boolean isSymmetric(FloatMatrix2D A) {
checkSquare(A);
return equals(A, A.viewDice());
}
/**
* A matrix A is triangular iff it is either upper or lower
* triangular. Matrix may but need not be square.
*/
public boolean isTriangular(FloatMatrix2D A) {
return isLowerTriangular(A) || isUpperTriangular(A);
}
/**
* A matrix A is tridiagonal if A[i,j]==0 whenever
* Math.abs(i-j) > 1. Matrix may but need not be square.
*/
public boolean isTridiagonal(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int row = rows; --row >= 0;) {
for (int column = columns; --column >= 0;) {
if (Math.abs(row - column) > 1) {
if (!(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
}
return true;
}
/**
* A matrix A is unit triangular if it is triangular and its
* diagonal elements all equal 1. Matrix may but need not be square.
*/
public boolean isUnitTriangular(FloatMatrix2D A) {
if (!isTriangular(A))
return false;
float epsilon = tolerance();
for (int i = Math.min(A.rows(), A.columns()); --i >= 0;) {
if (!(Math.abs(1 - A.getQuick(i, i)) <= epsilon))
return false;
}
return true;
}
/**
* A matrix A is upper bidiagonal if A[i,j]==0
* unless i==j || i==j-1. Matrix may but need not be square.
*/
public boolean isUpperBidiagonal(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int row = rows; --row >= 0;) {
for (int column = columns; --column >= 0;) {
if (!(row == column || row == column - 1)) {
if (!(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
}
return true;
}
/**
* A matrix A is upper triangular if A[i,j]==0
* whenever i > j. Matrix may but need not be square.
*/
public boolean isUpperTriangular(FloatMatrix2D A) {
float epsilon = tolerance();
int rows = A.rows();
int columns = A.columns();
for (int column = columns; --column >= 0;) {
for (int row = rows; --row > column;) {
if (!(Math.abs(A.getQuick(row, column)) <= epsilon))
return false;
}
}
return true;
}
/**
* A matrix A is zero if all its cells are zero.
*/
public boolean isZero(FloatMatrix2D A) {
return equals(A, 0);
}
/**
* The lower bandwidth of a square matrix A is the maximum
* i-j for which A[i,j] is nonzero and i > j.
* A banded matrix has a "band" about the diagonal. Diagonal,
* tridiagonal and triangular matrices are special cases.
*
* @param A
* the square matrix to analyze.
* @return the lower bandwith.
* @throws IllegalArgumentException
* if !isSquare(A).
* @see #semiBandwidth(FloatMatrix2D)
* @see #upperBandwidth(FloatMatrix2D)
*/
public int lowerBandwidth(FloatMatrix2D A) {
checkSquare(A);
float epsilon = tolerance();
int rows = A.rows();
for (int k = rows; --k >= 0;) {
for (int i = rows - k; --i >= 0;) {
int j = i + k;
if (!(Math.abs(A.getQuick(j, i)) <= epsilon))
return k;
}
}
return 0;
}
/**
* Returns the semi-bandwidth of the given square matrix A.
* A banded matrix has a "band" about the diagonal. It is a matrix
* with all cells equal to zero, with the possible exception of the cells
* along the diagonal line, the k diagonal lines above the
* diagonal, and the k diagonal lines below the diagonal. The
* semi-bandwith l is the number k+1. The bandwidth p
* is the number 2*k + 1. For example, a tridiagonal matrix
* corresponds to k=1, l=2, p=3, a diagonal or zero matrix
* corresponds to k=0, l=1, p=1,
*
* The upper bandwidth is the maximum j-i for which
* A[i,j] is nonzero and j > i. The lower
* bandwidth is the maximum i-j for which A[i,j] is
* nonzero and i > j. Diagonal, tridiagonal and triangular
* matrices are special cases.
*
* Examples:
*
*
* matrix
* 4 x 4
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
* 4 x 4
1 0 0 0
0 0 0 0
0 0 0 0
0 0 0 1
* 4 x 4
1 1 0 0
1 1 1 0
0 1 1 1
0 0 1 1
* 4 x 4
0 1 1 1
0 1 1 1
0 0 0 1
0 0 0 1
* 4 x 4
0 0 0 0
1 1 0 0
1 1 0 0
1 1 1 1
* 4 x 4
1 1 0 0
0 1 1 0
0 1 0 1
1 0 1 1
* 4 x 4
1 1 1 0
0 1 0 0
1 1 0 1
0 0 1 1
*
*
* upperBandwidth
* 0
* 0
* 1
* 3
* 0
* 1
* 2
*
*
* lowerBandwidth
* 0
* 0
* 1
* 0
* 3
* 3
* 2
*
*
* semiBandwidth
* 1
* 1
* 2
* 4
* 4
* 4
* 3
*
*
* description
* zero
* diagonal
* tridiagonal
* upper triangular
* lower triangular
*
* unstructured
*
* unstructured
*
*
*
* @param A
* the square matrix to analyze.
* @return the semi-bandwith l.
* @throws IllegalArgumentException
* if !isSquare(A).
* @see #lowerBandwidth(FloatMatrix2D)
* @see #upperBandwidth(FloatMatrix2D)
*/
public int semiBandwidth(FloatMatrix2D A) {
checkSquare(A);
float epsilon = tolerance();
int rows = A.rows();
for (int k = rows; --k >= 0;) {
for (int i = rows - k; --i >= 0;) {
int j = i + k;
if (!(Math.abs(A.getQuick(j, i)) <= epsilon))
return k + 1;
if (!(Math.abs(A.getQuick(i, j)) <= epsilon))
return k + 1;
}
}
return 1;
}
/**
* Sets the tolerance to Math.abs(newTolerance).
*
* @throws UnsupportedOperationException
* if this==DEFAULT || this==ZERO || this==TWELVE.
*/
public void setTolerance(float newTolerance) {
if (this == DEFAULT || this == ZERO || this == SEVEN) {
throw new IllegalArgumentException("Attempted to modify immutable object.");
}
tolerance = Math.abs(newTolerance);
}
/**
* Returns the current tolerance.
*/
public float tolerance() {
return tolerance;
}
/**
* Returns summary information about the given matrix A. That is a
* String with (propertyName, propertyValue) pairs. Useful for debugging or
* to quickly get the rough picture of a matrix. For example,
*
*
* density : 0.9
* isDiagonal : false
* isDiagonallyDominantByRow : false
* isDiagonallyDominantByColumn : false
* isIdentity : false
* isLowerBidiagonal : false
* isLowerTriangular : false
* isNonNegative : true
* isOrthogonal : Illegal operation or error: Matrix must be square.
* isPositive : true
* isSingular : Illegal operation or error: Matrix must be square.
* isSkewSymmetric : Illegal operation or error: Matrix must be square.
* isSquare : false
* isStrictlyLowerTriangular : false
* isStrictlyTriangular : false
* isStrictlyUpperTriangular : false
* isSymmetric : Illegal operation or error: Matrix must be square.
* isTriangular : false
* isTridiagonal : false
* isUnitTriangular : false
* isUpperBidiagonal : false
* isUpperTriangular : false
* isZero : false
* lowerBandwidth : Illegal operation or error: Matrix must be square.
* semiBandwidth : Illegal operation or error: Matrix must be square.
* upperBandwidth : Illegal operation or error: Matrix must be square.
*
*
*/
public String toString(FloatMatrix2D A) {
final cern.colt.list.tobject.ObjectArrayList names = new cern.colt.list.tobject.ObjectArrayList();
final cern.colt.list.tobject.ObjectArrayList values = new cern.colt.list.tobject.ObjectArrayList();
String unknown = "Illegal operation or error: ";
// determine properties
names.add("density");
try {
values.add(String.valueOf(density(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
// determine properties
names.add("isDiagonal");
try {
values.add(String.valueOf(isDiagonal(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
// determine properties
names.add("isDiagonallyDominantByRow");
try {
values.add(String.valueOf(isDiagonallyDominantByRow(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
// determine properties
names.add("isDiagonallyDominantByColumn");
try {
values.add(String.valueOf(isDiagonallyDominantByColumn(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isIdentity");
try {
values.add(String.valueOf(isIdentity(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isLowerBidiagonal");
try {
values.add(String.valueOf(isLowerBidiagonal(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isLowerTriangular");
try {
values.add(String.valueOf(isLowerTriangular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isNonNegative");
try {
values.add(String.valueOf(isNonNegative(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isOrthogonal");
try {
values.add(String.valueOf(isOrthogonal(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isPositive");
try {
values.add(String.valueOf(isPositive(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isSingular");
try {
values.add(String.valueOf(isSingular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isSkewSymmetric");
try {
values.add(String.valueOf(isSkewSymmetric(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isSquare");
try {
values.add(String.valueOf(isSquare(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isStrictlyLowerTriangular");
try {
values.add(String.valueOf(isStrictlyLowerTriangular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isStrictlyTriangular");
try {
values.add(String.valueOf(isStrictlyTriangular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isStrictlyUpperTriangular");
try {
values.add(String.valueOf(isStrictlyUpperTriangular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isSymmetric");
try {
values.add(String.valueOf(isSymmetric(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isTriangular");
try {
values.add(String.valueOf(isTriangular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isTridiagonal");
try {
values.add(String.valueOf(isTridiagonal(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isUnitTriangular");
try {
values.add(String.valueOf(isUnitTriangular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isUpperBidiagonal");
try {
values.add(String.valueOf(isUpperBidiagonal(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isUpperTriangular");
try {
values.add(String.valueOf(isUpperTriangular(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("isZero");
try {
values.add(String.valueOf(isZero(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("lowerBandwidth");
try {
values.add(String.valueOf(lowerBandwidth(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("semiBandwidth");
try {
values.add(String.valueOf(semiBandwidth(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
names.add("upperBandwidth");
try {
values.add(String.valueOf(upperBandwidth(A)));
} catch (IllegalArgumentException exc) {
values.add(unknown + exc.getMessage());
}
// sort ascending by property name
cern.colt.function.tint.IntComparator comp = new cern.colt.function.tint.IntComparator() {
public int compare(int a, int b) {
return get(names, a).compareTo(get(names, b));
}
};
cern.colt.Swapper swapper = new cern.colt.Swapper() {
public void swap(int a, int b) {
Object tmp;
tmp = names.get(a);
names.set(a, names.get(b));
names.set(b, tmp);
tmp = values.get(a);
values.set(a, values.get(b));
values.set(b, tmp);
}
};
cern.colt.GenericSorting.quickSort(0, names.size(), comp, swapper);
// determine padding for nice formatting
int maxLength = 0;
for (int i = 0; i < names.size(); i++) {
int length = ((String) names.get(i)).length();
maxLength = Math.max(length, maxLength);
}
// finally, format properties
StringBuffer buf = new StringBuffer();
for (int i = 0; i < names.size(); i++) {
String name = ((String) names.get(i));
buf.append(name);
buf.append(blanks(maxLength - name.length()));
buf.append(" : ");
buf.append(values.get(i));
if (i < names.size() - 1)
buf.append('\n');
}
return buf.toString();
}
/**
* The upper bandwidth of a square matrix A is the maximum
* j-i for which A[i,j] is nonzero and j > i.
* A banded matrix has a "band" about the diagonal. Diagonal,
* tridiagonal and triangular matrices are special cases.
*
* @param A
* the square matrix to analyze.
* @return the upper bandwith.
* @throws IllegalArgumentException
* if !isSquare(A).
* @see #semiBandwidth(FloatMatrix2D)
* @see #lowerBandwidth(FloatMatrix2D)
*/
public int upperBandwidth(FloatMatrix2D A) {
checkSquare(A);
float epsilon = tolerance();
int rows = A.rows();
for (int k = rows; --k >= 0;) {
for (int i = rows - k; --i >= 0;) {
int j = i + k;
if (!(Math.abs(A.getQuick(i, j)) <= epsilon))
return k;
}
}
return 0;
}
}
Ads help maintain this website.