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

Example 11 with IllegalDimension

use of dr.math.matrixAlgebra.IllegalDimension 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;
                }
            }
        }
    } 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) {
                double[] temp = transformations[i].inverse(transformedX, currentIndex, currentIndex + transformationSizes[i] - 1);
                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 12 with IllegalDimension

use of dr.math.matrixAlgebra.IllegalDimension 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 13 with IllegalDimension

use of dr.math.matrixAlgebra.IllegalDimension 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;
}