cern.colt.matrix.DoubleMatrix2D Java Examples

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

protected void run(final DoubleMatrix2D[] blocksA, final DoubleMatrix2D[] blocksB, final double[] results, final Matrix2DMatrix2DFunction function) {
	final RecursiveAction[] subTasks = new RecursiveAction[blocksA.length];
	for (int i=0; i<blocksA.length; i++) {
		final int k = i;
		subTasks[i] = new RecursiveAction() { 
			@Override
			protected void compute() {
				double result = function.apply(blocksA[k],blocksB != null ? blocksB[k] : null);
				if (results!=null) results[k] = result; 
				//System.out.print("."); 
			}
		};
	}

	// run tasks and wait for completion
		this.taskGroup.invoke(
			new RecursiveAction() {
				@Override
				protected void compute() {
					invokeAll(subTasks);	
				}
			}
                 );
}
 

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

public PepPlaneLinModel(Residue res1, Residue res2){
    //construct from the current conformations
    //of the residue involved in this peptide plane
    
    if(res1.template.name.equalsIgnoreCase("PRO") || res2.template.name.equalsIgnoreCase("PRO"))
        throw new RuntimeException("ERROR: CATS doesn't handle proline");
    
    double CA1[] = res1.getCoordsByAtomName("CA");
    double CA2[] = res2.getCoordsByAtomName("CA");
    double N[] = res2.getCoordsByAtomName("N");
    double C[] = res1.getCoordsByAtomName("C");
    double O[] = res1.getCoordsByAtomName("O");
    double H[] = res2.getCoordsByAtomName("H");
    
    DoubleMatrix2D M = makeAxisMatrix(CA1,N,CA2);
    DoubleMatrix2D vec = DoubleFactory2D.dense.make(3,1);
    
    vec.viewColumn(0).assign(VectorAlgebra.subtract(C, CA1));
    CCoeffs = Algebra.DEFAULT.solve(M, vec).viewColumn(0);
    
    vec.viewColumn(0).assign(VectorAlgebra.subtract(O, CA1));
    OCoeffs = Algebra.DEFAULT.solve(M, vec).viewColumn(0);
    
    vec.viewColumn(0).assign(VectorAlgebra.subtract(H, CA1));
    HCoeffs = Algebra.DEFAULT.solve(M, vec).viewColumn(0);
}
 

/**
 */
public static void doubleTest34() {
	double[][] data = 
	{ 
		{ 3, 0, 0, 0 },
		{ 0, 4, 2, 0 },
		{ 0, 0, 0, 0 },
		{ 0, 0, 0, 0 },
	};
	
	DoubleMatrix2D A = new DenseDoubleMatrix2D(data);
	Property.DEFAULT.generateNonSingular(A);
	DoubleMatrix2D inv = Algebra.DEFAULT.inverse(A);


	System.out.println("\n\n\n"+A);
	System.out.println("\n"+inv);
	DoubleMatrix2D B = A.zMult(inv,null);
	System.out.println(B);
	if (!(B.equals(DoubleFactory2D.dense.identity(A.rows)))) {
		throw new InternalError();
	}
}
 

public DoubleMatrix2D forEachNonZero(final cern.colt.function.IntIntDoubleFunction function) {
	if (this.isNoView) {
		this.elements.forEachPair(
			new cern.colt.function.IntDoubleProcedure() {
				public boolean apply(int key, double value) {
					int i = key/columns;
					int j = key%columns;
					double r = function.apply(i,j,value);
					if (r!=value) elements.put(key,r);
					return true;
				}
			}
		);
	}
	else {
		super.forEachNonZero(function);
	}
	return this;
}
 

/**
 * Replaces all cell values of the receiver with the values of another matrix.
 * Both matrices must have the same number of rows and columns.
 * If both matrices share the same cells (as is the case if they are views derived from the same matrix) and intersect in an ambiguous way, then replaces <i>as if</i> using an intermediate auxiliary deep copy of <tt>other</tt>.
 *
 * @param     source   the source matrix to copy from (may be identical to the receiver).
 * @return <tt>this</tt> (for convenience only).
 * @throws	IllegalArgumentException if <tt>columns() != source.columns() || rows() != source.rows()</tt>
 */
