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

Example 21 with StrictClockBranchRates

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

the class BeagleBranchLikelihood method main.

// END: LikelihoodColumn class
// ////////////
// ---TEST---//
// ////////////
public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("((SimSeq1:22.0,SimSeq2:22.0):12.0,(SimSeq3:23.1,SimSeq4:23.1):10.899999999999999);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 1.000);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogeneousBranchModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        Alignment alignment = simulator.simulate(false, false);
        System.out.println(alignment);
        BeagleTreeLikelihood btl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("BTL(homogeneous) = " + btl.getLogLikelihood());
        BeagleBranchLikelihood bbl = new BeagleBranchLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, freqModel, branchRateModel);
        int branchIndex = 4;
        System.out.println(bbl.getBranchLogLikelihood(branchIndex));
        bbl.finalizeBeagle();
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    } catch (Throwable e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 22 with StrictClockBranchRates

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

the class PartitionData method createClockRateModel.

public BranchRateModel createClockRateModel() {
    BranchRateModel branchRateModel = null;
    if (this.clockModelIndex == 0) {
        // Strict Clock
        Parameter rateParameter = new Parameter.Default(1, clockParameterValues[0]);
        branchRateModel = new StrictClockBranchRates(rateParameter);
    } else if (this.clockModelIndex == LRC_INDEX) {
        // Lognormal relaxed clock
        double numberOfBranches = 2 * (createTreeModel().getTaxonCount() - 1);
        Parameter rateCategoryParameter = new Parameter.Default(numberOfBranches);
        Parameter mean = new Parameter.Default(LogNormalDistributionModelParser.MEAN, 1, clockParameterValues[1]);
        Parameter stdev = new Parameter.Default(LogNormalDistributionModelParser.STDEV, 1, clockParameterValues[2]);
        // TODO: choose between log scale / real scale
        ParametricDistributionModel distributionModel = new LogNormalDistributionModel(mean, stdev, clockParameterValues[3], lrcParametersInRealSpace);
        branchRateModel = new // 
        DiscretizedBranchRates(// 
        createTreeModel(), // 
        rateCategoryParameter, // 
        distributionModel, // 
        1, // 
        false, // 
        Double.NaN, // randomizeRates
        true, // keepRates
        false, // cacheRates
        false);
    } else if (this.clockModelIndex == 2) {
        // Exponential relaxed clock
        double numberOfBranches = 2 * (createTreeModel().getTaxonCount() - 1);
        Parameter rateCategoryParameter = new Parameter.Default(numberOfBranches);
        Parameter mean = new Parameter.Default(DistributionModelParser.MEAN, 1, clockParameterValues[4]);
        ParametricDistributionModel distributionModel = new ExponentialDistributionModel(mean, clockParameterValues[5]);
        // branchRateModel = new DiscretizedBranchRates(createTreeModel(), rateCategoryParameter,
        // distributionModel, 1, false, Double.NaN);
        branchRateModel = new // 
        DiscretizedBranchRates(// 
        createTreeModel(), // 
        rateCategoryParameter, // 
        distributionModel, // 
        1, // 
        false, // 
        Double.NaN, // randomizeRates
        true, // keepRates
        false, // cacheRates
        false);
    } else if (this.clockModelIndex == 3) {
        // Inverse Gaussian
        double numberOfBranches = 2 * (createTreeModel().getTaxonCount() - 1);
        Parameter rateCategoryParameter = new Parameter.Default(numberOfBranches);
        Parameter mean = new Parameter.Default(InverseGaussianDistributionModelParser.MEAN, 1, clockParameterValues[6]);
        Parameter stdev = new Parameter.Default(InverseGaussianDistributionModelParser.STDEV, 1, clockParameterValues[7]);
        ParametricDistributionModel distributionModel = new InverseGaussianDistributionModel(mean, stdev, clockParameterValues[8], false);
        branchRateModel = new // 
        DiscretizedBranchRates(// 
        createTreeModel(), // 
        rateCategoryParameter, // 
        distributionModel, // 
        1, // 
        false, // 
        Double.NaN, // randomizeRates
        true, // keepRates
        false, // cacheRates
        false);
    } else {
        System.out.println("Not yet implemented");
    }
    return branchRateModel;
}

Example 23 with StrictClockBranchRates

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

the class ApproximatePoissonTreeLikelihoodTest method testOtherChildUpdates.

public void testOtherChildUpdates() throws IOException, Importer.ImportException, TreeUtils.MissingTaxonException {
    branchRateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    NewickImporter importer = new NewickImporter("((1:1.0,2:0.5,3:2.0):0.1,4:0.3,(7:0.2,8:0.1,9:0.3):0.2,6:0.01)");
    NewickImporter importer2 = new NewickImporter("((6:1,(4:1,((9:1,7:1):2,8:2):3.0):6):4,(1:1,(3:3,2:2):1):1);");
    tree = importer.importTree(null);
    treeModel = new BigFastTreeModel(importer2.importTree(null));
    CladeNodeModel cladeModel = new CladeNodeModel(tree, treeModel);
    BranchLengthProvider constrainedBranchLengthProvider = new ConstrainedTreeBranchLengthProvider(cladeModel);
    approximatePoissonTreeLikelihood = new ApproximatePoissonTreeLikelihood("approximateTreeLikelihood", 1, treeModel, branchRateModel, constrainedBranchLengthProvider);
    approximatePoissonTreeLikelihood.getLogLikelihood();
    NodeRef node = treeModel.getNode(1);
    NodeRef sibling = treeModel.getNode(9);
    NodeRef parent = treeModel.getParent(node);
    NodeRef grandparent = treeModel.getParent(parent);
    NodeRef root = treeModel.getRoot();
    CladeRef clade = cladeModel.getClade(parent);
    treeModel.beginTreeEdit();
    treeModel.removeChild(grandparent, parent);
    treeModel.removeChild(parent, sibling);
    treeModel.setNodeHeight(parent, treeModel.getNodeHeight(root) + 1);
    treeModel.setRoot(parent);
    // cladeModel.setRootNode(clade, parent);
    treeModel.addChild(parent, root);
    treeModel.addChild(grandparent, sibling);
    treeModel.endTreeEdit();
    double LL = approximatePoissonTreeLikelihood.getLogLikelihood();
    approximatePoissonTreeLikelihood.makeDirty();
    double newLL = approximatePoissonTreeLikelihood.getLogLikelihood();
    assertEquals(LL, newLL);
}

Example 24 with StrictClockBranchRates

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

the class ApproximatePoissonTreeLikelihoodTest method testRootAndRootChildUpdates.

// This tests updating when the root and rootchild2  change
// +- 6
// |
// +------| gp       +- 4 *
// |      |          |
// |      +----------| p       + 9
// |                 |     +---|
// |                 +-----|   + 7
// | r                     |
// |                       +--- 8
// |
// |+- 1
// +|
// | +----- 3
// +-|
// +--- 2
// To
// +-------------------- 4 *
// |
// |       +- 6
// |p      |
// |+------| gp                + 9
// ||      |               +---|
// ||      +---------------|   + 7
// +|                      |
// | r                    +--- 8
// |
// | +- 1
// +-|
// | +---- 3
// +-|
// +-- 2
// 
// 
public void testRootAndRootChildUpdates() throws IOException, Importer.ImportException, TreeUtils.MissingTaxonException {
    branchRateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    NewickImporter importer = new NewickImporter("((1:1.0,2:0.5,3:2.0):0.1,4:0.3,(7:0.2,8:0.1,9:0.3):0.2,6:0.01)");
    NewickImporter importer2 = new NewickImporter("((6:1,(4:1,((9:1,7:1):2,8:2):3.0):6):4,(1:1,(3:3,2:2):1):1);");
    tree = importer.importTree(null);
    treeModel = new BigFastTreeModel(importer2.importTree(null));
    CladeNodeModel cladeModel = new CladeNodeModel(tree, treeModel);
    BranchLengthProvider constrainedBranchLengthProvider = new ConstrainedTreeBranchLengthProvider(cladeModel);
    approximatePoissonTreeLikelihood = new ApproximatePoissonTreeLikelihood("approximateTreeLikelihood", 1, treeModel, branchRateModel, constrainedBranchLengthProvider);
    approximatePoissonTreeLikelihood.getLogLikelihood();
    NodeRef node = treeModel.getNode(9);
    NodeRef sibling = treeModel.getNode(1);
    NodeRef parent = treeModel.getParent(node);
    NodeRef grandparent = treeModel.getParent(parent);
    NodeRef root = treeModel.getRoot();
    CladeRef clade = cladeModel.getClade(node);
    treeModel.beginTreeEdit();
    treeModel.removeChild(grandparent, parent);
    treeModel.removeChild(parent, sibling);
    treeModel.setNodeHeight(parent, treeModel.getNodeHeight(root) + 1);
    treeModel.setRoot(parent);
    // cladeModel.setRootNode(clade, parent);
    treeModel.addChild(parent, root);
    treeModel.addChild(grandparent, sibling);
    treeModel.endTreeEdit();
    double LL = approximatePoissonTreeLikelihood.getLogLikelihood();
    approximatePoissonTreeLikelihood.makeDirty();
    double newLL = approximatePoissonTreeLikelihood.getLogLikelihood();
    assertEquals(LL, newLL);
}

Example 25 with StrictClockBranchRates

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

the class ApproximatePoissonTreeLikelihoodTest method testRootPolytomy.

public void testRootPolytomy() throws IOException, Importer.ImportException, TreeUtils.MissingTaxonException {
    branchRateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    NewickImporter importer = new NewickImporter("((1:1.0,2:1.0):1.0,3:1.0,4:1.0);");
    NewickImporter importer2 = new NewickImporter("(((1:1.0,2:1.0):1.0,3:1.1):0.1,4:1.0);");
    tree = importer.importTree(null);
    treeModel = new BigFastTreeModel(importer2.importTree(null));
    CladeNodeModel cladeModel = new CladeNodeModel(tree, treeModel);
    BranchLengthProvider constrainedBranchLengthProvider = new ConstrainedTreeBranchLengthProvider(cladeModel);
    approximatePoissonTreeLikelihood = new ApproximatePoissonTreeLikelihood("approximateTreeLikelihood", 1, treeModel, branchRateModel, constrainedBranchLengthProvider);
    approximatePoissonTreeLikelihood.getLogLikelihood();
    expectedLL = 0;
    double[] expectations = { 1d, 1d, 1d, 1.1, 1.1 };
    // time
    double[] mutations = { 1d, 1d, 1d, 1.0, 1 };
    for (int i = 0; i < expectations.length; i++) {
        PoissonDistribution p = new PoissonDistribution(expectations[i]);
        expectedLL += p.logPdf(mutations[i]);
    }
    NodeRef rootNode = treeModel.getRoot();
    NodeRef rootNode1 = treeModel.getNode(5);
    NodeRef tip3 = treeModel.getNode(2);
    CladeRef clade = cladeModel.getClade(rootNode);
    treeModel.beginTreeEdit();
    treeModel.removeChild(rootNode, rootNode1);
    treeModel.setRoot(rootNode1);
    treeModel.removeChild(rootNode1, tip3);
    treeModel.addChild(rootNode, tip3);
    treeModel.addChild(rootNode1, rootNode);
    treeModel.setNodeHeight(rootNode1, treeModel.getNodeHeight(rootNode) + 1);
    treeModel.endTreeEdit();
    // cladeModel.setRootNode(clade,rootNode1);
    double LL = approximatePoissonTreeLikelihood.getLogLikelihood();
    approximatePoissonTreeLikelihood.makeDirty();
    double newLL = approximatePoissonTreeLikelihood.getLogLikelihood();
    assertEquals(LL, newLL, 1E-13);
}