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

Example 6 with IllegalDimension

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

the class SemiConjugateMultivariateTraitLikelihood method integrateLogLikelihoodAtRoot.

protected double integrateLogLikelihoodAtRoot(double[] y, double[] Ay, double[][] AplusB, double[][] treePrecision, double rootPrecision) {
    double detAplusB = 0;
    double square = 0;
    if (dimTrait > 1) {
        for (int i = 0; i < dimTrait; i++) {
            // Ay is filled with sum, and original value is destroyed
            Ay[i] += Bz[i];
            for (int j = 0; j < dimTrait; j++) {
                AplusB[i][j] = treePrecision[i][j] * rootPrecision + rootPriorPrecision[i][j];
            }
        }
        Matrix mat = new Matrix(AplusB);
        try {
            detAplusB = mat.determinant();
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        double[][] invAplusB = mat.inverse().toComponents();
        for (int i = 0; i < dimTrait; i++) {
            for (int j = 0; j < dimTrait; j++) square += Ay[i] * invAplusB[i][j] * Ay[j];
        }
    } else {
        // 1D is very simple
        detAplusB = treePrecision[0][0] * rootPrecision + rootPriorPrecision[0][0];
        Ay[0] += Bz[0];
        square = Ay[0] * Ay[0] / detAplusB;
    }
    double retValue = 0.5 * (logRootPriorPrecisionDeterminant - Math.log(detAplusB) - zBz + square);
    if (DEBUG) {
        System.err.println("(Ay+Bz)(A+B)^{-1}(Ay+Bz) = " + square);
        System.err.println("density = " + retValue);
        System.err.println("zBz = " + zBz);
    }
    return retValue;
}

Example 7 with IllegalDimension

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

the class SemiConjugateMultivariateTraitLikelihood method setRootPrior.

private void setRootPrior(MultivariateNormalDistribution rootPrior) {
    rootPriorMean = rootPrior.getMean();
    rootPriorPrecision = rootPrior.getScaleMatrix();
    try {
        logRootPriorPrecisionDeterminant = Math.log(new Matrix(rootPriorPrecision).determinant());
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    setRootPriorSumOfSquares();
}

Example 8 with IllegalDimension

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

the class DeterminentalPointProcessPrior method computeLogLikelihood.

public double computeLogLikelihood() {
    int newSize = data.getColumnDimension();
    while (!changedList.isEmpty()) {
        int index = changedList.remove(0);
        int row = index % data.getRowDimension();
        int col = index / data.getRowDimension();
        for (int i = 0; i < data.getColumnDimension(); i++) {
            if (col != i) {
                if (data.getParameterValue(row, col) == data.getParameterValue(row, i)) {
                    relationshipList[col][i] *= Math.exp(1 / (theta * theta));
                    relationshipList[i][col] = relationshipList[col][i];
                } else {
                    relationshipList[col][i] *= Math.exp(-1 / (theta * theta));
                    relationshipList[i][col] = relationshipList[col][i];
                }
            }
        }
    }
    if (newSize != size) {
        size = newSize;
        relationshipList = new double[size][size];
        reset();
    }
    CholeskyDecomposition chol = null;
    try {
        chol = new CholeskyDecomposition(relationshipList);
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    double product = 0;
    for (int i = 0; i < newSize; i++) {
        product += Math.log(chol.getL()[i][i]);
    }
    product *= 2;
    return product;
}

Example 9 with IllegalDimension

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

the class GaussianProcessFromTree method nextRandomFast.

//    boolean firstTime=true;
public double[] nextRandomFast() {
    double[] random = new double[traitModel.getTreeModel().getExternalNodeCount() * traitModel.getDimTrait()];
    NodeRef root = traitModel.getTreeModel().getRoot();
    double[] traitStart = traitModel.getPriorMean();
    double[][] varianceCholesky = null;
    double[][] temp = new SymmetricMatrix(traitModel.getDiffusionModel().getPrecisionmatrix()).inverse().toComponents();
    try {
        varianceCholesky = (new CholeskyDecomposition(temp).getL());
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    //        }
    if (USE_BUFFER) {
        final int length = traitModel.getDimTrait();
        final int nodeCount = traitModel.getTreeModel().getNodeCount();
        double[] currentValue = new double[(nodeCount + 1) * length];
        double[] epsilon = new double[length];
        final int priorOffset = nodeCount * length;
        System.arraycopy(traitStart, 0, currentValue, priorOffset, length);
        nextRandomFast2(currentValue, priorOffset, root, random, varianceCholesky, epsilon);
    } else {
        nextRandomFast(traitStart, root, random, varianceCholesky);
    }
    //        }
    return random;
}

Example 10 with IllegalDimension

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

the class FullyConjugateMultivariateTraitLikelihood method getReport.

@Override
public String getReport() {
    StringBuilder sb = new StringBuilder();
    //        sb.append(this.g)
    //        System.err.println("Hello");
    sb.append("Tree:\n");
    sb.append(getId()).append("\t");
    sb.append(treeModel.toString());
    sb.append("\n\n");
    double[][] treeVariance = computeTreeVariance(true);
    double[][] traitPrecision = getDiffusionModel().getPrecisionmatrix();
    Matrix traitVariance = new Matrix(traitPrecision).inverse();
    double[][] jointVariance = KroneckerOperation.product(treeVariance, traitVariance.toComponents());
    sb.append("Tree variance:\n");
    sb.append(new Matrix(treeVariance));
    sb.append(matrixMin(treeVariance)).append("\t").append(matrixMax(treeVariance)).append("\t").append(matrixSum(treeVariance));
    sb.append("\n\n");
    sb.append("Trait variance:\n");
    sb.append(traitVariance);
    sb.append("\n\n");
    //        sb.append("Joint variance:\n");
    //        sb.append(new Matrix(jointVariance));
    //        sb.append("\n\n");
    sb.append("Tree dim: " + treeVariance.length + "\n");
    sb.append("data dim: " + jointVariance.length);
    sb.append("\n\n");
    double[] data = new double[jointVariance.length];
    System.arraycopy(meanCache, 0, data, 0, jointVariance.length);
    if (nodeToClampMap != null) {
        int offset = treeModel.getExternalNodeCount() * getDimTrait();
        for (Map.Entry<NodeRef, RestrictedPartials> clamps : nodeToClampMap.entrySet()) {
            double[] partials = clamps.getValue().getPartials();
            for (int i = 0; i < partials.length; ++i) {
                data[offset] = partials[i];
                ++offset;
            }
        }
    }
    sb.append("Data:\n");
    sb.append(new Vector(data)).append("\n");
    sb.append(data.length).append("\t").append(vectorMin(data)).append("\t").append(vectorMax(data)).append("\t").append(vectorSum(data));
    sb.append(treeModel.getNodeTaxon(treeModel.getExternalNode(0)).getId());
    sb.append("\n\n");
    MultivariateNormalDistribution mvn = new MultivariateNormalDistribution(new double[data.length], new Matrix(jointVariance).inverse().toComponents());
    double logDensity = mvn.logPdf(data);
    sb.append("logLikelihood: " + getLogLikelihood() + " == " + logDensity + "\n\n");
    final WishartSufficientStatistics sufficientStatistics = getWishartStatistics();
    final double[] outerProducts = sufficientStatistics.getScaleMatrix();
    sb.append("Outer-products (DP):\n");
    sb.append(new Vector(outerProducts));
    sb.append(sufficientStatistics.getDf() + "\n");
    Matrix treePrecision = new Matrix(treeVariance).inverse();
    final int n = data.length / traitPrecision.length;
    final int p = traitPrecision.length;
    double[][] tmp = new double[n][p];
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < p; ++j) {
            tmp[i][j] = data[i * p + j];
        }
    }
    Matrix y = new Matrix(tmp);
    Matrix S = null;
    try {
        // Using Matrix-Normal form
        S = y.transpose().product(treePrecision).product(y);
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    sb.append("Outer-products (from tree variance:\n");
    sb.append(S);
    sb.append("\n\n");
    return sb.toString();
}