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

Example 26 with Vector

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

the class WishartStatisticsWrapper method simulateMissingTraits.

private void simulateMissingTraits() {
    // Force new sample!
    likelihoodDelegate.fireModelChanged();
    //        ProcessSimulationDelegate.MeanAndVariance mv =
    //                (ProcessSimulationDelegate.MeanAndVariance) tipFullConditionalTrait.getTrait(
    //                        dataLikelihood.getTree(), dataLikelihood.getTree().getExternalNode(1));
    //
    //        System.err.println("DONE");
    //        System.exit(-1);
    double[] sample = (double[]) tipSampleTrait.getTrait(dataLikelihood.getTree(), null);
    if (DEBUG) {
        System.err.println("Attempting to simulate missing traits");
        System.err.println(new Vector(sample));
    }
    final ContinuousDiffusionIntegrator cdi = outerProductDelegate.getIntegrator();
    assert (cdi instanceof ContinuousDiffusionIntegrator.Basic);
    double[] buffer = new double[dimPartial * numTrait];
    for (int trait = 0; trait < numTrait; ++trait) {
        buffer[trait * dimPartial + dimTrait] = Double.POSITIVE_INFINITY;
    }
    for (int tip = 0; tip < tipCount; ++tip) {
        int sampleOffset = tip * dimTrait * numTrait;
        int bufferOffset = 0;
        for (int trait = 0; trait < numTrait; ++trait) {
            System.arraycopy(sample, sampleOffset, buffer, bufferOffset, dimTrait);
            sampleOffset += dimTrait;
            bufferOffset += dimPartial;
        }
        outerProductDelegate.setTipDataDirectly(tip, buffer);
    }
    if (DEBUG) {
        System.err.println("Finished draw");
    }
}

Example 27 with Vector

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

the class KernelDensityEstimator2D method main.

public static void main(String[] arg) {
    double[] x = { 3.4, 1.2, 5.6, 2.2, 3.1 };
    double[] y = { 1.0, 2.0, 1.0, 2.0, 1.0 };
    KernelDensityEstimator2D kde = new KernelDensityEstimator2D(x, y, 4);
    System.out.println(new Vector(kde.getXGrid()));
    System.out.println(new Vector(kde.getYGrid()));
    System.out.println(new Matrix(kde.getKDE()));
    System.exit(-1);
}

Example 28 with Vector

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

the class SVSComplexSubstitutionModel method main.

public static void main(String[] arg) {
    double[] rates = new double[] { 0.097, 0.515, 3.346, 0.623, 0.389, 0.631, 0.362, 1.127, 4.262, 0.424, 0.758, 1.297, 0.728, 0.228, 0.003, 0.075, 0.312, 0.356, 0.927, 4.420, 2.719, 0.264, 4.267, 0.741, 1.106, 1.568, 1.215, 0.172, 0.204, 1.493, 0.592, 0.105, 1.583, 1.201, 0.783, 2.224, 0.888, 1.401, 0.137, 0.259, 1.197, 0.056, 1.939, 1.385, 0.690, 0.815, 0.279, 1.100, 1.715, 0.011, 1.509, 0.961, 0.112, 1.305, 2.797, 0.578, 1.177, 1.009, 0.316, 1.143, 1.861, 0.176, 0.140, 0.104, 0.571, 0.521, 0.761, 1.795, 1.065, 1.563, 2.972, 2.295, 0.075, 1.690, 1.011, 0.128, 0.484, 0.355, 1.668, 1.052, 0.089, 0.104, 0.056, 1.591, 0.054, 0.487, 1.034, 1.145, 0.403, 0.254, 0.474, 0.181, 0.124, 2.067, 2.208, 0.120, 2.638, 0.195, 1.897, 1.955, 1.113, 1.399, 4.901, 3.218, 0.361, 1.934, 1.681, 3.572, 0.806, 0.077, 0.042, 0.310, 0.243, 1.526, 3.572, 4.173, 0.740, 3.086, 0.645, 2.755, 0.280, 2.476, 0.476, 5.174, 0.057, 0.225, 1.310, 0.201, 0.491, 1.507, 0.604, 0.404, 0.907, 0.048, 2.439, 0.676, 0.382, 0.062, 1.173, 2.026, 2.813, 0.655, 1.511, 0.452, 0.137, 0.435, 0.685, 0.826, 1.735, 0.935, 0.697, 1.230, 5.416, 1.043, 0.747, 0.945 };
    double[] indicators = new double[] { 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 1.000, 0.000, 0.000, 0.000, 1.000, 1.000, 1.000, 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000 };
    Parameter ratesP = new Parameter.Default(rates);
    Parameter indicatorsP = new Parameter.Default(indicators);
    DataType dataType = new DataType() {

        public String getDescription() {
            return null;
        }

        public int getType() {
            return 0;
        }

        @Override
        public char[] getValidChars() {
            return null;
        }

        public int getStateCount() {
            return 13;
        }
    };
    FrequencyModel freqModel = new FrequencyModel(dataType, new Parameter.Default(400, 1.0 / 400.0));
    SVSComplexSubstitutionModel substModel = new SVSComplexSubstitutionModel("test", dataType, freqModel, ratesP, indicatorsP);
    double logL = substModel.getLogLikelihood();
    System.out.println("Prior = " + logL);
    if (!Double.isInfinite(logL)) {
        double[] finiteTimeProbs = new double[substModel.getDataType().getStateCount() * substModel.getDataType().getStateCount()];
        double time = 1.0;
        substModel.getTransitionProbabilities(time, finiteTimeProbs);
        System.out.println("Probs = " + new Vector(finiteTimeProbs));
    }
}

