Examples with LUDecomposition cern.colt.matrix.linalg.LUDecomposition used on opensource projects

Example 1 with LUDecomposition

use of cern.colt.matrix.linalg.LUDecomposition 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));              
}