Examples with TreeDataLikelihood dr.evomodel.treedatalikelihood.TreeDataLikelihood used on opensource projects

Example 36 with TreeDataLikelihood

use of dr.evomodel.treedatalikelihood.TreeDataLikelihood in project beast-mcmc by beast-dev.

the class ContinuousDataLikelihoodDelegateTest method testLikelihoodDriftRelaxed.

public void testLikelihoodDriftRelaxed() {
    System.out.println("\nTest Likelihood using Drifted relaxed BM:");
    // Diffusion
    List<BranchRateModel> driftModels = new ArrayList<BranchRateModel>();
    ArbitraryBranchRates.BranchRateTransform transform = make(false, false, false);
    driftModels.add(new ArbitraryBranchRates(treeModel, new Parameter.Default("rate.1", new double[] { 0, 100, 200, 300, 400, 500, 600, 700, 800, 900 }), transform, false));
    driftModels.add(new ArbitraryBranchRates(treeModel, new Parameter.Default("rate.2", new double[] { 0, -100, 200, -300, 400, -500, 600, -700, 800, -900 }), transform, false));
    driftModels.add(new StrictClockBranchRates(new Parameter.Default("rate.3", new double[] { -2.0 })));
    DiffusionProcessDelegate diffusionProcessDelegate = new DriftDiffusionModelDelegate(treeModel, diffusionModel, driftModels);
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegate = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModel, rootPrior, rateTransformation, rateModel, false);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihood = new TreeDataLikelihood(likelihoodDelegate, treeModel, rateModel);
    testLikelihood("likelihoodDriftRelaxed", dataLikelihood);
    // Conditional moments (preorder)
    testConditionalMoments(dataLikelihood, likelihoodDelegate);
    // Conditional simulations
    MathUtils.setSeed(17890826);
    double[] expectedTraits = new double[] { -1.0, 2.0, 0.0, 2.843948876154644, 10.866053719140933, 3.467579698926694, 0.5, 12.000214659757933, 5.5, 2.0, 5.0, -8.0, 11.0, 1.0, -1.5, 1.0, 2.5, 4.0 };
    testConditionalSimulations(dataLikelihood, likelihoodDelegate, diffusionModel, dataModel, rootPrior, expectedTraits);
    // Fixed Root
    ContinuousDataLikelihoodDelegate likelihoodDelegateInf = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModel, rootPriorInf, rateTransformation, rateModel, true);
    TreeDataLikelihood dataLikelihoodInf = new TreeDataLikelihood(likelihoodDelegateInf, treeModel, rateModel);
    testLikelihood("likelihoodDriftRelaxedInf", dataLikelihoodInf);
    testConditionalMoments(dataLikelihoodInf, likelihoodDelegateInf);
}

Example 37 with TreeDataLikelihood

use of dr.evomodel.treedatalikelihood.TreeDataLikelihood in project beast-mcmc by beast-dev.

the class RepeatedMeasureFactorTest method testLikelihoodOU.

