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

Example 1 with Vector

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

the class BeagleOperationReport method calculateLogLikelihood.

protected double calculateLogLikelihood() {
    if (matrixUpdateIndices == null) {
        matrixUpdateIndices = new int[eigenCount][nodeCount];
        branchLengths = new double[eigenCount][nodeCount];
        branchUpdateCount = new int[eigenCount];
    //            scaleBufferIndices = new int[internalNodeCount];
    //            storedScaleBufferIndices = new int[internalNodeCount];
    }
    if (operations == null) {
        operations = new int[numRestrictedPartials + 1][internalNodeCount * Beagle.OPERATION_TUPLE_SIZE];
        operationCount = new int[numRestrictedPartials + 1];
    }
    recomputeScaleFactors = false;
    for (int i = 0; i < eigenCount; i++) {
        branchUpdateCount[i] = 0;
    }
    operationListCount = 0;
    if (hasRestrictedPartials) {
        for (int i = 0; i <= numRestrictedPartials; i++) {
            operationCount[i] = 0;
        }
    } else {
        operationCount[0] = 0;
    }
    System.out.println(alignmentString.toString());
    final NodeRef root = treeModel.getRoot();
    // Do not flip buffers
    traverse(treeModel, root, null, false);
    for (int i = 0; i < eigenCount; i++) {
        if (branchUpdateCount[i] > 0) {
            if (DEBUG_BEAGLE_OPERATIONS) {
                StringBuilder sb = new StringBuilder();
                sb.append("eval = ").append(new Vector(substitutionModel.getEigenDecomposition().getEigenValues())).append("\n");
                sb.append("evec = ").append(new Vector(substitutionModel.getEigenDecomposition().getEigenVectors())).append("\n");
                sb.append("ivec = ").append(new Vector(substitutionModel.getEigenDecomposition().getInverseEigenVectors())).append("\n");
                sb.append("Branch count: ").append(branchUpdateCount[i]);
                sb.append("\nNode indices:\n");
                if (SINGLE_LINE) {
                    sb.append("int n[] = {");
                }
                for (int k = 0; k < branchUpdateCount[i]; ++k) {
                    if (SINGLE_LINE) {
                        sb.append(" ").append(matrixUpdateIndices[i][k]);
                        if (k < (branchUpdateCount[i] - 1)) {
                            sb.append(",");
                        }
                    } else {
                        sb.append(matrixUpdateIndices[i][k]).append("\n");
                    }
                }
                if (SINGLE_LINE) {
                    sb.append(" };\n");
                }
                sb.append("\nBranch lengths:\n");
                if (SINGLE_LINE) {
                    sb.append("double b[] = {");
                }
                for (int k = 0; k < branchUpdateCount[i]; ++k) {
                    if (SINGLE_LINE) {
                        sb.append(" ").append(branchLengths[i][k]);
                        if (k < (branchUpdateCount[i] - 1)) {
                            sb.append(",");
                        }
                    } else {
                        sb.append(branchLengths[i][k]).append("\n");
                    }
                }
                if (SINGLE_LINE) {
                    sb.append(" };\n");
                }
                System.out.println(sb.toString());
            }
        }
    }
    if (DEBUG_BEAGLE_OPERATIONS) {
        StringBuilder sb = new StringBuilder();
        sb.append("Operation count: ").append(operationCount[0]);
        sb.append("\nOperations:\n");
        if (SINGLE_LINE) {
            sb.append("BeagleOperation o[] = {");
        }
        for (int k = 0; k < operationCount[0] * Beagle.OPERATION_TUPLE_SIZE; ++k) {
            if (SINGLE_LINE) {
                sb.append(" ").append(operations[0][k]);
                if (k < (operationCount[0] * Beagle.OPERATION_TUPLE_SIZE - 1)) {
                    sb.append(",");
                }
            } else {
                sb.append(operations[0][k]).append("\n");
            }
        }
        if (SINGLE_LINE) {
            sb.append(" };\n");
        }
        sb.append("Use scale factors: ").append(useScaleFactors).append("\n");
        System.out.println(sb.toString());
    }
    int rootIndex = partialBufferHelper.getOffsetIndex(root.getNumber());
    System.out.println("Root node: " + rootIndex);
    return 0.0;
}

Example 2 with Vector

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

the class ApproximateFactorAnalysisPrecisionMatrix method computeValuesImp.

