Examples with MCMCOptions dr.inference.mcmc.MCMCOptions used on opensource projects

Example 6 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions in project beast-mcmc by beast-dev.

the class GeneralSubstitutionModelTest method testGeneralSubstitutionModel.

public void testGeneralSubstitutionModel() {
    // Sub model
    FrequencyModel freqModel = new FrequencyModel(dataType, alignment.getStateFrequencies());
    // dimension="5" value="1.0"
    Parameter ratesPara = new Parameter.Default(GeneralSubstitutionModelParser.RATES, 5, 1.0);
    // relativeTo="5"
    GeneralSubstitutionModel generalSubstitutionModel = new GeneralSubstitutionModel(dataType, freqModel, ratesPara, 4);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(generalSubstitutionModel);
    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, null, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(ratesPara, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 10.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 1, 1, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(ratesPara);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(ratesPara);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(10000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, treeLikelihood, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("GeneralSubstitutionModelTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //      <expectation name="likelihood" value="-1815.75"/>
    //		<expectation name="treeModel.rootHeight" value="6.42048E-2"/>
    //		<expectation name="rateAC" value="6.08986E-2"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.75);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 0.0640787258170083);
    TraceCorrelation rateACStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(GeneralSubstitutionModelParser.RATES + "1"));
    assertExpectation(GeneralSubstitutionModelParser.RATES + "1", rateACStats, 0.061071756742081366);
}

Example 7 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions in project beast-mcmc by beast-dev.

the class ARGAddRemoveOperatorTest method flatPriorTester.

private void flatPriorTester(ARGModel arg, int chainLength, int sampleTreeEvery, double nodeCountSetting, double rootHeightAlpha, double rootHeightBeta, int maxCount) throws IOException, Importer.ImportException {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(chainLength);
    //        double nodeCountSetting = 2.0;
    //        double rootHeightAlpha = 100;
    //        double rootHeightBeta = 0.5;
    OperatorSchedule schedule = getSchedule(arg);
    ARGUniformPrior uniformPrior = new ARGUniformPrior(arg, maxCount, arg.getExternalNodeCount());
    PoissonDistribution poisson = new PoissonDistribution(nodeCountSetting);
    DistributionLikelihood nodeCountPrior = new DistributionLikelihood(poisson, 0.0);
    ARGReassortmentNodeCountStatistic nodeCountStatistic = new ARGReassortmentNodeCountStatistic("nodeCount", arg);
    nodeCountPrior.addData(nodeCountStatistic);
    DistributionLikelihood rootPrior = new DistributionLikelihood(new GammaDistribution(rootHeightAlpha, rootHeightBeta), 0.0);
    CompoundParameter rootHeight = (CompoundParameter) arg.createNodeHeightsParameter(true, false, false);
    rootPrior.addData(rootHeight);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(uniformPrior);
    likelihoods.add(rootPrior);
    likelihoods.add(nodeCountPrior);
    CompoundLikelihood compoundLikelihood = new CompoundLikelihood(1, likelihoods);
    compoundLikelihood.setId("likelihood1");
    MCLogger[] loggers = new MCLogger[3];
    loggers[0] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[0].add(compoundLikelihood);
    loggers[0].add(arg);
    File file = new File("test.args");
    file.deleteOnExit();
    FileOutputStream out = new FileOutputStream(file);
    loggers[1] = new ARGLogger(arg, new TabDelimitedFormatter(out), sampleTreeEvery, "test");
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    loggers[2] = new MCLogger(formatter, sampleTreeEvery, false);
    loggers[2].add(arg);
    arg.getRootHeightParameter().setId("root");
    loggers[2].add(arg.getRootHeightParameter());
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, compoundLikelihood, schedule, loggers);
    mcmc.run();
    out.flush();
    out.close();
    List<Trace> traces = formatter.getTraces();
    //        Set<String> uniqueTrees = new HashSet<String>();
    //
    //        NexusImporter importer = new NexusImporter(new FileReader(file));
    //        while (importer.hasTree()) {
    //            Tree t = importer.importNextTree();
    //            uniqueTrees.add(Tree.Utils.uniqueNewick(t, t.getRoot()));
    //        }
    //
    //        TestCase.assertEquals(numTopologies, uniqueTrees.size());            List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("ARGTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    TraceCorrelation nodeCountStats = traceList.getCorrelationStatistics(1);
    TraceCorrelation rootHeightStats = traceList.getCorrelationStatistics(4);
    assertExpectation("nodeCount", nodeCountStats, poisson.truncatedMean(maxCount));
    assertExpectation(TreeModelParser.ROOT_HEIGHT, rootHeightStats, rootHeightAlpha * rootHeightBeta);
}

Example 8 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions 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);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(treeModel, null, new ArrayList<TaxonList>(), 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);
    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 = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 15.0, 0.0077, true, false, false, false, CoercionMode.COERCION_ON);
    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 9 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions in project beast-mcmc by beast-dev.