public DoubleMatrix2D assign(DoubleMatrix2D source) {
	// overriden for performance only
	if (! (source instanceof SparseDoubleMatrix2D)) {
		return super.assign(source);
	}
	SparseDoubleMatrix2D other = (SparseDoubleMatrix2D) source;
	if (other==this) return this; // nothing to do
	checkShape(other);
	
	if (this.isNoView && other.isNoView) { // quickest
		this.elements.assign(other.elements);
		return this;
	}
	return super.assign(source);
}
 

/**
 * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
 * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
 * <tt>! (Math.abs(value - A[row,col]) > tolerance())</tt> holds for all coordinates.
 * @param   A   the first matrix to compare.
 * @param   value   the value to compare against.
 * @return  <tt>true</tt> if the matrix is equal to the value;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix2D A, double value) {
	if (A==null) return false;
	int rows = A.rows();
	int columns = A.columns();

	double epsilon = tolerance();
	for (int row=rows; --row >= 0;) {
		for (int column=columns; --column >= 0;) {
			//if (!(A.getQuick(row,column) == value)) return false;
			//if (Math.abs(value - A.getQuick(row,column)) > epsilon) return false;
			double x = A.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;
}
 

/**
 * Copies the rows of the indicated columns into a new sub matrix.
 * <tt>sub[0..rowTo-rowFrom,0..columnIndexes.length-1] = A[rowFrom..rowTo,columnIndexes(:)]</tt>;
 * The returned matrix is <i>not backed</i> by this matrix, so changes in the returned matrix are <i>not reflected</i> in this matrix, and vice-versa.
 *
 * @param   A   the source matrix to copy from.
 * @param   rowFrom the index of the first row to copy (inclusive).
 * @param   rowTo the index of the last row to copy (inclusive).
 * @param   columnIndexes the indexes of the columns to copy. May be unsorted.
 * @return  a new sub matrix; with <tt>sub.rows()==rowTo-rowFrom+1; sub.columns()==columnIndexes.length</tt>.
 * @throws	IndexOutOfBoundsException if <tt>rowFrom<0 || rowTo-rowFrom+1<0 || rowTo+1>matrix.rows() || for any col=columnIndexes[i]: col < 0 || col >= matrix.columns()</tt>.
 */
private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int rowFrom, int rowTo, int[] columnIndexes) {
	if (rowTo-rowFrom >= A.rows()) throw new IndexOutOfBoundsException("Too many rows");
	int height = rowTo-rowFrom+1;
	int columns = A.columns();
	A = A.viewPart(rowFrom,0,height,columns);
	DoubleMatrix2D sub = A.like(height, columnIndexes.length);
	
	for (int c = columnIndexes.length; --c >= 0; ) {
		int column = columnIndexes[c];
		if (column < 0 || column >= columns)
			throw new IndexOutOfBoundsException("Illegal Index");
		sub.viewColumn(c).assign(A.viewColumn(column));
	}
	return sub;
}
 

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

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

/**
 * Replaces all cell values of the receiver with the values of another matrix.
 * Both matrices must have the same number of rows and columns.
 * If both matrices share the same cells (as is the case if they are views derived from the same matrix) and intersect in an ambiguous way, then replaces <i>as if</i> using an intermediate auxiliary deep copy of <tt>other</tt>.
 *
 * @param     source   the source matrix to copy from (may be identical to the receiver).
 * @return <tt>this</tt> (for convenience only).
 * @throws	IllegalArgumentException if <tt>columns() != source.columns() || rows() != source.rows()</tt>
 */
public DoubleMatrix2D assign(DoubleMatrix2D source) {
	// overriden for performance only
	if (! (source instanceof SparseDoubleMatrix2D)) {
		return super.assign(source);
	}
	SparseDoubleMatrix2D other = (SparseDoubleMatrix2D) source;
	if (other==this) return this; // nothing to do
	checkShape(other);
	
	if (this.isNoView && other.isNoView) { // quickest
		this.elements.assign(other.elements);
		return this;
	}
	return super.assign(source);
}
 

