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

/**
 * Linear algebraic matrix-vector multiplication; <tt>z = A * y</tt>.
 * <tt>z[i] = alpha*Sum(A[i,j] * y[j]) + beta*z[i], i=0..A.rows()-1, j=0..y.size()-1</tt>.
 * Where <tt>A == this</tt>.
 * @param y the source vector.
 * @param z the vector where results are to be stored.
 * 
 * @throws IllegalArgumentException if <tt>A.columns() != y.size() || A.rows() > z.size())</tt>.
 */
protected void zMult(final DoubleMatrix1D y, final DoubleMatrix1D z, cern.colt.list.IntArrayList nonZeroIndexes, DoubleMatrix1D[] allRows, final double alpha, final double beta) {
	if (columns != y.size() || rows > z.size())	
		throw new IllegalArgumentException("Incompatible args: "+toStringShort()+", "+y.toStringShort()+", "+z.toStringShort());

	z.assign(cern.jet.math.Functions.mult(beta/alpha));
	for (int i = indexes.length; --i >= 0; ) {
		if (indexes[i] != null) {
			for (int k = indexes[i].size(); --k >= 0; ) {
				int j = indexes[i].getQuick(k);
				double value = values[i].getQuick(k);
				z.setQuick(i,z.getQuick(i) + value * y.getQuick(j));
			}
		}
	}
		
	z.assign(cern.jet.math.Functions.mult(alpha));
}
 

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

/**
 * Linear algebraic matrix-vector multiplication; <tt>z = A * y</tt>.
 * <tt>z[i] = alpha*Sum(A[i,j] * y[j]) + beta*z[i], i=0..A.rows()-1, j=0..y.size()-1</tt>.
 * Where <tt>A == this</tt>.
 * @param y the source vector.
 * @param z the vector where results are to be stored.
 * 
 * @throws IllegalArgumentException if <tt>A.columns() != y.size() || A.rows() > z.size())</tt>.
 */
protected void zMult(final DoubleMatrix1D y, final DoubleMatrix1D z, cern.colt.list.IntArrayList nonZeroIndexes, DoubleMatrix1D[] allRows, final double alpha, final double beta) {
	if (columns != y.size() || rows > z.size())	
		throw new IllegalArgumentException("Incompatible args: "+toStringShort()+", "+y.toStringShort()+", "+z.toStringShort());

	z.assign(cern.jet.math.Functions.mult(beta/alpha));
	for (int i = indexes.length; --i >= 0; ) {
		if (indexes[i] != null) {
			for (int k = indexes[i].size(); --k >= 0; ) {
				int j = indexes[i].getQuick(k);
				double value = values[i].getQuick(k);
				z.setQuick(i,z.getQuick(i) + value * y.getQuick(j));
			}
		}
	}
		
	z.assign(cern.jet.math.Functions.mult(alpha));
}
 

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

private static void doRR1(int L, DoubleMatrix1D w, double[][] x, double[] y, double[] c, double lambda) {
    int N = x.length;
    int K = x[0].length;
    
    double[] e = new double[N];
    for (int i = 0; i < N; i++) {
        double pred = 0.0;
        for (int k = 0; k < K; k++) {
            pred += w.getQuick(k) * x[i][k];
        }
        e[i] = y[i] - pred;
    }

    for (int l = 0; l < L; l++) {
        for (int k = 0; k < K; k++) {
            for (int i = 0; i < N; i++) {
                e[i] += w.getQuick(k) * x[i][k];
            }
            double a = 0.0;
            double d = 0.0;
            for (int i = 0; i < N; i++) {
                a += c[i] * x[i][k] * x[i][k];
                d += c[i] * x[i][k] * e[i];
            }
            w.setQuick(k, d / (lambda + a));
            for (int i = 0; i < N; i++) {
                e[i] -= w.getQuick(k) * x[i][k];
            }
        }
    }

}
 

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

