Examples with ArbitraryBranchRates dr.evomodel.branchratemodel.ArbitraryBranchRates used on opensource projects

Example 11 with ArbitraryBranchRates

use of dr.evomodel.branchratemodel.ArbitraryBranchRates in project beast-mcmc by beast-dev.

the class ArbitraryBranchRatesParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
    XMLObject cxo = xo.getChild(RATES);
    Parameter rateCategoryParameter = (Parameter) cxo.getChild(Parameter.class);
    boolean centerAtOne = xo.getAttribute(CENTER_AT_ONE, true);
    boolean randomizeRates = xo.getAttribute(RANDOMIZE_RATES, false);
    if (centerAtOne && randomizeRates) {
        throw new XMLParseException("Cannot centerAtOne and randomize the starting rates");
    }
    final int numBranches = tree.getNodeCount() - 1;
    if (rateCategoryParameter.getDimension() > 1 && (rateCategoryParameter.getDimension() != numBranches)) {
        throw new XMLParseException("Incorrect number of rate parameters");
    }
    if (rateCategoryParameter.getDimension() != numBranches) {
        rateCategoryParameter.setDimension(numBranches);
    }
    Logger.getLogger("dr.evomodel").info("\nUsing an scaled mixture of normals model.");
    Logger.getLogger("dr.evomodel").info("  rates = " + rateCategoryParameter.getDimension());
    Logger.getLogger("dr.evomodel").info("  NB: Make sure you have a prior on " + rateCategoryParameter.getId());
    ArbitraryBranchRates.BranchRateTransform transform = parseTransform(xo);
    if (randomizeRates) {
        for (int i = 0; i < rateCategoryParameter.getDimension(); i++) {
            rateCategoryParameter.setValue(i, MathUtils.uniform(0, 10));
        }
    }
    return new ArbitraryBranchRates(tree, rateCategoryParameter, transform, centerAtOne);
}

Example 12 with ArbitraryBranchRates

use of dr.evomodel.branchratemodel.ArbitraryBranchRates in project beast-mcmc by beast-dev.

the class AutoCorrelatedBranchRatesDistributionParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    ArbitraryBranchRates branchRates = (ArbitraryBranchRates) xo.getChild(ArbitraryBranchRates.class);
    ParametricMultivariateDistributionModel distribution = (ParametricMultivariateDistributionModel) xo.getChild(ParametricMultivariateDistributionModel.class);
    AutoCorrelatedBranchRatesDistribution.BranchVarianceScaling scaling = parseScaling(xo);
    boolean log = xo.getAttribute(LOG, false);
    return new AutoCorrelatedBranchRatesDistribution(xo.getId(), branchRates, distribution, scaling, log);
}

Example 13 with ArbitraryBranchRates

use of dr.evomodel.branchratemodel.ArbitraryBranchRates in project beast-mcmc by beast-dev.

the class ContinuousDataLikelihoodDelegateTest method testLikelihoodDiagonalOURelaxed.

public void testLikelihoodDiagonalOURelaxed() {
    System.out.println("\nTest Likelihood using Diagonal OU Relaxed:");
    // Diffusion
    List<BranchRateModel> optimalTraitsModels = new ArrayList<BranchRateModel>();
    ArbitraryBranchRates.BranchRateTransform transform = make(false, false, false);
    optimalTraitsModels.add(new ArbitraryBranchRates(treeModel, new Parameter.Default("rate.1", new double[] { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }), transform, false));
    optimalTraitsModels.add(new ArbitraryBranchRates(treeModel, new Parameter.Default("rate.2", new double[] { 0, -1, 2, -3, 4, -5, 6, -7, 8, -9 }), transform, false));
    optimalTraitsModels.add(new StrictClockBranchRates(new Parameter.Default("rate.3", new double[] { -2.0 })));
    DiagonalMatrix strengthOfSelectionMatrixParam = new DiagonalMatrix(new Parameter.Default(new double[] { 1.0, 100.0, 100.0 }));
    DiffusionProcessDelegate diffusionProcessDelegate = new OUDiffusionModelDelegate(treeModel, diffusionModel, optimalTraitsModels, new MultivariateElasticModel(strengthOfSelectionMatrixParam));
    // CDL
    ContinuousDataLikelihoodDelegate likelihoodDelegate = new ContinuousDataLikelihoodDelegate(treeModel, diffusionProcessDelegate, dataModel, rootPrior, rateTransformation, rateModel, false);
    // Likelihood Computation
    TreeDataLikelihood dataLikelihood = new TreeDataLikelihood(likelihoodDelegate, treeModel, rateModel);
    testLikelihood("likelihoodDiagonalOURelaxed", dataLikelihood);
    // Conditional moments (preorder)
    testConditionalMoments(dataLikelihood, likelihoodDelegate);
    // Conditional simulations
    MathUtils.setSeed(17890826);
    double[] expectedTraits = new double[] { -1.0, 2.0, 0.0, 1.811803424441062, 0.6837595819961084, -1.0607909328094163, 0.5, 3.8623525502275142, 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("likelihoodDiagonalOURelaxedInf", dataLikelihoodInf);
    testConditionalMoments(dataLikelihoodInf, likelihoodDelegateInf);
}

Example 14 with ArbitraryBranchRates

use of dr.evomodel.branchratemodel.ArbitraryBranchRates 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 15 with ArbitraryBranchRates

use of dr.evomodel.branchratemodel.ArbitraryBranchRates 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]));
    }
}