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