Java Code Examples for cern.colt.matrix.DoubleMatrix2D#getQuick()

private static DoubleMatrix2D blockDiagonalExponential(double distance, DoubleMatrix2D mat) {
    for (int i = 0; i < mat.rows(); i++) {
        if ((i + 1) < mat.rows() && mat.getQuick(i, i + 1) != 0) {
            double a = mat.getQuick(i, i);
            double b = mat.getQuick(i, i + 1);
            double expat = Math.exp(distance * a);
            double cosbt = Math.cos(distance * b);
            double sinbt = Math.sin(distance * b);
            mat.setQuick(i, i, expat * cosbt);
            mat.setQuick(i + 1, i + 1, expat * cosbt);
            mat.setQuick(i, i + 1, expat * sinbt);
            mat.setQuick(i + 1, i, -expat * sinbt);
            i++; // processed two entries in loop
        } else
            mat.setQuick(i, i, Math.exp(distance * mat.getQuick(i, i))); // 1x1 block
    }
    return mat;
}
 

/**
 * Modifies the given covariance matrix to be a correlation matrix (in-place).
 * The correlation matrix is a square, symmetric matrix consisting of nothing but correlation coefficients.
 * The rows and the columns represent the variables, the cells represent correlation coefficients. 
 * The diagonal cells (i.e. the correlation between a variable and itself) will equal 1, for the simple reason that the correlation coefficient of a variable with itself equals 1. 
 * The correlation of two column vectors x and y is given by <tt>corr(x,y) = cov(x,y) / (stdDev(x)*stdDev(y))</tt> (Pearson's correlation coefficient).
 * A correlation coefficient varies between -1 (for a perfect negative relationship) to +1 (for a perfect positive relationship). 
 * See the <A HREF="http://www.cquest.utoronto.ca/geog/ggr270y/notes/not05efg.html"> math definition</A>
 * and <A HREF="http://www.stat.berkeley.edu/users/stark/SticiGui/Text/gloss.htm#correlation_coef"> another def</A>.
 * Compares two column vectors at a time. Use dice views to compare two row vectors at a time.
 * 
 * @param covariance a covariance matrix, as, for example, returned by method {@link #covariance(DoubleMatrix2D)}.
 * @return the modified covariance, now correlation matrix (for convenience only).
 */
public static DoubleMatrix2D correlation(DoubleMatrix2D covariance) {
	for (int i=covariance.columns(); --i >= 0; ) {
		for (int j=i; --j >= 0; ) {
			double stdDev1 = Math.sqrt(covariance.getQuick(i,i));
			double stdDev2 = Math.sqrt(covariance.getQuick(j,j));
			double cov = covariance.getQuick(i,j);
			double corr = cov / (stdDev1*stdDev2);
			
			covariance.setQuick(i,j,corr);
			covariance.setQuick(j,i,corr); // symmetric
		}
	}
	for (int i=covariance.columns(); --i >= 0; ) covariance.setQuick(i,i,1);

	return covariance;	
}
 

/**
 * Modifies the given covariance matrix to be a correlation matrix (in-place).
 * The correlation matrix is a square, symmetric matrix consisting of nothing but correlation coefficients.
 * The rows and the columns represent the variables, the cells represent correlation coefficients. 
 * The diagonal cells (i.e. the correlation between a variable and itself) will equal 1, for the simple reason that the correlation coefficient of a variable with itself equals 1. 
 * The correlation of two column vectors x and y is given by <tt>corr(x,y) = cov(x,y) / (stdDev(x)*stdDev(y))</tt> (Pearson's correlation coefficient).
 * A correlation coefficient varies between -1 (for a perfect negative relationship) to +1 (for a perfect positive relationship). 
 * See the <A HREF="http://www.cquest.utoronto.ca/geog/ggr270y/notes/not05efg.html"> math definition</A>
 * and <A HREF="http://www.stat.berkeley.edu/users/stark/SticiGui/Text/gloss.htm#correlation_coef"> another def</A>.
 * Compares two column vectors at a time. Use dice views to compare two row vectors at a time.
 * 
 * @param covariance a covariance matrix, as, for example, returned by method {@link #covariance(DoubleMatrix2D)}.
 * @return the modified covariance, now correlation matrix (for convenience only).
 */
public static DoubleMatrix2D correlation(DoubleMatrix2D covariance) {
	for (int i=covariance.columns(); --i >= 0; ) {
		for (int j=i; --j >= 0; ) {
			double stdDev1 = Math.sqrt(covariance.getQuick(i,i));
			double stdDev2 = Math.sqrt(covariance.getQuick(j,j));
			double cov = covariance.getQuick(i,j);
			double corr = cov / (stdDev1*stdDev2);
			
			covariance.setQuick(i,j,corr);
			covariance.setQuick(j,i,corr); // symmetric
		}
	}
	for (int i=covariance.columns(); --i >= 0; ) covariance.setQuick(i,i,1);

	return covariance;	
}
 