public void testLikelihoodOU() {
    System.out.println("\nTest Likelihood using full OU:");
    // Diffusion
    List<BranchRateModel> optimalTraitsModels = new ArrayList<BranchRateModel>();
    optimalTraitsModels.add(new StrictClockBranchRates(new Parameter.Default("rate.1", new double[] { 1.0 })));
    optimalTraitsModels.add(new StrictClockBranchRates(new Parameter.Default("rate.2", new double[] { 2.0 })));
    optimalTraitsModels.add(new StrictClockBranchRates(new Parameter.Default("rate.3", new double[] { -2.0 })));
    optimalTraitsModels.add(new StrictClockBranchRates(new Parameter.Default("rate.4", new double[] { 10.0 })));
    optimalTraitsModels.add(new StrictClockBranchRates(new Parameter.Default("rate.5", new double[] { 20.0 })));
    optimalTraitsModels.add(new StrictClockBranchRates(new Parameter.Default("rate.6", new double[] { -20.0 })));
    Parameter[] strengthOfSelectionParameters = new Parameter[6];
    strengthOfSelectionParameters[0] = new Parameter.Default(new double[] { 1.0, 0.1, 0.0, 0.0, 0.5, 2.0 });
    strengthOfSelectionParameters[1] = new Parameter.Default(new double[] { 0.1, 10., 0.0, 0.0, 0.0, 0.0 });
    strengthOfSelectionParameters[2] = new Parameter.Default(new double[] { 0.0, 0.0, 20., 0.3, 0.0, 0.0 });
    strengthOfSelectionParameters[3] = new Parameter.Default(new double[] { 0.0, 0.0, 0.3, 30., 3.0, 0.0 });
    strengthOfSelectionParameters[4] = new Parameter.Default(new double[] { 1.0, 0.0, 0.0, 3.0, 40., 0.0 });
    strengthOfSelectionParameters[5] = new Parameter.Default(new double[] { 0.0, 0.0, 0.5, 0.0, 0.0, 50. });
    MatrixParameter strengthOfSelectionMatrixParam = new MatrixParameter("strengthOfSelectionMatrix", strengthOfSelectionParameters);
    DiffusionProcessDelegate diffusionProcessDelegate = new OUDiffusionModelDelegate(treeModel, diffusionModel, optimalTraitsModels, new MultivariateElasticModel(strengthOfSelectionMatrixParam));
    // Rates
    ContinuousRateTransformation rateTransformation = new ContinuousRateTransformation.Default(treeModel, false, false);
    BranchRateModel rateModel = new DefaultBranchRateModel();
    // // Factor Model //// *****************************************************************************************
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegateFactors = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModelFactor, rootPrior, rateTransformation, rateModel, true);
    dataModelFactor.setLikelihoodDelegate(likelihoodDelegateFactors);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihoodFactors = new TreeDataLikelihood(likelihoodDelegateFactors, treeModel, rateModel);
    double logDatumLikelihoodFactor = getLogDatumLikelihood(dataModelFactor);
    double likelihoodFactorData = dataLikelihoodFactors.getLogLikelihood();
    double likelihoodFactorDiffusion = dataModelFactor.getLogLikelihood();
    assertEquals("likelihoodOUFactor", format.format(logDatumLikelihoodFactor), format.format(likelihoodFactorData + likelihoodFactorDiffusion));
    System.out.println("likelihoodOUFactor: " + format.format(logDatumLikelihoodFactor));
    // Simulation
    MathUtils.setSeed(17890826);
    double[] traitsFactors = getConditionalSimulations(dataLikelihoodFactors, likelihoodDelegateFactors, diffusionModel, dataModelFactor, rootPrior, treeModel, rateTransformation);
    System.err.println(new Vector(traitsFactors));
    // // Repeated Measures //// ************************************************************************************
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegateRepMea = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModelRepeatedMeasures, rootPrior, rateTransformation, rateModel, true);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihoodRepMea = new TreeDataLikelihood(likelihoodDelegateRepMea, treeModel, rateModel);
    double logDatumLikelihoodRepMea = getLogDatumLikelihood(dataLikelihoodRepMea);
    double likelihoodRepMeaDiffusion = dataLikelihoodRepMea.getLogLikelihood();
    assertEquals("likelihoodOURepMea", format.format(logDatumLikelihoodRepMea), format.format(likelihoodRepMeaDiffusion));
    System.out.println("likelihoodOURepMea: " + format.format(logDatumLikelihoodRepMea));
    // Simulation
    MathUtils.setSeed(17890826);
    double[] traitsRepMea = getConditionalSimulations(dataLikelihoodRepMea, likelihoodDelegateRepMea, diffusionModel, dataModelRepeatedMeasures, rootPrior, treeModel, rateTransformation);
    System.err.println(new Vector(traitsRepMea));
    // // Equal ? //// **********************************************************************************************
    assertEquals("likelihoodOURepFactor", format.format(likelihoodFactorData + likelihoodFactorDiffusion), format.format(likelihoodRepMeaDiffusion));
    for (int i = 0; i < traitsFactors.length; i++) {
        assertEquals(format.format(traitsRepMea[i]), format.format(traitsFactors[i]));
    }
    // // Repeated Measures Full //// *******************************************************************************
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegateRepMeaFull = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModelRepeatedMeasuresFull, rootPrior, rateTransformation, rateModel, true);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihoodRepMeaFull = new TreeDataLikelihood(likelihoodDelegateRepMeaFull, treeModel, rateModel);
    double logDatumLikelihoodRepMeaFull = getLogDatumLikelihood(dataLikelihoodRepMeaFull);
    double likelihoodRepMeaDiffusionFull = dataLikelihoodRepMeaFull.getLogLikelihood();
    assertEquals("likelihoodBMRepMea", format.format(logDatumLikelihoodRepMeaFull), format.format(likelihoodRepMeaDiffusionFull));
    System.out.println("likelihoodBMRepMeaFull: " + format.format(logDatumLikelihoodRepMeaFull));
}

Example 38 with TreeDataLikelihood