/**
Modifies the given vector <tt>A</tt> such that it is permuted as specified; Useful for pivoting.
Cell <tt>A[i]</tt> will go into cell <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 vector to permute.
@param   indexes the permutation indexes, must satisfy <tt>indexes.length==A.size() && indexes[i] >= 0 && indexes[i] < A.size()</tt>;
@param   work the working storage, must satisfy <tt>work.length >= A.size()</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.size()</tt>.
*/
public DoubleMatrix1D permute(DoubleMatrix1D A, int[] indexes, double[] work) {
	// check validity
	int size = A.size();
	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
	*/

	if (work==null || size > work.length) {
		work = A.toArray();
	}
	else {
		A.toArray(work);
	}
	for (int i=size; --i >= 0; ) A.setQuick(i, work[indexes[i]]);
	return A;
}
 

/**
Decomposes 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>.
Uses a "left-looking", dot-product, Crout/Doolittle algorithm.
@param  A   any matrix.
*/	
public void decompose(DoubleMatrix2D A) {
	final int CUT_OFF = 10;
	// 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(LU);
		return; // nothing to do
	}
	
	//precompute and cache some views to avoid regenerating them time and again
	DoubleMatrix1D[] LUrows = new DoubleMatrix1D[m];
	for (int i = 0; i < m; i++) LUrows[i] = LU.viewRow(i);
	
	cern.colt.list.IntArrayList nonZeroIndexes = new cern.colt.list.IntArrayList(); // sparsity
	DoubleMatrix1D LUcolj = LU.viewColumn(0).like();  // blocked column j
	cern.jet.math.Mult multFunction = cern.jet.math.Mult.mult(0); 

	// Outer loop.
	for (int j = 0; j < n; j++) {
		// blocking (make copy of j-th column to localize references)
		LUcolj.assign(LU.viewColumn(j));
		
		// sparsity detection
		int maxCardinality = m/CUT_OFF; // == heuristic depending on speedup
		LUcolj.getNonZeros(nonZeroIndexes,null,maxCardinality);
		int cardinality = nonZeroIndexes.size(); 
		boolean sparse = (cardinality < maxCardinality);

		// Apply previous transformations.
		for (int i = 0; i < m; i++) {
			int kmax = Math.min(i,j);
			double s;
			if (sparse) {
				s = LUrows[i].zDotProduct(LUcolj,0,kmax,nonZeroIndexes);
			}
			else {
				s = LUrows[i].zDotProduct(LUcolj,0,kmax);
			}
			double before = LUcolj.getQuick(i);
			double after = before -s;
			LUcolj.setQuick(i, after); // LUcolj is a copy
			LU.setQuick(i,j, after);   // this is the original
			if (sparse) {
				if (before==0 && after!=0) { // nasty bug fixed!
					int pos = nonZeroIndexes.binarySearch(i);
					pos = -pos -1;
					nonZeroIndexes.beforeInsert(pos,i);
				}
				if (before!=0 && after==0) {
					nonZeroIndexes.remove(nonZeroIndexes.binarySearch(i));
				}
			}
		}
	
		// Find pivot and exchange if necessary.
		int p = j;
		if (p < m) {
			double max = Math.abs(LUcolj.getQuick(p));
			for (int i = j+1; i < m; i++) {
				double v = Math.abs(LUcolj.getQuick(i));
				if (v > max) {
					p = i;
					max = v;
				}
			}
		}
		if (p != j) {
			LUrows[p].swap(LUrows[j]);
			int k = piv[p]; piv[p] = piv[j]; piv[j] = k;
			pivsign = -pivsign;
		}
		
		// Compute multipliers.
		double jj;
		if (j < m && (jj=LU.getQuick(j,j)) != 0.0) {
			multFunction.multiplicator = 1 / jj;
			LU.viewColumn(j).viewPart(j+1,m-(j+1)).assign(multFunction);
		}
		
	}
	setLU(LU);
}
 

