Search in sources :

Example 41 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ConvertedRegionBoundaryShift method proposal.

@Override
public double proposal() {
    if (acg.getTotalConvCount() < 1 || acg.wholeLocusModeOn())
        return Double.NEGATIVE_INFINITY;
    // Select random conversion and region edge:
    Conversion conv = chooseConversion();
    boolean moveStart = Randomizer.nextBoolean();
    int currentLocus, minLocus, maxLocus;
    if (moveStart) {
        currentLocus = conv.getStartSite();
        maxLocus = conv.getEndSite();
        minLocus = 0;
    } else {
        currentLocus = conv.getEndSite();
        minLocus = conv.getStartSite();
        maxLocus = conv.getLocus().getSiteCount() - 1;
    }
    int radius = (int) Math.round(conv.getLocus().getSiteCount() * apertureSizeInput.get()) / 2;
    int newLocus = currentLocus + Randomizer.nextInt(2 * radius + 1) - radius;
    if (newLocus < minLocus || newLocus > maxLocus)
        return Double.NEGATIVE_INFINITY;
    if (moveStart)
        conv.setStartSite(newLocus);
    else
        conv.setEndSite(newLocus);
    return 0;
}
Also used : Conversion(bacter.Conversion)

Example 42 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ConvertedRegionShift method proposal.

@Override
public double proposal() {
    if (acg.getTotalConvCount() < 1 || acg.wholeLocusModeOn())
        return Double.NEGATIVE_INFINITY;
    Conversion conv = chooseConversion();
    int radius = (int) Math.round(conv.getLocus().getSiteCount() * apertureSizeInput.get()) / 2;
    int delta = Randomizer.nextInt(radius * 2 + 1) - radius;
    if (conv.getEndSite() + delta > conv.getLocus().getSiteCount() - 1)
        return Double.NEGATIVE_INFINITY;
    if (conv.getStartSite() + delta < 0)
        return Double.NEGATIVE_INFINITY;
    conv.setStartSite(conv.getStartSite() + delta);
    conv.setEndSite(conv.getEndSite() + delta);
    return 0.0;
}
Also used : Conversion(bacter.Conversion)

Example 43 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class MergeSplitConversion method mergeProposal.

/**
 * Perform merge portion of merge/split move.
 *
 * @param locus locus on which to apply move
 * @return log of move HR
 */
private double mergeProposal(Locus locus) {
    double logHGF = 0.0;
    if (acg.getConvCount(locus) < 2)
        return Double.NEGATIVE_INFINITY;
    Conversion conv1 = acg.getConversions(locus).get(Randomizer.nextInt(acg.getConvCount(locus)));
    Conversion conv2;
    do {
        conv2 = acg.getConversions(locus).get(Randomizer.nextInt(acg.getConvCount(locus)));
    } while (conv2 == conv1);
    if (conv2.getNode1() != conv1.getNode1() || conv2.getNode2() != conv1.getNode2())
        return Double.NEGATIVE_INFINITY;
    logHGF -= Math.log(1.0 / (acg.getConvCount(locus) * (acg.getConvCount(locus) - 1)));
    int minStart = conv1.getStartSite() < conv2.getStartSite() ? conv1.getStartSite() : conv2.getStartSite();
    int maxEnd = conv1.getEndSite() > conv2.getEndSite() ? conv1.getEndSite() : conv2.getEndSite();
    logHGF += 2.0 * Math.log(0.5 / (maxEnd - minStart + 1));
    logHGF += Math.log(1.0 / conv1.getNode1().getLength());
    if (conv1.getNode2().isRoot()) {
        double lambda = 1.0 / (conv1.getHeight2() - conv1.getNode2().getHeight());
        logHGF += -lambda * (conv2.getHeight2() - conv1.getNode2().getHeight()) + Math.log(lambda);
    } else {
        logHGF += Math.log(1.0 / conv1.getNode2().getLength());
    }
    acg.deleteConversion(conv2);
    conv1.setStartSite(minStart);
    conv1.setEndSite(maxEnd);
    logHGF += Math.log(1.0 / acg.getConvCount(locus));
    return logHGF;
}
Also used : Conversion(bacter.Conversion)

Example 44 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ReplaceConversion method proposal.