use of dr.evomodel.treedatalikelihood.TreeDataLikelihood in project beast-mcmc by beast-dev.

the class DiffusionGradientTest method testGradient.

private void testGradient(MultivariateDiffusionModel diffusionModel, DiffusionProcessDelegate diffusionProcessDelegate, ContinuousTraitPartialsProvider dataModel, ConjugateRootTraitPrior rootPrior, Parameter meanRoot, MatrixParameterInterface precision, Boolean wishart, MatrixParameterInterface attenuation, Parameter drift, MatrixParameterInterface samplingPrecision) {
    int dimLocal = rootPrior.getMean().length;
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegate = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModel, rootPrior, rateTransformation, rateModel, true);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihood = new TreeDataLikelihood(likelihoodDelegate, treeModel, rateModel);
    ProcessSimulationDelegate simulationDelegate = likelihoodDelegate.getPrecisionType() == PrecisionType.SCALAR ? new ConditionalOnTipsRealizedDelegate("trait", treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, likelihoodDelegate) : new MultivariateConditionalOnTipsRealizedDelegate("trait", treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, likelihoodDelegate);
    TreeTraitProvider traitProvider = new ProcessSimulation(dataLikelihood, simulationDelegate);
    dataLikelihood.addTraits(traitProvider.getTreeTraits());
    ProcessSimulationDelegate fullConditionalDelegate = new TipRealizedValuesViaFullConditionalDelegate("trait", treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, likelihoodDelegate);
    dataLikelihood.addTraits(new ProcessSimulation(dataLikelihood, fullConditionalDelegate).getTreeTraits());
    // Variance
    ContinuousDataLikelihoodDelegate cdld = (ContinuousDataLikelihoodDelegate) dataLikelihood.getDataLikelihoodDelegate();
    if (precision != null) {
        // Branch Specific
        ContinuousProcessParameterGradient traitGradient = new ContinuousProcessParameterGradient(rootPrior.getMean().length, treeModel, cdld, new ArrayList<>(Arrays.asList(DerivationParameter.WRT_VARIANCE)));
        BranchSpecificGradient branchSpecificGradient = new BranchSpecificGradient("trait", dataLikelihood, cdld, traitGradient, precision);
        GradientWrtPrecisionProvider gPPBranchSpecific = new GradientWrtPrecisionProvider.BranchSpecificGradientWrtPrecisionProvider(branchSpecificGradient);
        // Correlation Gradient Branch Specific
        CorrelationPrecisionGradient gradientProviderBranchSpecific = new CorrelationPrecisionGradient(gPPBranchSpecific, dataLikelihood, precision);
        double[] gradientAnalyticalBS = testOneGradient(gradientProviderBranchSpecific);
        // Diagonal Gradient Branch Specific
        DiagonalPrecisionGradient gradientDiagonalProviderBS = new DiagonalPrecisionGradient(gPPBranchSpecific, dataLikelihood, precision);
        double[] gradientDiagonalAnalyticalBS = testOneGradient(gradientDiagonalProviderBS);
        if (wishart) {
            // Wishart Statistic
            WishartStatisticsWrapper wishartStatistics = new WishartStatisticsWrapper("wishart", "trait", dataLikelihood, cdld);
            GradientWrtPrecisionProvider gPPWiwhart = new GradientWrtPrecisionProvider.WishartGradientWrtPrecisionProvider(wishartStatistics);
            // Correlation Gradient
            CorrelationPrecisionGradient gradientProviderWishart = new CorrelationPrecisionGradient(gPPWiwhart, dataLikelihood, precision);
            String sW = gradientProviderWishart.getReport();
            System.err.println(sW);
            double[] gradientAnalyticalW = parseGradient(sW, "analytic");
            assertEquals("Sizes", gradientAnalyticalW.length, gradientAnalyticalBS.length);
            for (int k = 0; k < gradientAnalyticalW.length; k++) {
                assertEquals("gradient correlation k=" + k, gradientAnalyticalW[k], gradientAnalyticalBS[k], delta);
            }
            // Diagonal Gradient
            DiagonalPrecisionGradient gradientDiagonalProviderW = new DiagonalPrecisionGradient(gPPWiwhart, dataLikelihood, precision);
            String sDiagW = gradientDiagonalProviderW.getReport();
            System.err.println(sDiagW);
            double[] gradientDiagonalAnalyticalW = parseGradient(sDiagW, "analytic");
            assertEquals("Sizes", gradientDiagonalAnalyticalW.length, gradientDiagonalAnalyticalBS.length);
            for (int k = 0; k < gradientDiagonalAnalyticalW.length; k++) {
                assertEquals("gradient diagonal k=" + k, gradientDiagonalAnalyticalW[k], gradientDiagonalAnalyticalBS[k], delta);
            }
        }
    }
    // Diagonal Attenuation Gradient Branch Specific
    if (attenuation != null) {
        ContinuousProcessParameterGradient traitGradientAtt = new ContinuousProcessParameterGradient(dimLocal, treeModel, cdld, new ArrayList<>(Arrays.asList(DerivationParameter.WRT_DIAGONAL_SELECTION_STRENGTH)));
        BranchSpecificGradient branchSpecificGradientAtt = new BranchSpecificGradient("trait", dataLikelihood, cdld, traitGradientAtt, attenuation);
        AbstractDiffusionGradient.ParameterDiffusionGradient gABranchSpecific = createDiagonalAttenuationGradient(branchSpecificGradientAtt, dataLikelihood, attenuation);
        testOneGradient(gABranchSpecific);
    }
    // WRT root mean
    boolean sameRoot = (drift == meanRoot);
    ContinuousProcessParameterGradient traitGradientRoot = new ContinuousProcessParameterGradient(dimLocal, treeModel, cdld, new ArrayList<>(Arrays.asList(sameRoot ? DerivationParameter.WRT_CONSTANT_DRIFT_AND_ROOT_MEAN : DerivationParameter.WRT_ROOT_MEAN)));
    BranchSpecificGradient branchSpecificGradientRoot = new BranchSpecificGradient("trait", dataLikelihood, cdld, traitGradientRoot, meanRoot);
    AbstractDiffusionGradient.ParameterDiffusionGradient gRootBranchSpecific = createDriftGradient(branchSpecificGradientRoot, dataLikelihood, meanRoot);
    testOneGradient(gRootBranchSpecific);
    // Drift Gradient Branch Specific
    if (drift != null && !sameRoot) {
        ContinuousProcessParameterGradient traitGradientDrift = new ContinuousProcessParameterGradient(dimLocal, treeModel, cdld, new ArrayList<>(Arrays.asList(DerivationParameter.WRT_CONSTANT_DRIFT)));
        BranchSpecificGradient branchSpecificGradientDrift = new BranchSpecificGradient("trait", dataLikelihood, cdld, traitGradientDrift, drift);
        AbstractDiffusionGradient.ParameterDiffusionGradient gDriftBranchSpecific = createDriftGradient(branchSpecificGradientDrift, dataLikelihood, drift);
        testOneGradient(gDriftBranchSpecific);
    }
    // Sampling Precision
    if (samplingPrecision != null) {
        ContinuousTraitGradientForBranch.SamplingVarianceGradient traitGradientSampling = new ContinuousTraitGradientForBranch.SamplingVarianceGradient(dimLocal, treeModel, likelihoodDelegate, (ModelExtensionProvider.NormalExtensionProvider) dataModel);
        BranchSpecificGradient branchSpecificGradientSampling = new BranchSpecificGradient("trait", dataLikelihood, cdld, traitGradientSampling, samplingPrecision);
        GradientWrtPrecisionProvider gPPBranchSpecificSampling = new GradientWrtPrecisionProvider.BranchSpecificGradientWrtPrecisionProvider(branchSpecificGradientSampling);
        // Correlation Gradient Branch Specific
        // CorrelationPrecisionGradient gradientProviderBranchSpecificSampling = new CorrelationPrecisionGradient(gPPBranchSpecificSampling, dataLikelihood, samplingPrecision);
        // 
        // testOneGradient(gradientProviderBranchSpecificSampling);
        // Diagonal Gradient Branch Specific
        DiagonalPrecisionGradient gradientDiagonalProviderBSSampling = new DiagonalPrecisionGradient(gPPBranchSpecificSampling, dataLikelihood, samplingPrecision);
        testOneGradient(gradientDiagonalProviderBSSampling);
    }
}

