Examples with Conversion bacter.Conversion used on opensource projects

Example 36 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class AddRemoveDetour method addDetour.

/**
 * Detour addition variant of move.
 *
 * @return log HGF
 */
private double addDetour() {
    double logHGF = 0.0;
    if (acg.getTotalConvCount() == 0)
        return Double.NEGATIVE_INFINITY;
    // Select conversion at random
    Conversion conv = chooseConversion();
    logHGF -= Math.log(1.0 / acg.getTotalConvCount());
    // Select detour times:
    double t1 = conv.getHeight1() + Randomizer.nextDouble() * (conv.getHeight2() - conv.getHeight1());
    double t2 = conv.getHeight1() + Randomizer.nextDouble() * (conv.getHeight2() - conv.getHeight1());
    double tLower, tUpper;
    if (t1 < t2) {
        tLower = t1;
        tUpper = t2;
    } else {
        tLower = t2;
        tUpper = t1;
    }
    logHGF -= Math.log(2.0) - 2.0 * Math.log(conv.getHeight2() - conv.getHeight1());
    // Select non-root node below detour edge at random
    Node detour = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
    logHGF -= Math.log(1.0 / (acg.getNodeCount() - 1));
    // Abort if selected detour edge does not contain tLower and tUpper
    if (detour.getHeight() > tLower || detour.getParent().getHeight() < tUpper)
        return Double.NEGATIVE_INFINITY;
    // Abort if conv is already attached to selected detour edge:
    if (detour == conv.getNode1() || detour == conv.getNode2())
        return Double.NEGATIVE_INFINITY;
    Conversion convA = new Conversion();
    convA.setLocus(conv.getLocus());
    convA.setNode1(conv.getNode1());
    convA.setHeight1(conv.getHeight1());
    convA.setNode2(detour);
    convA.setHeight2(tLower);
    convA.setStartSite(conv.getStartSite());
    convA.setEndSite(conv.getEndSite());
    Conversion convB = new Conversion();
    convB.setNode1(detour);
    convB.setHeight1(tUpper);
    convB.setNode2(conv.getNode2());
    convB.setHeight2(conv.getHeight2());
    logHGF -= drawAffectedRegion(convB);
    acg.deleteConversion(conv);
    acg.addConversion(convA);
    acg.addConversion(convB);
    // Count number of node1s and node2s attached to detour edge
    int node1Count = 0;
    int node2Count = 0;
    for (Locus locus : acg.getLoci()) {
        for (Conversion thisConv : acg.getConversions(locus)) {
            if (thisConv.getNode1() == detour && thisConv.getNode2() != detour)
                node1Count += 1;
            if (thisConv.getNode2() == detour && thisConv.getNode1() != detour)
                node2Count += 1;
        }
    }
    // Incorporate probability of reverse move:
    logHGF += Math.log(1.0 / ((acg.getNodeCount() - 1) * node1Count * node2Count));
    return logHGF;
}

Example 37 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ACGAnnotator method summarizeConversions.

/**
 * Add summarized conversions to given ACG.
 *
 * @param cladeSystem information summarizing ACG posterior
 * @param acg conversion graph
 * @param threshold significance threshold
 * @param summaryStrategy strategy used when summarizing event ages/heights
 */