the class TestCalibratedYuleModel method yuleTester.

private void yuleTester(TreeModel treeModel, OperatorSchedule schedule, Parameter brParameter, double S, int chainLength) throws IOException, TreeUtils.MissingTaxonException {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(chainLength);
    TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
    TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model("yule", brParameter, null, null, BirthDeathGernhard08Model.TreeType.UNSCALED, Units.Type.SUBSTITUTIONS, false);
    Likelihood speciationLikelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
    Taxa halfTaxa = new Taxa();
    for (int i = 0; i < taxa.getTaxonCount() / 2; i++) {
        halfTaxa.addTaxon(new Taxon("T" + Integer.toString(i)));
    }
    TMRCAStatistic tmrca = new TMRCAStatistic("tmrca(halfTaxa)", treeModel, halfTaxa, false, false);
    DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(M, S), // meanInRealSpace="false"
    0);
    logNormalLikelihood.addData(tmrca);
    MonophylyStatistic monophylyStatistic = new MonophylyStatistic("monophyly(halfTaxa)", treeModel, halfTaxa, null);
    BooleanLikelihood booleanLikelihood = new BooleanLikelihood();
    booleanLikelihood.addData(monophylyStatistic);
    //CompoundLikelihood
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(speciationLikelihood);
    likelihoods.add(logNormalLikelihood);
    likelihoods.add(booleanLikelihood);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    ArrayLogFormatter logformatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[1];
    loggers[0] = new MCLogger(logformatter, (int) (options.getChainLength() / 10000), false);
    loggers[0].add(speciationLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(tmrca);
    loggers[0].add(tls);
    loggers[0].add(brParameter);
    mcmc.setShowOperatorAnalysis(false);
    mcmc.init(options, prior, schedule, loggers);
    mcmc.run();
    List<Trace> traces = logformatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 1000);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    NumberFormatter formatter = new NumberFormatter(8);
    TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TL));
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("tmrca(halfTaxa)"));
    //        out.write("tmrcaHeight = \t");
    out.write(formatter.format(treeHeightStats.getMean()));
    out.write("\t");
    double expectedNodeHeight = Math.pow(Math.E, (M + (Math.pow(S, 2) / 2)));
    //        out.write("expectation = \t");
    out.write(formatter.format(expectedNodeHeight));
    out.write("\t");
    double error = Math.abs((treeHeightStats.getMean() - expectedNodeHeight) / expectedNodeHeight);
    NumberFormat percentFormatter = NumberFormat.getNumberInstance();
    percentFormatter.setMinimumFractionDigits(5);
    percentFormatter.setMinimumFractionDigits(5);
    //        out.write("error = \t");
    out.write(percentFormatter.format(error));
    out.write("\t");
    //        out.write("tl.ess = \t");
    out.write(Double.toString(tlStats.getESS()));
    System.out.println("tmrcaHeight = " + formatter.format(treeHeightStats.getMean()) + ";  expectation = " + formatter.format(expectedNodeHeight) + ";  error = " + percentFormatter.format(error) + ";  tl.ess = " + tlStats.getESS());
}

