Examples with DefaultTreeModel dr.evomodel.tree.DefaultTreeModel used on opensource projects

Example 31 with DefaultTreeModel

use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.

the class RateExchangeOperatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    DefaultTreeModel treeModel = (DefaultTreeModel) xo.getChild(DefaultTreeModel.class);
    double weight = xo.getDoubleAttribute(MCMCOperator.WEIGHT);
    boolean swapRates = xo.getBooleanAttribute(SWAP_RATES);
    boolean swapTraits = xo.getBooleanAttribute(SWAP_TRAITS);
    boolean swapAtRoot = xo.getBooleanAttribute(SWAP_AT_ROOT);
    boolean moveHeight = xo.getBooleanAttribute(MOVE_HEIGHT);
    return new RateExchangeOperator(treeModel, weight, swapRates, swapTraits, swapAtRoot, moveHeight);
}

Example 32 with DefaultTreeModel

use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.

the class RandomLocalClockTestProblem method testRandomLocalClock.

public void testRandomLocalClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 0.077, 0, Double.POSITIVE_INFINITY);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    TreeIntervals intervalList = new TreeIntervals(treeModel, null, null);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(intervalList, constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter ratesParameter = new Parameter.Default(RandomLocalClockModelParser.RATES, 10, 1);
    Parameter rateIndicatorParameter = new Parameter.Default(RandomLocalClockModelParser.RATE_INDICATORS, 10, 1);
    Parameter meanRateParameter = new Parameter.Default(RandomLocalClockModelParser.CLOCK_RATE, 1, 1.0);
    RandomLocalClockModel branchRateModel = new RandomLocalClockModel(treeModel, meanRateParameter, rateIndicatorParameter, ratesParameter, false, 0.5);
    SumStatistic rateChanges = new SumStatistic("rateChangeCount", true, null);
    rateChanges.addStatistic(rateIndicatorParameter);
    RateStatistic meanRate = new RateStatistic("meanRate", treeModel, branchRateModel, true, true, RateStatisticParser.MEAN);
    RateStatistic coefficientOfVariation = new RateStatistic(RateStatisticParser.COEFFICIENT_OF_VARIATION, treeModel, branchRateModel, true, true, RateStatisticParser.COEFFICIENT_OF_VARIATION);
    RateCovarianceStatistic covariance = new RateCovarianceStatistic("covariance", treeModel, branchRateModel);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, Double.POSITIVE_INFINITY);
    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(ratesParameter, 0.75);
    operator.setWeight(10.0);
    schedule.addOperator(operator);
    operator = new BitFlipOperator(rateIndicatorParameter, 15.0, true);
    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, 0.0077, 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
    OneOnXPrior likelihood1 = new OneOnXPrior();
    likelihood1.addData(popSize);
    OneOnXPrior likelihood2 = new OneOnXPrior();
    likelihood2.addData(kappa);
    DistributionLikelihood likelihood3 = new DistributionLikelihood(new GammaDistribution(0.5, 2.0), 0.0);
    likelihood3.addData(ratesParameter);
    DistributionLikelihood likelihood4 = new DistributionLikelihood(new PoissonDistribution(1.0), 0.0);
    likelihood4.addData(rateChanges);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(likelihood1);
    likelihoods.add(likelihood2);
    likelihoods.add(likelihood3);
    likelihoods.add(likelihood4);
    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, 1000, false);
    loggers[0].add(posterior);
    loggers[0].add(prior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(kappa);
    // loggers[0].add(meanRate);
    loggers[0].add(rateChanges);
    loggers[0].add(coefficientOfVariation);
    loggers[0].add(covariance);
    loggers[0].add(popSize);
    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(meanRate);
    loggers[1].add(rateChanges);
    // 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="-1818.26"/>
    // <expectation name="prior" value="-2.70143"/>
    // <expectation name="likelihood" value="-1815.56"/>
    // <expectation name="treeModel.rootHeight" value="6.363E-2"/>
    // <expectation name="constant.popSize" value="9.67405E-2"/>
    // <expectation name="hky.kappa" value="30.0394"/>
    // <expectation name="coefficientOfVariation" value="7.02408E-2"/>
    // covariance 0.47952
    // <expectation name="rateChangeCount" value="0.40786"/>
    // <expectation name="coalescent" value="7.29521"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -1818.26);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.PRIOR));
    assertExpectation(CompoundLikelihoodParser.PRIOR, likelihoodStats, -2.70143);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.56);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 6.363E-2);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 30.0394);
    TraceCorrelation rateChangeStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("rateChangeCount"));
    assertExpectation("rateChangeCount", rateChangeStats, 0.40786);
    TraceCorrelation coefficientOfVariationStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(RateStatisticParser.COEFFICIENT_OF_VARIATION));
    assertExpectation(RateStatisticParser.COEFFICIENT_OF_VARIATION, coefficientOfVariationStats, 7.02408E-2);
    TraceCorrelation covarianceStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("covariance"));
    assertExpectation("covariance", covarianceStats, 0.47952);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 9.67405E-2);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, 7.29521);
}

