Examples with SimpleTree dr.evolution.tree.SimpleTree used on opensource projects

Example 6 with SimpleTree

use of dr.evolution.tree.SimpleTree 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;
}

Example 7 with SimpleTree

use of dr.evolution.tree.SimpleTree in project beast-mcmc by beast-dev.

the class WeightedMultiplicativeBinary method report.

/**
     * Creates the report. The estimated posterior of the given tree is printed.
     *
     * @throws IOException if general I/O error occurs
     */
public void report(Tree tree) throws IOException {
    System.err.println("making report");
    SimpleTree sTree = new SimpleTree(tree);
    System.out.println("Estimated marginal posterior by condiational clade frequencies:");
    System.out.println(getTreeProbability(sTree));
    System.out.flush();
}

Example 8 with SimpleTree

use of dr.evolution.tree.SimpleTree 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());
    }
}

Example 9 with SimpleTree

use of dr.evolution.tree.SimpleTree in project beast-mcmc by beast-dev.

the class StructuredCoalescentSimulator method simulateTree.

/**
     * Simulates a coalescent tree, given a taxon list.
     *
     * @param taxa          the set of taxa to simulate a coalescent tree between
     * @param demoFunctions the demographic function to use
     */
public Tree simulateTree(TaxonList[] taxa, DemographicFunction[] demoFunctions, ColourChangeMatrix colourChangeMatrix) {
    SimpleNode[][] nodes = new SimpleNode[taxa.length][];
    for (int i = 0; i < taxa.length; i++) {
        nodes[i] = new SimpleNode[taxa[i].getTaxonCount()];
        for (int j = 0; j < taxa[i].getTaxonCount(); j++) {
            nodes[i][j] = new SimpleNode();
            nodes[i][j].setTaxon(taxa[i].getTaxon(j));
        }
    }
    dr.evolution.util.Date mostRecent = null;
    boolean usingDates = false;
    for (int i = 0; i < taxa.length; i++) {
        for (int j = 0; j < taxa[i].getTaxonCount(); j++) {
            if (TaxonList.Utils.hasAttribute(taxa[i], j, dr.evolution.util.Date.DATE)) {
                usingDates = true;
                dr.evolution.util.Date date = (dr.evolution.util.Date) taxa[i].getTaxonAttribute(j, dr.evolution.util.Date.DATE);
                if ((date != null) && (mostRecent == null || date.after(mostRecent))) {
                    mostRecent = date;
                }
            } else {
                // assume contemporaneous tips
                nodes[i][j].setHeight(0.0);
            }
        }
    }
    if (usingDates) {
        assert mostRecent != null;
        TimeScale timeScale = new TimeScale(mostRecent.getUnits(), true, mostRecent.getAbsoluteTimeValue());
        for (int i = 0; i < taxa.length; i++) {
            for (int j = 0; j < taxa[i].getTaxonCount(); j++) {
                dr.evolution.util.Date date = (dr.evolution.util.Date) taxa[i].getTaxonAttribute(j, dr.evolution.util.Date.DATE);
                if (date == null) {
                    throw new IllegalArgumentException("Taxon, " + taxa[i].getTaxonId(j) + ", is missing its date");
                }
                nodes[i][j].setHeight(timeScale.convertTime(date.getTimeValue(), date));
            }
            if (demoFunctions[0].getUnits() != mostRecent.getUnits()) {
            //throw new IllegalArgumentException("The units of the demographic model and the most recent date must match!");
            }
        }
    }
    return new SimpleTree(simulateCoalescent(nodes, demoFunctions, colourChangeMatrix));
}

Example 10 with SimpleTree

use of dr.evolution.tree.SimpleTree in project beast-mcmc by beast-dev.

the class AlloppDiploidHistory method makesimpletree.

// for testing
private SimpleTree makesimpletree() {
    SimpleNode[] snodes = new SimpleNode[dhnodes.length];
    for (int n = 0; n < dhnodes.length; n++) {
        snodes[n] = new SimpleNode();
        // I use taxon==null to identify joined leg node when removing hybtips
        snodes[n].setTaxon(null);
    }
    makesimplesubtree(snodes, 0, dhnodes[rootn]);
    return new SimpleTree(snodes[dhnodes.length - 1]);
}