Example 29 with Vector

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

the class IntegratedMultivariateTraitLikelihood method calculateLogLikelihood.

public double calculateLogLikelihood() {
    if (updateRestrictedNodePartials) {
        if (clampList != null) {
            setupClamps();
        }
        updateRestrictedNodePartials = false;
    }
    double logLikelihood = 0;
    double[][] traitPrecision = diffusionModel.getPrecisionmatrix();
    double logDetTraitPrecision = Math.log(diffusionModel.getDeterminantPrecisionMatrix());
    double[] conditionalRootMean = tmp2;
    final boolean computeWishartStatistics = getComputeWishartSufficientStatistics();
    if (computeWishartStatistics) {
        wishartStatistics = new WishartSufficientStatistics(dimTrait);
    }
    // Use dynamic programming to compute conditional likelihoods at each internal node
    postOrderTraverse(treeModel, treeModel.getRoot(), traitPrecision, logDetTraitPrecision, computeWishartStatistics);
    if (DEBUG) {
        System.err.println("mean: " + new Vector(cacheHelper.getMeanCache()));
        System.err.println("correctedMean: " + new Vector(cacheHelper.getCorrectedMeanCache()));
        System.err.println("upre: " + new Vector(upperPrecisionCache));
        System.err.println("lpre: " + new Vector(lowerPrecisionCache));
        System.err.println("cach: " + new Vector(logRemainderDensityCache));
    }
    // Compute the contribution of each datum at the root
    final int rootIndex = treeModel.getRoot().getNumber();
    // Precision scalar of datum conditional on root
    double conditionalRootPrecision = lowerPrecisionCache[rootIndex];
    for (int datum = 0; datum < numData; datum++) {
        double thisLogLikelihood = 0;
        // Get conditional mean of datum conditional on root
        // System.arraycopy(meanCache, rootIndex * dim + datum * dimTrait, conditionalRootMean, 0, dimTrait);
        System.arraycopy(cacheHelper.getMeanCache(), rootIndex * dim + datum * dimTrait, conditionalRootMean, 0, dimTrait);
        if (DEBUG) {
            System.err.println("Datum #" + datum);
            System.err.println("root mean: " + new Vector(conditionalRootMean));
            System.err.println("root prec: " + conditionalRootPrecision);
            System.err.println("diffusion prec: " + new Matrix(traitPrecision));
        }
        // B = root prior precision
        // z = root prior mean
        // A = likelihood precision
        // y = likelihood mean
        // y'Ay
        double yAy = computeWeightedAverageAndSumOfSquares(conditionalRootMean, Ay, traitPrecision, dimTrait, // Also fills in Ay
        conditionalRootPrecision);
        if (conditionalRootPrecision != 0) {
            thisLogLikelihood += -LOG_SQRT_2_PI * dimTrait + 0.5 * (logDetTraitPrecision + dimTrait * Math.log(conditionalRootPrecision) - yAy);
        }
        if (DEBUG) {
            double[][] T = new double[dimTrait][dimTrait];
            for (int i = 0; i < dimTrait; i++) {
                for (int j = 0; j < dimTrait; j++) {
                    T[i][j] = traitPrecision[i][j] * conditionalRootPrecision;
                }
            }
            System.err.println("Conditional root MVN precision = \n" + new Matrix(T));
            System.err.println("Conditional root MVN density = " + MultivariateNormalDistribution.logPdf(conditionalRootMean, new double[dimTrait], T, Math.log(MultivariateNormalDistribution.calculatePrecisionMatrixDeterminate(T)), 1.0));
        }
        if (integrateRoot) {
            // Integrate root trait out against rootPrior
            thisLogLikelihood += integrateLogLikelihoodAtRoot(conditionalRootMean, Ay, tmpM, traitPrecision, // Ay is destroyed
            conditionalRootPrecision);
        }
        if (DEBUG) {
            System.err.println("yAy = " + yAy);
            System.err.println("logLikelihood (before remainders) = " + thisLogLikelihood + " (should match conditional root MVN density when root not integrated out)");
        }
        logLikelihood += thisLogLikelihood;
    }
    logLikelihood += sumLogRemainders();
    if (DEBUG) {
        System.out.println("logLikelihood is " + logLikelihood);
    }
    if (DEBUG) {
        // Root trait is univariate!!!
        System.err.println("logLikelihood (final) = " + logLikelihood);
    //            checkViaLargeMatrixInversion();
    }
    if (DEBUG_PNAS) {
        checkLogLikelihood(logLikelihood, sumLogRemainders(), conditionalRootMean, conditionalRootPrecision, traitPrecision);
        for (int i = 0; i < logRemainderDensityCache.length; ++i) {
            if (logRemainderDensityCache[i] < -1E10) {
                System.err.println(logRemainderDensityCache[i] + " @ " + i);
            }
        }
    }
    // Should redraw internal node states when needed
    areStatesRedrawn = false;
    return logLikelihood;
}

