use of dr.evomodel.coalescent.DemographicModel in project beast-mcmc by beast-dev.
the class TwoEpochDemographicModelParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
Units.Type units = XMLUnits.Utils.getUnitsAttr(xo);
XMLObject cxo = xo.getChild(EPOCH_1);
DemographicModel demo1 = (DemographicModel) cxo.getChild(DemographicModel.class);
cxo = xo.getChild(EPOCH_2);
DemographicModel demo2 = (DemographicModel) cxo.getChild(DemographicModel.class);
cxo = xo.getChild(TRANSITION_TIME);
Parameter timeParameter = (Parameter) cxo.getChild(Parameter.class);
return new TwoEpochDemographicModel(demo1, demo2, timeParameter, units);
}
use of dr.evomodel.coalescent.DemographicModel in project beast-mcmc by beast-dev.
the class CoalescentLikelihoodParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
XMLObject cxo = xo.getChild(MODEL);
DemographicModel demoModel = (DemographicModel) cxo.getChild(DemographicModel.class);
List<TreeModel> trees = new ArrayList<TreeModel>();
List<Double> popFactors = new ArrayList<Double>();
MultiLociTreeSet treesSet = demoModel instanceof MultiLociTreeSet ? (MultiLociTreeSet) demoModel : null;
for (int k = 0; k < xo.getChildCount(); ++k) {
final Object child = xo.getChild(k);
if (child instanceof XMLObject) {
cxo = (XMLObject) child;
if (cxo.getName().equals(POPULATION_TREE)) {
final TreeModel t = (TreeModel) cxo.getChild(TreeModel.class);
assert t != null;
trees.add(t);
popFactors.add(cxo.getAttribute(POPULATION_FACTOR, 1.0));
}
}
// in the future we may have arbitrary multi-loci element
// else if( child instanceof MultiLociTreeSet ) {
// treesSet = (MultiLociTreeSet)child;
// }
}
TreeModel treeModel = null;
if (trees.size() == 1 && popFactors.get(0) == 1.0) {
treeModel = trees.get(0);
} else if (trees.size() > 1) {
treesSet = new MultiLociTreeSet.Default(trees, popFactors);
} else if (!(trees.size() == 0 && treesSet != null)) {
throw new XMLParseException("Incorrectly constructed likelihood element");
}
TaxonList includeSubtree = null;
if (xo.hasChildNamed(INCLUDE)) {
includeSubtree = (TaxonList) xo.getElementFirstChild(INCLUDE);
}
List<TaxonList> excludeSubtrees = new ArrayList<TaxonList>();
if (xo.hasChildNamed(EXCLUDE)) {
cxo = xo.getChild(EXCLUDE);
for (int i = 0; i < cxo.getChildCount(); i++) {
excludeSubtrees.add((TaxonList) cxo.getChild(i));
}
}
if (treeModel != null) {
try {
return new CoalescentLikelihood(treeModel, includeSubtree, excludeSubtrees, demoModel);
} catch (TreeUtils.MissingTaxonException mte) {
throw new XMLParseException("treeModel missing a taxon from taxon list in " + getParserName() + " element");
}
} else {
if (includeSubtree != null || excludeSubtrees.size() > 0) {
throw new XMLParseException("Include/Exclude taxa not supported for multi locus sets");
}
// a base - and modifing it will probsbly result in a bigger mess.
return new OldAbstractCoalescentLikelihood(treesSet, demoModel);
}
}
use of dr.evomodel.coalescent.DemographicModel 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.evomodel.coalescent.DemographicModel 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