protected void summarizeConversions(ACGCladeSystem cladeSystem, ConversionGraph acg, int nACGs, double threshold, SummaryStrategy summaryStrategy) {
    BitSet[] bitSets = cladeSystem.getBitSets(acg);
    for (int fromNr = 0; fromNr < acg.getNodeCount(); fromNr++) {
        BitSet from = bitSets[fromNr];
        for (int toNr = 0; toNr < acg.getNodeCount(); toNr++) {
            BitSet to = bitSets[toNr];
            for (Locus locus : acg.getLoci()) {
                List<ACGCladeSystem.ConversionSummary> conversionSummaries = cladeSystem.getConversionSummaries(from, to, locus, nACGs, threshold);
                for (ACGCladeSystem.ConversionSummary conversionSummary : conversionSummaries) {
                    Conversion conv = new Conversion();
                    conv.setLocus(locus);
                    conv.setNode1(acg.getNode(fromNr));
                    conv.setNode2(acg.getNode(toNr));
                    double posteriorSupport = conversionSummary.nIncludedACGs / (double) nACGs;
                    double[] height1s = new double[conversionSummary.summarizedConvCount()];
                    double[] height2s = new double[conversionSummary.summarizedConvCount()];
                    double[] startSites = new double[conversionSummary.summarizedConvCount()];
                    double[] endSites = new double[conversionSummary.summarizedConvCount()];
                    for (int i = 0; i < conversionSummary.summarizedConvCount(); i++) {
                        height1s[i] = conversionSummary.height1s.get(i);
                        height2s[i] = conversionSummary.height2s.get(i);
                        startSites[i] = conversionSummary.startSites.get(i);
                        endSites[i] = conversionSummary.ends.get(i);
                    }
                    if (summaryStrategy == SummaryStrategy.MEAN) {
                        conv.setHeight1(DiscreteStatistics.mean(height1s));
                        conv.setHeight2(DiscreteStatistics.mean(height2s));
                        conv.setStartSite((int) Math.round(DiscreteStatistics.mean(startSites)));
                        conv.setEndSite((int) Math.round(DiscreteStatistics.mean(endSites)));
                    } else {
                        conv.setHeight1(DiscreteStatistics.median(height1s));
                        conv.setHeight2(DiscreteStatistics.median(height2s));
                        conv.setStartSite((int) Math.round(DiscreteStatistics.median(startSites)));
                        conv.setEndSite((int) Math.round(DiscreteStatistics.median(endSites)));
                    }
                    Arrays.sort(height1s);
                    double minHeight1HPD = height1s[(int) (0.025 * height1s.length)];
                    double maxHeight1HPD = height1s[(int) (0.975 * height1s.length)];
                    Arrays.sort(height2s);
                    double minHeight2HPD = height2s[(int) (0.025 * height2s.length)];
                    double maxHeight2HPD = height2s[(int) (0.975 * height2s.length)];
                    Arrays.sort(startSites);
                    int minStartHPD = (int) startSites[(int) (0.025 * startSites.length)];
                    int maxStartHPD = (int) startSites[(int) (0.975 * startSites.length)];
                    Arrays.sort(endSites);
                    int minEndHPD = (int) endSites[(int) (0.025 * endSites.length)];
                    int maxEndHPD = (int) endSites[(int) (0.975 * endSites.length)];
                    conv.newickMetaDataBottom = "height_95%_HPD={" + minHeight1HPD + "," + maxHeight1HPD + "}";
                    conv.newickMetaDataMiddle = "posterior=" + String.format("%2.2f", posteriorSupport) + ", startSite_95%_HPD={" + minStartHPD + "," + maxStartHPD + "}" + ", endSite_95%_HPD={" + minEndHPD + "," + maxEndHPD + "}";
                    conv.newickMetaDataTop = "height_95%_HPD={" + minHeight2HPD + "," + maxHeight2HPD + "}";
                    acg.addConversion(conv);
                }
            }
        }
    }
}

Example 38 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ACGCladeSystem method collectConversions.

/**
 * Add conversions described on provided acg to the internal list
 * for later summary.
 *
 * @param acg conversion graph from which to extract conversions
 */
public void collectConversions(ConversionGraph acg) {
    getBitSets(acg);
    Map<BitSet, Map<BitSet, Long>> geneFlowTemp = new HashMap<>();
    // Assemble list of conversions for each pair of clades on each locus
    for (Locus locus : acg.getLoci()) {
        conversionListsTemp.clear();
        for (Conversion conv : acg.getConversions(locus)) {
            conv.acgIndex = acgIndex;
            BitSetPair bsPair = new BitSetPair(conv);
            if (!conversionListsTemp.containsKey(bsPair))
                conversionListsTemp.put(bsPair, new ArrayList<>());
            conversionListsTemp.get(bsPair).add(conv);
            // Record gene flow
            if (!geneFlowTemp.containsKey(bsPair.from))
                geneFlowTemp.put(bsPair.from, new HashMap<>());
            long oldFlow = 0;
            if (geneFlowTemp.get(bsPair.from).containsKey(bsPair.to))
                oldFlow = geneFlowTemp.get(bsPair.from).get(bsPair.to);
            geneFlowTemp.get(bsPair.from).put(bsPair.to, oldFlow + conv.getSiteCount());
        }
        // Merge overlapping conversions:
        for (BitSetPair bsPair : conversionListsTemp.keySet()) {
            List<Conversion> merged = mergeOverlappingConvs(conversionListsTemp.get(bsPair));
            if (!conversionLists.containsKey(bsPair))
                conversionLists.put(bsPair, new HashMap<>());
            if (!conversionLists.get(bsPair).containsKey(locus))
                conversionLists.get(bsPair).put(locus, new ArrayList<>());
            conversionLists.get(bsPair).get(locus).addAll(merged);
        }
    }
    geneFlow.add(geneFlowTemp);
    acgIndex += 1;
}

