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

Example 36 with Matrix

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

the class VarianceProportionStatistic method getMatrixSqrt.

// TODO: Move method below to a different class
// TODO: implement this for DenseMatrix64F rather than Matrix
private static Matrix getMatrixSqrt(Matrix M, Boolean invert) {
    DoubleMatrix2D S = new DenseDoubleMatrix2D(M.toComponents());
    RobustEigenDecomposition eigenDecomp = new RobustEigenDecomposition(S, 100);
    DoubleMatrix1D eigenValues = eigenDecomp.getRealEigenvalues();
    int dim = eigenValues.size();
    for (int i = 0; i < dim; i++) {
        double value = sqrt(eigenValues.get(i));
        if (invert) {
            value = 1 / value;
        }
        eigenValues.set(i, value);
    }
    DoubleMatrix2D eigenVectors = eigenDecomp.getV();
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            eigenVectors.set(i, j, eigenVectors.get(i, j) * eigenValues.get(j));
        }
    }
    DoubleMatrix2D storageMatrix = new DenseDoubleMatrix2D(dim, dim);
    eigenVectors.zMult(eigenDecomp.getV(), storageMatrix, 1, 0, false, true);
    return new Matrix(storageMatrix.toArray());
}

Example 37 with Matrix

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

the class MultivariateNormalGibbsOperator method getPrecision.

private Matrix getPrecision() throws IllegalDimension {
    Matrix currentPrecision = new Matrix(likelihoodPrecision.getParameterAsMatrix());
    currentPrecision = currentPrecision.product(likelihood.getDataList().size());
    /*
for(int i=0; i<currentPrecision.columns(); i++){
for(int j=0; j<currentPrecision.rows(); j++){
System.err.print(currentPrecision.toComponents()[i][j]);
System.err.print(" ");}
System.err.print("\n"); }
*/
    return priorPrecision.add(currentPrecision);
}

Example 38 with Matrix

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

the class MultivariateOUModel method calculateConditionPrecision.

private void calculateConditionPrecision() {
    int index = 0;
    double[] tempW = new double[Ksquared];
    double[][] G = gamma.getParameterAsMatrix();
    for (double deltaTime : deltaTimeList) {
        Q.getTransitionProbabilities(-deltaTime, tempW);
        System.arraycopy(tempW, 0, Wt, Ksquared * index, Ksquared);
        // needs to start with zeros
        double[][] WG = new double[K][K];
        for (int i = 0; i < K; i++) {
            for (int j = 0; j < K; j++) {
                for (int k = 0; k < K; k++) WG[i][j] += tempW[i * K + k] * G[k][j];
            }
        }
        // needs to start with zeros
        double[][] WGWt = new double[K][K];
        for (int i = 0; i < K; i++) {
            for (int j = 0; j < K; j++) {
                for (int k = 0; k < K; k++) WGWt[i][j] += WG[i][k] * tempW[j * K + k];
            }
        }
        for (int i = 0; i < K; i++) {
            for (int j = 0; j < K; j++) WGWt[i][j] = G[i][j] - WGWt[i][j];
        }
        WGWt = new Matrix(WGWt).inverse().toComponents();
        for (int i = 0; i < K; i++) System.arraycopy(WGWt[i], 0, conditionPrecisionVector, Ksquared * index + K * i, K);
        index++;
    }
    conditionalPrecisionKnown = true;
}

Example 39 with Matrix

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

the class MultivariateOUModel method calculateLogLikelihood.

public double calculateLogLikelihood() {
    double logLikelihood = 0;
    double[] previous = new double[K];
    double[] current = new double[K];
    double[] tmpHolder;
    double[][] G = gamma.getParameterAsMatrix();
    double[] theta = dependentParam.getParameterValues();
    double[] Xbeta = null;
    boolean hasEffects = getNumberOfFixedEffects() > 0;
    if (!conditionalPrecisionKnown)
        calculateConditionPrecision();
    try {
        if (new Matrix(G).determinant() < 0.01)
            return Double.NEGATIVE_INFINITY;
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    int index = 0;
    if (!hasEffects) {
        for (int i = 0; i < K; i++) previous[i] = theta[index++];
    } else {
        Xbeta = getXBeta();
        for (int i = 0; i < K; i++) {
            previous[i] = theta[index] - Xbeta[index];
            index++;
        }
    }
    initialPrior = new MultivariateNormalDistribution(initialPriorMean, new Matrix(G).inverse().toComponents());
    logLikelihood += initialPrior.logPdf(previous);
    double save = logLikelihood;
    double save2 = 0;
    int oldMapTime = -1;
    double[][] conditionalPrecision = new double[K][K];
    for (int timeStep = 0; timeStep < numTimeSteps; timeStep++) {
        int thisMapTime = mapTime[timeStep];
        if (thisMapTime != oldMapTime) {
            System.arraycopy(Wt, Ksquared * thisMapTime, W, 0, Ksquared);
            for (int i = 0; i < K; i++) System.arraycopy(conditionPrecisionVector, Ksquared * thisMapTime + K * i, conditionalPrecision[i], 0, K);
            oldMapTime = thisMapTime;
        }
        double[] mean = new double[K];
        int u = 0;
        for (int i = 0; i < K; i++) {
            for (int j = 0; j < K; j++) mean[i] += W[u++] * previous[j];
        }
        if (!hasEffects) {
            for (int i = 0; i < K; i++) current[i] = theta[index++];
        } else {
            for (int i = 0; i < K; i++) {
                current[i] = theta[index] - Xbeta[index];
                index++;
            }
        }
        MultivariateNormalDistribution density = new MultivariateNormalDistribution(mean, conditionalPrecision);
        double partialLogLikelihood = density.logPdf(current);
        if (partialLogLikelihood > 10) {
            return Double.NEGATIVE_INFINITY;
        }
        logLikelihood += partialLogLikelihood;
        // move to next point
        tmpHolder = previous;
        previous = current;
        current = tmpHolder;
    }
    if (logLikelihood > 100) {
        System.err.println("got here end");
        System.err.println("save1 = " + save);
        System.err.println("save2 = " + save2);
        System.exit(-1);
    }
    likelihoodKnown = true;
    return logLikelihood;
}

Example 40 with Matrix

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

the class PositiveDefiniteSubstitutionModel method setupMatrix.

/**
 * setup substitution matrix
 */
protected void setupMatrix() {
    if (!eigenInitialised)
        initialiseEigen();
    int i, j, k = 0;
    amat = new Matrix(rates.getParameterAsMatrix()).inverse().toComponents();
    // copy q matrix for unit testing
    for (i = 0; i < amat.length; i++) {
        System.arraycopy(amat[i], 0, q[i], 0, amat[i].length);
    }
    // compute eigenvalues and eigenvectors
    elmhes(amat, ordr, stateCount);
    eltran(amat, Evec, ordr, stateCount);
    hqr2(stateCount, 1, stateCount, amat, Evec, Eval, evali);
    luinverse(Evec, Ievc, stateCount);
    updateMatrix = false;
}