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

Example 6 with FixedBitSet

use of jebl.util.FixedBitSet 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 7 with FixedBitSet

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

the class MulTreeNodeSlide method operateOneNode.

public void operateOneNode(final double factor) {
    //            #print "operate: tree", ut.treerep(t)
    //   if( verbose)  System.out.println("  Mau at start: " + tree.getSimpleTree());
    final int count = multree.getExternalNodeCount();
    assert count == species.nSpSeqs();
    NodeRef[] order = new NodeRef[2 * count - 1];
    boolean[] swapped = new boolean[count - 1];
    mauCanonical(multree, order, swapped);
    // internal node to change
    // count-1 - number of internal nodes
    int which = MathUtils.nextInt(count - 1);
    FixedBitSet left = new FixedBitSet(count);
    FixedBitSet right = new FixedBitSet(count);
    for (int k = 0; k < 2 * which + 1; k += 2) {
        left.set(multree.speciesIndex(order[k]));
    }
    for (int k = 2 * (which + 1); k < 2 * count; k += 2) {
        right.set(multree.speciesIndex(order[k]));
    }
    double newHeight;
    if (factor > 0) {
        newHeight = multree.getNodeHeight(order[2 * which + 1]) * factor;
    } else {
        final double limit = species.speciationUpperBound(left, right);
        newHeight = MathUtils.nextDouble() * limit;
    }
    multree.beginTreeEdit();
    multree.setPreorderIndices(preOrderIndexBefore);
    final NodeRef node = order[2 * which + 1];
    multree.setNodeHeight(node, newHeight);
    mauReconstruct(multree, order, swapped);
    // restore pre-order of pops -
    {
        multree.setPreorderIndices(preOrderIndexAfter);
        double[] splitPopValues = null;
        for (int k = 0; k < preOrderIndexBefore.length; ++k) {
            final int b = preOrderIndexBefore[k];
            if (b >= 0) {
                final int a = preOrderIndexAfter[k];
                if (a != b) {
                    //if( verbose)  System.out.println("pops: " + a + " <- " + b);
                    final Parameter p1 = multree.sppSplitPopulations;
                    if (splitPopValues == null) {
                        splitPopValues = p1.getParameterValues();
                    }
                    if (multree.constPopulation()) {
                        p1.setParameterValue(count + a, splitPopValues[count + b]);
                    } else {
                        for (int i = 0; i < 2; ++i) {
                            p1.setParameterValue(count + 2 * a + i, splitPopValues[count + 2 * b + i]);
                        }
                    }
                }
            }
        }
    }
    multree.endTreeEdit();
}

Example 8 with FixedBitSet

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

the class SpeciesBindings method collectCoalInfo.

