Examples with Conversion bacter.Conversion used on opensource projects

Example 16 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class COACGLogFileReader method iterator.

@Override
public Iterator<ConversionGraph> iterator() {
    try {
        reset();
        skipBurnin();
    } catch (XMLStreamException | IOException e) {
        throw new IllegalStateException(e.getMessage());
    }
    return new Iterator<ConversionGraph>() {

        int current = 0;

        @Override
        public boolean hasNext() {
            return current < nACGs;
        }

        @Override
        public ConversionGraph next() {
            current += 1;
            ConversionGraph acg = new ConversionGraph();
            for (Locus locus : loci) acg.lociInput.setValue(locus, acg);
            try {
                acg.initAndValidate();
            } catch (Exception e) {
                throw new IllegalStateException(e.getMessage());
            }
            String newick = null;
            List<Integer> rStarts = new ArrayList<>();
            List<Integer> rEnds = new ArrayList<>();
            List<Integer> rEFroms = new ArrayList<>();
            List<Integer> rETos = new ArrayList<>();
            List<Double> rAFroms = new ArrayList<>();
            List<Double> rATos = new ArrayList<>();
            try {
                skipUntil("iteration");
                while (xmlStreamReader.hasNext()) {
                    int type = xmlStreamReader.next();
                    if (type == XMLStreamConstants.END_ELEMENT && xmlStreamReader.getLocalName().toLowerCase().equals("iteration"))
                        break;
                    if (type != XMLStreamConstants.START_ELEMENT)
                        continue;
                    switch(xmlStreamReader.getLocalName().toLowerCase()) {
                        case "tree":
                            xmlStreamReader.next();
                            newick = xmlStreamReader.getText().trim();
                            break;
                        case "recedge":
                            while ((type = xmlStreamReader.next()) != XMLStreamConstants.END_ELEMENT || !xmlStreamReader.getLocalName().toLowerCase().equals("recedge")) {
                                if (type != XMLStreamConstants.START_ELEMENT)
                                    continue;
                                switch(xmlStreamReader.getLocalName().toLowerCase()) {
                                    case "start":
                                        xmlStreamReader.next();
                                        rStarts.add(Integer.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "end":
                                        xmlStreamReader.next();
                                        rEnds.add(Integer.valueOf(xmlStreamReader.getText().trim()) - 1);
                                        break;
                                    case "efrom":
                                        xmlStreamReader.next();
                                        rEFroms.add(Integer.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "eto":
                                        xmlStreamReader.next();
                                        rETos.add(Integer.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "afrom":
                                        xmlStreamReader.next();
                                        rAFroms.add(Double.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "ato":
                                        xmlStreamReader.next();
                                        rATos.add(Double.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    default:
                                }
                            }
                            break;
                        default:
                    }
                }
            } catch (XMLStreamException e) {
                throw new IllegalStateException(e.getMessage());
            }
            acg.fromExtendedNewick(newick, true, 0);
            for (int i = 0; i < rStarts.size(); i++) {
                Node fromNode = acg.getNode(rEFroms.get(i));
                Node toNode = acg.getNode(rETos.get(i));
                Conversion conv = new Conversion(toNode, toNode.getHeight() + rATos.get(i), fromNode, fromNode.getHeight() + rAFroms.get(i), rStarts.get(i), rEnds.get(i), acg, loci.get(0));
                acg.addConversion(conv);
            }
            current += 1;
            return acg;
        }
    };
}

Example 17 with Conversion

use of bacter.Conversion 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 18 with Conversion

use of bacter.Conversion 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 19 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ConvertedEdgeFlip method proposal.

@Override
public double proposal() {
    if (acg.getTotalConvCount() == 0)
        return Double.NEGATIVE_INFINITY;
    Conversion recomb = chooseConversion();
    Node node1 = recomb.getNode1();
    Node node2 = recomb.getNode2();
    double height1 = recomb.getHeight1();
    double height2 = recomb.getHeight2();
    if (node1 == node2 || node2.isRoot())
        return Double.NEGATIVE_INFINITY;
    if (height1 < node2.getHeight() || height1 > node2.getParent().getHeight() || height2 < node1.getHeight() || height2 > node1.getParent().getHeight())
        return Double.NEGATIVE_INFINITY;
    recomb.setNode1(node2);
    recomb.setNode2(node1);
    return 0.0;
}

Example 20 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ConvertedEdgeHop method proposal.

@Override
public double proposal() {
    if (acg.getTotalConvCount() == 0)
        return Double.NEGATIVE_INFINITY;
    // Select recombination at random
    Conversion recomb = chooseConversion();
    // Choose whether to move departure or arrival point
    boolean moveDeparture;
    moveDeparture = recomb.getNode2().isRoot() || Randomizer.nextBoolean();
    // Select new attachment point:
    double u = Randomizer.nextDouble() * acg.getClonalFrameLength();
    Node nodeBelow = null;
    double newHeight = -1;
    for (Node node : acg.getNodesAsArray()) {
        if (node.isRoot())
            continue;
        if (u < node.getLength()) {
            newHeight = node.getHeight() + u;
            nodeBelow = node;
            break;
        }
        u -= node.getLength();
    }
    if (newHeight < 0.0)
        throw new IllegalStateException("Problem with recombinant edge " + "hop operator!  This is a bug.");
    // Check that new height does not lie out of bounds
    if (moveDeparture) {
        if (newHeight > recomb.getHeight2())
            return Double.NEGATIVE_INFINITY;
        else {
            recomb.setHeight1(newHeight);
            recomb.setNode1(nodeBelow);
        }
    } else {
        if (newHeight < recomb.getHeight1())
            return Double.NEGATIVE_INFINITY;
        else {
            recomb.setHeight2(newHeight);
            recomb.setNode2(nodeBelow);
        }
    }
    return 0.0;
}