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

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

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

DoubleMatrix1D getStDVec(DoubleMatrix1D sum, DoubleMatrix1D sumsq, int nsamp){
    //for n samples with given sum and sum of squares (component-wise)
    //get standard deviations for each component
    DoubleMatrix1D ans = sumsq.copy().assign(Functions.mult(1./nsamp));
    ans.assign( sum.copy().assign(Functions.mult(1./nsamp)).assign(Functions.square), Functions.minus  );
    ans.assign(Functions.sqrt);
    return ans;
}
 

/**
 * Normalizes matrix of p(z|u) such that \forall_u: \sum_z p(z|u) = 1.
 *
 * @param pu_z normalized matrix of p(z|u)
 */
@Override
protected void normalizePuz(DoubleMatrix2D pu_z) {
    for (int u = 0; u < pu_z.rows(); u++) {
        DoubleMatrix1D tmp = pu_z.viewRow(u);
        double norm = tmp.aggregate(Functions.plus, Functions.identity);
        if (norm != 0.0) {
            tmp.assign(Functions.mult(1 / norm));
        }
    }
}
 

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

/**
 * Demonstrates advanced sorting.
 * Sorts by sinus of cell values.
 */
public static void zdemo3() {
	Sorting sort = quickSort;
	double[] values = {0.5, 1.5, 2.5, 3.5};
	DoubleMatrix1D matrix = new DenseDoubleMatrix1D(values);
	cern.colt.function.DoubleComparator comp = new cern.colt.function.DoubleComparator() {
		public int compare(double a, double b) {
			double as = Math.sin(a); double bs = Math.sin(b);
			return as < bs ? -1 : as == bs ? 0 : 1;
		}
	};
	System.out.println("unsorted:"+matrix);
	
	DoubleMatrix1D sorted = sort.sort(matrix,comp);
	System.out.println("sorted  :"+sorted);

	// check whether it is really sorted
	sorted.assign(cern.jet.math.Functions.sin);
	/*
	sorted.assign(
		new cern.colt.function.DoubleFunction() {
			public double apply(double arg) { return Math.sin(arg); }
		}
	);
	*/
	System.out.println("sined  :"+sorted);
}
 

public static DoubleMatrix1D getInteriorPt(ArrayList<LinearConstraint> polytope){
    //Find a point in the interior of a polytope,
    //First attempt:
    //maximizing and minimizing a random function, and taking the average
    //Then if still any active constraints, maximize to get to corner opposite constraint,
    //will 
    //average these corners to get an interior pt--just need to have no constraint active
    //at all the corners
    
    
    if(polytope.isEmpty())
        return DoubleFactory1D.dense.make(0);
    
    /*ArrayList<DoubleMatrix1D> corners = new ArrayList<>();
    SimplexSolver ss = new SimplexSolver();
    LinearConstraintSet polytopeConstr = new LinearConstraintSet(polytope);
    DoubleMatrix1D rand = DoubleFactory1D.dense.random(polytope.get(0).getCoefficients().getDimension());
    LinearObjectiveFunction func = new LinearObjectiveFunction(rand.toArray(), 0);
    PointValuePair max = ss.optimize(func, polytopeConstr, GoalType.MAXIMIZE);
    PointValuePair min = ss.optimize(func, polytopeConstr, GoalType.MINIMIZE);
    
    DoubleMatrix1D ans = DoubleFactory1D.dense.make(max.getPoint());
    ans.assign(DoubleFactory1D.dense.make(min.getPoint()),Functions.plus);
    ans.assign(Functions.mult(0.5));
    return ans;*/
    
    //DEBUG!!  how about this just maximize all the constraints and average
    //slow but should be in the interior
    ArrayList<DoubleMatrix1D> corners = new ArrayList<>();
    SimplexSolver ss = new SimplexSolver();
    LinearConstraintSet polytopeConstr = new LinearConstraintSet(polytope);
    for(LinearConstraint constr : polytope){
        LinearObjectiveFunction func = constrAsFunction(constr);
        PointValuePair otherCorner;
        try{
            otherCorner = ss.optimize(func, polytopeConstr, GoalType.MINIMIZE);
        }
        catch(NoFeasibleSolutionException e){
            return null;//there are no interior points
        }
        corners.add(DoubleFactory1D.dense.make(otherCorner.getPoint()));
    }
    
    //do the average
    DoubleMatrix1D ans = DoubleFactory1D.dense.make(polytope.get(0).getCoefficients().getDimension());
    for(DoubleMatrix1D corner : corners){
        ans.assign(corner,Functions.plus);
    }
    ans.assign(Functions.mult(1./corners.size()));
    return ans;
}
 

public void dscal(double alpha, DoubleMatrix1D x) {
	x.assign(F.mult(alpha));
}
 

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

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

public void dscal(double alpha, DoubleMatrix1D x) {
	x.assign(F.mult(alpha));
}
 

/**
 * <tt>A = A * B <=> A[i] = A[i] * B[i]</tt>.
 * @param A the matrix to modify.
 * @param B the matrix to stay unaffected.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D mult(DoubleMatrix1D A, DoubleMatrix1D B) {
	return A.assign(B,F.mult);
}
 

/**
 * <tt>A = A - B <=> A[i] = A[i] - B[i]</tt>.
 * @param A the matrix to modify.
 * @param B the matrix to stay unaffected.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D minus(DoubleMatrix1D A, DoubleMatrix1D B) {
	return A.assign(B,F.minus);
}
 

/**
 * <tt>A = A * B <=> A[i] = A[i] * B[i]</tt>.
 * @param A the matrix to modify.
 * @param B the matrix to stay unaffected.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D mult(DoubleMatrix1D A, DoubleMatrix1D B) {
	return A.assign(B,F.mult);
}
 

/**
 * <tt>A = A * s <=> A[i] = A[i] * s</tt>.
 * @param A the matrix to modify.
 * @param s the scalar; can have any value.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D mult(DoubleMatrix1D A, double s) {
	return A.assign(F.mult(s));
}
 

/**
 * <tt>A = A - s <=> A[i] = A[i] - s</tt>.
 * @param A the matrix to modify.
 * @param s the scalar; can have any value.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D minus(DoubleMatrix1D A, double s) {
	return A.assign(F.minus(s));
}
 

/**
 * <tt>A = A / s <=> A[i] = A[i] / s</tt>.
 * @param A the matrix to modify.
 * @param s the scalar; can have any value.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D div(DoubleMatrix1D A, double s) {
	return A.assign(F.div(s));
}
 

/**
 * <tt>A = A<sup>s</sup> <=> A[i] = Math.pow(A[i], s)</tt>.
 * @param A the matrix to modify.
 * @param s the scalar; can have any value.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D pow(DoubleMatrix1D A, double s) {
	return A.assign(F.pow(s));
}
 

/**
 * <tt>A[i] = Math.abs(A[i])</tt>.
 * @param A the matrix to modify.
 * @return <tt>A</tt> (for convenience only).
 */
public static DoubleMatrix1D abs(DoubleMatrix1D A) {
	return A.assign(F.abs);
}