Example 39 with TreeDataLikelihood

use of dr.evomodel.treedatalikelihood.TreeDataLikelihood in project beast-mcmc by beast-dev.

the class BranchSpecificGradientTest method testRateGradient.

public void testRateGradient() {
    System.out.println("\nTest Likelihood using vanilla BM:");
    // Diffusion
    List<BranchRateModel> driftModels = new ArrayList<BranchRateModel>();
    driftModels.add(new StrictClockBranchRates(new Parameter.Default("rate.1", new double[] { 1.0 })));
    driftModels.add(new StrictClockBranchRates(new Parameter.Default("rate.2", new double[] { 2.0 })));
    driftModels.add(new StrictClockBranchRates(new Parameter.Default("rate.3", new double[] { -2.0 })));
    DiffusionProcessDelegate diffusionProcessDelegate = new DriftDiffusionModelDelegate(treeModel, diffusionModel, driftModels);
    // Rates
    ArbitraryBranchRates.BranchRateTransform transform = make(false, true, false, null, null);
    Parameter branchRates = new Parameter.Default(new double[] { 10, 9, 8, 7, 6, 5, 4, 3, 2, 1 });
    ArbitraryBranchRates rateModel = new ArbitraryBranchRates(treeModel, branchRates, transform, false);
    ContinuousRateTransformation rateTransformation = new ContinuousRateTransformation.Default(treeModel, false, false);
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegate = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModel, rootPrior, rateTransformation, rateModel, false);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihood = new TreeDataLikelihood(likelihoodDelegate, treeModel, rateModel);
    // Gradient (Rates)
    BranchRateGradient branchGradient1 = new BranchRateGradient("trait", dataLikelihood, likelihoodDelegate, branchRates);
    double[] gradient1 = branchGradient1.getGradientLogDensity();
    // Gradient (Specific)
    ContinuousTraitGradientForBranch.RateGradient traitGradient = new ContinuousTraitGradientForBranch.RateGradient(dimTrait, treeModel, rateModel);
    BranchSpecificGradient branchGradient2 = new BranchSpecificGradient("trait", dataLikelihood, likelihoodDelegate, traitGradient, branchRates);
    double[] gradient2 = branchGradient2.getGradientLogDensity();
    double[] numericalGradient = branchGradient1.getNumericalGradient();
    System.err.println("\tGradient with rate method    = " + new dr.math.matrixAlgebra.Vector(gradient1));
    System.err.println("\tGradient with general method = " + new dr.math.matrixAlgebra.Vector(gradient2));
    System.err.println("\tNumerical gradient           = " + new dr.math.matrixAlgebra.Vector(numericalGradient));
    assertEquals("length", gradient1.length, gradient2.length);
    for (int i = 0; i < gradient1.length; i++) {
        assertEquals("numeric " + i, gradient1[i], numericalGradient[i], 1E-4);
    }
    for (int i = 0; i < gradient1.length; i++) {
        assertEquals("gradient " + i, format.format(gradient1[i]), format.format(gradient2[i]));
    }
}

