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

Example 6 with Node

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

the class TipDatesRandomWalker method proposal.

@Override
public double proposal() {
    // randomly select leaf node
    int i = Randomizer.nextInt(taxonIndices.length);
    Node node = treeInput.get().getNode(taxonIndices[i]);
    double value = node.getHeight();
    double newValue = value;
    if (useGaussian) {
        newValue += Randomizer.nextGaussian() * windowSize;
    } else {
        newValue += Randomizer.nextDouble() * 2 * windowSize - windowSize;
    }
    if (newValue > node.getParent().getHeight()) {
        // || newValue < 0.0) {
        if (reflectValue) {
            newValue = reflectValue(newValue, 0.0, node.getParent().getHeight());
        } else {
            return Double.NEGATIVE_INFINITY;
        }
    }
    if (newValue == value) {
        // this saves calculating the posterior
        return Double.NEGATIVE_INFINITY;
    }
    node.setHeight(newValue);
    return 0.0;
}

Example 7 with Node

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

the class WilsonBalding method proposal.

/**
 * WARNING: Assumes strictly bifurcating beast.tree.
 */
/**
 * override this for proposals,
 *
 * @return log of Hastings Ratio, or Double.NEGATIVE_INFINITY if proposal should not be accepted *
 */
@Override
public double proposal() {
    Tree tree = treeInput.get(this);
    double oldMinAge, newMinAge, newRange, oldRange, newAge, hastingsRatio;
    // choose a random node avoiding root
    final int nodeCount = tree.getNodeCount();
    Node i;
    do {
        i = tree.getNode(Randomizer.nextInt(nodeCount));
    } while (i.isRoot());
    final Node p = i.getParent();
    // choose another random node to insert i above
    Node j;
    Node jP;
    // make sure that the target branch <k, j> is above the subtree being moved
    do {
        j = tree.getNode(Randomizer.nextInt(nodeCount));
        jP = j.getParent();
    } while ((jP != null && jP.getHeight() <= i.getHeight()) || (i.getNr() == j.getNr()));
    // disallow moves that change the root.
    if (j.isRoot() || p.isRoot()) {
        return Double.NEGATIVE_INFINITY;
    }
    // j != root tested above
    assert jP != null;
    final int pnr = p.getNr();
    final int jPnr = jP.getNr();
    if (jPnr == pnr || j.getNr() == pnr || jPnr == i.getNr())
        return Double.NEGATIVE_INFINITY;
    final Node CiP = getOtherChild(p, i);
    final Node PiP = p.getParent();
    newMinAge = Math.max(i.getHeight(), j.getHeight());
    newRange = jP.getHeight() - newMinAge;
    newAge = newMinAge + (Randomizer.nextDouble() * newRange);
    oldMinAge = Math.max(i.getHeight(), CiP.getHeight());
    oldRange = PiP.getHeight() - oldMinAge;
    hastingsRatio = newRange / Math.abs(oldRange);
    if (oldRange == 0 || newRange == 0) {
        // For symmetry, newRange = 0 should therefore be ruled out as well
        return Double.NEGATIVE_INFINITY;
    }
    // root changing moves disallowed earlier. Comment out just like in BEAST 1
    // update
    // if (j.isRoot()) {
    // // 1. remove edges <p, CiP>
    // // 2. add edges <k, p>, <p, j>, <PiP, CiP>
    // 
    // replace(p, CiP, j);
    // replace(PiP, p, CiP);
    // 
    // // p is the new root
    // p.setParent(null);
    // tree.setRoot(p);
    // 
    // } else if (p.isRoot()) {
    // // 1. remove edges <k, j>, <p, CiP>, <PiP, p>
    // // 2. add edges <k, p>, <p, j>, <PiP, CiP>
    // 
    // replace(jP, j, p);
    // //replace(p, CiP, p);
    // replace(p, CiP, j);
    // 
    // // CiP is the new root
    // CiP.setParent(null);
    // tree.setRoot(CiP);
    // 
    // } else {
    // // 1. remove edges <k, j>, <p, CiP>, <PiP, p>
    // 2. add edges <k, p>, <p, j>, <PiP, CiP>
    // disconnect p
    final Node pP = p.getParent();
    replace(pP, p, CiP);
    // re-attach, first child node to p
    replace(p, CiP, j);
    // then parent node of j to p
    replace(jP, j, p);
    // mark paths to common ancestor as changed
    if (markCladesInput.get()) {
        Node iup = pP;
        Node jup = p;
        while (iup != jup) {
            if (iup.getHeight() < jup.getHeight()) {
                assert !iup.isRoot();
                iup = iup.getParent();
                iup.makeDirty(Tree.IS_FILTHY);
            } else {
                assert !jup.isRoot();
                jup = jup.getParent();
                jup.makeDirty(Tree.IS_FILTHY);
            }
        }
    }
    // }
    p.setHeight(newAge);
    return Math.log(hastingsRatio);
}

