Examples with Node beast.evolution.tree.Node used on opensource projects

Example 51 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ACGLikelihoodBeagle method setStates.

/**
 * Set leaf states in a Beagle instance
 *
 * @param beagle beagle instance object
 * @param patterns leaf state patterns
 */
void setStates(Beagle beagle, Multiset<int[]> patterns) {
    for (Node node : acg.getExternalNodes()) {
        int[] states = new int[patterns.size()];
        int taxon = alignment.getTaxonIndex(node.getID());
        int i = 0;
        for (int[] pattern : patterns.elementSet()) {
            // int code = pattern[taxon];
            // int[] statesForCode = alignment.getDataType().getStatesForCode(code);
            // if (statesForCode.length==1)
            // states[i] = statesForCode[0];
            // else
            // states[i] = code; // Causes ambiguous states to be ignored.
            states[i] = pattern[taxon];
            i += 1;
        }
        beagle.setTipStates(node.getNr(), states);
    }
}

Example 52 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ConvertedRegionSwap method proposal.

@Override
public double proposal() {
    if (acg.getTotalConvCount() < 2)
        return Double.NEGATIVE_INFINITY;
    // Select a random pair of recombinations
    Conversion conv1 = chooseConversion();
    Conversion conv2;
    do {
        conv2 = chooseConversion();
    } while (conv2 == conv1);
    // Switch edges corresponding to recombinations.  (Can't switch
    // loci, as this would break ordering constraint.)
    double tmpHeight1 = conv1.getHeight1();
    double tmpHeight2 = conv1.getHeight2();
    Node tmpNode1 = conv1.getNode1();
    Node tmpNode2 = conv1.getNode2();
    conv1.setHeight1(conv2.getHeight1());
    conv1.setHeight2(conv2.getHeight2());
    conv1.setNode1(conv2.getNode1());
    conv1.setNode2(conv2.getNode2());
    conv2.setHeight1(tmpHeight1);
    conv2.setHeight2(tmpHeight2);
    conv2.setNode1(tmpNode1);
    conv2.setNode2(tmpNode2);
    return 0.0;
}

Example 53 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class EdgeCreationOperator method coalesceEdge.

/**
 * Take a recombination with an existing departure point and determine
 * the arrival point by allowing it to coalesce with the clonal frame.
 *
 * @param conv recombination to modify
 * @return log probability density of coalescent point chosen.
 */
public double coalesceEdge(Conversion conv) {
    double logP = 0.0;
    List<CFEventList.Event> events = acg.getCFEvents();
    // Locate event immediately below departure point
    int startIdx = 0;
    while (events.get(startIdx + 1).getHeight() < conv.getHeight1()) startIdx += 1;
    // Choose edge length in dimensionless time.
    double u = Randomizer.nextExponential(1.0);
    // Determine arrival point in real time
    for (int i = startIdx; i < events.size(); i++) {
        CFEventList.Event event = events.get(i);
        double t = Math.max(conv.getHeight1(), event.getHeight());
        // Determine length of interval in dimensionless time
        double intervalArea;
        if (i < events.size() - 1)
            intervalArea = popFunc.getIntegral(t, events.get(i + 1).getHeight());
        else
            intervalArea = Double.POSITIVE_INFINITY;
        // Check whether arrival falls within this interval
        if (u < intervalArea * event.getLineageCount()) {
            // Set arrival point in real time
            conv.setHeight2(popFunc.getInverseIntensity(popFunc.getIntensity(t) + u / event.getLineageCount()));
            // Attach to random clonal frame lineage extant at this time
            int z = Randomizer.nextInt(event.getLineageCount());
            for (Node node : acg.getNodesAsArray()) {
                if (conv.getHeight2() > node.getHeight() && (node.isRoot() || conv.getHeight2() < node.getParent().getHeight())) {
                    if (z == 0) {
                        conv.setNode2(node);
                        break;
                    } else
                        z -= 1;
                }
            }
            logP += -u + Math.log(1.0 / popFunc.getPopSize(conv.getHeight2()));
            break;
        } else {
            u -= intervalArea * event.getLineageCount();
            logP += -intervalArea * event.getLineageCount();
        }
    }
    return logP;
}

Example 54 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class EdgeCreationOperator method attachEdge.

/**
 * Attach chosen recombination to the clonal frame.  Note that only the
 * attachment points (nodes and heights) are set, the affected region of
 * the alignment is not modified.
 *
 * @param conv conversion
 * @return log probability density of chosen attachment.
 */
public double attachEdge(Conversion conv) {
    double logP = 0.0;
    // Select departure point
    double u = Randomizer.nextDouble() * acg.getClonalFrameLength();
    logP += Math.log(1.0 / acg.getClonalFrameLength());
    for (Node node : acg.getNodesAsArray()) {
        if (node.isRoot())
            continue;
        if (u < node.getLength()) {
            conv.setNode1(node);
            conv.setHeight1(node.getHeight() + u);
            break;
        } else
            u -= node.getLength();
    }
    // Select arrival point
    logP += coalesceEdge(conv);
    return logP;
}

Example 55 with Node

use of beast.evolution.tree.Node 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;
        }
    };
}