Examples with Locus bacter.Locus used on opensource projects

Example 11 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class CFOperator method collapseConversions.

/**
 * Take conversions which connect to edge above srcNode at times greater than
 * destTime and attach them instead to the lineage above destNode.
 *
 * Assumes topology has not yet been altered.
 *
 * @param srcNode   source node for move
 * @param destNode  dest node for move
 * @param destTime  new time of attachment of edge above srcNode to edge
 *                  above destNode
 * @return log probability of the collapsed attachments.
 */
protected double collapseConversions(Node srcNode, Node destNode, double destTime) {
    double logP = 0.0;
    boolean reverseRootMove = srcNode.getParent().isRoot();
    Node srcNodeP = srcNode.getParent();
    Node srcNodeS = getSibling(srcNode);
    double maxChildHeight = getMaxRootChildHeight();
    double volatileHeight = Math.max(maxChildHeight, destTime);
    // Collapse non-root conversions
    Node node = destNode;
    while (!node.isRoot() && node.getHeight() < srcNodeP.getHeight()) {
        double lowerBound = Math.max(destTime, node.getHeight());
        double upperBound = Math.min(node.getParent().getHeight(), srcNodeP.getHeight());
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if (conv.getHeight1() > lowerBound && conv.getHeight1() < upperBound) {
                    if (conv.getNode1() == srcNode)
                        conv.setNode1(node);
                    if (conv.getNode1() == node && (!reverseRootMove || conv.getHeight1() < volatileHeight))
                        logP += Math.log(0.5);
                }
                if (conv.getHeight2() > lowerBound && conv.getHeight2() < upperBound) {
                    if (conv.getNode2() == srcNode)
                        conv.setNode2(node);
                    if (conv.getNode2() == node && (!reverseRootMove || conv.getNode1() != node || conv.getHeight1() < volatileHeight))
                        logP += Math.log(0.5);
                }
            }
        }
        node = node.getParent();
    }
    if (reverseRootMove) {
        double L = 2.0 * (srcNode.getParent().getHeight() - volatileHeight);
        double Nexp = L * rhoInput.get().getValue() * (acg.getTotalSequenceLength() + acg.getLoci().size() * (deltaInput.get().getValue() - 1.0));
        List<Conversion> toRemove = new ArrayList<>();
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if (conv.getHeight1() > volatileHeight)
                    toRemove.add(conv);
            }
        }
        logP += -Nexp + toRemove.size() * Math.log(Nexp);
        for (Conversion conv : toRemove) {
            logP += Math.log(1.0 / L) + getAffectedRegionProb(conv) + getEdgeCoalescenceProb(conv);
            acg.deleteConversion(conv);
        }
    }
    // Apply topology modifications.
    disconnectEdge(srcNode);
    connectEdge(srcNode, destNode, destTime);
    if (reverseRootMove && destTime < maxChildHeight)
        acg.setRoot(srcNodeS);
    return logP;
}

Example 12 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class CFUniform method proposal.

@Override
public double proposal() {
    double logHGF = 0.0;
    double logHalf = Math.log(0.5);
    // Select internal non-root node at random.
    Node node = acg.getNode(acg.getLeafNodeCount() + Randomizer.nextInt(acg.getInternalNodeCount()));
    Node leftChild = node.getLeft();
    Node rightChild = node.getRight();
    double oldHeight = node.getHeight();
    double maxChildHeight = Math.max(leftChild.getHeight(), rightChild.getHeight());
    // Choose new height:
    double newHeight;
    if (node.isRoot()) {
        double fMin = Math.min(scaleFactorInput.get(), 1.0 / scaleFactorInput.get());
        double fMax = 1.0 / fMin;
        double f = fMin + (fMax - fMin) * Randomizer.nextDouble();
        newHeight = node.getHeight() * f;
        logHGF += Math.log(1.0 / f);
        if (newHeight < maxChildHeight)
            return Double.NEGATIVE_INFINITY;
    } else {
        Node parent = node.getParent();
        newHeight = maxChildHeight + Randomizer.nextDouble() * (parent.getHeight() - maxChildHeight);
    }
    if (newHeight > oldHeight) {
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if (conv.getNode1() == node && conv.getHeight1() < newHeight) {
                    conv.setNode1(Randomizer.nextBoolean() ? leftChild : rightChild);
                    logHGF -= logHalf;
                }
                if (conv.getNode2() == node && conv.getHeight2() < newHeight) {
                    conv.setNode2(Randomizer.nextBoolean() ? leftChild : rightChild);
                    logHGF -= logHalf;
                }
            }
        }
        node.setHeight(newHeight);
        if (node.isRoot()) {
            // Draw a number of conversions
            double L = 2.0 * (newHeight - oldHeight);
            double Nexp = L * rhoInput.get().getValue() * (acg.getTotalSequenceLength() + acg.getLoci().size() * deltaInput.get().getValue());
            int N = (int) Randomizer.nextPoisson(Nexp);
            // N! cancels
            logHGF -= -Nexp + N * Math.log(Nexp);
            for (int i = 0; i < N; i++) {
                Conversion conv = new Conversion();
                double u = L * Randomizer.nextDouble();
                if (u < 0.5 * L) {
                    conv.setNode1(leftChild);
                    conv.setHeight1(oldHeight + u);
                } else {
                    conv.setNode1(rightChild);
                    conv.setHeight1(oldHeight + u - 0.5 * L);
                }
                logHGF -= Math.log(1.0 / L) + coalesceEdge(conv) + drawAffectedRegion(conv);
                acg.addConversion(conv);
            }
        }
    } else {
        List<Conversion> toRemove = new ArrayList<>();
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if ((conv.getNode1() == leftChild || conv.getNode1() == rightChild) && conv.getHeight1() > newHeight) {
                    if (node.isRoot()) {
                        toRemove.add(conv);
                        continue;
                    } else {
                        conv.setNode1(node);
                        logHGF += logHalf;
                    }
                }
                if ((conv.getNode2() == leftChild || conv.getNode2() == rightChild) && conv.getHeight2() > newHeight) {
                    conv.setNode2(node);
                    logHGF += logHalf;
                }
            }
        }
        if (node.isRoot()) {
            double L = 2.0 * (oldHeight - newHeight);
            double Nexp = L * rhoInput.get().getValue() * (acg.getTotalSequenceLength() + acg.getLoci().size() * deltaInput.get().getValue());
            // N! cancels
            logHGF += -Nexp + toRemove.size() * Math.log(Nexp);
            for (Conversion conv : toRemove) {
                logHGF += Math.log(1.0 / L) + getEdgeCoalescenceProb(conv) + getAffectedRegionProb(conv);
                acg.deleteConversion(conv);
            }
        }
        node.setHeight(newHeight);
    }
    if (acg.isInvalid())
        throw new IllegalStateException("Something is broken: CFUniform proposed illegal state!");
    return logHGF;
}