Example 33 with DefaultTreeModel

use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.

the class BigFastTreeTreeIntervalsTest method testCompareIntervals.

public void testCompareIntervals() throws TreeUtils.MissingTaxonException, IOException, Importer.ImportException {
    NewickImporter importer = new NewickImporter("(Lishui/LS557/2020:0,((Netherlands/Utrecht_10015/2020:0.00006795400000000001,USA/IL-NM073/2020:0.00006799599999999999):0.000033976,England/LOND-D604F/2020:0.000101963):0.000033968,Guangdong/2020XN4459-P0041/2020:0.000000005,(Portugal/PT0063/2020:0,(Spain/Zaragoza2486/2020:0.000102605,Scotland/CVR746/2020:0.000000005,Spain/COV000882/2020:0.000067956,Colombia/INS-79253/2020:0.000101944,Uruguay/UY-4/2020:0.000031515):0.000033979,(Spain/CastillaLaMancha201329/2020:0.000000005,Netherlands/NoordHolland_10011/2020:0.000033987):0.00006799,England/LIVE-9CE87/2020:0.00013727299999999998,Spain/Granada-COV002916/2020:0.000033979999999999997):0.000033968,((USA/VI-CDC-3705/2020:0.000000005,Australia/VIC229/2020:0,USA/MA-MGH-00063/2020:0,(USA/WA-S41/2020:0.000068895,USA/WA-UW114/2020:0.000067978,USA/WA-UW17/2020:0.000000005,(USA/WA-S582/2020:0,USA/WA-UW-1682/2020:0.000000005,USA/WA-S994/2020:0.000101934):0.000033955,USA/WA-S121/2020:0.000000005,USA/WA-S154/2020:0.000067982,USA/WA-UW37/2020:0,USA/WA-S321/2020:0,USA/WA-S445/2020:0,USA/WA-S512/2020:0,USA/WA-S33/2020:0.000033979,Canada/BC_6981299/2020:0.000033972,USA/WA-UW-1294/2020:0.000033972,USA/WA-UW-2247/2020:0.000033988,Australia/VIC140/2020:0.000033984,USA/WA-UW61/2020:0.000033972,Canada/BC_8606204/2020:0.000166157,(USA/WA-S734/2020:0,USA/WA-S844/2020:0.000033983):0.000067965,(USA/WA-S1191/2020:0.000067947,USA/WA-S951/2020:0.000101914):0.000095803,Australia/NSW99/2020:0.000101953,(USA/WA-S317/2020:0.000000005,USA/WA-S721/2020:0.00003397):0.00010195700000000001,USA/WA-UW139/2020:0.000135916,USA/WA-S572/2020:0.000033979999999999997,USA/WA-S279/2020:0.000033972,USA/WA-UW28/2020:0.000034002,USA/WA-S114/2020:0.000033969,(USA/WA-S852/2020:0.000203899,(USA/WA-S568/2020:0,USA/WA-S791/2020:0.00006794599999999999):0.000033983):0.000101964,USA/WA-S842/2020:0.000067951):0.000033986,Singapore/302/2020:0.000101947):0.00016677,(((USA/IL-NM0112/2020:0.00003397,USA/IL-NM053/2020:0.000034229,USA/IL-NM059/2020:0.000101967):0.00003397,USA/WI-UW-218/2020:0.000033995):0.000030539,(USA/UT-QDX-63/2020:0,USA/CA-QDX-111/2020:0,USA/TX-HMH0427/2020:0.000203861):0.000101955):0.00023787300000000002):0.000033959,(((Scotland/CVR3203/2020:0.000000005,Scotland/CVR2246/2020:0.000000005,Scotland/GCVR-1714B2/2020:0.000033975999999999995,Scotland/CVR3514/2020:0.000068628):0.000067954,Australia/NT08/2020:0.000034000999999999995):0.00003397,Spain/COV001440/2020:0,Spain/Alcaniz2449/2020:0.000068985,Spain/COV001548/2020:0,USA/WI-WSLH-200057/2020:0.000000005,Spain/Valencia6/2020:0.0000343,Spain/Granada-COV002944/2020:0.000000005,Spain/COV001929/2020:0.000000005,Spain/COV002049/2020:0.000000005,(Spain/Valencia59/2020:0,Spain/Valencia306/2020:0.000000005):0.000033996,(Spain/COV001117/2020:0.00010265,Spain/COV002055/2020:0.000000005,England/20126000104/2020:0.00006758400000000001):0.000033997,Spain/COV001576/2020:0.000000005,Chile/Santiago-1/2020:0.000000005,Spain/COV000721/2020:0.000000005,(Spain/COV001575/2020:0,Spain/COV001505/2020:0):0.000067968,Spain/Madrid_H12_28/2020:0.000067957,Spain/COV001568/2020:0.000033975,England/CAMB-83357/2020:0.000068619,(Spain/Almeria-COV002842/2020:0.000000005,Spain/Malaga-COV002841/2020:0.000000005):0.000067851):0.000169854,Spain/Madrid_LP16_6193/2020:0.00006795299999999999,Singapore/51/2020:0.000044697,(Thailand/Nonthaburi_193/2020:0,Thailand/Bangkok_237/2020:0,Thailand/Bangkok_238/2020:0,((Thailand/Bangkok-0034/2020:0.000000005,Thailand/Bangkok_2295/2020:0,Thailand/Bangkok-0065/2020:0.000047826,Thailand/Bangkok-CONI-0147/2020:0.000033997):0.000033983,Thailand/SI202769-NT/2020:0.000203899):0.000101951):0.00006797500000000001,Shenzhen/SZTH-002/2020:0.000033999);");
    // NewickImporter importer = new NewickImporter("(((0:0.5,(1:1.0,2:1.0)n6:1.0)n7:1.0,3:1.5)n8:1.0,(4:2.0,5:1.51)n9:1.5)n10;");
    // NewickImporter constraintsImporter = new NewickImporter("(((0:0.5,(1:1.0,2:1.0)n6:1.0)n7:1.0,3:1.5)n8:1.0,(4:2.0,5:1.51)n9:1.5)n10;");
    tree = new DefaultTreeModel(importer.importTree(null));
    IntervalList intervals = new TreeIntervals(tree, null, null);
    BigFastTreeIntervals bigFastTreeIntervals = new BigFastTreeIntervals(tree);
    SubtreeLeapOperator op = new SubtreeLeapOperator(tree, 1, 0.0001, SubtreeLeapOperator.DistanceKernelType.NORMAL, AdaptationMode.ADAPTATION_OFF, 0.2);
    NodeHeightOperator nh = new NodeHeightOperator(tree, 1, 1, NodeHeightOperator.OperatorType.UNIFORM, AdaptationMode.ADAPTATION_OFF, 0.25);
    NodeHeightOperator root = new NodeHeightOperator(tree, 1, 0.75, NodeHeightOperator.OperatorType.SCALEROOT, AdaptationMode.ADAPTATION_OFF, 0.25);
    boolean pass = true;
    MathUtils.setSeed(2);
    for (int i = 0; i < 100000; i++) {
        op.doOperation();
        intervals.calculateIntervals();
        // bigFastIntervals.makeDirty();
        bigFastTreeIntervals.calculateIntervals();
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getInterval(j) != bigFastTreeIntervals.getInterval(j)) {
                System.out.println(i);
                System.out.println("interval wrong");
                pass = false;
                break;
            }
        }
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getLineageCount(j) != bigFastTreeIntervals.getLineageCount(j)) {
                System.out.println(i);
                System.out.println("lineage Counts wrong: " + j);
                System.out.println("expected: " + intervals.getLineageCount(j));
                System.out.println("got " + bigFastTreeIntervals.getLineageCount(j));
                pass = false;
                break;
            }
        }
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getIntervalTime(j) != bigFastTreeIntervals.getIntervalTime(j)) {
                System.out.println(i);
                System.out.println("times wrong");
                pass = false;
                break;
            }
        }
        if (!pass) {
            break;
        }
    }
    assertTrue(pass);
}