/**
 */
public static void doubleTest3() {
int rows = 4;
int columns = 5; // make a 4*5 matrix
DoubleMatrix2D master = new DenseDoubleMatrix2D(rows,columns);
System.out.println(master);
master.assign(1); // set all cells to 1
System.out.println("\n"+master);
master.viewPart(2,0,2,3).assign(2); // set [2,1] .. [3,3] to 2
System.out.println("\n"+master);

DoubleMatrix2D flip1 = master.viewColumnFlip();
System.out.println("flip around columns="+flip1);
DoubleMatrix2D flip2 = flip1.viewRowFlip();
System.out.println("further flip around rows="+flip2);

flip2.viewPart(0,0,2,2).assign(3);
System.out.println("master replaced"+master);
System.out.println("flip1 replaced"+flip1);
System.out.println("flip2 replaced"+flip2);


/*
DoubleMatrix2D copyPart = master.copyPart(2,1,2,3);
copyPart.assign(3); // modify an independent copy
copyPart.set(0,0,4);
System.out.println("\n"+copyPart); // has changed
System.out.println("\n"+master); // master has not changed

DoubleMatrix2D view1 = master.viewPart(0,3,4,2); // [0,3] .. [3,4]
DoubleMatrix2D view2 = view1.viewPart(0,0,4,1); // a view from a view 
System.out.println("\n"+view1);
System.out.println("\n"+view2);
*/
}
 

/**
 * Copies the columns of the indicated rows into a new sub matrix.
 * <tt>sub[0..rowIndexes.length-1,0..columnTo-columnFrom] = A[rowIndexes(:),columnFrom..columnTo]</tt>;
 * The returned matrix is <i>not backed</i> by this matrix, so changes in the returned matrix are <i>not reflected</i> in this matrix, and vice-versa.
 *
 * @param   A   the source matrix to copy from.
 * @param   rowIndexes the indexes of the rows to copy. May be unsorted.
 * @param   columnFrom the index of the first column to copy (inclusive).
 * @param   columnTo the index of the last column to copy (inclusive).
 * @return  a new sub matrix; with <tt>sub.rows()==rowIndexes.length; sub.columns()==columnTo-columnFrom+1</tt>.
 * @throws	IndexOutOfBoundsException if <tt>columnFrom<0 || columnTo-columnFrom+1<0 || columnTo+1>matrix.columns() || for any row=rowIndexes[i]: row < 0 || row >= matrix.rows()</tt>.
 */
private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int[] rowIndexes, int columnFrom, int columnTo) {
	int width = columnTo-columnFrom+1;
	int rows = A.rows();
	A = A.viewPart(0,columnFrom,rows,width);
	DoubleMatrix2D sub = A.like(rowIndexes.length, width);
	
	for (int r = rowIndexes.length; --r >= 0; ) {
		int row = rowIndexes[r];
		if (row < 0 || row >= rows) 
			throw new IndexOutOfBoundsException("Illegal Index");
		sub.viewRow(r).assign(A.viewRow(row));
	}
	return sub;
}
 

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

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

/**
 */
public static void doubleTest29(int size,DoubleFactory2D f) {
	
	DoubleMatrix2D x = new DenseDoubleMatrix2D(size,size).assign(0.5);
	DoubleMatrix2D matrix = f.sample(size,size,0.5,0.001);
	
	cern.colt.Timer timer = new cern.colt.Timer().start();
	DoubleMatrix2D res = matrix.zMult(x,null);
	timer.stop().display();
	
	//System.out.println(res);
}
 