public void dsymv(boolean isUpperTriangular, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
	if (isUpperTriangular) A = A.viewDice();
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size() || size!=y.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort() + ", " + y.toStringShort());
	}
	DoubleMatrix1D tmp = x.like();
	for (int i = 0; i < size; i++) {
		double sum = 0;
		for (int j = 0; j <= i; j++) {
			sum += A.getQuick(i,j) * x.getQuick(j);
		}
		for (int j = i + 1; j < size; j++) {
			sum += A.getQuick(j,i) * x.getQuick(j);
		}
		tmp.setQuick(i, alpha * sum + beta * y.getQuick(i));
	}
	y.assign(tmp);
}
 

static DoubleMatrix1D diagMult(DoubleMatrix2D A, DoubleMatrix2D B){
    //Return the diagonal of the product of two matrices
    //A^T and B should have the same dimensions
    int m = A.rows();
    int n = A.columns();
    if(B.rows()!=n || B.columns()!=m){
        throw new Error("Size mismatch in diagMult: A is "+m+"x"+n+
                ", B is "+B.rows()+"x"+B.columns());
    }

    DoubleMatrix1D ans = DoubleFactory1D.dense.make(m);

    for(int i=0; i<m; i++){
        double s = 0;
        for(int k=0; k<n; k++)
            s += A.getQuick(i, k)*B.getQuick(k, i);
        ans.setQuick(i,s);
    }

    return ans;
}
 

static DoubleMatrix2D QuQt(DoubleMatrix2D Q, DoubleMatrix1D u){
    //Return matrix product of Q * diagonal version of u * Q^T
    //answer = \sum_i ( u_i * q_i * q_i')  is a (d+1)x(d+1) matrix

    int m = Q.rows();
    int n = Q.columns();
    if(u.size()!=n){
        throw new Error("Size mismatch in QuQt: "+n+" columns in Q, u length="+u.size());
    }

    DoubleMatrix2D ans = DoubleFactory2D.dense.make(m,m);

    for(int i=0; i<m; i++){
        for(int j=0; j<m; j++){
            double s = 0;
            for(int k=0; k<n; k++)
                s += Q.getQuick(i,k)*Q.getQuick(j,k)*u.get(k);
            ans.setQuick(i, j, s);
        }
    }

    return ans;
}
 

/**
 * Returns whether both given matrices <tt>A</tt> and <tt>B</tt> are equal.
 * The result is <tt>true</tt> if <tt>A==B</tt>. 
 * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>, 
 * have the same number of columns and rows and 
 * <tt>! (Math.abs(A[row,col] - B[row,col]) > tolerance())</tt> holds for all coordinates.
 * @param   A   the first matrix to compare.
 * @param   B   the second matrix to compare.
 * @return  <tt>true</tt> if both matrices are equal;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix2D A, DoubleMatrix2D B) {
	if (A==B) return true;
	if (! (A != null && B != null)) return false;
	int rows = A.rows();
	int columns = A.columns();
	if (columns != B.columns() || rows != B.rows()) return false;

	double 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;
			double x = A.getQuick(row,column);
			double value = B.getQuick(row,column);
			double 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 <tt>A</tt> and <tt>B</tt> are equal.
 * The result is <tt>true</tt> if <tt>A==B</tt>. 
 * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>, 
 * have the same number of columns and rows and 
 * <tt>! (Math.abs(A[row,col] - B[row,col]) > tolerance())</tt> holds for all coordinates.
 * @param   A   the first matrix to compare.
 * @param   B   the second matrix to compare.
 * @return  <tt>true</tt> if both matrices are equal;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix2D A, DoubleMatrix2D B) {
	if (A==B) return true;
	if (! (A != null && B != null)) return false;
	int rows = A.rows();
	int columns = A.columns();
	if (columns != B.columns() || rows != B.rows()) return false;

	double 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;
			double x = A.getQuick(row,column);
			double value = B.getQuick(row,column);
			double diff = Math.abs(value - x);
			if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
			if (!(diff <= epsilon)) return false;
		}
	}
	return true;
}
 

public void dsymv(boolean isUpperTriangular, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
	if (isUpperTriangular) A = A.viewDice();
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size() || size!=y.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort() + ", " + y.toStringShort());
	}
	DoubleMatrix1D tmp = x.like();
	for (int i = 0; i < size; i++) {
		double sum = 0;
		for (int j = 0; j <= i; j++) {
			sum += A.getQuick(i,j) * x.getQuick(j);
		}
		for (int j = i + 1; j < size; j++) {
			sum += A.getQuick(j,i) * x.getQuick(j);
		}
		tmp.setQuick(i, alpha * sum + beta * y.getQuick(i));
	}
	y.assign(tmp);
}
 

/**
 * A matrix <tt>A</tt> is <i>non-negative</i> if <tt>A[i,j] &gt;= 0</tt> holds for all cells.
 * <p>
 * Note: Ignores tolerance.
 */