Example 30 with Vector

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

the class SampledMultivariateTraitLikelihood method traitLogLikelihood.

private double traitLogLikelihood(double[] parentTrait, NodeRef node) {
    double logL = 0.0;
    double[] childTrait = treeModel.getMultivariateNodeTrait(node, traitName);
    if (parentTrait != null) {
        double time = getRescaledBranchLengthForPrecision(node);
        logL = diffusionModel.getLogLikelihood(parentTrait, childTrait, time);
        if (new Double(logL).isNaN()) {
            System.err.println("AbstractMultivariateTraitLikelihood: likelihood is undefined");
            System.err.println("time = " + time);
            System.err.println("parent trait value = " + new Vector(parentTrait));
            System.err.println("child trait value = " + new Vector(childTrait));
            double[][] precisionMatrix = diffusionModel.getPrecisionmatrix();
            if (precisionMatrix != null) {
                System.err.println("precision matrix = " + new Matrix(diffusionModel.getPrecisionmatrix()));
                if (diffusionModel.getPrecisionParameter() instanceof CompoundSymmetricMatrix) {
                    CompoundSymmetricMatrix csMatrix = (CompoundSymmetricMatrix) diffusionModel.getPrecisionParameter();
                    System.err.println("diagonals = " + new Vector(csMatrix.getDiagonals()));
                    System.err.println("off diagonal = " + csMatrix.getOffDiagonal());
                }
            }
        }
    }
    int childCount = treeModel.getChildCount(node);
    for (int i = 0; i < childCount; i++) {
        logL += traitLogLikelihood(childTrait, treeModel.getChild(node, i));
    }
    if (new Double(logL).isNaN()) {
        System.err.println("logL = " + logL);
        //            System.err.println(new Matrix(diffusionModel.getPrecisionmatrix()));
        System.exit(-1);
    }
    return logL;
}