Java Code Examples for org.apache.commons.math3.linear.RealMatrix#transpose()

/**
 * Reads a matrix from the underlying PoN file and check its dimensions.
 * @param fullPath the target data-set full path within the HDF5 file.
 * @param expectedRowCount a predicate that returns true iff its argument is an expected number of rows.
 * @param expectedColumnCount a predicate that returns true iff its argument is an expected number of columns.
 * @return GATKException if the result matrix dimensions do not match the expectations or
 *  any other cause as described in {@link #readMatrix(String)}.
 */
private RealMatrix readMatrixAndCheckDimensions(final String fullPath, final IntPredicate expectedRowCount, final IntPredicate expectedColumnCount) {
    final RealMatrix result = readMatrix(fullPath);
    if (expectedRowCount.test(result.getRowDimension())
            && expectedColumnCount.test(result.getColumnDimension())) {
        return result;
    }
    final RealMatrix transpose = result.transpose();
    if (!expectedRowCount.test(transpose.getRowDimension())) {
        throw new GATKException(String.format("wrong number of rows in '%s' matrix from file '%s': %d",
                fullPath, file.getFile(), result.getRowDimension()));
    }
    if (!expectedColumnCount.test(transpose.getColumnDimension())) {
        throw new GATKException(String.format("wrong number of columns in '%s' from file '%s': %d",
                fullPath, file.getFile(), result.getColumnDimension()));
    }
    return transpose;
}
 

/**
 *
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBeta(RealVector y, RealMatrix x) {
    if (solver == Solver.QR) {
        return computeBetaQR(y, x);
    } else {
        final int n = x.getRowDimension();
        final int p = x.getColumnDimension();
        final int offset = hasIntercept() ? 1 : 0;
        final RealMatrix xT = x.transpose();
        final RealMatrix xTxInv = new LUDecomposition(xT.multiply(x)).getSolver().getInverse();
        final RealVector betaVector = xTxInv.multiply(xT).operate(y);
        final RealVector residuals = y.subtract(x.operate(betaVector));
        this.rss = residuals.dotProduct(residuals);
        this.errorVariance = rss / (n - p);
        this.stdError = Math.sqrt(errorVariance);
        this.residuals = createResidualsFrame(residuals);
        final RealMatrix covMatrix = xTxInv.scalarMultiply(errorVariance);
        final RealMatrix result = new Array2DRowRealMatrix(p, 2);
        if (hasIntercept()) {
            result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
            result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
        }
        for (int i = 0; i < getRegressors().size(); i++) {
            final int index = i + offset;
            final double variance = covMatrix.getEntry(index, index);
            result.setEntry(index, 1, variance);
            result.setEntry(index, 0, betaVector.getEntry(index));
        }
        return result;
    }
}
 

/**
 * Find eigenvalues and eigenvectors of given tridiagonal matrix T.
 *
 * @see http://web.csulb.edu/~tgao/math423/s94.pdf
 * @see http://stats.stackexchange.com/questions/20643/finding-matrix-eigenvectors-using-qr-
 *      decomposition
 * @param T target tridiagonal matrix
 * @param nIter number of iterations for the QR method
 * @param eigvals eigenvalues are stored here
 * @param eigvecs eigenvectors are stored here
 */
