Examples with AlloppDiploidHistory dr.evomodel.alloppnet.speciation.AlloppDiploidHistory used on opensource projects

Example 1 with AlloppDiploidHistory

use of dr.evomodel.alloppnet.speciation.AlloppDiploidHistory in project beast-mcmc by beast-dev.

the class AlloppNetworkNodeSlide method operateOneNodeInTetraTree.

private void operateOneNodeInTetraTree(AlloppLeggedTree tettree, int which, double factor) {
    // As TreeNodeSlide(). Randomly flip children at each node,
    // keeping track of node order (in-order order, left to right).
    NodeRef[] order = SlidableTree.Utils.mnlCanonical(tettree);
    // Find the time of the most recent gene coalescence which
    // has (species,sequence)'s to left and right of this node.
    FixedBitSet left = apsp.speciesseqEmptyUnion();
    FixedBitSet right = apsp.speciesseqEmptyUnion();
    for (int k = 0; k < 2 * which + 1; k += 2) {
        FixedBitSet left0 = apsp.taxonseqToTipUnion(tettree.getSlidableNodeTaxon(order[k]), 0);
        FixedBitSet left1 = apsp.taxonseqToTipUnion(tettree.getSlidableNodeTaxon(order[k]), 1);
        left.union(left0);
        left.union(left1);
    }
    for (int k = 2 * (which + 1); k < order.length; k += 2) {
        FixedBitSet right0 = apsp.taxonseqToTipUnion(tettree.getSlidableNodeTaxon(order[k]), 0);
        FixedBitSet right1 = apsp.taxonseqToTipUnion(tettree.getSlidableNodeTaxon(order[k]), 1);
        right.union(right0);
        right.union(right1);
    }
    double genelimit = apsp.spseqUpperBound(left, right);
    // also keep this node more recent than the hybridization event that led to this tree.
    AlloppDiploidHistory diphist = apspnet.getDiploidHistory();
    double hybridheight = diphist.getHybHeight(tettree);
    final double limit = Math.min(genelimit, hybridheight);
    // On direct call, factor==0.0 and use limit. Else use passed in scaling factor
    double newHeight = -1.0;
    if (factor > 0) {
        newHeight = tettree.getSlidableNodeHeight(order[2 * which + 1]) * factor;
    } else {
        newHeight = MathUtils.nextDouble() * limit;
    }
    apspnet.beginNetworkEdit();
    final NodeRef node = order[2 * which + 1];
    tettree.setSlidableNodeHeight(node, newHeight);
    SlidableTree.Utils.mnlReconstruct(tettree, order);
    apspnet.endNetworkEdit();
}

Example 2 with AlloppDiploidHistory

use of dr.evomodel.alloppnet.speciation.AlloppDiploidHistory in project beast-mcmc by beast-dev.

the class AlloppNetworkNodeSlide method operateHybridHeight.

private void operateHybridHeight(int footindex) {
    AlloppDiploidHistory diphist = apspnet.getDiploidHistory();
    ArrayList<Integer> feet = diphist.collectFeet();
    assert footindex < feet.size();
    int foot = feet.get(footindex);
    int tt = diphist.getNodeTettree(foot);
    AlloppLeggedTree tettree = apspnet.getTetraploidTree(tt);
    double minh = tettree.getRootHeight();
    int f1 = tettree.getDiphistLftLeg();
    int f2 = tettree.getDiphistRgtLeg();
    assert (foot == f1) || (foot == f2);
    apspnet.beginNetworkEdit();
    diphist.moveHybridHeight(f1, f2, minh);
    apspnet.endNetworkEdit();
}

Example 3 with AlloppDiploidHistory

use of dr.evomodel.alloppnet.speciation.AlloppDiploidHistory in project beast-mcmc by beast-dev.

the class AlloppNetworkNodeSlide method operateOneNodeInDiploidHistory.

private void operateOneNodeInDiploidHistory(int which, double factor) {
    apspnet.beginNetworkEdit();
    int slidingn = 2 * which + 1;
    AlloppDiploidHistory diphist = apspnet.getDiploidHistory();
    NodeRef[] order = SlidableTree.Utils.mnlCanonical(diphist);
    // Find the time of the most recent gene coalescence which
    // has (species,sequence)'s to left and right of this node.
    FixedBitSet left = apsp.speciesseqEmptyUnion();
    FixedBitSet right = apsp.speciesseqEmptyUnion();
    for (int k = 0; k < slidingn; k += 2) {
        FixedBitSet u = apspnet.calculateDipHistTipUnion(order[k]);
        left.union(u);
    }
    for (int k = slidingn + 1; k < order.length; k += 2) {
        FixedBitSet u = apspnet.calculateDipHistTipUnion(order[k]);
        right.union(u);
    }
    double genelimit = apsp.spseqUpperBound(left, right);
    // find limit due to hyb-tips - must be bigger than adjacent heights
    // Note that adjacent nodes in order[] are tips; if they are not hyb-tips
    // they have height zero anyway.
    double hybtiplimit = 0.0;
    if (slidingn - 1 >= 0) {
        hybtiplimit = Math.max(hybtiplimit, diphist.getSlidableNodeHeight(order[slidingn - 1]));
    }
    if (slidingn + 1 < order.length) {
        hybtiplimit = Math.max(hybtiplimit, diphist.getSlidableNodeHeight(order[slidingn + 1]));
    }
    RootHeightRange rootrange = new RootHeightRange(0.0, Double.MAX_VALUE);
    if (apspnet.getDiploidRootIsRoot()) {
        rootrange = findRootRangeForDiploidRootIsRoot(diphist, order, slidingn);
    }
    final double upperlimit = Math.min(genelimit, rootrange.upperlimit);
    final double lowerlimit = Math.max(hybtiplimit, rootrange.lowerlimit);
    // On direct call, factor==0.0 and use limit. Else use passed in scaling factor
    double newHeight = -1.0;
    if (factor > 0) {
        newHeight = diphist.getSlidableNodeHeight(order[slidingn]) * factor;
    } else {
        newHeight = MathUtils.uniform(lowerlimit, upperlimit);
    }
    assert diphist.diphistOK(apspnet.getDiploidRootIsRoot());
    final NodeRef node = order[slidingn];
    diphist.setSlidableNodeHeight(node, newHeight);
    SlidableTree.Utils.mnlReconstruct(diphist, order);
    if (!diphist.diphistOK(apspnet.getDiploidRootIsRoot())) {
        System.out.println("BUG in operateOneNodeInDiploidHistory()");
    }
    assert diphist.diphistOK(apspnet.getDiploidRootIsRoot());
    apspnet.endNetworkEdit();
}