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

Example 1 with CholeskyDecomposition

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

the class DeterminentalPointProcessPrior method computeLogLikelihood.

public double computeLogLikelihood() {
    if (pathSampling && notZero != null && sum != data.getColumnDimension()) {
        return Double.NEGATIVE_INFINITY;
    }
    CholeskyDecomposition chol = null;
    try {
        chol = new CholeskyDecomposition(relationshipList);
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    double product = 0;
    for (int i = 0; i < relationshipList.length; i++) {
        product += Math.log(chol.getL()[i][i]);
    }
    product *= 2;
    return product;
}

Example 2 with CholeskyDecomposition

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

the class LKJTransformConstrained method transform.

@Override
protected double[] transform(double[] values) {
    SymmetricMatrix R = compoundCorrelationSymmetricMatrix(values, dimVector);
    double[] L;
    try {
        L = (new CholeskyDecomposition(R)).getStrictlyUpperTriangular();
    } catch (IllegalDimension illegalDimension) {
        throw new RuntimeException("Unable to decompose matrix in LKJ inverse transform.");
    }
    double[] results = super.transform(L);
    if (DEBUG) {
        System.err.println("R: " + compoundCorrelationSymmetricMatrix(values, dimVector));
        System.err.println("L: " + new WrappedMatrix.WrappedStrictlyUpperTriangularMatrix(L, dimVector));
        System.err.println("Z: " + compoundCorrelationSymmetricMatrix(results, dimVector));
    }
    return results;
}

Example 3 with CholeskyDecomposition

use of dr.math.matrixAlgebra.CholeskyDecomposition 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 4 with CholeskyDecomposition

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

the class AdaptableVarianceMultivariateNormalOperator method doOperation.

public double doOperation() {
    iterations++;
    if (DEBUG) {
        System.err.println("\nAVMVN Iteration: " + iterations);
        System.err.println("Using AdaptableVarianceMultivariateNormalOperator: " + iterations + " for " + parameter.getParameterName());
        System.err.println("Old parameter values:");
        for (int i = 0; i < dim; i++) {
            System.err.println(parameter.getParameterValue(i));
        }
    }
    double[] x = parameter.getParameterValues();
    // transform to the appropriate scale
    double[] transformedX = new double[dim];
    /*for (int i = 0; i < dim; i++) {
            transformedX[i] = transformations[i].transform(x[i]);
        }*/
    // iterate over transformation sizes rather than number of parameters
    // as a transformation might impact multiple parameters
    int currentIndex = 0;
    for (int i = 0; i < transformationSizes.length; i++) {
        if (DEBUG) {
            System.err.println("currentIndex = " + currentIndex);
            System.err.println("transformationSizes[i] = " + transformationSizes[i]);
        }
        if (transformationSizes[i] > 1) {
            System.arraycopy(transformations[i].transform(x, currentIndex, currentIndex + transformationSizes[i] - 1), 0, transformedX, currentIndex, transformationSizes[i]);
        } else {
            transformedX[currentIndex] = transformations[i].transform(x[currentIndex]);
            if (DEBUG) {
                System.err.println("x[" + currentIndex + "] = " + x[currentIndex] + " -> " + transformedX[currentIndex]);
            }
        }
        currentIndex += transformationSizes[i];
    }
    if (DEBUG) {
        System.err.println("Old transformed parameter values:");
        for (int i = 0; i < dim; i++) {
            System.err.println(transformedX[i]);
        }
    }
    // store MH-ratio in logq
    double logJacobian = 0.0;
    // change this: make a rule for when iterations == burnin
    if (iterations > 1 && iterations > burnin) {
        if (DEBUG) {
            System.err.println("  AVMVN iterations > burnin");
        }
        if (iterations > (burnin + 1)) {
            if (iterations % every == 0) {
                updates++;
                if (DEBUG) {
                    System.err.println("updates = " + updates);
                }
                // first recalculate the means using recursion
                for (int i = 0; i < dim; i++) {
                    newMeans[i] = ((oldMeans[i] * (updates - 1)) + transformedX[i]) / updates;
                }
                if (updates > 1) {
                    // here we can simply use the double[][] matrix
                    for (int i = 0; i < dim; i++) {
                        for (int j = i; j < dim; j++) {
                            empirical[i][j] = calculateCovariance(updates, empirical[i][j], transformedX, i, j);
                            empirical[j][i] = empirical[i][j];
                        }
                    }
                }
                if (DEBUG) {
                    System.err.println("Old means:");
                    for (int i = 0; i < dim; i++) {
                        System.err.println(oldMeans[i]);
                    }
                    System.err.println("New means:");
                    for (int i = 0; i < dim; i++) {
                        System.err.println(newMeans[i]);
                    }
                    System.err.println("Empirical covariance matrix:");
                    for (int i = 0; i < dim; i++) {
                        for (int j = 0; j < dim; j++) {
                            System.err.print(empirical[i][j] + " ");
                        }
                        System.err.println();
                    }
                }
            }
        /*if (iterations == 17) {
                    System.exit(0);
                }*/
        } else if (iterations == (burnin + 1)) {
            // this will not be reached when burnin is set to 0
            for (int i = 0; i < dim; i++) {
                // oldMeans[i] = transformedX[i];
                // newMeans[i] = transformedX[i];
                oldMeans[i] = 0.0;
                newMeans[i] = 0.0;
            }
            for (int i = 0; i < dim; i++) {
                for (int j = 0; j < dim; j++) {
                    empirical[i][j] = 0.0;
                }
            }
        }
    // TODO End of adaptable covariance -- move into separate class
    } else if (iterations == 1) {
        if (DEBUG) {
            System.err.println("\niterations == 1");
        }
        // iterations == 1
        for (int i = 0; i < dim; i++) {
            // oldMeans[i] = transformedX[i];
            // newMeans[i] = transformedX[i];
            oldMeans[i] = 0.0;
            newMeans[i] = 0.0;
        }
        for (int i = 0; i < dim; i++) {
            for (int j = 0; j < dim; j++) {
                empirical[i][j] = 0.0;
                proposal[i][j] = matrix[i][j];
            }
        }
    }
    for (int i = 0; i < dim; i++) {
        epsilon[i] = scaleFactor * MathUtils.nextGaussian();
    }
    if (iterations > initial) {
        if (DEBUG) {
            System.err.println("  iterations > initial");
        }
        if (iterations % every == 0) {
            // double[][] proposal = new double[dim][dim];
            for (int i = 0; i < dim; i++) {
                for (int j = i; j < dim; j++) {
                    // symmetric matrix
                    proposal[j][i] = proposal[i][j] = // constantFactor *  /* auto-tuning using scaleFactor */
                    (1 - beta) * empirical[i][j] + beta * matrix[i][j];
                }
            }
            // not necessary for first test phase, but will need to be performed when covariance matrix is being updated
            try {
                cholesky = (new CholeskyDecomposition(proposal)).getL();
            } catch (IllegalDimension illegalDimension) {
                throw new RuntimeException("Unable to decompose matrix in AdaptableVarianceMultivariateNormalOperator");
            }
        // double end = System.nanoTime();
        // double baseResult = end - start;
        // System.err.println("Cholesky decomposition took: " + baseResult);
        }
    }
    if (DEBUG) {
        System.err.println("  Drawing new values");
    }
    for (int i = 0; i < dim; i++) {
        for (int j = i; j < dim; j++) {
            transformedX[i] += cholesky[j][i] * epsilon[j];
        // caution: decomposition returns lower triangular
        }
    }
    if (DEBUG) {
        System.err.println("\nTransformed X values:");
        for (int i = 0; i < dim; i++) {
            System.err.println(transformedX[i]);
        }
        System.err.println();
    }
    // iterate over transformation sizes rather than number of parameters
    // as a transformation might impact multiple parameters
    currentIndex = 0;
    for (int i = 0; i < transformationSizes.length; i++) {
        if (DEBUG) {
            System.err.println("currentIndex = " + currentIndex);
            System.err.println("transformationSizes[i] = " + transformationSizes[i]);
        }
        if (MULTI) {
            if (transformationSizes[i] > 1) {
                if (DEBUG) {
                    System.err.println("Current transformation sum = " + transformationSums[i]);
                }
                double[] temp = transformations[i].inverse(transformedX, currentIndex, currentIndex + transformationSizes[i] - 1, transformationSums[i]);
                for (int k = 0; k < temp.length; k++) {
                    parameter.setParameterValueQuietly(currentIndex + k, temp[k]);
                }
                logJacobian += transformations[i].getLogJacobian(x, currentIndex, currentIndex + transformationSizes[i] - 1) - transformations[i].getLogJacobian(temp, 0, transformationSizes[i] - 1);
            } else {
                parameter.setParameterValueQuietly(currentIndex, transformations[i].inverse(transformedX[currentIndex]));
                logJacobian += transformations[i].getLogJacobian(x[currentIndex]) - transformations[i].getLogJacobian(parameter.getParameterValue(currentIndex));
            }
            if (DEBUG) {
                System.err.println("Current logJacobian = " + logJacobian);
            }
        } else {
            if (transformationSizes[i] > 1) {
                // TODO: figure out if this is really a problem ...
                throw new RuntimeException("Transformations on more than 1 parameter value should be set quietly");
            } else {
                parameter.setParameterValue(currentIndex, transformations[i].inverse(transformedX[currentIndex]));
                logJacobian += transformations[i].getLogJacobian(x[currentIndex]) - transformations[i].getLogJacobian(parameter.getParameterValue(currentIndex));
            }
            if (DEBUG) {
                System.err.println("Current logJacobian = " + logJacobian);
            }
        }
        currentIndex += transformationSizes[i];
    }
    if (DEBUG) {
        System.err.println("Proposed parameter values:");
        for (int i = 0; i < dim; i++) {
            System.err.println(x[i] + " -> " + parameter.getValue(i));
        }
        System.err.println("LogJacobian: " + logJacobian);
    }
    if (MULTI) {
        // Signal once.
        parameter.fireParameterChangedEvent();
    }
    if (iterations % every == 0) {
        if (DEBUG) {
            System.err.println("  Copying means");
        }
        // copy new means to old means for next update iteration
        // System.arraycopy(newMeans, 0, oldMeans, 0, dim);
        double[] tmp = oldMeans;
        oldMeans = newMeans;
        // faster to swap pointers
        newMeans = tmp;
    }
    // return 0.0;
    return logJacobian;
}

Example 5 with CholeskyDecomposition

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

the class WishartDistribution method nextWishart.

/**
 * Generate a random draw from a Wishart distribution
 * Follows Odell and Feiveson (1996) JASA 61, 199-203
 * <p/>
 * Returns a random variable with expectation = df * scaleMatrix
 *
 * @param df          degrees of freedom
 * @param scaleMatrix scaleMatrix
 * @return a random draw
 */
public static double[][] nextWishart(double df, double[][] scaleMatrix) {
    int dim = scaleMatrix.length;
    double[][] draw = new double[dim][dim];
    double[][] z = new double[dim][dim];
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < i; j++) {
            z[i][j] = MathUtils.nextGaussian();
        }
    }
    for (int i = 0; i < dim; i++) // sqrt of chisq with df-i dfs
    z[i][i] = Math.sqrt(MathUtils.nextGamma((df - i) * 0.5, 0.5));
    double[][] cholesky = new double[dim][dim];
    for (int i = 0; i < dim; i++) {
        for (int j = i; j < dim; j++) cholesky[i][j] = cholesky[j][i] = scaleMatrix[i][j];
    }
    try {
        cholesky = (new CholeskyDecomposition(cholesky)).getL();
    // caution: this returns the lower triangular form
    } catch (IllegalDimension illegalDimension) {
        throw new RuntimeException("Numerical exception in WishartDistribution");
    }
    double[][] result = new double[dim][dim];
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            // lower triangular
            for (// can also be shortened
            int k = 0; // can also be shortened
            k < dim; // can also be shortened
            k++) result[i][j] += cholesky[i][k] * z[k][j];
        }
    }
    for (int i = 0; i < dim; i++) {
        // lower triangular, so more efficiency is possible
        for (int j = 0; j < dim; j++) {
            for (int k = 0; k < dim; k++) // transpose of 2nd element
            draw[i][j] += result[i][k] * result[j][k];
        }
    }
    return draw;
}