public static void tridiagonalEigen(@Nonnull final RealMatrix T, @Nonnull final int nIter,
        @Nonnull final double[] eigvals, @Nonnull final RealMatrix eigvecs) {
    Preconditions.checkArgument(Arrays.deepEquals(T.getData(), T.transpose().getData()),
        "Target matrix T must be a symmetric (tridiagonal) matrix");
    Preconditions.checkArgument(eigvecs.getRowDimension() == eigvecs.getColumnDimension(),
        "eigvecs must be a square matrix");
    Preconditions.checkArgument(T.getRowDimension() == eigvecs.getRowDimension(),
        "T and eigvecs must be the same shape");
    Preconditions.checkArgument(eigvals.length == eigvecs.getRowDimension(),
        "Number of eigenvalues and eigenvectors must be same");

    int nEig = eigvals.length;

    // initialize eigvecs as an identity matrix
    eigvecs.setSubMatrix(eye(nEig), 0, 0);

    RealMatrix T_ = T.copy();

    for (int i = 0; i < nIter; i++) {
        // QR decomposition for the tridiagonal matrix T
        RealMatrix R = new Array2DRowRealMatrix(nEig, nEig);
        RealMatrix Qt = new Array2DRowRealMatrix(nEig, nEig);
        tridiagonalQR(T_, R, Qt);

        RealMatrix Q = Qt.transpose();
        T_ = R.multiply(Q);
        eigvecs.setSubMatrix(eigvecs.multiply(Q).getData(), 0, 0);
    }

    // diagonal elements correspond to the eigenvalues
    for (int i = 0; i < nEig; i++) {
        eigvals[i] = T_.getEntry(i, i);
    }
}
 

@Test
public void testCombinedMatrices2D_Toeplitz() {
    RealMatrix e1 = new Array2DRowRealMatrix(new double[] {1, 1.1});
    RealMatrix e2 = new Array2DRowRealMatrix(new double[] {2, 2.2});
    RealMatrix e2T = e2.transpose();
    RealMatrix e3 = new Array2DRowRealMatrix(new double[] {3, 3.3});
    RealMatrix e3T = e3.transpose();
    RealMatrix[] m1 = new RealMatrix[] {e1, e2, e3};
    // {1.0,1.1}
    // {2.0,2.2}
    // {3.0,3.3}
    RealMatrix[][] toeplitz1 = MatrixUtils.toeplitz(m1, 3);
    Assert.assertEquals(3, toeplitz1.length);
    Assert.assertEquals(3, toeplitz1[0].length);
    Assert.assertEquals(3, toeplitz1[1].length);
    Assert.assertEquals(3, toeplitz1[2].length);
    Assert.assertEquals(e1, toeplitz1[0][0]);
    Assert.assertEquals(e1, toeplitz1[1][1]);
    Assert.assertEquals(e1, toeplitz1[2][2]);
    Assert.assertEquals(e2, toeplitz1[1][0]);
    Assert.assertEquals(e2, toeplitz1[2][1]);
    Assert.assertEquals(e3, toeplitz1[2][0]);
    Assert.assertEquals(e2T, toeplitz1[0][1]);
    Assert.assertEquals(e2T, toeplitz1[1][2]);
    Assert.assertEquals(e3T, toeplitz1[0][2]);
    // {1.0,1.1}  {2.0,2.2}' {3.0,3.3}'
    // {2.0,2.2}  {1.0,1.1}  {2.0,2.2}'
    // {3.0,3.3}  {2.0,2.2}  {1.0,1.1}
    RealMatrix flatten1 = MatrixUtils.combinedMatrices(toeplitz1, 2);
    // 1.0 0.0 2.0 2.2 3.0 3.3
    // 1.1 0.0 0.0 0.0 0.0 0.0
    // 2.0 0.0 1.0 0.0 2.0 2.2
    // 2.2 0.0 1.1 0.0 0.0 0.0
    // 3.0 0.0 2.0 0.0 1.0 0.0
    // 3.3 0.0 2.2 0.0 1.1 0.0
    Assert.assertEquals(6, flatten1.getRowDimension());
    Assert.assertEquals(6, flatten1.getColumnDimension());
}
 

/**
 * Tangent normalize given the raw PoN data using Spark:  the code here is a little more complex for optimization purposes.
 *
 *  Please see notes in docs/PoN ...
 *
 *  Ahat^T = (C^T P^T) A^T
 *  Therefore, C^T is the RowMatrix
 *
 *  pinv: P
 *  panel: A
 *  projection: Ahat
 *  cases: C
 *  betahat: C^T P^T
 *  tangentNormalizedCounts: C - Ahat
 */