/**
Modifies the given matrix <tt>A</tt> such that it's rows are permuted as specified; Useful for pivoting.
Row <tt>A[i]</tt> will go into row <tt>A[indexes[i]]</tt>.
<p>
<b>Example:</b>
<pre>
Reordering
[A,B,C,D,E] with indexes [0,4,2,3,1] yields 
[A,E,C,D,B]
In other words A[0]<--A[0], A[1]<--A[4], A[2]<--A[2], A[3]<--A[3], A[4]<--A[1].

Reordering
[A,B,C,D,E] with indexes [0,4,1,2,3] yields 
[A,E,B,C,D]
In other words A[0]<--A[0], A[1]<--A[4], A[2]<--A[1], A[3]<--A[2], A[4]<--A[3].
</pre>

@param   A   the matrix to permute.
@param   indexes the permutation indexes, must satisfy <tt>indexes.length==A.rows() && indexes[i] >= 0 && indexes[i] < A.rows()</tt>;
@param   work the working storage, must satisfy <tt>work.length >= A.rows()</tt>; set <tt>work==null</tt> if you don't care about performance.
@return the modified <tt>A</tt> (for convenience only).
@throws	IndexOutOfBoundsException if <tt>indexes.length != A.rows()</tt>.
*/
public DoubleMatrix2D permuteRows(final DoubleMatrix2D A, int[] indexes, int[] work) {
	// check validity
	int size = A.rows();
	if (indexes.length != size) throw new IndexOutOfBoundsException("invalid permutation");

	/*
	int i=size;
	int a;
	while (--i >= 0 && (a=indexes[i])==i) if (a < 0 || a >= size) throw new IndexOutOfBoundsException("invalid permutation");
	if (i<0) return; // nothing to permute
	*/

	int columns = A.columns();
	if (columns < size/10) { // quicker
		double[] doubleWork = new double[size];
		for (int j=A.columns(); --j >= 0; ) permute(A.viewColumn(j), indexes, doubleWork);
		return A;
	}

	cern.colt.Swapper swapper = new cern.colt.Swapper() {
		public void swap(int a, int b) {
			A.viewRow(a).swap(A.viewRow(b));
		}
	};

	cern.colt.GenericPermuting.permute(indexes, swapper, work, null);
	return A;
}
 

/**
 * Sets all cells to the state specified by <tt>value</tt>.
 * @param    value the value to be filled into the cells.
 * @return <tt>this</tt> (for convenience only).
 */
public DoubleMatrix2D assign(double value) {
	// overriden for performance only
	if (this.isNoView && value==0) this.elements.clear();
	else super.assign(value);
	return this;
}
 

public DoubleMatrix2D assign(final cern.colt.function.DoubleFunction function) {
	if (function instanceof cern.jet.math.Mult) { // x[i] = mult*x[i]
		final double alpha = ((cern.jet.math.Mult) function).multiplicator;
		if (alpha==1) return this;
		if (alpha==0) return assign(0);
		if (alpha!=alpha) return assign(alpha); // the funny definition of isNaN(). This should better not happen.

		/*
		double[] vals = values.elements();
		for (int j=values.size(); --j >= 0; ) {
			vals[j] *= alpha;
		}
		*/

		forEachNonZero(
			new cern.colt.function.IntIntDoubleFunction() {
				public double apply(int i, int j, double value) {
					return function.apply(value);
				}
			}
		);
	}
	else {
		super.assign(function);
	}
	return this;
}
 

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

public double[] getCartesianCoords(double[] polars) {
	if (polars.length==0) { return new double[0]; }
	
	EigenvalueDecomposition evd = new EigenvalueDecomposition(A);
	DoubleMatrix2D Q = evd.getV();
	DoubleMatrix2D L = evd.getD();
	DoubleMatrix2D qT = Q.viewDice().copy();
	Algebra alg = new Algebra();		
	
	int n = polars.length;		
	double radius = polars[0];
	double[] phis = new double[n-1];
	for (int i=1; i<n; i++) { phis[i-1] = polars[i]; }
	
	double[] cartCoords = new double[n];
	for (int i=0; i<n; i++) {
		double prod = 1;
		for (int j=0; j<i; j++) { prod = prod * Math.sin(phis[j]); }
		if (i<n-1) { prod = prod * Math.cos(phis[i]); } 			
		cartCoords[i] = radius*prod;
	}
	
	// transform cartesian coordinates back!
	
	double[] u = new double[cartCoords.length];
	for (int i=0; i<u.length; i++) {
		u[i] = cartCoords[i]*Math.sqrt(L.get(i, i));
	}
	double[] s = alg.mult(Q, DoubleFactory1D.dense.make(u)).toArray();
	double[] x = new double[s.length];
	for (int i=0; i<x.length; i++) { x[i] = s[i] + c.get(i); }
	
	return x;
}
 

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

