Examples with SingularValueDecomposition org.apache.commons.math3.linear.SingularValueDecomposition used on opensource projects

Example 1 with SingularValueDecomposition

use of org.apache.commons.math3.linear.SingularValueDecomposition in project jstructure by JonStargaryen.

the class SVDSuperimposer method align.

@Override
public StructureAlignmentResult align(AtomContainer reference, AtomContainer query) {
    AtomContainer originalReference = reference;
    AtomContainer originalCandidate = query;
    Pair<GroupContainer, GroupContainer> atomContainerPair = AbstractAlignmentAlgorithm.comparableGroupContainerPair(reference, query, minimalSetOfAtomNames, maximalSetOfAtomNames);
    reference = atomContainerPair.getLeft();
    query = atomContainerPair.getRight();
    // calculate centroids
    double[] centroid1 = reference.calculate().centroid().getValue();
    double[] centroid2 = query.calculate().centroid().getValue();
    // center atoms
    reference.calculate().center();
    query.calculate().center();
    // compose covariance matrix and calculate SVD
    RealMatrix matrix1 = convertToMatrix(reference);
    RealMatrix matrix2 = convertToMatrix(query);
    RealMatrix covariance = matrix2.transpose().multiply(matrix1);
    SingularValueDecomposition svd = new SingularValueDecomposition(covariance);
    // R = (V * U')'
    RealMatrix ut = svd.getU().transpose();
    RealMatrix rotationMatrix = svd.getV().multiply(ut).transpose();
    // check if reflection
    if (new LUDecomposition(rotationMatrix).getDeterminant() < 0) {
        RealMatrix v = svd.getV().transpose();
        v.setEntry(2, 0, (0 - v.getEntry(2, 0)));
        v.setEntry(2, 1, (0 - v.getEntry(2, 1)));
        v.setEntry(2, 2, (0 - v.getEntry(2, 2)));
        rotationMatrix = v.transpose().multiply(ut).transpose();
    }
    double[][] rotation = rotationMatrix.getData();
    // calculate translation
    double[] translation = LinearAlgebra.on(centroid1).subtract(LinearAlgebra.on(centroid2).multiply(rotation)).getValue();
    logger.trace("rotation matrix\n{}\ntranslation vector\n{}", Arrays.deepToString(rotationMatrix.getData()), Arrays.toString(translation));
    /* transform 2nd atom select - employ neutral translation (3D vector of zeros), because the atoms are already
        * centered and calculate RMSD */
    query.calculate().transform(new Transformation(rotation));
    double rmsd = calculateRmsd(reference, query);
    // return alignment
    return new StructureAlignmentResult(originalReference, originalCandidate, query, rmsd, translation, rotation);
}

Example 2 with SingularValueDecomposition

use of org.apache.commons.math3.linear.SingularValueDecomposition in project gatk by broadinstitute.

the class SparkSingularValueDecomposer method createSVD.

/**
     *  Create a SVD of the given matrix using the given Java Spark Context.
     *
     * @param realMat the matrix target.  Not {@code null}
     * @return never {@code null}
     */
@Override
public SVD createSVD(final RealMatrix realMat) {
    Utils.nonNull(realMat, "Cannot perform Spark MLLib SVD on a null matrix.");
    final RowMatrix mat = SparkConverter.convertRealMatrixToSparkRowMatrix(sc, realMat, NUM_SLICES);
    // Compute all of the singular values and corresponding singular vectors.
    final SingularValueDecomposition<RowMatrix, Matrix> svd = mat.computeSVD((int) mat.numCols(), true, 1.0E-9d);
    // Get our distributed results
    final RowMatrix u = svd.U();
    final Vector s = svd.s();
    final Matrix v = svd.V().transpose();
    // Move the matrices from Spark/distributed space to Apache Commons space
    logger.info("Converting distributed Spark matrix to local matrix...");
    final RealMatrix uReal = SparkConverter.convertSparkRowMatrixToRealMatrix(u, realMat.getRowDimension());
    logger.info("Done converting distributed Spark matrix to local matrix...");
    logger.info("Converting Spark matrix to local matrix...");
    final RealMatrix vReal = SparkConverter.convertSparkMatrixToRealMatrix(v);
    logger.info("Done converting Spark matrix to local matrix...");
    final double[] singularValues = s.toArray();
    logger.info("Calculating the pseudoinverse...");
    logger.info("Pinv: calculating tolerance...");
    // Note that the pinv of realMat is V * invS * U'
    final double tolerance = Math.max(realMat.getColumnDimension(), realMat.getRowDimension()) * realMat.getNorm() * EPS;
    logger.info("Pinv: inverting the singular values (with tolerance) and creating a diagonal matrix...");
    final double[] invS = Arrays.stream(singularValues).map(sv -> invertSVWithTolerance(sv, tolerance)).toArray();
    final Matrix invSMat = Matrices.diag(Vectors.dense(invS));
    logger.info("Pinv: Multiplying V * invS * U' to get the pinv (using pinv transpose = U * invS' * V') ...");
    final RowMatrix pinvT = u.multiply(invSMat).multiply(v);
    logger.info("Pinv: Converting back to local matrix ...");
    final RealMatrix pinv = SparkConverter.convertSparkRowMatrixToRealMatrix(pinvT, realMat.getRowDimension()).transpose();
    logger.info("Done calculating the pseudoinverse and converting it...");
    return new SimpleSVD(uReal, s.toArray(), vReal, pinv);
}

Example 3 with SingularValueDecomposition

use of org.apache.commons.math3.linear.SingularValueDecomposition in project gatk by broadinstitute.

the class ApacheSingularValueDecomposer method createSVD.

/** Create a SVD instance using Apache Commons Math.
     *
     * @param m matrix that is not {@code null}
     * @return SVD instance that is never {@code null}
     */
@Override
public SVD createSVD(final RealMatrix m) {
    Utils.nonNull(m, "Cannot create SVD on a null matrix.");
    final SingularValueDecomposition svd = new SingularValueDecomposition(m);
    final RealMatrix pinv = svd.getSolver().getInverse();
    return new SimpleSVD(svd.getU(), svd.getSingularValues(), svd.getV(), pinv);
}