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

/**
 * 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;	
}
 

DoubleMatrix2D homMatrix(){
    //Return M, where this ellipse is represented as [x^T 1]M[x;1]=0
    int n = c.size();
    DoubleMatrix2D ans = DoubleFactory2D.dense.make(n+1,n+1);
    DoubleMatrix1D Qnc = A.zMult(nc, null);
    for(int a=0; a<n; a++){
        for(int b=0; b<n; b++)
            ans.setQuick( a, b, A.get(a,b) );
        
        ans.setQuick(a, n, Qnc.get(a));
        ans.setQuick(n, a, Qnc.get(a));
    }
    
    ans.setQuick(n, n, Qnc.zDotProduct(nc)-1 );
    
    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;
}
 

/** 
Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>.
@return <tt>Q</tt>
*/
public DoubleMatrix2D getQ () {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	DoubleMatrix2D Q = QR.like();
	//double[][] Q = X.getArray();
	for (int k = n-1; k >= 0; k--) {
		DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
		Q.setQuick(k,k, 1);
		for (int j = k; j < n; j++) {
			if (QR.getQuick(k,k) != 0) {
				DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k);
				double s = QRcolk.zDotProduct(Qcolj);
				s = -s / QR.getQuick(k,k);
				Qcolj.assign(QRcolk, F.plusMult(s));
			}
		}
	}
	return Q;
}
 

/** 
Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>.
@return <tt>Q</tt>
*/
public DoubleMatrix2D getQ () {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	DoubleMatrix2D Q = QR.like();
	//double[][] Q = X.getArray();
	for (int k = n-1; k >= 0; k--) {
		DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
		Q.setQuick(k,k, 1);
		for (int j = k; j < n; j++) {
			if (QR.getQuick(k,k) != 0) {
				DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k);
				double s = QRcolk.zDotProduct(Qcolj);
				s = -s / QR.getQuick(k,k);
				Qcolj.assign(QRcolk, F.plusMult(s));
			}
		}
	}
	return Q;
}
 

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;
}
 

/** 
Least squares solution of <tt>A*X = B</tt>; <tt>returns X</tt>.
@param B    A matrix with as many rows as <tt>A</tt> and any number of columns.
@return     <tt>X</tt> that minimizes the two norm of <tt>Q*R*X - B</tt>.
@exception  IllegalArgumentException  if <tt>B.rows() != A.rows()</tt>.
@exception  IllegalArgumentException  if <tt>!this.hasFullRank()</tt> (<tt>A</tt> is rank deficient).
*/
public DoubleMatrix2D solve(DoubleMatrix2D B) {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	if (B.rows() != m) {
		throw new IllegalArgumentException("Matrix row dimensions must agree.");
	}
	if (!this.hasFullRank()) {
		throw new IllegalArgumentException("Matrix is rank deficient.");
	}
	
	// Copy right hand side
	int nx = B.columns();
	DoubleMatrix2D X = B.copy();
	
	// Compute Y = transpose(Q)*B
	for (int k = 0; k < n; k++) {
		for (int j = 0; j < nx; j++) {
			double s = 0.0; 
			for (int i = k; i < m; i++) {
				s += QR.getQuick(i,k)*X.getQuick(i,j);
			}
			s = -s / QR.getQuick(k,k);
			for (int i = k; i < m; i++) {
				X.setQuick(i,j, X.getQuick(i,j) + s*QR.getQuick(i,k));
			}
		}
	}
	// Solve R*X = Y;
	for (int k = n-1; k >= 0; k--) {
		for (int j = 0; j < nx; j++) {
			X.setQuick(k,j, X.getQuick(k,j) / Rdiag.getQuick(k));
		}
		for (int i = 0; i < k; i++) {
			for (int j = 0; j < nx; j++) {
				X.setQuick(i,j, X.getQuick(i,j) - X.getQuick(k,j)*QR.getQuick(i,k));
			}
		}
	}
	return X.viewPart(0,0,n,nx);
}
 

/**
Modifies the matrix to be an upper triangular matrix.
@return <tt>A</tt> (for convenience only).
@see #triangulateLower(DoubleMatrix2D)
*/
protected DoubleMatrix2D upperTriangular(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	int min = Math.min(rows,columns);
	for (int r = min; --r >= 0; ) {
		for (int c = min; --c >= 0; ) {
			if (r > c) A.setQuick(r,c, 0);
		}
	}
	if (columns<rows) A.viewPart(min,0,rows-min,columns).assign(0);

	return A;
}
 

