use of dr.evomodel.operators.SubtreeSlideOperator in project beast-mcmc by beast-dev.
the class SubtreeSlideOperatorParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
boolean swapRates = xo.getAttribute(SWAP_RATES, false);
boolean swapTraits = xo.getAttribute(SWAP_TRAITS, false);
boolean scaledDirichletBranches = xo.getAttribute(DIRICHLET_BRANCHES, false);
AdaptationMode mode = AdaptationMode.DEFAULT;
if (xo.hasAttribute(AdaptableMCMCOperator.AUTO_OPTIMIZE)) {
if (xo.getBooleanAttribute(AdaptableMCMCOperator.AUTO_OPTIMIZE)) {
mode = AdaptationMode.ADAPTATION_ON;
} else {
mode = AdaptationMode.ADAPTATION_OFF;
}
}
DefaultTreeModel treeModel = (DefaultTreeModel) xo.getChild(DefaultTreeModel.class);
final double weight = xo.getDoubleAttribute(MCMCOperator.WEIGHT);
final double targetAcceptance = xo.getAttribute(TARGET_ACCEPTANCE, 0.234);
final double size = xo.getAttribute("size", 1.0);
if (Double.isInfinite(size) || size <= 0.0) {
throw new XMLParseException("size attribute must be positive and not infinite. was " + size + " for tree " + treeModel.getId());
}
final boolean gaussian = xo.getBooleanAttribute("gaussian");
SubtreeSlideOperator operator = new SubtreeSlideOperator(treeModel, weight, size, gaussian, swapRates, swapTraits, scaledDirichletBranches, mode, targetAcceptance);
return operator;
}
use of dr.evomodel.operators.SubtreeSlideOperator in project beast-mcmc by beast-dev.
the class YuleModelTest method testYuleWithSubtreeSlide.
public void testYuleWithSubtreeSlide() {
DefaultTreeModel treeModel = new DefaultTreeModel("treeModel", tree);
OperatorSchedule schedule = new SimpleOperatorSchedule();
MCMCOperator operator = new SubtreeSlideOperator(treeModel, 1, 1, true, false, false, false, AdaptationMode.ADAPTATION_ON, AdaptableMCMCOperator.DEFAULT_ADAPTATION_TARGET);
schedule.addOperator(operator);
yuleTester(treeModel, schedule);
}
use of dr.evomodel.operators.SubtreeSlideOperator 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);
}
use of dr.evomodel.operators.SubtreeSlideOperator 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 = ((DefaultTreeModel) treeModel).getRootHeightParameter();
rootHeight.setId(TREE_HEIGHT);
operator = new ScaleOperator(rootHeight, 0.5);
operator.setWeight(1.0);
schedule.addOperator(operator);
Parameter internalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(false, true, false);
operator = new UniformOperator(internalHeights, 10.0);
schedule.addOperator(operator);
operator = new SubtreeSlideOperator(((DefaultTreeModel) treeModel), 1, 1, true, false, false, false, AdaptationMode.ADAPTATION_ON, AdaptableMCMCOperator.DEFAULT_ADAPTATION_TARGET);
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);
}
use of dr.evomodel.operators.SubtreeSlideOperator 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);
}
Aggregations