private static PCATangentNormalizationResult tangentNormalizeSpark(final ReadCountCollection targetFactorNormalizedCounts,
                                                                   final RealMatrix reducedPanelCounts,
                                                                   final RealMatrix reducedPanelPInvCounts,
                                                                   final CaseToPoNTargetMapper targetMapper,
                                                                   final RealMatrix tangentNormalizationInputCounts,
                                                                   final JavaSparkContext ctx) {
    // Make the C^T a distributed matrix (RowMatrix)
    final RowMatrix caseTDistMat = SparkConverter.convertRealMatrixToSparkRowMatrix(
            ctx, tangentNormalizationInputCounts.transpose(), TN_NUM_SLICES_SPARK);

    // Spark local matrices (transposed)
    final Matrix pinvTLocalMat = new DenseMatrix(
            reducedPanelPInvCounts.getRowDimension(), reducedPanelPInvCounts.getColumnDimension(),
            Doubles.concat(reducedPanelPInvCounts.getData()), true).transpose();
    final Matrix panelTLocalMat = new DenseMatrix(
            reducedPanelCounts.getRowDimension(), reducedPanelCounts.getColumnDimension(),
            Doubles.concat(reducedPanelCounts.getData()), true).transpose();

    // Calculate the projection transpose in a distributed matrix, then convert to Apache Commons matrix (not transposed)
    final RowMatrix betahatDistMat = caseTDistMat.multiply(pinvTLocalMat);
    final RowMatrix projectionTDistMat = betahatDistMat.multiply(panelTLocalMat);
    final RealMatrix projection = SparkConverter.convertSparkRowMatrixToRealMatrix(
            projectionTDistMat, tangentNormalizationInputCounts.transpose().getRowDimension()).transpose();

    // Subtract the projection from the cases
    final RealMatrix tangentNormalizedCounts = tangentNormalizationInputCounts.subtract(projection);

    // Construct the result object and return it with the correct targets.
    final ReadCountCollection tangentNormalized = targetMapper.fromPoNtoCaseCountCollection(
            tangentNormalizedCounts, targetFactorNormalizedCounts.columnNames());
    final ReadCountCollection preTangentNormalized = targetMapper.fromPoNtoCaseCountCollection(
            tangentNormalizationInputCounts, targetFactorNormalizedCounts.columnNames());
    final RealMatrix tangentBetaHats = SparkConverter.convertSparkRowMatrixToRealMatrix(
            betahatDistMat, tangentNormalizedCounts.getColumnDimension());

    return new PCATangentNormalizationResult(tangentNormalized, preTangentNormalized, tangentBetaHats.transpose(), targetFactorNormalizedCounts);
}
 

/**
 * Asserts that a collection of beta-hats corresponds to the expected value given
 * the input pre-tangent normalization matrix and the PoN file.
 */
