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

Example 1 with StrictClockBranchRates

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

the class BeagleTreeLikelihood method main.

public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 1: simulateOnePartition");
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        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);
        Parameter kappa2 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        HKY hky2 = new HKY(kappa2, freqModel);
        HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        substitutionModels.add(hky2);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        Parameter epochTimes = new Parameter.Default(1, 20);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 0.001);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogenousBranchSubstitutionModel, //
        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);
        BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
        nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 2 with StrictClockBranchRates

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

the class AncestralSequenceAnnotator method processTree.

private Tree processTree(Tree tree) {
    // Remake tree to fix node ordering - Marc
    GammaSiteRateModel siteModel = loadSiteModel(tree);
    SimpleAlignment alignment = new SimpleAlignment();
    alignment.setDataType(siteModel.getSubstitutionModel().getDataType());
    if (siteModel.getSubstitutionModel().getDataType().getClass().equals(Codons.class)) {
        // System.out.println("trololo");
        alignment.setDataType(Nucleotides.INSTANCE);
    }
    // System.out.println("BOO BOO " + siteModel.getSubstitutionModel().getDataType().getClass().getName()+"\t" + Codons.UNIVERSAL.getClass().getName() + "\t" + alignment.getDataType().getClass().getName());
    // Get sequences
    String[] sequence = new String[tree.getNodeCount()];
    for (int i = 0; i < tree.getNodeCount(); i++) {
        NodeRef node = tree.getNode(i);
        sequence[i] = (String) tree.getNodeAttribute(node, SEQ_STRING);
        if (tree.isExternal(node)) {
            Taxon taxon = tree.getNodeTaxon(node);
            alignment.addSequence(new Sequence(taxon, sequence[i]));
        // System.out.println("seq " + sequence[i]);
        }
    }
    // Make evolutionary model
    BranchRateModel rateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    FlexibleTree flexTree;
    if (siteModel.getSubstitutionModel().getDataType().getClass().equals(Codons.class)) {
        ConvertAlignment convertAlignment = new ConvertAlignment(siteModel.getSubstitutionModel().getDataType(), ((Codons) siteModel.getSubstitutionModel().getDataType()).getGeneticCode(), alignment);
        flexTree = sampleTree(tree, convertAlignment, siteModel, rateModel);
    // flexTree = sampleTree(tree, alignment, siteModel, rateModel);
    } else {
        flexTree = sampleTree(tree, alignment, siteModel, rateModel);
    }
    introduceGaps(flexTree, tree);
    return flexTree;
}

Example 3 with StrictClockBranchRates

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

the class MicrosatelliteSamplerTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    MicrosatelliteSamplerTreeModel mstModel = (MicrosatelliteSamplerTreeModel) xo.getChild(MicrosatelliteSamplerTreeModel.class);
    MicrosatelliteModel microsatelliteModel = (MicrosatelliteModel) xo.getChild(MicrosatelliteModel.class);
    BranchRateModel branchRateModel;
    Object cxo = xo.getChild(BranchRateModel.class);
    if (xo.getChild(BranchRateModel.class) != null) {
        branchRateModel = (BranchRateModel) cxo;
        System.out.println("BranchRateModel provided to MicrosatelliteSamplerTreeLikelihood");
    } else if (xo.hasChildNamed(MUTATION_RATE)) {
        Parameter muRate = (Parameter) xo.getElementFirstChild(MUTATION_RATE);
        branchRateModel = new StrictClockBranchRates(muRate);
        System.out.println("mutation rate provided to MicrosatelliteSamplerTreeLikelihood");
    } else {
        branchRateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    }
    return new MicrosatelliteSamplerTreeLikelihood(mstModel, microsatelliteModel, branchRateModel);
}

Example 4 with StrictClockBranchRates

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

the class StrictClockTest method testStrictClock.

public void testStrictClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 380.0, 0, 38000.0);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    TreeIntervals intervalList = new TreeIntervals(treeModel, null, null);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(intervalList, constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter rateParameter = new Parameter.Default(StrictClockBranchRates.RATE, 2.3E-5, 0, 100.0);
    StrictClockBranchRates branchRateModel = new StrictClockBranchRates(rateParameter);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 100.0);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    // siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    // treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(rateParameter, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter allInternalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(true, true, false);
    operator = new UpDownOperator(new Scalable[] { new Scalable.Default(rateParameter) }, new Scalable[] { new Scalable.Default(allInternalHeights) }, 0.75, 3.0, AdaptationMode.ADAPTATION_ON);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = ((DefaultTreeModel) treeModel).getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(((DefaultTreeModel) treeModel), 15.0, 1.0, true, false, false, false, AdaptationMode.ADAPTATION_ON, AdaptableMCMCOperator.DEFAULT_ADAPTATION_TARGET);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 15.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 3.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 3.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    // CompoundLikelihood
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(coalescent);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    likelihoods.clear();
    likelihoods.add(treeLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    posterior.setId(CompoundLikelihoodParser.POSTERIOR);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 500, false);
    loggers[0].add(posterior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(rateParameter);
    loggers[0].add(popSize);
    loggers[0].add(kappa);
    loggers[0].add(coalescent);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[1].add(posterior);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(rateParameter);
    loggers[1].add(coalescent);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("RandomLocalClockTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    // <expectation name="posterior" value="-3928.71"/>
    // <expectation name="clock.rate" value="8.04835E-4"/>
    // <expectation name="constant.popSize" value="37.3762"/>
    // <expectation name="hky.kappa" value="18.2782"/>
    // <expectation name="treeModel.rootHeight" value="69.0580"/>
    // <expectation name="treeLikelihood" value="-3856.59"/>
    // <expectation name="coalescent" value="-72.1285"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -3928.71);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -3856.59);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 69.0580);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 18.2782);
    TraceCorrelation rateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(StrictClockBranchRates.RATE));
    assertExpectation(StrictClockBranchRates.RATE, rateStats, 8.04835E-4);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 37.3762);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, -72.1285);
}

Example 5 with StrictClockBranchRates

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

the class ApproximatePoissonTreeLikelihoodTest method setUp.

public void setUp() throws Exception {
    super.setUp();
    branchRateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    NewickImporter importer = new NewickImporter("((1:1.0,2:1.0,3:1.0):1.0,4:1.0);");
    NewickImporter importer2 = new NewickImporter("(((1:1.0,2:1.0):0.1,3:1.1):1.0,4:1.0);");
    tree = importer.importTree(null);
    treeModel = new BigFastTreeModel(importer2.importTree(null));
    cladeModel = new CladeNodeModel(tree, treeModel);
    BranchLengthProvider constrainedBranchLengthProvider = new ConstrainedTreeBranchLengthProvider(cladeModel);
    approximatePoissonTreeLikelihood = new ApproximatePoissonTreeLikelihood("approximateTreeLikelihood", 1, treeModel, branchRateModel, constrainedBranchLengthProvider);
    expectedLL = 0;
    double[] expectations = { 1d, 1d, 1.1, 2d, 0.1 };
    // time
    double[] mutations = { 1d, 1d, 1.0, 2d, 0 };
    for (int i = 0; i < expectations.length; i++) {
        PoissonDistribution p = new PoissonDistribution(expectations[i]);
        expectedLL += p.logPdf(mutations[i]);
    }
    approximatePoissonTreeLikelihood.getLogLikelihood();
}