/**
     * 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, 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(nSpecies());
        if (tree.isExternal(child)) {
            info[loc].sinfo[i].set(taxon2Species.get(tree.getNodeTaxon(child)));
            assert tree.getNodeHeight(child) == 0;
        } else {
            final int used = collectCoalInfo(tree, child, 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 9 with FixedBitSet

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

the class SpeciesTreeModel method setSPsets.

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

Example 10 with FixedBitSet

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

the class TreeNodeSlide method operateOneNode.

public void operateOneNode(final double factor) {
    //            #print "operate: tree", ut.treerep(t)
    //   if( verbose)  System.out.println("  Mau at start: " + tree.getSimpleTree());
    final int count = tree.getExternalNodeCount();
    assert count == species.nSpecies();
    NodeRef[] order = new NodeRef[2 * count - 1];
    boolean[] swapped = new boolean[count - 1];
    mauCanonical(tree, order, swapped);
    // internal node to change
    // count-1 - number of internal nodes
    int which = MathUtils.nextInt(count - 1);
    //        if( outgroupOnly ) {
    //           if( order[1] == tree.getRoot() ) {
    //               which = 0;
    //           } else if( order[2*count - 3] == tree.getRoot() ) {
    //               which = count - 2;
    //           }
    //        }
    //        if( verbose)  {
    //            System.out.print("order:");
    //            for(int k = 0 ; k < 2*count; k += 2) {
    //                System.out.print(tree.getNodeTaxon(order[k]).getId() + ((k == 2*which) ? "*" : " "));
    //            }
    //            System.out.println();
    //        }
    FixedBitSet left = new FixedBitSet(count);
    FixedBitSet right = new FixedBitSet(count);
    for (int k = 0; k < 2 * which + 1; k += 2) {
        left.set(tree.speciesIndex(order[k]));
    }
    for (int k = 2 * (which + 1); k < 2 * count; k += 2) {
        right.set(tree.speciesIndex(order[k]));
    }
    double newHeight;
    if (factor > 0) {
        newHeight = tree.getNodeHeight(order[2 * which + 1]) * factor;
    } else {
        final double limit = species.speciationUpperBound(left, right);
        newHeight = MathUtils.nextDouble() * limit;
    }
    if (false) {
        double totChange = 0;
        double tot = 0;
        double totChangeUp = 0.0;
        int topChange = 0;
        for (int which1 = 0; which1 < count - 1; which1++) {
            FixedBitSet left1 = new FixedBitSet(count);
            FixedBitSet right1 = new FixedBitSet(count);
            final NodeRef node = order[2 * which1 + 1];
            final double h = tree.getNodeHeight(node);
            for (int k = 0; k < 2 * which1 + 1; k += 2) {
                final NodeRef ref = order[k];
                left1.set(tree.speciesIndex(ref));
            }
            for (int k = 2 * (which1 + 1); k < 2 * count; k += 2) {
                right1.set(tree.speciesIndex(order[k]));
            }
            final double limit = species.speciationUpperBound(left1, right1);
            final NodeRef parent = tree.getParent(node);
            final double ph = parent != null ? tree.getNodeHeight(parent) : limit;
            final double h1 = tree.getNodeHeight(tree.getChild(node, 0));
            final double h2 = tree.getNodeHeight(tree.getChild(node, 1));
            final double chDown = Math.max(h1, h2);
            final double chup = Math.max(limit - ph, 0);
            totChangeUp += chup;
            double change = chDown + chup;
            totChange += change;
            if (change > limit) {
                assert change <= limit;
            }
            tot += limit;
            if (which == which1 && (newHeight < chDown || newHeight > limit - chup)) {
                topChange = 1;
            }
        //               final double h1 = tree.getNodeHeight(order[2 * which1]);
        //               final double h2 = tree.getNodeHeight(order[2 * which1 + 2]);
        //               double lowerLimit = 0;
        //               double upperLimit = limit;
        //               if( h > h1 ) {
        //                   assert tree.getParent(order[2 * which]) == node;
        //                   lowerLimit = h1;
        //               } else {
        //                   assert tree.getParent() == node;
        //                   upperLimit = Math.min(upperLimit, h1);
        //               }
        //               if( h > h2 && h2 > lowerLimit ) {
        //                   assert tree.getParent(order[2 * which+2]) == node;
        //                   lowerLimit = h2;
        //               }
        //
        }
        final double p = totChange / tot;
        System.err.println("top% " + p + " " + totChangeUp / tot + (topChange > 0 ? "++" : ""));
    }
    tree.beginTreeEdit();
    tree.setPreorderIndices(preOrderIndexBefore);
    final NodeRef node = order[2 * which + 1];
    if (false) {
        System.err.println("Before " + node.getNumber() + " " + newHeight);
        ptree(tree, order, preOrderIndexBefore);
        if (newHeight > tree.getNodeHeight(node)) {
            NodeRef parent = tree.getParent(node);
            while (parent != null && tree.getNodeHeight(parent) < newHeight) {
                // swap
                final int pi = getIndexOf(parent, order);
                final int side = pi < which ? 0 : 1;
                // 0 for right, 1 for left
                //final int side = (1+getSide(which, parent, order, swapped))/2;
                swapPops(node, swapped[which] ? 1 - side : side, parent, swapped[pi] ? side : 1 - side, preOrderIndexBefore);
                parent = tree.getParent(parent);
            }
        } else {
            // NodeRef child = tree.getChild(node, 0);
            if (which > 0) {
                NodeRef nb = order[2 * which - 1];
                if (tree.getNodeHeight(node) > tree.getNodeHeight(nb)) {
                    while (nb != node && tree.getNodeHeight(nb) < newHeight) {
                        nb = tree.getParent(nb);
                    }
                    if (nb != node) {
                        final int side = swapped[which] ? 1 : 0;
                        swapPops(node, side, nb, 1 - side, preOrderIndexBefore);
                        NodeRef pnb = tree.getParent(nb);
                        while (pnb != node) {
                            swapPops(nb, 1 - side, pnb, 1 - side, preOrderIndexBefore);
                            nb = pnb;
                            pnb = tree.getParent(nb);
                        }
                    }
                }
            }
            if (which < count - 2) {
                NodeRef nb = order[2 * which + 3];
                if (tree.getNodeHeight(node) > tree.getNodeHeight(nb)) {
                    while (nb != node && tree.getNodeHeight(nb) < newHeight) {
                        nb = tree.getParent(nb);
                    }
                    if (nb != node) {
                        final int side = swapped[which] ? 0 : 1;
                        swapPops(node, side, nb, 1 - side, preOrderIndexBefore);
                        NodeRef pnb = tree.getParent(nb);
                        while (pnb != node) {
                            swapPops(nb, 1 - side, pnb, 1 - side, preOrderIndexBefore);
                            nb = pnb;
                            pnb = tree.getParent(nb);
                        }
                    }
                }
            }
        }
    }
    tree.setNodeHeight(node, newHeight);
    mauReconstruct(tree, order, swapped);
    // restore pre-order of pops -
    {
        tree.setPreorderIndices(preOrderIndexAfter);
        double[] splitPopValues = null;
        for (int k = 0; k < preOrderIndexBefore.length; ++k) {
            final int b = preOrderIndexBefore[k];
            if (b >= 0) {
                final int a = preOrderIndexAfter[k];
                if (a != b) {
                    //if( verbose)  System.out.println("pops: " + a + " <- " + b);
                    final Parameter p1 = tree.sppSplitPopulations;
                    if (splitPopValues == null) {
                        splitPopValues = p1.getParameterValues();
                    }
                    if (tree.constPopulation()) {
                        p1.setParameterValue(count + a, splitPopValues[count + b]);
                    } else {
                        for (int i = 0; i < 2; ++i) {
                            p1.setParameterValue(count + 2 * a + i, splitPopValues[count + 2 * b + i]);
                        }
                    }
                }
            }
        }
    }
    // System.err.println("After");
    //ptree(tree, order, preOrderIndexAfter);
    // if( verbose)  System.out.println("  Mau after: " + tree.getSimpleTree());
    // }
    tree.endTreeEdit();
}