Examples with FixedBitSet jebl.util.FixedBitSet used on opensource projects

Example 11 with FixedBitSet

use of jebl.util.FixedBitSet in project beast-mcmc by beast-dev.

the class SpeciesTreeModel method isCompatible.

// Not very efficient, should do something better, based on traversing the cList once
private int isCompatible(NodeRef node, SpeciesBindings.CoalInfo[] cList, int loc) {
    if (!isExternal(node)) {
        int l = -1;
        for (int nc = 0; nc < getChildCount(node); ++nc) {
            int l1 = isCompatible(getChild(node, nc), cList, loc);
            if (l1 < 0) {
                return -1;
            }
            assert l == -1 || l1 == l;
            l = l1;
        }
        loc = l;
        assert cList[loc].ctime >= getNodeHeight(node);
    }
    if (node == getRoot()) {
        return cList.length;
    }
    // spSet guaranteed to be ready by caller
    final FixedBitSet nodeSps = props.get(node).spSet;
    final double limit = getNodeHeight(getParent(node));
    while (loc < cList.length) {
        final SpeciesBindings.CoalInfo ci = cList[loc];
        if (ci.ctime >= limit) {
            break;
        }
        boolean allIn = true, noneIn = true;
        for (int i = 0; i < 2; ++i) {
            final FixedBitSet s = ci.sinfo[i];
            final int in1 = s.intersectCardinality(nodeSps);
            if (in1 > 0) {
                noneIn = false;
            }
            if (s.cardinality() != in1) {
                allIn = false;
            }
        }
        if (!(allIn || noneIn)) {
            return -1;
        }
        ++loc;
    }
    return loc;
}

Example 12 with FixedBitSet

use of jebl.util.FixedBitSet 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();
}

Example 13 with FixedBitSet

use of jebl.util.FixedBitSet in project beast-mcmc by beast-dev.

the class AlloppChangeNumHybridizations method splitTettree.

private double splitTettree(int tt, AlloppNode root1, AlloppNode root2) {
    double hr = 0.0;
    // collect info from old TetraTree
    AlloppLeggedTree tetTree = apspnet.getTetraploidTree(tt);
    AlloppDiploidHistory adhist = apspnet.getDiploidHistory();
    double rooth = tetTree.getRootHeight();
    int lftleg = tetTree.getDiphistLftLeg();
    int rgtleg = tetTree.getDiphistRgtLeg();
    double lftanchgt = adhist.getAncHeight(lftleg);
    double rgtanchgt = adhist.getAncHeight(rgtleg);
    // account for the hybhgt that will be lost
    hr += Math.log(uniformpdf(rooth, Math.min(lftanchgt, rgtanchgt)));
    // make two new trees
    AlloppLeggedTree tetTree1 = new AlloppLeggedTree(tetTree, root1);
    AlloppLeggedTree tetTree2 = new AlloppLeggedTree(tetTree, root2);
    // tetree2 gets old one's legs, with new height
    tetTree2.setDiphistLftLeg(lftleg);
    tetTree2.setDiphistRgtLeg(rgtleg);
    double hybhgt2 = MathUtils.uniform(tetTree2.getRootHeight(), rooth);
    hr -= Math.log(uniformpdf(tetTree2.getRootHeight(), rooth));
    adhist.setHybridHeight(tetTree2, hybhgt2);
    // remove old and add new ones to list.
    // tetTree2 replaces tetTree, that is, same index, so dip tips stay consistent
    apspnet.setTetTree(tt, tetTree2);
    int tt2 = tt;
    int tt1 = apspnet.addTetTree(tetTree1);
    // new hybhgt for tree1
    double hybhgt1 = MathUtils.uniform(tetTree1.getRootHeight(), rooth);
    hr -= Math.log(uniformpdf(tetTree1.getRootHeight(), rooth));
    // it is constrained by gene trees and existing node height
    if (MathUtils.nextBoolean()) {
        FixedBitSet tt1leg0 = apspnet.unionOfWholeTetTree(tt1, 0);
        FixedBitSet tt2leg0 = apspnet.unionOfWholeTetTree(tt2, 0);
        double genelimit = apsp.spseqUpperBound(tt1leg0, tt2leg0);
        double maxfootanchgt = Math.min(genelimit, lftanchgt);
        double footanchgt = MathUtils.uniform(rooth, maxfootanchgt);
        hr -= Math.log(uniformpdf(rooth, maxfootanchgt));
        adhist.addTwoDipTips(apspnet, tt1, tt2, footanchgt, rooth, hybhgt1);
    } else {
        FixedBitSet tt1leg1 = apspnet.unionOfWholeTetTree(tt1, 1);
        FixedBitSet tt2leg1 = apspnet.unionOfWholeTetTree(tt2, 1);
        double genelimit = apsp.spseqUpperBound(tt1leg1, tt2leg1);
        double maxfootanchgt = Math.min(genelimit, rgtanchgt);
        double footanchgt = MathUtils.uniform(rooth, maxfootanchgt);
        hr -= Math.log(uniformpdf(rooth, maxfootanchgt));
        adhist.addTwoDipTips(apspnet, tt1, tt2, rooth, footanchgt, hybhgt1);
    }
    // grjtodo-soon The only difference between the two states is the time-order of the nodes.
    // Should topologies or histories be counted?
    // Account for left/right choice. This says histories
    hr += Math.log(2.0);
    return hr;
}

