Examples with DenseDoubleMatrix2D cern.colt.matrix.impl.DenseDoubleMatrix2D used on opensource projects

Example 1 with DenseDoubleMatrix2D

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

the class ComplexSubstitutionModel method setupMatrix.

public void setupMatrix() {
    if (!eigenInitialised) {
        initialiseEigen();
        storedEvalImag = new double[stateCount];
    }
    int i = 0;
    storeIntoAmat();
    makeValid(amat, stateCount);
    // compute eigenvalues and eigenvectors
    //        EigenvalueDecomposition eigenDecomp = new EigenvalueDecomposition(new DenseDoubleMatrix2D(amat));
    RobustEigenDecomposition eigenDecomp;
    try {
        eigenDecomp = new RobustEigenDecomposition(new DenseDoubleMatrix2D(amat), maxIterations);
    } catch (ArithmeticException ae) {
        System.err.println(ae.getMessage());
        wellConditioned = false;
        System.err.println("amat = \n" + new Matrix(amat));
        return;
    }
    DoubleMatrix2D eigenV = eigenDecomp.getV();
    DoubleMatrix1D eigenVReal = eigenDecomp.getRealEigenvalues();
    DoubleMatrix1D eigenVImag = eigenDecomp.getImagEigenvalues();
    DoubleMatrix2D eigenVInv;
    if (checkConditioning) {
        RobustSingularValueDecomposition svd;
        try {
            svd = new RobustSingularValueDecomposition(eigenV, maxIterations);
        } catch (ArithmeticException ae) {
            System.err.println(ae.getMessage());
            wellConditioned = false;
            return;
        }
        if (svd.cond() > maxConditionNumber) {
            wellConditioned = false;
            return;
        }
    }
    try {
        eigenVInv = alegbra.inverse(eigenV);
    } catch (IllegalArgumentException e) {
        wellConditioned = false;
        return;
    }
    Ievc = eigenVInv.toArray();
    Evec = eigenV.toArray();
    Eval = eigenVReal.toArray();
    EvalImag = eigenVImag.toArray();
    // Check for valid decomposition
    for (i = 0; i < stateCount; i++) {
        if (Double.isNaN(Eval[i]) || Double.isNaN(EvalImag[i]) || Double.isInfinite(Eval[i]) || Double.isInfinite(EvalImag[i])) {
            wellConditioned = false;
            return;
        } else if (Math.abs(Eval[i]) < 1e-10) {
            Eval[i] = 0.0;
        }
    }
    updateMatrix = false;
    wellConditioned = true;
    // compute normalization and rescale eigenvalues
    computeStationaryDistribution();
    if (doNormalization) {
        double subst = 0.0;
        for (i = 0; i < stateCount; i++) subst += -amat[i][i] * stationaryDistribution[i];
        for (i = 0; i < stateCount; i++) {
            Eval[i] /= subst;
            EvalImag[i] /= subst;
        }
    }
}

Example 2 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);
    }
}

Example 3 with DenseDoubleMatrix2D

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

the class GeneralizedLinearModel method getAllIndependentVariablesIdentifiable.

public boolean getAllIndependentVariablesIdentifiable() {
    int totalColDim = 0;
    for (DesignMatrix mat : designMatrix) {
        totalColDim += mat.getColumnDimension();
    }
    double[][] grandDesignMatrix = new double[N][totalColDim];
    int offset = 0;
    for (DesignMatrix mat : designMatrix) {
        final int length = mat.getColumnDimension();
        for (int i = 0; i < N; ++i) {
            for (int j = 0; j < length; ++j) {
                grandDesignMatrix[i][offset + j] = mat.getParameterValue(i, j);
            }
        }
        offset += length;
    }
    double[][] mat = grandDesignMatrix;
    if (grandDesignMatrix.length < grandDesignMatrix[0].length) {
        mat = new double[grandDesignMatrix[0].length][grandDesignMatrix.length];
        for (int i = 0; i < grandDesignMatrix.length; ++i) {
            for (int j = 0; j < grandDesignMatrix[i].length; ++j) {
                mat[j][i] = grandDesignMatrix[i][j];
            }
        }
    }
    SingularValueDecomposition svd = new SingularValueDecomposition(new DenseDoubleMatrix2D(mat));
    int rank = svd.rank();
    boolean isFullRank = (totalColDim == rank);
    Logger.getLogger("dr.inference").info("\tTotal # of predictors = " + totalColDim + " and rank = " + rank);
    return isFullRank;
}

Example 4 with DenseDoubleMatrix2D

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

the class ColtEigenSystem method decomposeMatrix.

public EigenDecomposition decomposeMatrix(double[][] matrix) {
    final int stateCount = matrix.length;
    RobustEigenDecomposition eigenDecomp = new RobustEigenDecomposition(new DenseDoubleMatrix2D(matrix), maxIterations);
    DoubleMatrix2D eigenV = eigenDecomp.getV();
    DoubleMatrix2D eigenVInv;
    if (checkConditioning) {
        RobustSingularValueDecomposition svd;
        try {
            svd = new RobustSingularValueDecomposition(eigenV, maxIterations);
        } catch (ArithmeticException ae) {
            System.err.println(ae.getMessage());
            return getEmptyDecomposition(stateCount);
        }
        if (svd.cond() > maxConditionNumber) {
            return getEmptyDecomposition(stateCount);
        }
    }
    try {
        eigenVInv = alegbra.inverse(eigenV);
    } catch (IllegalArgumentException e) {
        return getEmptyDecomposition(stateCount);
    }
    double[][] Evec = eigenV.toArray();
    double[][] Ievc = eigenVInv.toArray();
    double[] Eval = getAllEigenValues(eigenDecomp);
    if (checkConditioning) {
        for (int i = 0; i < Eval.length; i++) {
            if (Double.isNaN(Eval[i]) || Double.isInfinite(Eval[i])) {
                return getEmptyDecomposition(stateCount);
            } else if (Math.abs(Eval[i]) < 1e-10) {
                Eval[i] = 0.0;
            }
        }
    }
    double[] flatEvec = new double[stateCount * stateCount];
    double[] flatIevc = new double[stateCount * stateCount];
    for (int i = 0; i < Evec.length; i++) {
        System.arraycopy(Evec[i], 0, flatEvec, i * stateCount, stateCount);
        System.arraycopy(Ievc[i], 0, flatIevc, i * stateCount, stateCount);
    }
    return new EigenDecomposition(flatEvec, flatIevc, Eval);
}

Example 5 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));              
}