Examples with TreeInterface beast.evolution.tree.TreeInterface used on opensource projects

Example 11 with TreeInterface

use of beast.evolution.tree.TreeInterface in project beast2 by CompEvol.

the class CalibratedBirthDeathModel method compatibleInitialTree.

Tree compatibleInitialTree() throws MathException {
    final int calCount = orderedCalibrations.length;
    final double[] lowBound = new double[calCount];
    final double[] cladeHeight = new double[calCount];
    // get lower  bound: max(lower bound of dist , bounds of nested clades)
    for (int k = 0; k < calCount; ++k) {
        final CalibrationPoint cal = orderedCalibrations[k];
        lowBound[k] = cal.dist().inverseCumulativeProbability(0);
        // those are node heights
        if (lowBound[k] < 0) {
            lowBound[k] = 0;
        }
        for (final int i : taxaPartialOrder[k]) {
            lowBound[k] = Math.max(lowBound[k], lowBound[i]);
        }
        cladeHeight[k] = cal.dist().inverseCumulativeProbability(1);
    }
    for (int k = calCount - 1; k >= 0; --k) {
        // cladeHeight[k] should be the upper bound of k
        double upper = cladeHeight[k];
        if (Double.isInfinite(upper)) {
            upper = lowBound[k] + 1;
        }
        cladeHeight[k] = (upper + lowBound[k]) / 2.0;
        for (final int i : taxaPartialOrder[k]) {
            cladeHeight[i] = Math.min(cladeHeight[i], cladeHeight[k]);
        }
    }
    final TreeInterface tree = treeInput.get();
    final int nodeCount = tree.getLeafNodeCount();
    final boolean[] used = new boolean[nodeCount];
    int curLeaf = -1;
    int curInternal = nodeCount - 1;
    final Node[] subTree = new Node[calCount];
    for (int k = 0; k < calCount; ++k) {
        final List<Integer> freeTaxa = new ArrayList<>();
        for (final int ti : xclades[k]) {
            freeTaxa.add(ti);
        }
        for (final int i : taxaPartialOrder[k]) {
            for (final int u : xclades[i]) {
                freeTaxa.remove(new Integer(u));
            }
        }
        final List<Node> sbs = new ArrayList<>();
        for (final int i : freeTaxa) {
            final Node n = new Node(tree.getNode(i).getID());
            n.setNr(++curLeaf);
            n.setHeight(0.0);
            sbs.add(n);
            used[i] = true;
        }
        for (final int i : taxaPartialOrder[k]) {
            sbs.add(subTree[i]);
        }
        final double base = sbs.get(sbs.size() - 1).getHeight();
        final double step = (cladeHeight[k] - base) / (sbs.size() - 1);
        Node tr = sbs.get(0);
        for (int i = 1; i < sbs.size(); ++i) {
            tr = Node.connect(tr, sbs.get(i), base + i * step);
            tr.setNr(++curInternal);
        }
        subTree[k] = tr;
    }
    Node finalTree = subTree[calCount - 1];
    double h = cladeHeight[calCount - 1];
    for (int k = 0; k < calCount - 1; ++k) {
        final Node s = subTree[k];
        h = Math.max(h, cladeHeight[k]) + 1;
        finalTree = Node.connect(finalTree, s, h);
        finalTree.setNr(++curInternal);
    }
    for (int k = 0; k < used.length; ++k) {
        if (!used[k]) {
            final String tx = tree.getNode(k).getID();
            final Node n = new Node(tx);
            n.setHeight(0.0);
            n.setNr(++curLeaf);
            finalTree = Node.connect(finalTree, n, h + 1);
            finalTree.setNr(++curInternal);
            h += 1;
        }
    }
    final Tree t = new Tree();
    t.setRoot(finalTree);
    t.initAndValidate();
    return t;
}

Example 12 with TreeInterface

use of beast.evolution.tree.TreeInterface in project beast2 by CompEvol.

the class CalibratedYuleModel method log.

@Override
public void log(final long sample, final PrintStream out) {
    out.print(getCurrentLogP() + "\t");
    if (calcCalibrations) {
        final TreeInterface tree = treeInput.get();
        for (int k = 0; k < orderedCalibrations.length; ++k) {
            final CalibrationPoint cal = orderedCalibrations[k];
            Node c;
            final int[] taxk = xclades[k];
            if (taxk.length > 1) {
                // find MRCA of taxa
                c = getCommonAncestor(tree, taxk);
            } else {
                c = tree.getNode(taxk[0]);
            }
            if (cal.forParent()) {
                c = c.getParent();
            }
            final double h = c.getHeight();
            out.print(h + "\t");
        }
    }
}

Example 13 with TreeInterface

use of beast.evolution.tree.TreeInterface in project beast2 by CompEvol.

the class CalibratedYuleModel method compatibleInitialTree.