private void assertBetaHats(final ReadCountCollection preTangentNormalized,
                            final RealMatrix actual, final File ponFile) {
    Assert.assertEquals(actual.getColumnDimension(), preTangentNormalized.columnNames().size());
    final double epsilon = PCATangentNormalizationUtils.EPSILON;

    try (final HDF5File ponReader = new HDF5File(ponFile)) {
        final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
        final List<String> ponTargets = pon.getPanelTargetNames();
        final RealMatrix inCounts = reorderTargetsToPoNOrder(preTangentNormalized, ponTargets);

        // obtain subset of relevant targets to calculate the beta-hats;
        final int[][] betaHatTargets = new int[inCounts.getColumnDimension()][];
        for (int i = 0; i < inCounts.getColumnDimension(); i++) {
            final List<Integer> relevantTargets = new ArrayList<>();
            for (int j = 0; j < inCounts.getRowDimension(); j++) {
                if (inCounts.getEntry(j, i) > 1  +  (Math.log(epsilon) / Math.log(2))) {
                    relevantTargets.add(j);
                }
            }
            betaHatTargets[i] = relevantTargets.stream().mapToInt(Integer::intValue).toArray();
        }
        // calculate beta-hats per column and check with actual values.
        final RealMatrix normalsInv = pon.getReducedPanelPInverseCounts();
        Assert.assertEquals(actual.getRowDimension(), normalsInv.getRowDimension());
        final RealMatrix normalsInvT = normalsInv.transpose();
        for (int i = 0; i < inCounts.getColumnDimension(); i++) {
            final RealMatrix inValues = inCounts.getColumnMatrix(i).transpose().getSubMatrix(new int[] { 0 }, betaHatTargets[i]);
            final RealMatrix normalValues = normalsInvT.getSubMatrix(betaHatTargets[i], IntStream.range(0, normalsInvT.getColumnDimension()).toArray());
            final RealMatrix betaHats = inValues.multiply(normalValues);
            for (int j = 0; j < actual.getRowDimension(); j++) {
                Assert.assertEquals(actual.getEntry(j, i), betaHats.getEntry(0, j),0.000001,"Col " + i + " row " + j);
            }
        }
    }
}
 

/**
 * Assert that the given matrix is unitary.
 * @param m
 */
public static void assertUnitaryMatrix(final RealMatrix m){
    final RealMatrix mInv = new QRDecomposition(m).getSolver().getInverse();
    final RealMatrix mT = m.transpose();

    for (int i = 0; i < mInv.getRowDimension(); i ++) {
        for (int j = 0; j < mInv.getColumnDimension(); j ++) {
            Assert.assertEquals(mInv.getEntry(i, j), mT.getEntry(i, j), 1e-7);
        }
    }
}
 

/**
 * Check that the given matrix is unitary.
 */
public static boolean isUnitaryMatrix(final RealMatrix m){
    //Note can't use MatrixUtils.inverse because m may not be square
    final RealMatrix mInv = new QRDecomposition(m).getSolver().getInverse();
    final RealMatrix mT = m.transpose();

    for (int i = 0; i < mInv.getRowDimension(); i ++) {
        for (int j = 0; j < mInv.getColumnDimension(); j ++) {
            if (Math.abs(mInv.getEntry(i, j) - mT.getEntry(i, j)) > 1.0e-7){
                return false;
            }
        }
    }
    return true;
}
 

/** Calculates the coefficients of linear ridge regression. */
public static double[] performRegression(Matrix a, Matrix b, double ridge) {
	RealMatrix x = MatrixUtils.createRealMatrix(a.getArray());
	RealMatrix y = MatrixUtils.createRealMatrix(b.getArray());
	int numberOfColumns = x.getColumnDimension();
	double[] coefficients = new double[numberOfColumns];
	RealMatrix xTransposed = x.transpose();
	Matrix result;
	boolean finished = false;
	while (!finished) {
		RealMatrix xTx = xTransposed.multiply(x);

		for (int i = 0; i < numberOfColumns; i++) {
			xTx.addToEntry(i, i, ridge);
		}

		RealMatrix xTy = xTransposed.multiply(y);
		coefficients = xTy.getColumn(0);

		try {
			// do not use Apache LUDecomposition for solve instead because it creates different
			// results
			result = new Matrix(xTx.getData()).solve(new Matrix(coefficients, coefficients.length));
			for (int i = 0; i < numberOfColumns; i++) {
				coefficients[i] = result.get(i, 0);
			}
			finished = true;
		} catch (Exception ex) {
			double ridgeOld = ridge;
			if (ridge > 0) {
				ridge *= 10;
			} else {
				ridge = 0.0000001;
			}
			finished = false;
			logger.warning("Error during calculation: " + ex.getMessage() + ": Increasing ridge factor from " + ridgeOld
					+ " to " + ridge);
		}
	}
	return coefficients;
}