/** 
Returns whether the matrix is nonsingular.
@return true if <tt>matrix</tt> is nonsingular; false otherwise.
*/
protected boolean isNonsingular(DoubleMatrix2D matrix) {
	int m = matrix.rows();
	int n = matrix.columns();
	double epsilon = algebra.property().tolerance(); // consider numerical instability
	for (int j = Math.min(n,m); --j >= 0;) {
		//if (matrix.getQuick(j,j) == 0) return false;
		if (Math.abs(matrix.getQuick(j,j)) <= epsilon) return false;
	}
	return true;
}
 

public void dger(double alpha, DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix2D A) {
	cern.jet.math.PlusMult fun = cern.jet.math.PlusMult.plusMult(0);
	for (int i=A.rows(); --i >= 0; ) {
		fun.multiplicator = alpha * x.getQuick(i);
 		A.viewRow(i).assign(y,fun);
		
	}
}
 

/**
 * Sets all cells to the state specified by <tt>value</tt>.
 * @param    value the value to be filled into the cells.
 * @return <tt>this</tt> (for convenience only).
 */
public DoubleMatrix2D assign(double value) {
	// overriden for performance only
	if (value==0) {
		indexes.clear();
		values.clear();
		for (int i=starts.length; --i >= 0; ) starts[i] = 0;
	}
	else super.assign(value);
	return this;
}
 

/**
 * A matrix <tt>A</tt> is <i>tridiagonal</i> if <tt>A[i,j]==0</tt> whenever <tt>Math.abs(i-j) > 1</tt>.
 * Matrix may but need not be square.
 */
public boolean isTridiagonal(DoubleMatrix2D A) {
	double 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 (A.getQuick(row,column) != 0) return false;
				if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
			}
		}
	}
	return true;
}
 

public void dcopy(DoubleMatrix2D A, DoubleMatrix2D B) {
	B.assign(A);
}
 

/**
Constructs and returns a new 2-dimensional <i>slice view</i> representing the slices and columns of the given row.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
<p>
To obtain a slice view on subranges, construct a sub-ranging view (<tt>view().part(...)</tt>), then apply this method to the sub-range view.
To obtain 1-dimensional views, apply this method, then apply another slice view (methods <tt>viewColumn</tt>, <tt>viewRow</tt>) on the intermediate 2-dimensional view.
To obtain 1-dimensional views on subranges, apply both steps.

@param row the index of the row to fix.
@return a new 2-dimensional slice view.
@throws IndexOutOfBoundsException if <tt>row < 0 || row >= row()</tt>.
@see #viewSlice(int)
@see #viewColumn(int)
*/
public DoubleMatrix2D viewRow(int row) {
	checkRow(row);
	
	int viewRows = this.slices;
	int viewColumns = this.columns;
	
	int viewRowZero = sliceZero;
	int viewColumnZero = columnZero;
	int viewOffset = this.offset + _rowOffset(_rowRank(row));

	int viewRowStride = this.sliceStride;
	int viewColumnStride = this.columnStride;

	int[] viewRowOffsets = this.sliceOffsets;
	int[] viewColumnOffsets = this.columnOffsets;

	return new SelectedDenseDoubleMatrix2D(viewRows,viewColumns,this.elements,viewRowZero,viewColumnZero,viewRowStride,viewColumnStride,viewRowOffsets,viewColumnOffsets,viewOffset);
}
 

/**
 * Returns the matrix's fraction of non-zero cells; <tt>A.cardinality() / A.size()</tt>.
 */
public double density(DoubleMatrix2D A) {
	return A.cardinality() / (double) A.size();
}