public Tree compatibleInitialTree() throws MathException {
    final int calCount = orderedCalibrations.length;
    final double[] lowBound = new double[calCount];
    final double[] cladeHeight = new double[calCount];
    // get lower  bound: max(lower bound of dist , bounds of nested clades)
    for (int k = 0; k < calCount; ++k) {
        final CalibrationPoint cal = orderedCalibrations[k];
        final ParametricDistribution dist = cal.dist();
        // final double offset = dist.getOffset();
        lowBound[k] = dist.inverseCumulativeProbability(0);
        // those are node heights
        if (lowBound[k] < 0) {
            lowBound[k] = 0;
        }
        for (final int i : taxaPartialOrder[k]) {
            lowBound[k] = Math.max(lowBound[k], lowBound[i]);
        }
        cladeHeight[k] = dist.inverseCumulativeProbability(1);
    // if (! Double.isInfinite(cladeHeight[k])) {
    // cladeHeight[k] += offset;
    // }
    }
    for (int k = calCount - 1; k >= 0; --k) {
        // cladeHeight[k] should be the upper bound of k
        double upper = cladeHeight[k];
        if (Double.isInfinite(upper)) {
            upper = lowBound[k] + 1;
        }
        cladeHeight[k] = (upper + lowBound[k]) / 2.0;
        for (final int i : taxaPartialOrder[k]) {
            cladeHeight[i] = Math.min(cladeHeight[i], cladeHeight[k]);
        }
    }
    final TreeInterface tree = treeInput.get();
    final int nodeCount = tree.getLeafNodeCount();
    final boolean[] used = new boolean[nodeCount];
    int curLeaf = -1;
    int curInternal = nodeCount - 1;
    final Node[] subTree = new Node[calCount];
    for (int k = 0; k < calCount; ++k) {
        final List<Integer> freeTaxa = new ArrayList<>();
        for (final int ti : xclades[k]) {
            freeTaxa.add(ti);
        }
        for (final int i : taxaPartialOrder[k]) {
            for (final int u : xclades[i]) {
                freeTaxa.remove(new Integer(u));
            }
        }
        final List<Node> sbs = new ArrayList<>();
        for (final int i : freeTaxa) {
            final Node n = new Node(tree.getNode(i).getID());
            n.setNr(++curLeaf);
            n.setHeight(0.0);
            sbs.add(n);
            used[i] = true;
        }
        for (final int i : taxaPartialOrder[k]) {
            sbs.add(subTree[i]);
            subTree[i] = null;
        }
        final double base = sbs.get(sbs.size() - 1).getHeight();
        final double step = (cladeHeight[k] - base) / (sbs.size() - 1);
        Node tr = sbs.get(0);
        for (int i = 1; i < sbs.size(); ++i) {
            tr = Node.connect(tr, sbs.get(i), base + i * step);
            tr.setNr(++curInternal);
        }
        subTree[k] = tr;
    }
    Node finalTree = subTree[calCount - 1];
    double h = cladeHeight[calCount - 1];
    for (int k = 0; k < calCount - 1; ++k) {
        final Node s = subTree[k];
        if (s != null) {
            h = Math.max(h, cladeHeight[k]) + 1;
            finalTree = Node.connect(finalTree, s, h);
            finalTree.setNr(++curInternal);
        }
    }
    for (int k = 0; k < used.length; ++k) {
        if (!used[k]) {
            final String tx = tree.getNode(k).getID();
            final Node n = new Node(tx);
            n.setHeight(0.0);
            n.setNr(++curLeaf);
            finalTree = Node.connect(finalTree, n, h + 1);
            finalTree.setNr(++curInternal);
            h += 1;
        }
    }
    final Tree t = new Tree();
    t.setRoot(finalTree);
    t.initAndValidate();
    return t;
}

Example 14 with TreeInterface

use of beast.evolution.tree.TreeInterface in project beast2 by CompEvol.

the class TreeLikelihood method calculateLogP.

@Override
public double calculateLogP() {
    if (beagle != null) {
        logP = beagle.calculateLogP();
        return logP;
    }
    final TreeInterface tree = treeInput.get();
    try {
        if (traverse(tree.getRoot()) != Tree.IS_CLEAN)
            calcLogP();
    } catch (ArithmeticException e) {
        return Double.NEGATIVE_INFINITY;
    }
    m_nScale++;
    if (logP > 0 || (likelihoodCore.getUseScaling() && m_nScale > X)) {
    // System.err.println("Switch off scaling");
    // m_likelihoodCore.setUseScaling(1.0);
    // m_likelihoodCore.unstore();
    // m_nHasDirt = Tree.IS_FILTHY;
    // X *= 2;
    // traverse(tree.getRoot());
    // calcLogP();
    // return logP;
    } else if (logP == Double.NEGATIVE_INFINITY && m_fScale < 10 && !scaling.get().equals(Scaling.none)) {
        // && !m_likelihoodCore.getUseScaling()) {
        m_nScale = 0;
        m_fScale *= 1.01;
        Log.warning.println("Turning on scaling to prevent numeric instability " + m_fScale);
        likelihoodCore.setUseScaling(m_fScale);
        likelihoodCore.unstore();
        hasDirt = Tree.IS_FILTHY;
        traverse(tree.getRoot());
        calcLogP();
        return logP;
    }
    return logP;
}

Example 15 with TreeInterface

use of beast.evolution.tree.TreeInterface in project beast2 by CompEvol.

the class SpeciesTreeDistribution method calculateLogP.

// SpeciesTreeDistribution extends TreeDistribution
/**
 * Calculates the log likelihood of this set of coalescent intervals,
 * given a demographic model.
 *
 * @return the log likelihood
 */
@Override
public double calculateLogP() {
    // (Q2R): what if tree intervals?
    // (Q2R): always the same tree, no? so why pass in argument
    final TreeInterface tree = treeInput.get();
    logP = calculateTreeLogLikelihood(tree);
    return logP;
}