@Override
public double proposal() {
    if (acg.getTotalConvCount() == 0)
        return Double.NEGATIVE_INFINITY;
    double logHGF = 0;
    // Select conversion to replace:
    Conversion conv = chooseConversion();
    // Incorporate conversion prior probability into HGF
    logHGF += getEdgeAttachmentProb(conv) + getAffectedRegionProb(conv);
    // Remove conversion
    acg.deleteConversion(conv);
    // Draw replacement conversion from prior, incoroporating
    // probability into HGF
    conv = new Conversion();
    logHGF -= attachEdge(conv) + drawAffectedRegion(conv);
    // Add replacement conversion
    acg.addConversion(conv);
    return logHGF;
}
Also used : Conversion(bacter.Conversion)

Example 45 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class CFOperator method expandConversions.

/**
 * Take length of new edge above srcNode that is greater than the
 * original height of srcNode.parent and shifts a random fraction of
 * conversion attachments to it from the lineage above destNode.
 *
 * In the case that destNode was the root, the conversions starting
 * above destNode are drawn from the prior.
 *
 * Assumes topology has not yet been altered.
 *
 * @param srcNode source node for the move
 * @param destNode dest node for the move
 * @param destTime new time drawn for srcNode.P.
 * @return log probability of new conversion configuration.
 */
protected double expandConversions(Node srcNode, Node destNode, double destTime) {
    double logP = 0.0;
    double volatileHeight = acg.getRoot().getHeight();
    boolean forwardRootMove = destTime > volatileHeight;
    Node node = srcNode.getParent();
    while (node != null) {
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if (conv.getNode1() == node && conv.getHeight1() < destTime) {
                    if (Randomizer.nextBoolean())
                        conv.setNode1(srcNode);
                    logP += Math.log(0.5);
                }
                if (conv.getNode2() == node && conv.getHeight2() < destTime) {
                    if (Randomizer.nextBoolean())
                        conv.setNode2(srcNode);
                    logP += Math.log(0.5);
                }
            }
        }
        node = node.getParent();
    }
    // Apply topology modifications.
    disconnectEdge(srcNode);
    connectEdge(srcNode, destNode, destTime);
    // this was a forward root move
    if (forwardRootMove) {
        acg.setRoot(srcNode.getParent());
        double L = 2.0 * (destTime - volatileHeight);
        double Nexp = L * rhoInput.get().getValue() * (acg.getTotalSequenceLength() + acg.getLoci().size() * (deltaInput.get().getValue() - 1.0));
        int N = (int) Randomizer.nextPoisson(Nexp);
        // Factorial cancels
        logP += -Nexp + N * Math.log(Nexp);
        for (int i = 0; i < N; i++) {
            Conversion conv = new Conversion();
            double u = Randomizer.nextDouble() * L;
            if (u < 0.5 * L) {
                conv.setNode1(destNode);
                conv.setHeight1(volatileHeight + u);
            } else {
                conv.setNode1(srcNode);
                conv.setHeight1(volatileHeight + u - 0.5 * L);
            }
            logP += Math.log(1.0 / L) + drawAffectedRegion(conv) + coalesceEdge(conv);
            acg.addConversion(conv);
        }
    }
    return logP;
}
Also used : Node(beast.evolution.tree.Node) Locus(bacter.Locus) Conversion(bacter.Conversion)

Aggregations

Conversion (bacter.Conversion)49 Locus (bacter.Locus)27 Node (beast.evolution.tree.Node)24 ConversionGraph (bacter.ConversionGraph)7 ArrayList (java.util.ArrayList)7 Test (org.junit.Test)6 RealParameter (beast.core.parameter.RealParameter)4 SiteModel (beast.evolution.sitemodel.SiteModel)4 JukesCantor (beast.evolution.substitutionmodel.JukesCantor)4 TaxonSet (beast.evolution.alignment.TaxonSet)3 ConstantPopulation (beast.evolution.tree.coalescent.ConstantPopulation)3 ClusterTree (beast.util.ClusterTree)3 PrintStream (java.io.PrintStream)3 SimulatedACG (bacter.model.SimulatedACG)2 CFEventList (bacter.CFEventList)1 ACGCoalescent (bacter.model.ACGCoalescent)1 ACGLogReader (bacter.util.ACGLogReader)1 BacterACGLogReader (bacter.util.BacterACGLogReader)1 COACGLogFileReader (bacter.util.COACGLogFileReader)1 State (beast.core.State)1