Example 14 with FixedBitSet

use of jebl.util.FixedBitSet in project beast-mcmc by beast-dev.

the class AlloppChangeNumHybridizations method mergeTettreePair.

private double mergeTettreePair(int tt1, int tt2) {
    double hr = 0.0;
    AlloppLeggedTree ttree1 = apspnet.getTetraploidTree(tt1);
    AlloppLeggedTree ttree2 = apspnet.getTetraploidTree(tt2);
    AlloppDiploidHistory adhist = apspnet.getDiploidHistory();
    // collect height info
    AlloppDiploidHistory.FootAncHeights lftleg2 = adhist.intervalOfFootAncestor(ttree2, AlloppDiploidHistory.LegLorR.left);
    AlloppDiploidHistory.FootAncHeights rgtleg2 = adhist.intervalOfFootAncestor(ttree2, AlloppDiploidHistory.LegLorR.right);
    // Choose most recent footanc height as root height of merged tree.
    // Account for loss of the other footanc height.
    // Use gene limit on the lost footanc height for hr calculation.
    // grjtodo-soon test the gene limit calculation somehow
    double rooth;
    if (lftleg2.anchgt < rgtleg2.anchgt) {
        rooth = lftleg2.anchgt;
        FixedBitSet tt1leg1 = apspnet.unionOfWholeTetTree(tt1, 1);
        FixedBitSet tt2leg1 = apspnet.unionOfWholeTetTree(tt2, 1);
        double genelimit = apsp.spseqUpperBound(tt1leg1, tt2leg1);
        double maxfootanchgt = Math.min(genelimit, rgtleg2.ancanchgt);
        hr += Math.log(uniformpdf(rooth, maxfootanchgt));
    } else {
        rooth = rgtleg2.anchgt;
        FixedBitSet tt1leg0 = apspnet.unionOfWholeTetTree(tt1, 0);
        FixedBitSet tt2leg0 = apspnet.unionOfWholeTetTree(tt2, 0);
        double genelimit = apsp.spseqUpperBound(tt1leg0, tt2leg0);
        double maxfootanchgt = Math.min(genelimit, lftleg2.ancanchgt);
        hr += Math.log(uniformpdf(rooth, maxfootanchgt));
    }
    // account for loss of two old hybhgts
    hr += Math.log(uniformpdf(ttree1.getRootHeight(), rooth));
    hr += Math.log(uniformpdf(ttree2.getRootHeight(), rooth));
    // merge the trees and replace tt2 with result
    AlloppLeggedTree merged = new AlloppLeggedTree(ttree1, ttree2, rooth);
    apspnet.setTetTree(tt2, merged);
    apspnet.removeTetree(tt1);
    // Fix up the links from diploid history.
    // Get rid of old links first, to enable later assertions
    adhist.clearAllNodeTettree();
    for (int i = 0; i < apspnet.getNumberOfTetraTrees(); i++) {
        AlloppLeggedTree ttree = apspnet.getTetraploidTree(i);
        int dhlftleg = ttree.getDiphistLftLeg();
        assert adhist.getNodeTettree(dhlftleg) == -1;
        adhist.setNodeTettree(dhlftleg, i);
        int dhrgtleg = ttree.getDiphistRgtLeg();
        assert adhist.getNodeTettree(dhrgtleg) == -1;
        adhist.setNodeTettree(dhrgtleg, i);
    }
    // new hybhgt for merged tree
    double maxhybhgt = Math.min(lftleg2.ancanchgt, rgtleg2.ancanchgt);
    double hybght = MathUtils.uniform(rooth, maxhybhgt);
    adhist.setHybridHeight(merged, hybght);
    hr -= Math.log(uniformpdf(rooth, maxhybhgt));
    adhist.removeFeet(apspnet, ttree1);
    return hr;
}

Example 15 with FixedBitSet

use of jebl.util.FixedBitSet in project beast-mcmc by beast-dev.

the class MulSpeciesTreeModel method getDemographicPoints.