Example 34 with DefaultTreeModel

use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.

the class AncestralStateBeagleTreeLikelihoodTest method testJointLikelihood.

public void testJointLikelihood() {
    TreeModel treeModel = new DefaultTreeModel("treeModel", tree);
    Sequence[] sequence = new Sequence[3];
    sequence[0] = new Sequence(new Taxon("0"), "A");
    sequence[1] = new Sequence(new Taxon("1"), "C");
    sequence[2] = new Sequence(new Taxon("2"), "C");
    Taxa taxa = new Taxa();
    for (Sequence s : sequence) {
        taxa.addTaxon(s.getTaxon());
    }
    SimpleAlignment alignment = new SimpleAlignment();
    for (Sequence s : sequence) {
        alignment.addSequence(s);
    }
    Parameter mu = new Parameter.Default(1, 1.0);
    Parameter kappa = new Parameter.Default(1, 1.0);
    double[] pi = { 0.25, 0.25, 0.25, 0.25 };
    Parameter freqs = new Parameter.Default(pi);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    GammaSiteRateModel siteRateModel = new GammaSiteRateModel("gammaModel", mu, null, -1, null);
    siteRateModel.setSubstitutionModel(hky);
    BranchModel branchModel = new HomogeneousBranchModel(siteRateModel.getSubstitutionModel());
    BranchRateModel branchRateModel = null;
    AncestralStateBeagleTreeLikelihood treeLikelihood = new AncestralStateBeagleTreeLikelihood(alignment, treeModel, branchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true, null, hky.getDataType(), "stateTag", // useMap = true
    true, false);
    double logLike = treeLikelihood.getLogLikelihood();
    StringBuffer buffer = new StringBuffer();
    // Tree.Utils.newick(treeModel, treeModel.getRoot(), false, Tree.BranchLengthType.LENGTHS_AS_TIME,
    // null, null, new NodeAttributeProvider[]{treeLikelihood}, null, null, buffer);
    TreeUtils.newick(treeModel, treeModel.getRoot(), false, TreeUtils.BranchLengthType.LENGTHS_AS_TIME, null, null, new TreeTraitProvider[] { treeLikelihood }, null, buffer);
    System.out.println(buffer);
    System.out.println("t_CA(2) = " + t(false, 2.0));
    System.out.println("t_CC(1) = " + t(true, 1.0));
    double trueValue = 0.25 * t(false, 2.0) * Math.pow(t(true, 1.0), 3.0);
    assertEquals(logLike, Math.log(trueValue), 1e-6);
}