/**
Modifies the given vector <tt>A</tt> such that it is permuted as specified; Useful for pivoting.
Cell <tt>A[i]</tt> will go into cell <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 vector to permute.
@param   indexes the permutation indexes, must satisfy <tt>indexes.length==A.size() && indexes[i] >= 0 && indexes[i] < A.size()</tt>;
@param   work the working storage, must satisfy <tt>work.length >= A.size()</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.size()</tt>.
*/
public DoubleMatrix1D permute(DoubleMatrix1D A, int[] indexes, double[] work) {
	// check validity
	int size = A.size();
	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
	*/

	if (work==null || size > work.length) {
		work = A.toArray();
	}
	else {
		A.toArray(work);
	}
	for (int i=size; --i >= 0; ) A.setQuick(i, work[indexes[i]]);
	return A;
}
 

/**
Decomposes 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>.
Uses a "left-looking", dot-product, Crout/Doolittle algorithm.
@param  A   any matrix.
*/	
public void decompose(DoubleMatrix2D A) {
	final int CUT_OFF = 10;
	// 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(LU);
		return; // nothing to do
	}
	
	//precompute and cache some views to avoid regenerating them time and again
	DoubleMatrix1D[] LUrows = new DoubleMatrix1D[m];
	for (int i = 0; i < m; i++) LUrows[i] = LU.viewRow(i);
	
	cern.colt.list.IntArrayList nonZeroIndexes = new cern.colt.list.IntArrayList(); // sparsity
	DoubleMatrix1D LUcolj = LU.viewColumn(0).like();  // blocked column j
	cern.jet.math.Mult multFunction = cern.jet.math.Mult.mult(0); 

	// Outer loop.
	for (int j = 0; j < n; j++) {
		// blocking (make copy of j-th column to localize references)
		LUcolj.assign(LU.viewColumn(j));
		
		// sparsity detection
		int maxCardinality = m/CUT_OFF; // == heuristic depending on speedup
		LUcolj.getNonZeros(nonZeroIndexes,null,maxCardinality);
		int cardinality = nonZeroIndexes.size(); 
		boolean sparse = (cardinality < maxCardinality);

		// Apply previous transformations.
		for (int i = 0; i < m; i++) {
			int kmax = Math.min(i,j);
			double s;
			if (sparse) {
				s = LUrows[i].zDotProduct(LUcolj,0,kmax,nonZeroIndexes);
			}
			else {
				s = LUrows[i].zDotProduct(LUcolj,0,kmax);
			}
			double before = LUcolj.getQuick(i);
			double after = before -s;
			LUcolj.setQuick(i, after); // LUcolj is a copy
			LU.setQuick(i,j, after);   // this is the original
			if (sparse) {
				if (before==0 && after!=0) { // nasty bug fixed!
					int pos = nonZeroIndexes.binarySearch(i);
					pos = -pos -1;
					nonZeroIndexes.beforeInsert(pos,i);
				}
				if (before!=0 && after==0) {
					nonZeroIndexes.remove(nonZeroIndexes.binarySearch(i));
				}
			}
		}
	
		// Find pivot and exchange if necessary.
		int p = j;
		if (p < m) {
			double max = Math.abs(LUcolj.getQuick(p));
			for (int i = j+1; i < m; i++) {
				double v = Math.abs(LUcolj.getQuick(i));
				if (v > max) {
					p = i;
					max = v;
				}
			}
		}
		if (p != j) {
			LUrows[p].swap(LUrows[j]);
			int k = piv[p]; piv[p] = piv[j]; piv[j] = k;
			pivsign = -pivsign;
		}
		
		// Compute multipliers.
		double jj;
		if (j < m && (jj=LU.getQuick(j,j)) != 0.0) {
			multFunction.multiplicator = 1 / jj;
			LU.viewColumn(j).viewPart(j+1,m-(j+1)).assign(multFunction);
		}
		
	}
	setLU(LU);
}