/**
Modifies the given matrix square matrix <tt>A</tt> 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 <tt>!isSquare(A)</tt>.
*/
public void generateNonSingular(DoubleMatrix2D A) {
	checkSquare(A);
	cern.jet.math.Functions F = cern.jet.math.Functions.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; ) {
		double rowSum = A.viewRow(i).aggregate(F.plus,F.abs);
		double colSum = A.viewColumn(i).aggregate(F.plus,F.abs);
		A.setQuick(i,i, Math.max(rowSum,colSum) + i+1);
	}
}
 

/**
Modifies the matrix to be a lower trapezoidal matrix.
@return <tt>A</tt> (for convenience only).
@see #triangulateLower(DoubleMatrix2D)
*/
protected DoubleMatrix2D trapezoidalLower(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	for (int r = rows; --r >= 0; ) {
		for (int c = columns; --c >= 0; ) {
			if (r < c) A.setQuick(r,c, 0);
		}
	}
	return A;
}
 

/** 
Returns the upper triangular factor, <tt>R</tt>.
@return <tt>R</tt>
*/
public DoubleMatrix2D getR() {
	DoubleMatrix2D R = QR.like(n,n);
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			if (i < j) 
				R.setQuick(i,j, QR.getQuick(i,j));
			else if (i == j) 
				R.setQuick(i,j, Rdiag.getQuick(i));
			else 
				R.setQuick(i,j, 0);
		}
	}
	return R;
}
 

/**
 * 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;
}
 

/**
 */
public static void doubleTest27() {

	
System.out.println("\n\n");
System.out.println("initializing...");

int rows=51;
int columns=10;
double[][] trainingSet = new double[columns][rows];
for (int i=columns; --i >= 0; ) trainingSet[i][i]=2.0;

int patternIndex        = 0;
int unitIndex           = 0;

DoubleMatrix2D patternMatrix       = null;
DoubleMatrix2D transposeMatrix     = null;
DoubleMatrix2D QMatrix             = null;
DoubleMatrix2D inverseQMatrix      = null;
DoubleMatrix2D pseudoInverseMatrix = null;
DoubleMatrix2D weightMatrix        = null;

// form a matrix with the columns as training vectors
patternMatrix = DoubleFactory2D.dense.make (rows, columns);

// copy the patterns into the matrix
for (patternIndex=0;patternIndex<columns;patternIndex++) {
	for (unitIndex=0;unitIndex<rows;unitIndex++) {
		patternMatrix.setQuick (unitIndex, patternIndex, trainingSet[patternIndex][unitIndex]);
	}
}

transposeMatrix     = Algebra.DEFAULT.transpose (patternMatrix);
QMatrix             = Algebra.DEFAULT.mult (transposeMatrix,patternMatrix);
inverseQMatrix      = Algebra.DEFAULT.inverse (QMatrix);
pseudoInverseMatrix = Algebra.DEFAULT.mult (inverseQMatrix,transposeMatrix);
weightMatrix        = Algebra.DEFAULT.mult (patternMatrix,pseudoInverseMatrix);
System.out.println("done.");
	
}
 

/**
Modifies the given matrix square matrix <tt>A</tt> 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 <tt>!isSquare(A)</tt>.
*/
public void generateNonSingular(DoubleMatrix2D A) {
	checkSquare(A);
	cern.jet.math.Functions F = cern.jet.math.Functions.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; ) {
		double rowSum = A.viewRow(i).aggregate(F.plus,F.abs);
		double colSum = A.viewColumn(i).aggregate(F.plus,F.abs);
		A.setQuick(i,i, Math.max(rowSum,colSum) + i+1);
	}
}
 

/**
 */
public static void doubleTest27() {

	
System.out.println("\n\n");
System.out.println("initializing...");

int rows=51;
int columns=10;
double[][] trainingSet = new double[columns][rows];
for (int i=columns; --i >= 0; ) trainingSet[i][i]=2.0;

int patternIndex        = 0;
int unitIndex           = 0;

DoubleMatrix2D patternMatrix       = null;
DoubleMatrix2D transposeMatrix     = null;
DoubleMatrix2D QMatrix             = null;
DoubleMatrix2D inverseQMatrix      = null;
DoubleMatrix2D pseudoInverseMatrix = null;
DoubleMatrix2D weightMatrix        = null;

// form a matrix with the columns as training vectors
patternMatrix = DoubleFactory2D.dense.make (rows, columns);

// copy the patterns into the matrix
for (patternIndex=0;patternIndex<columns;patternIndex++) {
	for (unitIndex=0;unitIndex<rows;unitIndex++) {
		patternMatrix.setQuick (unitIndex, patternIndex, trainingSet[patternIndex][unitIndex]);
	}
}

transposeMatrix     = Algebra.DEFAULT.transpose (patternMatrix);
QMatrix             = Algebra.DEFAULT.mult (transposeMatrix,patternMatrix);
inverseQMatrix      = Algebra.DEFAULT.inverse (QMatrix);
pseudoInverseMatrix = Algebra.DEFAULT.mult (inverseQMatrix,transposeMatrix);
weightMatrix        = Algebra.DEFAULT.mult (patternMatrix,pseudoInverseMatrix);
System.out.println("done.");
	
}
 

