use of dr.evolution.util.TaxonList 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);
}
use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.
the class IndependentCoalescentSampler method doOperation.
/**
* change the parameter and return the hastings ratio.
*/
public double doOperation() {
CoalescentSimulator simulator = new CoalescentSimulator();
List<TaxonList> taxonLists = new ArrayList<TaxonList>();
double rootHeight = -1.0;
double oldLikelihood = 0.0;
double newLikelihood = 0.0;
// should have one child that is node
for (int i = 0; i < xo.getChildCount(); i++) {
final Object child = xo.getChild(i);
//careful: Trees are TaxonLists ... (AER); see OldCoalescentSimulatorParser
if (child instanceof Tree) {
//do nothing
} else if (child instanceof TaxonList) {
//taxonLists.add((TaxonList) child);
taxonLists.add((Taxa) child);
//taxa added
break;
}
}
try {
Tree[] trees = new Tree[taxonLists.size()];
// simulate each taxonList separately
for (int i = 0; i < taxonLists.size(); i++) {
trees[i] = simulator.simulateTree(taxonLists.get(i), demoModel);
}
oldLikelihood = coalescent.getLogLikelihood();
SimpleTree simTree = simulator.simulateTree(trees, demoModel, rootHeight, trees.length != 1);
//this would be the normal way to do it
treeModel.beginTreeEdit();
//now it's allowed to adjust the tree structure
treeModel.adoptTreeStructure(simTree);
//endTreeEdit() would then fire the events
treeModel.endTreeEdit();
newLikelihood = coalescent.getLogLikelihood();
} catch (IllegalArgumentException iae) {
try {
throw new XMLParseException(iae.getMessage());
} catch (XMLParseException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
return oldLikelihood - newLikelihood;
}
use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.
the class HiddenLinkageTreeLogger method processTree.
/**
* Create a tree of the same topology as the input tree, but with
* the linkage groups replaced by their constituent reads
*/
protected static SimpleTree processTree(Tree tree, HiddenLinkageModel hlm) {
TaxonList reads = hlm.getData().getReadsTaxa();
TaxonList reference = hlm.getData().getReferenceTaxa();
// allocate space
int nodeCount = tree.getTaxonCount() + reads.getTaxonCount();
nodeCount = 2 * nodeCount - 1;
SimpleNode[] nodes = new SimpleNode[nodeCount];
for (int i = 0; i < nodes.length; i++) {
nodes[i] = new SimpleNode();
nodes[i].setNumber(i);
}
SimpleNode root = null;
// copy the tree structure
for (int i = 0; i < tree.getNodeCount(); i++) {
NodeRef n = tree.getNode(i);
for (int cI = 0; cI < tree.getChildCount(n); cI++) {
NodeRef c1 = tree.getChild(n, cI);
nodes[n.getNumber()].addChild(nodes[c1.getNumber()]);
}
nodes[n.getNumber()].setHeight(tree.getNodeHeight(n));
nodes[n.getNumber()].setRate(tree.getNodeRate(n));
nodes[n.getNumber()].setTaxon(tree.getNodeTaxon(n));
}
root = nodes[tree.getRoot().getNumber()];
// now replace linkage groups with their constituent reads
// first free up anything in the range of read leaf nodes
int nextFree = tree.getNodeCount();
int readI = 0;
for (int i = reference.getTaxonCount(); i < reference.getTaxonCount() + reads.getTaxonCount(); i++) {
SimpleNode tmp = nodes[nextFree];
nodes[nextFree] = nodes[i];
nodes[nextFree].setNumber(nextFree);
nodes[i] = tmp;
nodes[i].setNumber(i);
nodes[i].setTaxon(reads.getTaxon(readI));
readI++;
nextFree++;
}
// if the linkage group has many reads, build a ladder
for (int i = 0; i < nodes.length; i++) {
SimpleNode n = nodes[i];
if (n.getTaxon() == null)
continue;
if (reads.getTaxonIndex(n.getTaxon()) >= 0 || reference.getTaxonIndex(n.getTaxon()) >= 0)
// not a linkage group
continue;
int gid = hlm.getTaxonIndex(n.getTaxon()) - reference.getTaxonCount();
if (gid < 0) {
System.err.println("big trouble, little china");
}
Set<Taxon> group = hlm.getGroup(gid);
if (group.size() == 0) {
// remove the group completely
SimpleNode parent = n.getParent();
parent.removeChild(n);
if (parent.getChildCount() == 1) {
SimpleNode grandparent = parent.getParent();
SimpleNode child = parent.getChild(0);
parent.removeChild(child);
if (grandparent == null) {
// parent is root! other child should become root.
root = child;
} else {
grandparent.removeChild(parent);
grandparent.addChild(child);
}
}
} else if (group.size() == 1) {
// swap the group with the constituent read
Taxon[] tax = new Taxon[group.size()];
tax = (Taxon[]) group.toArray(tax);
int rI = getTaxonNode(tax[0], nodes);
SimpleNode parent = n.getParent();
parent.removeChild(n);
parent.addChild(nodes[rI]);
} else {
// create a star tree with the reads
Taxon[] tax = new Taxon[group.size()];
tax = (Taxon[]) group.toArray(tax);
SimpleNode parent = n.getParent();
parent.removeChild(n);
parent.addChild(nodes[nextFree]);
int tI = 0;
for (; tI < tax.length - 2; tI++) {
int rI = getTaxonNode(tax[tI], nodes);
nodes[nextFree].addChild(nodes[rI]);
nodes[nextFree].addChild(nodes[nextFree + 1]);
nextFree++;
}
int rI = getTaxonNode(tax[tI], nodes);
nodes[nextFree].addChild(nodes[rI]);
int rJ = getTaxonNode(tax[tI + 1], nodes);
nodes[nextFree].addChild(nodes[rJ]);
nextFree++;
}
}
SimpleTree st = new SimpleTree(root);
return st;
}
use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.
the class OldCoalescentSimulatorParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
CoalescentSimulator simulator = new CoalescentSimulator();
DemographicModel demoModel = (DemographicModel) xo.getChild(DemographicModel.class);
List<TaxonList> taxonLists = new ArrayList<TaxonList>();
List<Tree> subtrees = new ArrayList<Tree>();
double rootHeight = xo.getAttribute(ROOT_HEIGHT, -1.0);
if (xo.hasAttribute(RESCALE_HEIGHT)) {
rootHeight = xo.getDoubleAttribute(RESCALE_HEIGHT);
}
// should have one child that is node
for (int i = 0; i < xo.getChildCount(); i++) {
final Object child = xo.getChild(i);
// AER - swapped the order of these round because Trees are TaxonLists...
if (child instanceof Tree) {
subtrees.add((Tree) child);
} else if (child instanceof TaxonList) {
taxonLists.add((TaxonList) child);
} else if (xo.getChildName(i).equals(CONSTRAINED_TAXA)) {
XMLObject constrainedTaxa = (XMLObject) child;
// all taxa
final TaxonList taxa = (TaxonList) constrainedTaxa.getChild(TaxonList.class);
List<CoalescentSimulator.TaxaConstraint> constraints = new ArrayList<CoalescentSimulator.TaxaConstraint>();
final String setsNotCompatibleMessage = "taxa sets not compatible";
for (int nc = 0; nc < constrainedTaxa.getChildCount(); ++nc) {
final Object object = constrainedTaxa.getChild(nc);
if (object instanceof XMLObject) {
final XMLObject constraint = (XMLObject) object;
if (constraint.getName().equals(TMRCA_CONSTRAINT)) {
TaxonList taxaSubSet = (TaxonList) constraint.getChild(TaxonList.class);
ParametricDistributionModel dist = (ParametricDistributionModel) constraint.getChild(ParametricDistributionModel.class);
boolean isMono = constraint.getAttribute(IS_MONOPHYLETIC, true);
final CoalescentSimulator.TaxaConstraint taxaConstraint = new CoalescentSimulator.TaxaConstraint(taxaSubSet, dist, isMono);
int insertPoint;
for (insertPoint = 0; insertPoint < constraints.size(); ++insertPoint) {
// if new <= constraints[insertPoint] insert before insertPoint
final CoalescentSimulator.TaxaConstraint iConstraint = constraints.get(insertPoint);
if (iConstraint.isMonophyletic) {
if (!taxaConstraint.isMonophyletic) {
continue;
}
final TaxonList taxonsip = iConstraint.taxons;
final int nIn = simulator.sizeOfIntersection(taxonsip, taxaSubSet);
if (nIn == taxaSubSet.getTaxonCount()) {
break;
}
if (nIn > 0 && nIn != taxonsip.getTaxonCount()) {
throw new XMLParseException(setsNotCompatibleMessage);
}
} else {
// reached non mono area
if (!taxaConstraint.isMonophyletic) {
if (iConstraint.upper >= taxaConstraint.upper) {
break;
}
} else {
break;
}
}
}
constraints.add(insertPoint, taxaConstraint);
}
}
}
final int nConstraints = constraints.size();
if (nConstraints == 0) {
if (taxa != null) {
taxonLists.add(taxa);
}
} else {
for (int nc = 0; nc < nConstraints; ++nc) {
CoalescentSimulator.TaxaConstraint cnc = constraints.get(nc);
if (!cnc.isMonophyletic) {
for (int nc1 = nc - 1; nc1 >= 0; --nc1) {
CoalescentSimulator.TaxaConstraint cnc1 = constraints.get(nc1);
int x = simulator.sizeOfIntersection(cnc.taxons, cnc1.taxons);
if (x > 0) {
Taxa combinedTaxa = new Taxa(cnc.taxons);
combinedTaxa.addTaxa(cnc1.taxons);
cnc = new CoalescentSimulator.TaxaConstraint(combinedTaxa, cnc.lower, cnc.upper, cnc.isMonophyletic);
constraints.set(nc, cnc);
}
}
}
}
// determine upper bound for each set.
double[] upper = new double[nConstraints];
for (int nc = nConstraints - 1; nc >= 0; --nc) {
final CoalescentSimulator.TaxaConstraint cnc = constraints.get(nc);
if (cnc.realLimits()) {
upper[nc] = cnc.upper;
} else {
upper[nc] = Double.POSITIVE_INFINITY;
}
}
for (int nc = nConstraints - 1; nc >= 0; --nc) {
final CoalescentSimulator.TaxaConstraint cnc = constraints.get(nc);
if (upper[nc] < Double.POSITIVE_INFINITY) {
for (int nc1 = nc - 1; nc1 >= 0; --nc1) {
final CoalescentSimulator.TaxaConstraint cnc1 = constraints.get(nc1);
if (simulator.contained(cnc1.taxons, cnc.taxons)) {
upper[nc1] = Math.min(upper[nc1], upper[nc]);
if (cnc1.realLimits() && cnc1.lower > upper[nc1]) {
throw new XMLParseException(setsNotCompatibleMessage);
}
break;
}
}
}
}
// collect subtrees here
List<Tree> st = new ArrayList<Tree>();
for (int nc = 0; nc < constraints.size(); ++nc) {
final CoalescentSimulator.TaxaConstraint nxt = constraints.get(nc);
// collect all previously built subtrees which are a subset of taxa set to be added
List<Tree> subs = new ArrayList<Tree>();
Taxa newTaxons = new Taxa(nxt.taxons);
for (int k = 0; k < st.size(); ++k) {
final Tree stk = st.get(k);
int x = simulator.sizeOfIntersection(stk, nxt.taxons);
if (x == st.get(k).getTaxonCount()) {
final Tree tree = st.remove(k);
--k;
subs.add(tree);
newTaxons.removeTaxa(tree);
}
}
SimpleTree tree = simulator.simulateTree(newTaxons, demoModel);
final double lower = nxt.realLimits() ? nxt.lower : 0;
if (upper[nc] < Double.MAX_VALUE) {
simulator.attemptToScaleTree(tree, (lower + upper[nc]) / 2);
}
if (subs.size() > 0) {
if (tree.getTaxonCount() > 0)
subs.add(tree);
double h = -1;
if (upper[nc] < Double.MAX_VALUE) {
for (Tree t : subs) {
// protect against 1 taxa tree with height 0
final double rh = t.getNodeHeight(t.getRoot());
h = Math.max(h, rh > 0 ? rh : (lower + upper[nc]) / 2);
}
h = (h + upper[nc]) / 2;
}
tree = simulator.simulateTree(subs.toArray(new Tree[subs.size()]), demoModel, h, true);
}
st.add(tree);
}
// add a taxon list for remaining taxa
if (taxa != null) {
final Taxa list = new Taxa();
for (int j = 0; j < taxa.getTaxonCount(); ++j) {
Taxon taxonj = taxa.getTaxon(j);
for (Tree aSt : st) {
if (aSt.getTaxonIndex(taxonj) >= 0) {
taxonj = null;
break;
}
}
if (taxonj != null) {
list.addTaxon(taxonj);
}
}
if (list.getTaxonCount() > 0) {
taxonLists.add(list);
}
}
if (st.size() > 1) {
final Tree t = simulator.simulateTree(st.toArray(new Tree[st.size()]), demoModel, -1, false);
subtrees.add(t);
} else {
subtrees.add(st.get(0));
}
}
}
}
if (taxonLists.size() == 0) {
if (subtrees.size() == 1) {
return subtrees.get(0);
}
throw new XMLParseException("Expected at least one taxonList or two subtrees in " + getParserName() + " element.");
}
try {
Tree[] trees = new Tree[taxonLists.size() + subtrees.size()];
// simulate each taxonList separately
for (int i = 0; i < taxonLists.size(); i++) {
trees[i] = simulator.simulateTree(taxonLists.get(i), demoModel);
}
// add the preset trees
for (int i = 0; i < subtrees.size(); i++) {
trees[i + taxonLists.size()] = subtrees.get(i);
}
return simulator.simulateTree(trees, demoModel, rootHeight, trees.length != 1);
} catch (IllegalArgumentException iae) {
throw new XMLParseException(iae.getMessage());
}
}
use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.
the class CoalescentSimulatorParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
CoalescentSimulator simulator = new CoalescentSimulator();
DemographicModel demoModel = (DemographicModel) xo.getChild(DemographicModel.class);
List<TaxonList> taxonLists = new ArrayList<TaxonList>();
List<Tree> subtrees = new ArrayList<Tree>();
double height = xo.getAttribute(HEIGHT, Double.NaN);
// should have one child that is node
for (int i = 0; i < xo.getChildCount(); i++) {
final Object child = xo.getChild(i);
// AER - swapped the order of these round because Trees are TaxonLists...
if (child instanceof Tree) {
subtrees.add((Tree) child);
} else if (child instanceof TaxonList) {
taxonLists.add((TaxonList) child);
}
}
if (taxonLists.size() == 0) {
if (subtrees.size() == 1) {
return subtrees.get(0);
}
throw new XMLParseException("Expected at least one taxonList or two subtrees in " + getParserName() + " element.");
}
Taxa remainingTaxa = new Taxa();
for (int i = 0; i < taxonLists.size(); i++) {
remainingTaxa.addTaxa(taxonLists.get(i));
}
for (int i = 0; i < subtrees.size(); i++) {
remainingTaxa.removeTaxa(subtrees.get(i));
}
try {
Tree[] trees = new Tree[subtrees.size() + remainingTaxa.getTaxonCount()];
// add the preset trees
for (int i = 0; i < subtrees.size(); i++) {
trees[i] = subtrees.get(i);
}
// add all the remaining taxa in as single tip trees...
for (int i = 0; i < remainingTaxa.getTaxonCount(); i++) {
Taxa tip = new Taxa();
tip.addTaxon(remainingTaxa.getTaxon(i));
trees[i + subtrees.size()] = simulator.simulateTree(tip, demoModel);
}
return simulator.simulateTree(trees, demoModel, height, trees.length != 1);
} catch (IllegalArgumentException iae) {
throw new XMLParseException(iae.getMessage());
}
}
Aggregations