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

Example 41 with Matrix

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

the class CorrelationSymmetricMatrixTest method testGetParameterAsMatrix.

public void testGetParameterAsMatrix() {
    CorrelationSymmetricMatrix matrix = getMatrix(CorrelationSymmetricMatrix.Type.AS_IS);
    double[][] M = matrix.getParameterAsMatrix();
    System.out.println(new Matrix(M));
    double sum = 0.0;
    for (double x : M[0]) {
        sum += x;
    }
    assertEquals(sum, 7.0);
}

Example 42 with Matrix

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

the class LKJTransformTest method testTransformation.

public void testTransformation() {
    System.out.println("\nTest LKJ transform.");
    double[] transformedValue = transform.inverse(CPCs, 0, CPCs.length);
    double[] inverseTransformedValues = transform.transform(transformedValue, 0, transformedValue.length);
    SymmetricMatrix R = compoundCorrelationSymmetricMatrix(transformedValue, dim);
    System.out.println("transformedValue=" + R);
    try {
        assertTrue("Positive Definite", R.isPD());
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    System.out.println("iCPC=" + new Matrix(inverseTransformedValues, dim * (dim - 1) / 2, 1));
    assertEquals("size CPCs", format.format(CPCs.length), format.format(inverseTransformedValues.length));
    for (int k = 0; k < CPCs.length; k++) {
        assertEquals("inverse transform k=" + k, format.format(CPCs[k]), format.format(inverseTransformedValues[k]));
    }
}

Example 43 with Matrix

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

the class LKJTransformTest method testJacobianCholesky.

public void testJacobianCholesky() {
    System.out.println("\nTest LKJ Cholesky Jacobian.");
    // Matrix
    double[][] jacobianMat = transformChol.computeJacobianMatrixInverse(CPCs);
    Matrix Jac = new Matrix(jacobianMat);
    System.out.println("Jacobian Matrix=" + Jac.transpose());
    assertEquals("size Jacobian Matrix", format.format(dim * (dim - 1) / 2), format.format(Jac.rows()));
    assertEquals("size Jacobian Matrix", format.format(dim * (dim - 1) / 2), format.format(Jac.columns()));
    // Determinant
    double jacobianDet = (new Transform.InverseMultivariate(transformChol)).getLogJacobian(CPCs, 0, CPCs.length);
    double jacobianDetBis = 0;
    for (int i = 0; i < jacobianMat[0].length; i++) {
        jacobianDetBis += Math.log(jacobianMat[i][i]);
    }
    System.out.println("Log Jacobiant Det direct=" + jacobianDet);
    System.out.println("Log Jacobiant Det matrix=" + jacobianDetBis);
    assertEquals("jacobian log det", format.format(jacobianDet), format.format(jacobianDetBis));
}

Example 44 with Matrix

use of dr.math.matrixAlgebra.Matrix 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();
    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: ").append(treeVariance.length).append("\n");
    sb.append("data dim: ").append(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 (double partial : partials) {
                data[offset] = partial;
                ++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: ").append(getLogLikelihood()).append(" == ").append(logDensity).append("\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()).append("\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) {
        System.arraycopy(data, i * p, tmp[i], 0, p);
    }
    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();
}

Example 45 with Matrix

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

the class PrecisionMatrixGibbsOperator method doOperationDontFireChange.

public void doOperationDontFireChange() {
    if (wishartIsModel) {
        // TODO Deprecate
        setupWishartStatistics(priorModel);
        priorStatistics = setupStatistics(priorModel);
    }
    final double[][] scaleMatrix = getOperationScaleMatrixAndSetObservationCount();
    final double treeDf = numberObservations;
    final double df = priorDf + treeDf * pathWeight;
    double[][] draw = WishartDistribution.nextWishart(df, scaleMatrix);
    if (DEBUG) {
        System.err.println("draw = " + new Matrix(draw));
    }
    for (int i = 0; i < dim; i++) {
        Parameter column = precisionParam.getParameter(i);
        for (int j = 0; j < dim; j++) column.setParameterValueQuietly(j, draw[j][i]);
    }
}