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;
}
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;
}
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;
}
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;
}
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;
}
Aggregations