Example 13 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ConversionCreationOperator method drawAffectedRegionUnrestricted.

/**
 * Choose region to be affected by this conversion using the
 * full ClonalOrigin model.
 *
 * @param conv Conversion object where these sites are stored.
 * @return log probability density of chosen attachment.
 */
protected double drawAffectedRegionUnrestricted(Conversion conv) {
    double logP = 0.0;
    // Total effective number of possible start sites
    double alpha = acg.getTotalSequenceLength() + acg.getLoci().size() * (deltaInput.get().getValue() - 1.0);
    // Draw location of converted region.
    int startSite = -1;
    int endSite;
    Locus locus = null;
    double u = Randomizer.nextDouble() * alpha;
    for (Locus thisLocus : acg.getLoci()) {
        if (u < deltaInput.get().getValue() - 1.0 + thisLocus.getSiteCount()) {
            locus = thisLocus;
            if (u < deltaInput.get().getValue()) {
                startSite = 0;
                logP += Math.log(deltaInput.get().getValue() / alpha);
            } else {
                startSite = (int) Math.ceil(u - deltaInput.get().getValue());
                logP += Math.log(1.0 / alpha);
            }
            break;
        }
        u -= deltaInput.get().getValue() - 1.0 + thisLocus.getSiteCount();
    }
    if (locus == null)
        throw new IllegalStateException("Programmer error: " + "loop in drawAffectedRegion() fell through.");
    endSite = startSite + (int) Randomizer.nextGeometric(1.0 / deltaInput.get().getValue());
    endSite = Math.min(endSite, locus.getSiteCount() - 1);
    // Probability of end site:
    if (endSite == locus.getSiteCount() - 1) {
        logP += (locus.getSiteCount() - 1 - startSite) * Math.log(1.0 - 1.0 / deltaInput.get().getValue());
    } else {
        logP += (endSite - startSite) * Math.log(1.0 - 1.0 / deltaInput.get().getValue()) - Math.log(deltaInput.get().getValue());
    }
    conv.setLocus(locus);
    conv.setStartSite(startSite);
    conv.setEndSite(endSite);
    return logP;
}

Example 14 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class CFConvSwapExperiment method connectEdge.

public static void connectEdge(ConversionGraph acg, Node node, Node destEdgeBase, double destTime) {
    if (node.isRoot())
        throw new IllegalArgumentException("Programmer error: " + "root argument passed to connectEdge().");
    Node parent = node.getParent();
    if (destEdgeBase.isRoot()) {
        parent.addChild(destEdgeBase);
    } else {
        Node grandParent = destEdgeBase.getParent();
        grandParent.removeChild(destEdgeBase);
        grandParent.addChild(parent);
        parent.addChild(destEdgeBase);
    }
    parent.setHeight(destTime);
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (conv.getNode1() == destEdgeBase && conv.getHeight1() > destTime)
                conv.setNode1(parent);
            if (conv.getNode2() == destEdgeBase && conv.getHeight2() > destTime)
                conv.setNode2(parent);
        }
    }
}

Example 15 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ConvertedRegionLogger method log.

@Override
public void log(long nSample, PrintStream out) {
    for (Locus locus : acgInput.get().getLoci()) {
        if (acgInput.get().getConvCount(locus) == 0) {
            out.print("NA\t");
            return;
        }
        for (int r = 0; r < acgInput.get().getConvCount(locus); r++) {
            if (r > 0)
                out.print(",");
            Conversion recomb = acgInput.get().getConversions(locus).get(r);
            out.print(recomb.getStartSite() + ":" + recomb.getEndSite());
        }
        out.print("\t");
    }
}