Example 39 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class CFConvSwapExperiment method main.

public static void main(String[] args) throws Exception {
    // Load model
    XMLParser parser = new XMLParser();
    MCMC mcmc = (MCMC) parser.parseFile(new File("inferencePreSimulatedData.xml"));
    State state = mcmc.startStateInput.get();
    state.setStateFileName("problem.state");
    state.restoreFromFile();
    double oldPosterior = state.robustlyCalcPosterior(mcmc.posteriorInput.get());
    ConversionGraph acg = (ConversionGraph) state.getStateNode(0);
    PrintStream ps = new PrintStream("proposal.trees");
    ps.println(acg);
    Node srcNode = acg.getNode(3);
    Node srcNodeS = getSibling(srcNode);
    Node destNode = acg.getNode(6);
    double t_srcNodeP = srcNode.getParent().getHeight();
    disconnectEdge(acg, srcNode);
    Locus locus = acg.getLoci().get(0);
    Conversion convToReplace = acg.getConversions(locus).get(27);
    acg.deleteConversion(convToReplace);
    Node srcNodeP = srcNode.getParent();
    connectEdge(acg, srcNode, destNode, convToReplace.getHeight2());
    // Move Conversions
    for (Conversion conv : acg.getConversions(locus)) {
        boolean moved = false;
        if (conv.getNode1() == srcNode && conv.getHeight1() > srcNodeP.getHeight()) {
            conv.setNode1(destNode);
            moved = true;
        }
        if (conv.getNode2() == srcNode && conv.getHeight2() > srcNodeP.getHeight()) {
            conv.setNode2(destNode);
            moved = true;
        }
        if (moved) {
            while (conv.getHeight1() > conv.getNode1().getParent().getHeight()) conv.setNode1(conv.getNode1().getParent());
            while (conv.getHeight2() > conv.getNode2().getParent().getHeight()) conv.setNode2(conv.getNode2().getParent());
        }
    }
    Conversion convNew = new Conversion();
    convNew.setLocus(locus);
    convNew.setStartSite(0);
    convNew.setEndSite(4000);
    convNew.setNode1(srcNode);
    convNew.setNode2(srcNodeS);
    convNew.setHeight1(convToReplace.getHeight1());
    convNew.setHeight2(t_srcNodeP);
    acg.addConversion(convNew);
    ps.println(acg);
    double newPosterior = state.robustlyCalcPosterior(mcmc.posteriorInput.get());
    System.out.println(newPosterior - oldPosterior);
// state.setStateFileName("proposal.state");
// state.storeToFile(1);
// Open state file
}

Example 40 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class CFConvSwapExperiment method disconnectEdge.

public static void disconnectEdge(ConversionGraph acg, Node node) {
    if (node.isRoot())
        throw new IllegalArgumentException("Programmer error: " + "root argument passed to disconnectEdge().");
    Node parent = node.getParent();
    Node sister = getSibling(node);
    if (parent.isRoot()) {
        parent.removeChild(sister);
        sister.setParent(null);
    } else {
        Node grandParent = parent.getParent();
        grandParent.removeChild(parent);
        parent.setParent(null);
        parent.removeChild(sister);
        grandParent.addChild(sister);
    }
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (conv.getNode1() == parent)
                conv.setNode1(sister);
            if (conv.getNode2() == parent)
                conv.setNode2(sister);
        }
    }
}