Search in sources :

Example 6 with DenseDoubleMatrix2D

use of cern.colt.matrix.impl.DenseDoubleMatrix2D in project beast-mcmc by beast-dev.

the class GeneralizedLinearModelParser method checkFullRank.

private void checkFullRank(DesignMatrix designMatrix) throws XMLParseException {
    int fullRank = designMatrix.getColumnDimension();
    //        System.err.println("designMatrix getColumnDimension = "+fullRank);
    SingularValueDecomposition svd = new SingularValueDecomposition(new DenseDoubleMatrix2D(designMatrix.getParameterAsMatrix()));
    int realRank = svd.rank();
    if (realRank != fullRank) {
        throw new XMLParseException("rank(" + designMatrix.getId() + ") = " + realRank + ".\nMatrix is not of full rank as colDim(" + designMatrix.getId() + ") = " + fullRank);
    }
}
Also used : DenseDoubleMatrix2D(cern.colt.matrix.impl.DenseDoubleMatrix2D) SingularValueDecomposition(cern.colt.matrix.linalg.SingularValueDecomposition)

Example 7 with DenseDoubleMatrix2D

use of cern.colt.matrix.impl.DenseDoubleMatrix2D in project beast-mcmc by beast-dev.

the class MarkovModulatedFrequencyModel method computeStationaryDistribution.

private void computeStationaryDistribution(double[] statDistr) {
    if (allRatesAreZero(switchingRates)) {
        return;
    }
    // Uses an LU decomposition to solve Q^t \pi = 0 and \sum \pi_i = 1
    DoubleMatrix2D mat2 = new DenseDoubleMatrix2D(numBaseModel + 1, numBaseModel);
    int index2 = 0;
    for (int i = 0; i < numBaseModel; ++i) {
        for (int j = i + 1; j < numBaseModel; ++j) {
            // Transposed
            mat2.set(j, i, switchingRates.getParameterValue(index2));
            index2++;
        }
    }
    for (int j = 0; j < numBaseModel; ++j) {
        for (int i = j + 1; i < numBaseModel; ++i) {
            // Transposed
            mat2.set(j, i, switchingRates.getParameterValue(index2));
            index2++;
        }
    }
    for (int i = 0; i < numBaseModel; ++i) {
        double rowTotal = 0.0;
        for (int j = 0; j < numBaseModel; ++j) {
            if (i != j) {
                // Transposed
                rowTotal += mat2.get(j, i);
            }
        }
        mat2.set(i, i, -rowTotal);
    }
    // Add row for sum-to-one constraint
    for (int i = 0; i < numBaseModel; ++i) {
        mat2.set(numBaseModel, i, 1.0);
    }
    LUDecomposition decomp = new LUDecomposition(mat2);
    DoubleMatrix2D x = new DenseDoubleMatrix2D(numBaseModel + 1, 1);
    x.set(numBaseModel, 0, 1.0);
    DoubleMatrix2D y = decomp.solve(x);
    for (int i = 0; i < numBaseModel; ++i) {
        statDistr[i] = y.get(i, 0);
    }
//System.err.println(new Vector(statDistr));              
}
Also used : DenseDoubleMatrix2D(cern.colt.matrix.impl.DenseDoubleMatrix2D) DoubleMatrix2D(cern.colt.matrix.DoubleMatrix2D) LUDecomposition(cern.colt.matrix.linalg.LUDecomposition) DenseDoubleMatrix2D(cern.colt.matrix.impl.DenseDoubleMatrix2D)

Aggregations

DenseDoubleMatrix2D (cern.colt.matrix.impl.DenseDoubleMatrix2D)7 SingularValueDecomposition (cern.colt.matrix.linalg.SingularValueDecomposition)4 DoubleMatrix2D (cern.colt.matrix.DoubleMatrix2D)3 RobustEigenDecomposition (dr.math.matrixAlgebra.RobustEigenDecomposition)2 RobustSingularValueDecomposition (dr.math.matrixAlgebra.RobustSingularValueDecomposition)2 DoubleMatrix1D (cern.colt.matrix.DoubleMatrix1D)1 LUDecomposition (cern.colt.matrix.linalg.LUDecomposition)1 Matrix (dr.math.matrixAlgebra.Matrix)1