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

Example 1 with FixedBitSet

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

the class AlloppSpeciesBindings method taxonseqToTipUnion.

// Taxons vs species.
// Taxons may have a final "0", "1",... to distinguish sequences, while
// species do not. AlloppLeggedTree uses a SimpleTree, which only has
// Taxons, so same thing there. Multree needs distinguishable Taxons
// so has suffices.
public FixedBitSet taxonseqToTipUnion(Taxon tx, int seq) {
    FixedBitSet union = speciesseqEmptyUnion();
    int sp = apspeciesId2index(tx.getId());
    int spseq = spandseq2spseqindex(sp, seq);
    union.set(spseq);
    return union;
}

Example 2 with FixedBitSet

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

the class AlloppSpeciesNetworkModel method calculateDipHistTipUnion.

// finds the union of a tip (diploid or hyb-tip) in the diploid history.
// Used by move that slides node in diploid history to check gene-tree compatibility.
public FixedBitSet calculateDipHistTipUnion(NodeRef node) {
    if (node == null) {
        System.out.println("BUG in calculateDipHistTipUnion()");
    }
    assert node != null;
    int tt = adhist.getNodeTettree(node.getNumber());
    AlloppDiploidHistory.LegLorR leg = adhist.getNodeLeg(node.getNumber());
    int seq = (leg == AlloppDiploidHistory.LegLorR.left) ? 0 : 1;
    FixedBitSet union;
    if (tt < 0) {
        // ordinary tip
        union = apsp.taxonseqToTipUnion(adhist.getSlidableNodeTaxon(node), 0);
    } else {
        union = unionOfWholeTetTree(tt, seq);
    }
    return union;
}

Example 3 with FixedBitSet

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

the class MulMSCoalescent method treeLogLikelihood.

private double treeLogLikelihood(MulSpeciesBindings.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 MulSpeciesBindings.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) {
            System.err.println("ERROR in evomodel.speciation.MulMSCoalescent.treeLogLikelihood()");
        }
        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;
}

Example 4 with FixedBitSet

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

the class MulSpeciesBindings method collectCoalInfo.

// jh + grj
/**
     * Collect coalescence information for sub-tree rooted at 'node'.
     *
     * @param tree
     * @param node
     * @param loc  Place node data in loc, sub-tree info before that.
     * @param info array to fill
     * @return location of next available location
     */
private int collectCoalInfo(Tree tree, NodeRef node, GeneTreeInfo.SequenceAssignment[] seqassigns, int loc, CoalInfo[] info) {
    info[loc] = new CoalInfo(tree.getNodeHeight(node), tree.getChildCount(node));
    int newLoc = loc - 1;
    for (int i = 0; i < 2; i++) {
        NodeRef child = tree.getChild(node, i);
        info[loc].sinfo[i] = new FixedBitSet(nSpSeqs());
        if (tree.isExternal(child)) {
            Taxon taxon = tree.getNodeTaxon(child);
            int ti = taxon2index.get(taxon);
            int spseq = spsq[seqassigns[ti].spIndex][seqassigns[ti].seqIndex];
            info[loc].sinfo[i].set(spseq);
            assert tree.getNodeHeight(child) == 0;
        } else {
            final int used = collectCoalInfo(tree, child, seqassigns, newLoc, info);
            for (int j = 0; j < info[newLoc].sinfo.length; ++j) {
                info[loc].sinfo[i].union(info[newLoc].sinfo[j]);
            }
            newLoc = used;
        }
    }
    return newLoc;
}

Example 5 with FixedBitSet

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

the class MulSpeciesTreeModel method setSPsets.

private NodeProperties setSPsets(NodeRef nodeID) {
    final NodeProperties nprop = props.get(nodeID);
    if (!isExternal(nodeID)) {
        nprop.spSet = new FixedBitSet(mulspb.nSpSeqs());
        for (int nc = 0; nc < getChildCount(nodeID); ++nc) {
            NodeProperties p = setSPsets(getChild(nodeID, nc));
            nprop.spSet.union(p.spSet);
        }
    }
    return nprop;
}