private void computeValuesImp() {
    dim = L.getRowDimension();
    double[][] matrix = new double[dim][dim];
    for (int row = 0; row < L.getRowDimension(); ++row) {
        for (int col = 0; col < L.getRowDimension(); ++col) {
            double sum = 0;
            for (int k = 0; k < L.getColumnDimension(); ++k) {
                sum += L.getParameterValue(row, k) * L.getParameterValue(col, k);
            }
            matrix[row][col] = sum;
        }
    }
    for (int row = 0; row < dim; row++) {
        matrix[row][row] += 1 / gamma.getParameterValue(row);
    }
    if (DEBUG) {
        System.err.println("mult:");
        System.err.println(new Matrix(L.getParameterAsMatrix()));
        System.err.println(new Vector(gamma.getParameterValues()) + "\n");
        System.err.println(new Matrix(matrix));
    }
    matrix = new Matrix(matrix).inverse().toComponents();
    int index = 0;
    for (int row = 0; row < dim; ++row) {
        for (int col = 0; col < dim; ++col) {
            values[index] = matrix[row][col];
            ++index;
        }
    }
}

Example 3 with Vector

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

the class MultivariateNormalGibbsOperator method getMean.

private Vector getMean() throws IllegalDimension {
    Vector meanSum = new Vector(getMeanSum());
    Matrix workingPrecision = new Matrix(likelihoodPrecision.getParameterAsMatrix());
    Vector meanPart = workingPrecision.product(meanSum);
    meanPart = meanPart.add(priorPrecision.product(priorMean));
    Matrix varPart = getPrecision().inverse();
    Vector answer = varPart.product(meanPart);
    return answer;
}

Example 4 with Vector

use of dr.math.matrixAlgebra.Vector 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) {
    }
}

Example 5 with Vector

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

the class FullyConjugateMultivariateTraitLikelihood method calculateAscertainmentCorrection.

protected double calculateAscertainmentCorrection(int taxonIndex) {
    NodeRef tip = treeModel.getNode(taxonIndex);
    int nodeIndex = treeModel.getNode(taxonIndex).getNumber();
    if (ascertainedData == null) {
        // Assumes that ascertained data are fixed
        ascertainedData = new double[dimTrait];
    }
    //        diffusionModel.diffusionPrecisionMatrixParameter.setParameterValue(0,2); // For debugging non-1 values
    double[][] traitPrecision = diffusionModel.getPrecisionmatrix();
    double logDetTraitPrecision = Math.log(diffusionModel.getDeterminantPrecisionMatrix());
    double lengthToRoot = getRescaledLengthToRoot(tip);
    double marginalPrecisionScalar = 1.0 / lengthToRoot + rootPriorSampleSize;
    double logLikelihood = 0;
    for (int datum = 0; datum < numData; ++datum) {
        // Get observed trait value
        System.arraycopy(meanCache, nodeIndex * dim + datum * dimTrait, ascertainedData, 0, dimTrait);
        if (DEBUG_ASCERTAINMENT) {
            System.err.println("Datum #" + datum);
            System.err.println("Value: " + new Vector(ascertainedData));
            System.err.println("Cond : " + lengthToRoot);
            System.err.println("MargV: " + 1.0 / marginalPrecisionScalar);
            System.err.println("MargP: " + marginalPrecisionScalar);
            System.err.println("diffusion prec: " + new Matrix(traitPrecision));
        }
        double SSE;
        if (dimTrait > 1) {
            throw new RuntimeException("Still need to implement multivariate ascertainment correction");
        } else {
            double precision = traitPrecision[0][0] * marginalPrecisionScalar;
            SSE = ascertainedData[0] * precision * ascertainedData[0];
        }
        double thisLogLikelihood = -LOG_SQRT_2_PI * dimTrait + 0.5 * (logDetTraitPrecision + dimTrait * Math.log(marginalPrecisionScalar) - SSE);
        if (DEBUG_ASCERTAINMENT) {
            System.err.println("LogLik: " + thisLogLikelihood);
            dr.math.distributions.NormalDistribution normal = new dr.math.distributions.NormalDistribution(0, Math.sqrt(1.0 / (traitPrecision[0][0] * marginalPrecisionScalar)));
            System.err.println("TTTLik: " + normal.logPdf(ascertainedData[0]));
            if (datum >= 10) {
                System.exit(-1);
            }
        }
        logLikelihood += thisLogLikelihood;
    }
    return logLikelihood;
}