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

Example 16 with FixedBitSet

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

the class MulSpeciesTreeModel method isSubtreeCompatible.

// Not very efficient, should do something better, based on traversing the cList once
private int isSubtreeCompatible(NodeRef node, MulSpeciesBindings.CoalInfo[] cList, int loc) {
    if (!isExternal(node)) {
        int l = -1;
        for (int nc = 0; nc < getChildCount(node); ++nc) {
            int l1 = isSubtreeCompatible(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 MulSpeciesBindings.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 17 with FixedBitSet

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

the class AlloppMulLabTree method fillmlnodesforlhoodtest1.

private void fillmlnodesforlhoodtest1() {
    mlnodes = new MulLabNode[7];
    for (int i = 0; i < mlnodes.length; i++) {
        mlnodes[i] = new MulLabNode(i);
    }
    // 4 tips, 2 dips, 1 tet
    mlnodes[0].taxon = new Taxon("a");
    mlnodes[1].taxon = new Taxon("b");
    mlnodes[2].taxon = new Taxon("z0");
    mlnodes[3].taxon = new Taxon("z1");
    mlnodes[2].tetraancestor = true;
    mlnodes[3].tetraancestor = true;
    mlnodes[2].intetratree = true;
    mlnodes[3].intetratree = true;
    mlnodes[2].tetraroot = true;
    mlnodes[3].tetraroot = true;
    mlnodes[0].union = new FixedBitSet(4);
    mlnodes[0].union.set(0);
    mlnodes[1].union = new FixedBitSet(4);
    mlnodes[1].union.set(1);
    mlnodes[2].union = new FixedBitSet(4);
    mlnodes[2].union.set(2);
    mlnodes[3].union = new FixedBitSet(4);
    mlnodes[3].union.set(3);
    // toplogy and times
    mlnodes[4].addChildren(mlnodes[0], mlnodes[2]);
    mlnodes[5].addChildren(mlnodes[1], mlnodes[3]);
    mlnodes[6].addChildren(mlnodes[4], mlnodes[5]);
    rootn = 6;
    mlnodes[0].height = 0.0;
    mlnodes[1].height = 0.0;
    mlnodes[2].height = 0.0;
    mlnodes[3].height = 0.0;
    mlnodes[4].height = 0.01;
    mlnodes[5].height = 0.02;
    mlnodes[6].height = 0.03;
    mlnodes[2].hybridheight = 0.005;
    mlnodes[3].hybridheight = 0.005;
    // nlineages and coalheights
    mlnodes[0].nlineages = 1;
    mlnodes[1].nlineages = 2;
    mlnodes[2].nlineages = 1;
    mlnodes[3].nlineages = 1;
    mlnodes[4].nlineages = 2;
    mlnodes[4].coalheights.add(0.015);
    mlnodes[5].nlineages = 3;
    mlnodes[5].coalheights.add(0.025);
    mlnodes[6].nlineages = 3;
    mlnodes[6].coalheights.add(0.035);
    mlnodes[6].coalheights.add(0.045);
}

Example 18 with FixedBitSet

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

the class SpeciesTreeModel 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 nSpecies = species.nSpecies();
    // Species assignment from the tips never changes
    if (!isExternal(nodeID)) {
        nprop.spSet = new FixedBitSet(nSpecies);
        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 * nSpecies - 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;
}

Example 19 with FixedBitSet

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

the class MultiSpeciesCoalescent method treeLogLikelihood.

private double treeLogLikelihood(SpeciesBindings.GeneTreeInfo geneTree, NodeRef node, int[] info, double popFactor) {
    // number of lineages remaining at node
    int nLineages;
    // location in coalescent list (optimization)
    int indexInClist = 0;
    // accumulated log-likelihood inBranchh from node to it's parent
    double like = 0;
    final double t0 = spTree.getNodeHeight(node);
    final SpeciesBindings.CoalInfo[] cList = geneTree.getCoalInfo();
    if (verbose) {
        if (spTree.isRoot(node)) {
            System.err.println("gtree:" + geneTree.tree.getId());
            System.err.println("t0 " + t0);
            for (int k = 0; k < cList.length; ++k) {
                System.err.println(k + " " + cList[k].ctime + " " + cList[k].sinfo[0] + " " + cList[k].sinfo[1]);
            }
        }
    }
    if (spTree.isExternal(node)) {
        nLineages = geneTree.nLineages(spTree.speciesIndex(node));
        indexInClist = 0;
    } else {
        //assert spTree.getChildCount(node) == 2;
        nLineages = 0;
        for (int nc = 0; nc < 2; ++nc) {
            final NodeRef child = spTree.getChild(node, nc);
            like += treeLogLikelihood(geneTree, child, info, popFactor);
            nLineages += info[0];
            indexInClist = Math.max(indexInClist, info[1]);
        }
        // root of species tree
        if (indexInClist >= cList.length) {
            int k = 1;
        }
        assert indexInClist < cList.length;
        while (cList[indexInClist].ctime < t0) {
            ++indexInClist;
        }
    }
    final boolean isRoot = spTree.isRoot(node);
    // Upper limit
    // use of (t0 + spTree.getBranchLength(node)) caused problem since there was a tiny difference
    // between those (supposedly equal) values. we should track where the discrepancy comes from.
    final double stopTime = isRoot ? Double.MAX_VALUE : spTree.getNodeHeight(spTree.getParent(node));
    // demographic function is 0 based (relative to node height)
    // time away from node
    double lastTime = 0.0;
    // demographic function across branch
    DemographicFunction demog = spTree.getNodeDemographic(node);
    if (popFactor > 0) {
        demog = new ScaledDemographic(demog, popFactor);
    }
    //        if(false) {
    //            final double duration = isRoot ? cList[cList.length - 1].ctime - t0 : stopTime;
    //            double demographicAtCoalPoint = demog.getDemographic(duration);
    //            double intervalArea = demog.getIntegral(0, duration);
    //            if( demographicAtCoalPoint < 1e-12 * (duration/intervalArea) ) {
    //               return Double.NEGATIVE_INFINITY;
    //            }
    //        }
    // Species sharing this branch
    FixedBitSet subspeciesSet = spTree.spSet(node);
    if (verbose) {
        System.err.println(TreeUtils.uniqueNewick(spTree, node) + " nl " + nLineages + " " + subspeciesSet + " t0 - st " + t0 + " - " + stopTime);
    }
    while (nLineages > 1) {
        //            }
        assert (indexInClist < cList.length);
        final double nextT = cList[indexInClist].ctime;
        // while rare they can be equal
        if (nextT >= stopTime) {
            break;
        }
        if (nonEmptyIntersection(cList[indexInClist].sinfo, subspeciesSet)) {
            final double time = nextT - t0;
            if (time > 0) {
                final double interval = demog.getIntegral(lastTime, time);
                lastTime = time;
                final int nLineageOver2 = (nLineages * (nLineages - 1)) / 2;
                like -= nLineageOver2 * interval;
                final double pop = demog.getDemographic(time);
                assert (pop > 0);
                like -= Math.log(pop);
            }
            --nLineages;
        }
        ++indexInClist;
    }
    if (nLineages > 1) {
        // add term for No coalescent until root
        final double interval = demog.getIntegral(lastTime, stopTime - t0);
        final int nLineageOver2 = (nLineages * (nLineages - 1)) / 2;
        like -= nLineageOver2 * interval;
    }
    info[0] = nLineages;
    info[1] = indexInClist;
    if (verbose) {
        System.err.println(TreeUtils.uniqueNewick(spTree, node) + " stopTime " + stopTime + " nl " + nLineages + " icl " + indexInClist);
    }
    return like;
}