Example 8 with Node

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

the class CalibratedBirthDeathModel 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 9 with Node

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

the class CalibratedBirthDeathModel method getCorrection.

public double getCorrection(final TreeInterface tree, final double lam, final double a, final double rho) {
    double logL = 0.0;
    final int calCount = orderedCalibrations.length;
    final double[] hs = new double[calCount];
    for (int k = 0; k < calCount; ++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);
            // only monophyletics clades can be calibrated
            if (getLeafCount(c) != taxk.length) {
                return Double.NEGATIVE_INFINITY;
            }
        } else {
            c = tree.getNode(taxk[0]);
            assert cal.forParent();
        }
        if (cal.forParent()) {
            c = c.getParent();
        }
        final double h = c.getHeight();
        // add calibration density for point
        if (calcCalibrations) {
            logL += cal.logPdf(h);
        }
        hs[k] = h;
    }
    if (Double.isInfinite(logL)) {
        // some calibration points out of range
        return logL;
    }
    if (type == Type.NONE) {
        return logL;
    }
    if (userPDF == null) {
        final int taxonCount = tree.getLeafNodeCount();
        switch(type) {
            case OVER_ALL_TOPOS:
                {
                    if (calCount == 1 && isYule) {
                        logL -= logMarginalDensity(lam, taxonCount, hs[0], xclades[0].length, orderedCalibrations[0].forParent());
                    } else if (calCount == 2 && taxaPartialOrder[1].length == 1 && isYule) {
                        // assert !forParent[0] && !forParent[1];
                        logL -= logMarginalDensity(lam, taxonCount, hs[0], xclades[0].length, hs[1], xclades[1].length);
                    } else {
                        if (lastLam == lam) {
                            int k = 0;
                            for (; k < hs.length; ++k) {
                                if (hs[k] != lastHeights[k]) {
                                    break;
                                }
                            }
                            if (k == hs.length) {
                                return lastValue;
                            }
                        }
                        // the slow and painful way
                        final double[] hss = new double[hs.length];
                        final int[] ranks = new int[hs.length];
                        for (int k = 0; k < hs.length; ++k) {
                            int r = 0;
                            for (int j = 0; j < k; ++j) {
                                r += (hs[j] <= hs[k]) ? 1 : 0;
                            }
                            for (int j = k + 1; j < hs.length; ++j) {
                                r += (hs[j] < hs[k]) ? 1 : 0;
                            }
                            // for (final double h : hs) {
                            // r += (h < hs[k]) ? 1 : 0;
                            // }
                            ranks[k] = r + 1;
                            hss[r] = hs[k];
                        }
                        logL -= logMarginalDensity(lam, hss, ranks, linsIter);
                        lastLam = lam;
                        System.arraycopy(hs, 0, lastHeights, 0, lastHeights.length);
                        lastValue = logL;
                    }
                    break;
                }
            case OVER_RANKED_COUNTS:
                {
                    final int[] ranks = new int[hs.length];
                    for (int k = 0; k < hs.length; ++k) {
                        int r = 0;
                        for (int j = 0; j < k; ++j) {
                            r += (hs[j] <= hs[k]) ? 1 : 0;
                        }
                        for (int j = k + 1; j < hs.length; ++j) {
                            r += (hs[j] < hs[k]) ? 1 : 0;
                        }
                        ranks[k] = r + 1;
                    }
                    final double nt = countTrees(ranks, linsIter);
                    logL -= nt;
                    Arrays.sort(hs);
                    final int[] cs = new int[calCount + 1];
                    for (final Node n : tree.getInternalNodes()) {
                        final double nhk = n.getHeight();
                        int i = 0;
                        for (; /**/
                        i < hs.length; ++i) {
                            if (hs[i] >= nhk) {
                                break;
                            }
                        }
                        if (i == hs.length) {
                            cs[i]++;
                        } else {
                            if (nhk < hs[i]) {
                                cs[i]++;
                            }
                        }
                    }
                    if (isYule) {
                        double ll = 0;
                        ll += cs[0] * Math.log1p(-Math.exp(-lam * hs[0])) - lam * hs[0] - lfactorials[cs[0]];
                        for (int i = 1; i < cs.length - 1; ++i) {
                            final int c = cs[i];
                            ll += c * (Math.log1p(-Math.exp(-lam * (hs[i] - hs[i - 1]))) - lam * hs[i - 1]);
                            ll += -lam * hs[i] - lfactorials[c];
                        }
                        ll += -lam * (cs[calCount] + 1) * hs[calCount - 1] - lfactorials[cs[calCount] + 1];
                        ll += Math.log(lam) * calCount;
                        logL -= ll;
                    } else {
                        assert (a < 1);
                        // We try to minimize the effects of numerical unstability,
                        // so the evaluations look different than the nice equations in the MS.
                        final double r = lam;
                        final double ia = 1 - a;
                        final double ira = ia - rho;
                        final double v = (ia * ia) / r;
                        // accumulate terms here
                        double ll = 0;
                        for (int i = 0; i < cs.length - 1; ++i) {
                            final int c = cs[i];
                            final double mrh = -r * hs[i];
                            final double mrh1 = (i > 0) ? -r * hs[i - 1] : 0;
                            ll -= lfactorials[c];
                            // double xx =  Math.log(1 - Math.exp(mrh - mrh1));
                            // can we have accurate  log(1+e^x) ??
                            final double xxo = Math.log1p(-Math.exp(mrh - mrh1));
                            final double xx1 = rho + ira * Math.exp(mrh);
                            final double xx2 = rho + ira * Math.exp(mrh1);
                            ll += c * (mrh1 + xxo + Math.log(v / (xx1 * xx2)));
                            // contribution of  calibration point
                            final double z = Math.log(ia / (rho + ira * Math.exp(mrh)));
                            ll += 2 * z + mrh;
                        }
                        // from last calibration to oo contribution
                        final double xu = ia / (rho + ira * Math.exp(-r * hs[calCount - 1]));
                        // fixme
                        final int n2 = cs[calCount];
                        // fixme
                        assert n2 > 0;
                        final double rhoia = rho / ia;
                        ll += -lfactorials[n2 + 1] - (n2) * Math.log(r * rhoia) - (n2 + 1) * r * hs[calCount - 1] + (n2 + 1) * Math.log(xu);
                        // non-node terms
                        ll += lfactorials[taxonCount] + (taxonCount - 1) * Math.log(r * rhoia);
                        logL -= ll;
                    }
                    break;
                }
            default:
                break;
        }
    } else {
        final double value = userPDF.getArrayValue();
        if (Double.isNaN(value) || Double.isInfinite(value)) {
            logL = Double.NEGATIVE_INFINITY;
        } else {
            logL -= value;
        }
    }
    return logL;
}

Example 10 with Node

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

the class CalibratedYuleModel method calculateYuleLikelihood.

private static double calculateYuleLikelihood(final TreeInterface tree, final double lam) {
    final int taxonCount = tree.getLeafNodeCount();
    // add all lambda multipliers here
    // No normalization at the moment.  for n! use logGamma(taxonCount + 1);
    double logL = (taxonCount - 1) * Math.log(lam);
    final Node[] nodes = tree.getNodesAsArray();
    for (int i = taxonCount; i < nodes.length; i++) {
        final Node node = nodes[i];
        assert (!node.isLeaf());
        final double height = node.getHeight();
        final double mrh = -lam * height;
        logL += mrh + (node.isRoot() ? mrh : 0);
    }
    return logL;
}