Example 10 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions in project beast-mcmc by beast-dev.

the class MCMCParser method parseXMLObject.

/**
     * @return an mcmc object based on the XML element it was passed.
     */
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    MCMC mcmc = new MCMC(xo.getAttribute(NAME, "mcmc1"));
    long chainLength = xo.getLongIntegerAttribute(CHAIN_LENGTH);
    boolean useCoercion = xo.getAttribute(COERCION, true);
    long coercionDelay = chainLength / 100;
    if (xo.hasAttribute(PRE_BURNIN)) {
        coercionDelay = xo.getIntegerAttribute(PRE_BURNIN);
    }
    coercionDelay = xo.getAttribute(COERCION_DELAY, coercionDelay);
    double temperature = xo.getAttribute(TEMPERATURE, 1.0);
    long fullEvaluationCount = xo.getAttribute(FULL_EVALUATION, 2000);
    double evaluationTestThreshold = MarkovChain.EVALUATION_TEST_THRESHOLD;
    if (System.getProperty("mcmc.evaluation.threshold") != null) {
        evaluationTestThreshold = Double.parseDouble(System.getProperty("mcmc.evaluation.threshold"));
    }
    evaluationTestThreshold = xo.getAttribute(EVALUATION_THRESHOLD, evaluationTestThreshold);
    int minOperatorCountForFullEvaluation = xo.getAttribute(MIN_OPS_EVALUATIONS, 1);
    MCMCOptions options = new MCMCOptions(chainLength, fullEvaluationCount, minOperatorCountForFullEvaluation, evaluationTestThreshold, useCoercion, coercionDelay, temperature);
    OperatorSchedule opsched = (OperatorSchedule) xo.getChild(OperatorSchedule.class);
    Likelihood likelihood = (Likelihood) xo.getChild(Likelihood.class);
    likelihood.setUsed();
    if (Boolean.valueOf(System.getProperty("show_warnings", "false"))) {
        // check that all models, parameters and likelihoods are being used
        for (Likelihood l : Likelihood.FULL_LIKELIHOOD_SET) {
            if (!l.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Likelihood, " + l.getId() + ", of class " + l.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
        for (Model m : Model.FULL_MODEL_SET) {
            if (!m.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Model, " + m.getId() + ", of class " + m.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
        for (Parameter p : Parameter.FULL_PARAMETER_SET) {
            if (!p.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Parameter, " + p.getId() + ", of class " + p.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
    }
    ArrayList<Logger> loggers = new ArrayList<Logger>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        Object child = xo.getChild(i);
        if (child instanceof Logger) {
            loggers.add((Logger) child);
        }
    }
    mcmc.setShowOperatorAnalysis(true);
    if (xo.hasAttribute(OPERATOR_ANALYSIS)) {
        mcmc.setOperatorAnalysisFile(XMLParser.getLogFile(xo, OPERATOR_ANALYSIS));
    }
    Logger[] loggerArray = new Logger[loggers.size()];
    loggers.toArray(loggerArray);
    java.util.logging.Logger.getLogger("dr.inference").info("\nCreating the MCMC chain:" + "\n  chainLength=" + options.getChainLength() + "\n  autoOptimize=" + options.useCoercion() + (options.useCoercion() ? "\n  autoOptimize delayed for " + options.getCoercionDelay() + " steps" : "") + (options.getFullEvaluationCount() == 0 ? "\n  full evaluation test off" : ""));
    mcmc.init(options, likelihood, opsched, loggerArray);
    MarkovChain mc = mcmc.getMarkovChain();
    double initialScore = mc.getCurrentScore();
    if (initialScore == Double.NEGATIVE_INFINITY) {
        String message = "The initial posterior is zero";
        if (likelihood instanceof CompoundLikelihood) {
            message += ": " + ((CompoundLikelihood) likelihood).getDiagnosis(2);
        } else {
            message += "!";
        }
        throw new IllegalArgumentException(message);
    }
    if (!xo.getAttribute(SPAWN, true))
        mcmc.setSpawnable(false);
    return mcmc;
}