/**
 * 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;
}
 

/**
Modifies the matrix to be a lower trapezoidal matrix.
@return <tt>A</tt> (for convenience only).
@see #triangulateLower(DoubleMatrix2D)
*/
protected DoubleMatrix2D trapezoidalLower(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	for (int r = rows; --r >= 0; ) {
		for (int c = columns; --c >= 0; ) {
			if (r < c) A.setQuick(r,c, 0);
		}
	}
	return A;
}
 

/**
 */
public static void doubleTest15(int size, int runs) {
System.out.println("\n\n");
double[][] values = 
{
	{ 0, 5, 9 },
	{ 2, 6, 10 },
	{ 3, 7, 11 }
};

//DoubleMatrix2D A = Factory2D.make(values);
DoubleMatrix2D A = Factory2D.make(size,size);
double value = 5;
for (int i=size; --i >= 0; ) {
	A.setQuick(i,i, value);
}
A.viewRow(0).assign(value);

//DoubleMatrix2D A = Factory2D.makeIdentity(size,size);


//DoubleMatrix2D A = Factory2D.makeAscending(size,size).assign(new cern.jet.random.engine.MersenneTwister());
cern.colt.Timer timer = new cern.colt.Timer().start();
DoubleMatrix2D inv = null;
for (int run=0; run<runs; run++) {
	 inv = LinearAlgebra.inverse(A);
}
timer.stop().display();

/*
timer.reset().start();
for (int run=0; run<runs; run++) {
	new Jama.Matrix(A.toArray()).inverse();
}
timer.stop().display();
*/
//System.out.println("A="+A);
//System.out.println("inverse(A)="+inv);
//System.out.println("formatted inverse(A)="+ new Jama.Matrix(inv.toArray()));

/*
-1.0000000000000018, 2.000000000000007, -1.0000000000000047
2.000000000000007, -6.750000000000024, 4.500000000000016
-1.000000000000004, 3.7500000000000133, -2.500000000000009
*/
}
 

/**
Modifies the matrix to be a lower triangular matrix.
<p>
<b>Examples:</b> 
<table border="0">
  <tr nowrap> 
	<td valign="top">3 x 5 matrix:<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9 </td>
	<td align="center">triang.Upper<br>
	  ==></td>
	<td valign="top">3 x 5 matrix:<br>
	  9, 9, 9, 9, 9<br>
	  0, 9, 9, 9, 9<br>
	  0, 0, 9, 9, 9</td>
  </tr>
  <tr nowrap> 
	<td valign="top">5 x 3 matrix:<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9 </td>
	<td align="center">triang.Upper<br>
	  ==></td>
	<td valign="top">5 x 3 matrix:<br>
	  9, 9, 9<br>
	  0, 9, 9<br>
	  0, 0, 9<br>
	  0, 0, 0<br>
	  0, 0, 0</td>
  </tr>
  <tr nowrap> 
	<td valign="top">3 x 5 matrix:<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9<br>
	  9, 9, 9, 9, 9 </td>
	<td align="center">triang.Lower<br>
	  ==></td>
	<td valign="top">3 x 5 matrix:<br>
	  1, 0, 0, 0, 0<br>
	  9, 1, 0, 0, 0<br>
	  9, 9, 1, 0, 0</td>
  </tr>
  <tr nowrap> 
	<td valign="top">5 x 3 matrix:<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9<br>
	  9, 9, 9 </td>
	<td align="center">triang.Lower<br>
	  ==></td>
	<td valign="top">5 x 3 matrix:<br>
	  1, 0, 0<br>
	  9, 1, 0<br>
	  9, 9, 1<br>
	  9, 9, 9<br>
	  9, 9, 9</td>
  </tr>
</table>

@return <tt>A</tt> (for convenience only).
@see #triangulateUpper(DoubleMatrix2D)
*/
protected DoubleMatrix2D lowerTriangular(DoubleMatrix2D A) {
	int rows = A.rows();
	int columns = A.columns();
	int min = Math.min(rows,columns);
	for (int r = min; --r >= 0; ) {
		for (int c = min; --c >= 0; ) {
			if (r < c) A.setQuick(r,c, 0);
			else if (r == c) A.setQuick(r,c, 1);
		}
	}
	if (columns>rows) A.viewPart(0,min,rows,columns-min).assign(0);

	return A;
}