Example 40 with TreeDataLikelihood

use of dr.evomodel.treedatalikelihood.TreeDataLikelihood in project beast-mcmc by beast-dev.

the class ContinuousDataLikelihoodDelegateTest method testLikelihoodDriftRelaxedFactor.

public void testLikelihoodDriftRelaxedFactor() {
    System.out.println("\nTest Likelihood using drifted Relaxed BM and factor:");
    // Diffusion
    List<BranchRateModel> driftModels = new ArrayList<BranchRateModel>();
    ArbitraryBranchRates.BranchRateTransform transform = make(false, false, false);
    driftModels.add(new ArbitraryBranchRates(treeModel, new Parameter.Default("rate.1", new double[] { 0, 10, 20, 30, 40, 40, 30, 20, 10, 0 }), transform, false));
    driftModels.add(new StrictClockBranchRates(new Parameter.Default("rate.2", new double[] { 0.0 })));
    DiffusionProcessDelegate diffusionProcessDelegate = new DriftDiffusionModelDelegate(treeModel, diffusionModelFactor, driftModels);
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegateFactors = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModelFactor, rootPriorFactor, rateTransformation, rateModel, false);
    dataModelFactor.setLikelihoodDelegate(likelihoodDelegateFactors);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihoodFactors = new TreeDataLikelihood(likelihoodDelegateFactors, treeModel, rateModel);
    testLikelihood("likelihoodDriftRelaxedFactor", dataModelFactor, dataLikelihoodFactors);
    // Conditional simulations
    MathUtils.setSeed(17890826);
    double[] expectedTraits = new double[] { 0.21992781609528125, 1.271388273711557, 0.40761548539751596, 1.3682648770877144, 0.6599021787120436, 1.2830636141108613, 1.1488658943588324, 1.472103688153391, 0.8971632986744889, 1.20748933414854, 1.603739823726808, 1.4761482401796842 };
    testConditionalSimulations(dataLikelihoodFactors, likelihoodDelegateFactors, diffusionModelFactor, dataModelFactor, rootPriorFactor, expectedTraits);
}