Example 35 with DefaultTreeModel

use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.

the class PMDTestProblem method testPMD.

public void testPMD() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 496432.69917113904, 0, Double.POSITIVE_INFINITY);
    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, 4.0E-7, 0, 100.0);
    StrictClockBranchRates branchRateModel = new StrictClockBranchRates(rateParameter);
    // Sub model
    Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    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);
    // SequenceErrorModel
    Parameter ageRelatedRateParameter = new Parameter.Default(SequenceErrorModelParser.AGE_RELATED_RATE, 4.0E-7, 0, 100.0);
    TipStatesModel aDNADamageModel = new SequenceErrorModel(null, null, SequenceErrorModel.ErrorType.TRANSITIONS_ONLY, null, ageRelatedRateParameter, null);
    // treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, aDNADamageModel, 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);
    operator = new ScaleOperator(ageRelatedRateParameter, 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, 49643.2699171139, 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);
    // ??? correct?
    operator = new DeltaExchangeOperator(freqs, new int[] { 1, 1, 1, 1 }, 0.01, 1.0, false, AdaptationMode.ADAPTATION_ON);
    schedule.addOperator(operator);
    // CompoundLikelihood
    OneOnXPrior likelihood1 = new OneOnXPrior();
    likelihood1.addData(popSize);
    OneOnXPrior likelihood2 = new OneOnXPrior();
    likelihood2.addData(kappa);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(likelihood1);
    likelihoods.add(likelihood2);
    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, 1000, false);
    loggers[0].add(posterior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(rateParameter);
    loggers[0].add(ageRelatedRateParameter);
    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);
    // 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("PMDTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    // <expectation name="clock.rate" value="1.5E-7"/>
    // <expectation name="errorModel.ageRate" value="0.7E-7"/>
    // <expectation name="hky.kappa" value="10"/>
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 10);
    TraceCorrelation rateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(StrictClockBranchRates.RATE));
    assertExpectation(StrictClockBranchRates.RATE, rateStats, 1.5E-7);
    TraceCorrelation ageRateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(SequenceErrorModelParser.AGE_RELATED_RATE));
    assertExpectation(SequenceErrorModelParser.AGE_RELATED_RATE, ageRateStats, 0.7E-7);
}