public boolean isNonNegative(DoubleMatrix2D 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;
}
 

public void dtrmv(boolean isUpperTriangular, boolean transposeA, boolean isUnitTriangular, DoubleMatrix2D A, DoubleMatrix1D x) {
	if (transposeA) {
		A = A.viewDice();
		isUpperTriangular = !isUpperTriangular;
	}
	
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort());
	}
	    
	DoubleMatrix1D b = x.like();
	DoubleMatrix1D y = x.like();
	if (isUnitTriangular) {
		y.assign(1);
	}
	else {
		for (int i = 0; i < size; i++) {
			y.setQuick(i, A.getQuick(i,i));
		}
	}
	
	for (int i = 0; i < size; i++) {
		double sum = 0;
		if (!isUpperTriangular) {
			for (int j = 0; j < i; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);
			}
			sum += y.getQuick(i) * x.getQuick(i);
		}
		else {
			sum += y.getQuick(i) * x.getQuick(i);
			for (int j = i + 1; j < size; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);			}
		}
		b.setQuick(i,sum);
	}
	x.assign(b);
}
 

/**
Decomposes the banded and square matrix <tt>A</tt> into <tt>L</tt> and <tt>U</tt> (in-place).
Upon return <tt>A</tt> is overridden with the result <tt>LU</tt>, such that <tt>L*U = A</tt>.
Currently supports diagonal and tridiagonal matrices, all other cases fall through to {@link #decompose(DoubleMatrix2D)}.
@param semiBandwidth == 1 --> A is diagonal, == 2 --> A is tridiagonal.
@param  A   any matrix.
*/	
public void decompose(DoubleMatrix2D A, int semiBandwidth) {
	if (! algebra.property().isSquare(A) || semiBandwidth<0 || semiBandwidth>2) {
		decompose(A);
		return;
	}
	// setup
	LU = A;
	int m = A.rows();
	int n = A.columns();

	// setup pivot vector
	if (this.piv==null || this.piv.length != m) this.piv = new int[m];
	for (int i = m; --i >= 0; ) piv[i] = i;
	pivsign = 1;

	if (m*n == 0) {
		setLU(A);
		return; // nothing to do
	}
	
	//if (semiBandwidth == 1) { // A is diagonal; nothing to do
	if (semiBandwidth == 2) { // A is tridiagonal
		// currently no pivoting !
		if (n>1) A.setQuick(1,0, A.getQuick(1,0) / A.getQuick(0,0));

		for (int i=1; i<n; i++) {
			double ei = A.getQuick(i,i) - A.getQuick(i,i-1) * A.getQuick(i-1,i);
			A.setQuick(i,i, ei);
			if (i<n-1) A.setQuick(i+1,i, A.getQuick(i+1,i) / ei);
		}
	}
	setLU(A);
}
 

/**
 * A matrix <tt>A</tt> is <i>positive</i> if <tt>A[i,j] &gt; 0</tt> holds for all cells.
 * <p>
 * Note: Ignores tolerance.
 */
public boolean isPositive(DoubleMatrix2D 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;
}
 

/**
 * Returns the sum of the diagonal elements of matrix <tt>A</tt>; <tt>Sum(A[i,i])</tt>.
 */
public double trace(DoubleMatrix2D A) {
	double sum = 0;
	for (int i=Math.min(A.rows(),A.columns()); --i >= 0;) {
		sum += A.getQuick(i,i);
	}
	return sum;
}
 

/**
 * Modifies A to hold its inverse.
 * @param x the first vector.
 * @param y the second vector.
 * @return isNonSingular.
 * @throws IllegalArgumentException if <tt>x.size() != y.size()</tt>.
 */
public static boolean inverse(DoubleMatrix2D A) {
	Property.DEFAULT.checkSquare(A);
	boolean isNonSingular = true;
	for (int i=A.rows(); --i >= 0;) {
		double v = A.getQuick(i,i);
		isNonSingular &= (v!=0);
		A.setQuick(i,i, 1/v);
	}
	return isNonSingular;
}
 

/**
 * A matrix <tt>A</tt> is <i>non-negative</i> if <tt>A[i,j] &gt;= 0</tt> holds for all cells.
 * <p>
 * Note: Ignores tolerance.
 */
