Examples with Matrix dr.math.matrixAlgebra.Matrix used on opensource projects

Example 21 with Matrix

use of dr.math.matrixAlgebra.Matrix in project beast-mcmc by beast-dev.

the class InverseWishartDistribution method logPdf.

public double logPdf(double[] x) {
    Matrix W = new Matrix(x, dim, dim);
    double logDensity = 0;
    try {
        logDensity = Math.log(W.determinant());
        logDensity *= -0.5;
        logDensity *= df + dim + 1;
        Matrix product = S.product(W.inverse());
        for (int i = 0; i < dim; i++) logDensity -= 0.5 * product.component(i, i);
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    logDensity += logNormalizationConstant;
    return logDensity;
}

Example 22 with Matrix

use of dr.math.matrixAlgebra.Matrix in project beast-mcmc by beast-dev.

the class WishartDistribution method logPdf.

public static double logPdf(Matrix W, Matrix Sinv, double df, int dim, double logNormalizationConstant) {
    double logDensity = 0;
    try {
        // if (!W.isPD()) { // TODO isPD() does not appear to work
        // return Double.NEGATIVE_INFINITY;
        // }
        // Returns NaN is W is not positive-definite.
        logDensity = W.logDeterminant();
        if (Double.isInfinite(logDensity) || Double.isNaN(logDensity)) {
            return Double.NEGATIVE_INFINITY;
        }
        logDensity *= 0.5;
        logDensity *= df - dim - 1;
        // the whole matrix
        if (Sinv != null) {
            Matrix product = Sinv.product(W);
            for (int i = 0; i < dim; i++) logDensity -= 0.5 * product.component(i, i);
        }
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    logDensity += logNormalizationConstant;
    return logDensity;
}

Example 23 with Matrix

use of dr.math.matrixAlgebra.Matrix in project beast-mcmc by beast-dev.

the class ComplexSubstitutionModel method printDebugSetupMatrix.

private void printDebugSetupMatrix() {
    System.out.println("Normalized infinitesimal rate matrix:");
    System.out.println(new Matrix(amat));
    System.out.println(new Matrix(amat).toStringOctave());
    // System.out.println("Normalization = " + normalization);
    System.out.println("Values in setupMatrix():");
// System.out.println(eigenV);
// System.out.println(eigenVInv);
// System.out.println(eigenVReal);
}

Example 24 with Matrix

use of dr.math.matrixAlgebra.Matrix 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 25 with Matrix

use of dr.math.matrixAlgebra.Matrix in project beast-mcmc by beast-dev.

the class BinaryCovarionModelTest method testEquilibriumDistribution.

/**
 * Tests that pi*Q = 0
 */
public void testEquilibriumDistribution() {
    alpha.setParameterValue(0, 0.1);
    switchingRate.setParameterValue(0, 1.0);
    model.setupMatrix();
    double[] pi = model.getFrequencyModel().getFrequencies();
    try {
        Matrix m = new Matrix(model.getQ());
        Vector p = new Vector(pi);
        Vector y = m.product(p);
        assertEquals(0.0, y.norm(), 1e-14);
    } catch (IllegalDimension illegalDimension) {
    }
}