//  Assign positions in 'pointsList' for the sub-tree rooted at the ancestor of
//  nodeID.
//
//  pointsList is indexed by node-id. Every element is a list of internal
//  population points for the branch between nodeID and it's ancestor
//
private NodeProperties getDemographicPoints(final NodeRef nodeID, Args args, Points[][] pointsList) {
    final NodeProperties nprop = props.get(nodeID);
    final int nSpSeqs = mulspb.nSpSeqs();
    // Species assignment from the tips never changes
    if (!isExternal(nodeID)) {
        nprop.spSet = new FixedBitSet(nSpSeqs);
        for (int nc = 0; nc < getChildCount(nodeID); ++nc) {
            final NodeProperties p = getDemographicPoints(getChild(nodeID, nc), args, pointsList);
            nprop.spSet.union(p.spSet);
        }
    }
    if (args == null) {
        return nprop;
    }
    // parent height
    final double cHeight = nodeID != getRoot() ? getNodeHeight(getParent(nodeID)) : Double.MAX_VALUE;
    // points along branch
    // not sure what a good default size is?
    List<Points> allPoints = new ArrayList<Points>(5);
    if (bmp) {
        for (int isp = nprop.spSet.nextOnBit(0); isp >= 0; isp = nprop.spSet.nextOnBit(isp + 1)) {
            final double[] cp = args.cps[isp];
            final int upi = singleStartPoints[isp];
            int i = args.iSingle[isp];
            for (; /**/
            i < cp.length && cp[i] < cHeight; ++i) {
                allPoints.add(new Points(cp[i], args.indicators[upi + i] > 0));
            }
            args.iSingle[isp] = i;
        }
    } else {
        for (int isp = nprop.spSet.nextOnBit(0); isp >= 0; isp = nprop.spSet.nextOnBit(isp + 1)) {
            final double nodeHeight = spTree.getNodeHeight(nodeID);
            {
                double[] cp = args.cps[isp];
                final int upi = singleStartPoints[isp];
                int i = args.iSingle[isp];
                while (i < cp.length && cp[i] < cHeight) {
                    if (args.indicators[upi + i] > 0) {
                        //System.out.println("  popbit s");
                        args.iSingle[isp] = i;
                        double prev = args.findPrev(cp[i], nodeHeight);
                        double mid = (prev + cp[i]) / 2.0;
                        assert nodeHeight < mid;
                        allPoints.add(new Points(mid, args.pops[upi + i]));
                    }
                    ++i;
                }
                args.iSingle[isp] = i;
            }
            final int kx = (isp * (2 * nSpSeqs - isp - 3)) / 2 - 1;
            for (int y = nprop.spSet.nextOnBit(isp + 1); y >= 0; y = nprop.spSet.nextOnBit(y + 1)) {
                assert isp < y;
                int k = kx + y;
                double[] cp = args.cpp[k];
                int i = args.iPair[k];
                final int upi = pairStartPoints[k];
                while (i < cp.length && cp[i] < cHeight) {
                    if (args.indicators[upi + i] > 0) {
                        //System.out.println("  popbit p");
                        args.iPair[k] = i;
                        final double prev = args.findPrev(cp[i], nodeHeight);
                        double mid = (prev + cp[i]) / 2.0;
                        assert nodeHeight < mid;
                        allPoints.add(new Points(mid, args.pops[upi + i]));
                    }
                    ++i;
                }
                args.iPair[k] = i;
            }
        }
    }
    Points[] all = null;
    if (allPoints.size() > 0) {
        all = allPoints.toArray(new Points[allPoints.size()]);
        if (all.length > 1) {
            HeapSort.sort(all);
        }
        int len = all.length;
        if (bmp) {
            int k = 0;
            while (k + 1 < len) {
                final double t = all[k].time;
                if (t == all[k + 1].time) {
                    int j = k + 2;
                    boolean use = all[k].use || all[k + 1].use;
                    while (j < len && t == all[j].time) {
                        use = use || all[j].use;
                        j += 1;
                    }
                    int removed = (j - k - 1);
                    all[k] = new Points(t, use);
                    for (int i = k + 1; i < len - removed; ++i) {
                        all[i] = all[i + removed];
                    }
                    len -= removed;
                }
                ++k;
            }
        } else {
            // duplications
            int k = 0;
            while (k + 1 < len) {
                double t = all[k].time;
                if (t == all[k + 1].time) {
                    int j = k + 2;
                    double v = all[k].population + all[k + 1].population;
                    while (j < len && t == all[j].time) {
                        v += all[j].population;
                        j += 1;
                    }
                    int removed = (j - k - 1);
                    all[k] = new Points(t, v / (removed + 1));
                    for (int i = k + 1; i < len - removed; ++i) {
                        all[i] = all[i + removed];
                    }
                    //System.arraycopy(all, j, all, k + 1, all.length - j + 1);
                    len -= removed;
                }
                ++k;
            }
        }
        if (len != all.length) {
            Points[] a = new Points[len];
            System.arraycopy(all, 0, a, 0, len);
            all = a;
        }
        if (bmp) {
            for (Points p : all) {
                double t = p.time;
                assert p.population == 0;
                for (int isp = nprop.spSet.nextOnBit(0); isp >= 0; isp = nprop.spSet.nextOnBit(isp + 1)) {
                    SimpleDemographicFunction d = args.dms[isp];
                    if (t <= d.upperBound()) {
                        p.population += d.population(t);
                    }
                }
            }
        }
    }
    pointsList[nodeID.getNumber()] = all;
    return nprop;
}