public boolean isNonNegative(DoubleMatrix2D 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;
}
 

/**
Decomposes the banded and square matrix <tt>A</tt> into <tt>L</tt> and <tt>U</tt> (in-place).
Upon return <tt>A</tt> is overridden with the result <tt>LU</tt>, such that <tt>L*U = A</tt>.
Currently supports diagonal and tridiagonal matrices, all other cases fall through to {@link #decompose(DoubleMatrix2D)}.
@param semiBandwidth == 1 --> A is diagonal, == 2 --> A is tridiagonal.
@param  A   any matrix.
*/	
public void decompose(DoubleMatrix2D A, int semiBandwidth) {
	if (! algebra.property().isSquare(A) || semiBandwidth<0 || semiBandwidth>2) {
		decompose(A);
		return;
	}
	// setup
	LU = A;
	int m = A.rows();
	int n = A.columns();

	// setup pivot vector
	if (this.piv==null || this.piv.length != m) this.piv = new int[m];
	for (int i = m; --i >= 0; ) piv[i] = i;
	pivsign = 1;

	if (m*n == 0) {
		setLU(A);
		return; // nothing to do
	}
	
	//if (semiBandwidth == 1) { // A is diagonal; nothing to do
	if (semiBandwidth == 2) { // A is tridiagonal
		// currently no pivoting !
		if (n>1) A.setQuick(1,0, A.getQuick(1,0) / A.getQuick(0,0));

		for (int i=1; i<n; i++) {
			double ei = A.getQuick(i,i) - A.getQuick(i,i-1) * A.getQuick(i-1,i);
			A.setQuick(i,i, ei);
			if (i<n-1) A.setQuick(i+1,i, A.getQuick(i+1,i) / ei);
		}
	}
	setLU(A);
}
 

/**
 * Returns the sum of the diagonal elements of matrix <tt>A</tt>; <tt>Sum(A[i,i])</tt>.
 */
public double trace(DoubleMatrix2D A) {
	double sum = 0;
	for (int i=Math.min(A.rows(),A.columns()); --i >= 0;) {
		sum += A.getQuick(i,i);
	}
	return sum;
}
 

/**
Constructs and returns a new eigenvalue decomposition object; 
The decomposed matrices can be retrieved via instance methods of the returned decomposition object.
Checks for symmetry, then constructs the eigenvalue decomposition.
@param A    A square matrix.
@return     A decomposition object to access <tt>D</tt> and <tt>V</tt>.
@throws IllegalArgumentException if <tt>A</tt> is not square.
*/
public EigenvalueDecomposition(DoubleMatrix2D A) {
	Property.DEFAULT.checkSquare(A);
	
	n = A.columns();
	V = new double[n][n];
	d = new double[n];
	e = new double[n];
	
	issymmetric = Property.DEFAULT.isSymmetric(A);
	
	if (issymmetric) {
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				V[i][j] = A.getQuick(i,j);
			}
		}
	
		// Tridiagonalize.
		tred2();
		
		// Diagonalize.
		tql2();
		
	} 
	else {
		H = new double[n][n];
		ort = new double[n];
		 
		for (int j = 0; j < n; j++) {
			for (int i = 0; i < n; i++) {
				H[i][j] = A.getQuick(i,j);
			}
		}
		
		// Reduce to Hessenberg form.
		orthes();
		
		// Reduce Hessenberg to real Schur form.
		hqr2();
	}
}
 

/**
Constructs and returns a new eigenvalue decomposition object; 
The decomposed matrices can be retrieved via instance methods of the returned decomposition object.
Checks for symmetry, then constructs the eigenvalue decomposition.
@param A    A square matrix.
@return     A decomposition object to access <tt>D</tt> and <tt>V</tt>.
@throws IllegalArgumentException if <tt>A</tt> is not square.
*/
public EigenvalueDecomposition(DoubleMatrix2D A) {
	Property.DEFAULT.checkSquare(A);
	
	n = A.columns();
	V = new double[n][n];
	d = new double[n];
	e = new double[n];
	
	issymmetric = Property.DEFAULT.isSymmetric(A);
	
	if (issymmetric) {
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				V[i][j] = A.getQuick(i,j);
			}
		}
	
		// Tridiagonalize.
		tred2();
		
		// Diagonalize.
		tql2();
		
	} 
	else {
		H = new double[n][n];
		ort = new double[n];
		 
		for (int j = 0; j < n; j++) {
			for (int i = 0; i < n; i++) {
				H[i][j] = A.getQuick(i,j);
			}
		}
		
		// Reduce to Hessenberg form.
		orthes();
		
		// Reduce Hessenberg to real Schur form.
		hqr2();
	}
}