Examples with Conversion bacter.Conversion used on opensource projects

Example 31 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class DifferenceFromTrueACG method main.

public static void main(String[] args) throws IOException, XMLStreamException {
    Options options = processArguments(args);
    // Load true ARG
    BacterACGLogReader truthReader = new BacterACGLogReader(options.truthFile, 0);
    if (truthReader.getACGCount() != 1) {
        System.out.println("Expected exactly 1 ACG in truth file. Found " + truthReader.getACGCount());
        System.exit(1);
    }
    ConversionGraph trueACG = null;
    for (ConversionGraph acg : truthReader) trueACG = acg;
    // Determine clades present in truth
    Clade[] trueClades = new Clade[trueACG.getNodeCount()];
    getClades(trueClades, trueACG.getRoot());
    Set<Clade> trueCladeSet = new HashSet<>(Arrays.asList(trueClades));
    // Set up histograms
    Map<Clade, Integer> cladeHist = new HashMap<>();
    for (Clade clade : trueClades) cladeHist.put(clade, 0);
    Map<Conversion, Integer> convHist = new HashMap<>();
    for (Conversion conv : trueACG.getConversions(trueACG.getLoci().get(0))) convHist.put(conv, 0);
    // Set up ARG log file reader
    ACGLogReader logReader;
    if (options.useCOFormat) {
        logReader = new COACGLogFileReader(options.logFile, options.burninPerc);
    } else {
        logReader = new BacterACGLogReader(options.logFile, options.burninPerc);
    }
    try (PrintStream ps = new PrintStream(options.outFile)) {
        ps.println("trueCladeCount sampledTrueCladeCount trueConvCount sampledConvCount sampledTrueConvCount");
        for (ConversionGraph acg : logReader) {
            Clade[] clades = new Clade[acg.getNodeCount()];
            getClades(clades, acg.getRoot());
            List<Double> timeErrors = new ArrayList<>();
            int sampledTrueClades = countSampledTrueClades(trueClades, clades, options.ageTol, cladeHist);
            int sampledTrueConvs = countSampledTrueConversions(trueACG, trueClades, acg, clades, options.boundaryTol, options.ageTol, convHist);
            ps.println(trueACG.getNodeCount() + "\t" + sampledTrueClades + "\t" + trueACG.getConvCount(trueACG.getLoci().get(0)) + "\t" + acg.getConvCount(acg.getLoci().get(0)) + "\t" + sampledTrueConvs);
        }
    }
    try (PrintStream ps = new PrintStream(options.summaryFile)) {
        ps.println("trueCladeCount recoveredCladeCount trueConvCount recoveredConvCount");
        int recoveredTrueClades = countRecoveredTrueClades(cladeHist, logReader.getCorrectedACGCount(), 0.5);
        int recoveredTrueConvs = countRecoveredTrueConvs(convHist, logReader.getCorrectedACGCount(), 0.5);
        ps.println(trueACG.getNodeCount() + "\t" + recoveredTrueClades + "\t" + trueACG.getConvCount(trueACG.getLoci().get(0)) + "\t" + recoveredTrueConvs);
    }
}

Example 32 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ACGCoalescent method calculateLogP.

@Override
public double calculateLogP() {
    // Check whether conversion count exceeds bounds.
    if (acg.getTotalConvCount() < lowerCCBoundInput.get() || acg.getTotalConvCount() > upperCCBoundInput.get())
        return Double.NEGATIVE_INFINITY;
    logP = calculateClonalFrameLogP();
    double poissonMean = rhoInput.get().getValue() * acg.getClonalFrameLength() * (acg.getTotalSequenceLength() + acg.getLoci().size() * (deltaInput.get().getValue() - 1.0));
    // Probability of conversion count:
    if (poissonMean > 0.0) {
        logP += -poissonMean + acg.getTotalConvCount() * Math.log(poissonMean);
    // - GammaFunction.lnGamma(acg.getConvCount()+1);
    } else {
        if (acg.getTotalConvCount() > 0)
            logP = Double.NEGATIVE_INFINITY;
    }
    for (Locus locus : acg.getLoci()) for (Conversion conv : acg.getConversions(locus)) logP += calculateConversionLogP(conv);
    if (lowerCCBoundInput.get() > 0 || upperCCBoundInput.get() < Integer.MAX_VALUE) {
        try {
            logP -= new PoissonDistributionImpl(poissonMean).cumulativeProbability(lowerCCBoundInput.get(), upperCCBoundInput.get());
        } catch (MathException e) {
            throw new RuntimeException("Error computing modification to ARG " + "prior density required by conversion number constraint.");
        }
    }
    return logP;
}

Example 33 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class SimulatedACG method generateConversions.

private void generateConversions() {
    // Draw number of conversions:
    int Nconv = (int) Randomizer.nextPoisson(rho * getClonalFrameLength() * (getTotalSequenceLength() + (delta - 1.0) * getLoci().size()));
    // Generate conversions:
    for (int i = 0; i < Nconv; i++) {
        // Choose alignment
        double u = Randomizer.nextDouble() * (getTotalSequenceLength() + (delta - 1.0) * getLoci().size());
        Locus affectedLocus = null;
        for (Locus locus : getLoci()) {
            if (u < locus.getSiteCount() + delta - 1.0) {
                affectedLocus = locus;
                break;
            } else
                u -= locus.getSiteCount() + delta - 1.0;
        }
        if (affectedLocus == null)
            throw new IllegalStateException("Programmer error: " + "locus choice loop fell through.");
        int startSite, endSite;
        if (u < delta) {
            startSite = 0;
        } else {
            startSite = (int) Math.ceil(u - delta);
        }
        endSite = startSite + (int) Randomizer.nextGeometric(1.0 / delta);
        endSite = Math.min(endSite, affectedLocus.getSiteCount() - 1);
        Conversion conv = new Conversion();
        conv.setLocus(affectedLocus);
        conv.setStartSite(startSite);
        conv.setEndSite(endSite);
        associateConversionWithCF(conv);
        addConversion(conv);
    }
}

Example 34 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ACGOperator method disconnectEdge.

/**
 * Disconnect edge <node, node.parent> from its sister and
 * grandparent (if the grandparent exists), forming a new edge
 * <sister, grandparent>. All conversions originally above node.parent
 * are re-attached to sister.
 *
 * Conversions on edge above node are not modified.
 *
 * @param node base of edge to detach.
 */
protected void disconnectEdge(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);
        }
    }
}

Example 35 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class AddRemoveConversion method proposal.

@Override
public double proposal() {
    double logHGF = 0;
    if (Randomizer.nextBoolean()) {
        // Add
        logHGF += Math.log(1.0 / (acg.getTotalConvCount() + 1));
        logHGF -= drawNewConversion();
    } else {
        if (acg.getTotalConvCount() == 0)
            return Double.NEGATIVE_INFINITY;
        // Select conversion to remove:
        Conversion conv = chooseConversion();
        // Calculate HGF
        logHGF += getConversionProb(conv);
        logHGF -= Math.log(1.0 / acg.getTotalConvCount());
        // Remove conversion
        acg.deleteConversion(conv);
    }
    return logHGF;
}