use of bacter.Conversion in project bacter by tgvaughan.
the class CFConversionSwap method proposal.
@Override
public double proposal() {
double logHGF = 0.0;
// Determine whether we can apply this operator:
if (acg.getLeafNodeCount() < 3 || acg.getTotalConvCount() == 0)
return Double.NEGATIVE_INFINITY;
// Acquire list of conversions compatible with swap:
List<Conversion> compatible = getCompatibleConversions();
if (compatible.isEmpty())
return Double.NEGATIVE_INFINITY;
// Select conversion to replace
Conversion replaceConversion = compatible.get(Randomizer.nextInt(compatible.size()));
logHGF -= Math.log(1.0 / compatible.size());
acg.deleteConversion(replaceConversion);
Node srcNode = replaceConversion.getNode1();
Node destNode = replaceConversion.getNode2();
Node srcNodeP = srcNode.getParent();
Node srcNodeS = getSibling(srcNode);
double t_srcNodeP = srcNodeP.getHeight();
if (destNode == srcNode.getParent())
destNode = srcNodeS;
double newTime = replaceConversion.getHeight2();
// Conversion modification:
replaceConversion.setNode2(srcNodeS);
replaceConversion.setHeight2(t_srcNodeP);
logHGF += getAffectedRegionProb(replaceConversion);
logHGF -= drawAffectedRegion(replaceConversion);
acg.addConversion(replaceConversion);
if (!includeRootInput.get() && (srcNodeP.isRoot() || destNode.isRoot()))
return Double.NEGATIVE_INFINITY;
// Perform necessary conversion expansions/collapses:
if (newTime < t_srcNodeP) {
logHGF += collapseConversions(srcNode, destNode, newTime);
} else {
logHGF -= expandConversions(srcNode, destNode, newTime);
}
logHGF += Math.log(1.0 / getCompatibleConversions().size());
assert !acg.isInvalid() : "CFCS proposed invalid state.";
return logHGF;
}
use of bacter.Conversion in project bacter by tgvaughan.
the class DifferenceFromTrueACG method countSampledTrueConversions.
/**
* Count the number of true conversions which have correspondences on the
* provided sampled ARG.
*
* @param trueACG true arg
* @param trueClades clades in true arg
* @param acg sampled arg
* @param clades clades in sampled arg
* @param boundaryTol minimum relative error in region boundaries to allow.
* @param convHist
* @return number of found conversions
*/
public static int countSampledTrueConversions(ConversionGraph trueACG, Clade[] trueClades, ConversionGraph acg, Clade[] clades, double boundaryTol, double ageTol, Map<Conversion, Integer> convHist) {
int count = 0;
for (Conversion trueConv : trueACG.getConversions(trueACG.getLoci().get(0))) {
Clade trueFromClade = trueClades[trueConv.getNode1().getNr()];
Clade trueToClade = trueClades[trueConv.getNode1().getNr()];
for (Conversion conv : acg.getConversions(acg.getLoci().get(0))) {
Clade fromClade = clades[conv.getNode1().getNr()];
Clade toClade = clades[conv.getNode1().getNr()];
if (fromClade.equals(trueFromClade) && toClade.equals(trueToClade)) {
if (Math.abs(trueConv.getStartSite() - conv.getStartSite()) / trueConv.getSiteCount() <= boundaryTol && Math.abs(trueConv.getEndSite() - conv.getEndSite()) / trueConv.getSiteCount() <= boundaryTol) // && Math.abs(trueConv.getHeight1()-conv.getHeight1())/trueConv.getHeight1() <= ageTol
// && Math.abs(trueConv.getHeight2()-conv.getHeight2())/trueConv.getHeight2() <= ageTol
{
count += 1;
convHist.put(trueConv, convHist.get(trueConv) + 1);
break;
}
}
}
}
return count;
}
use of bacter.Conversion in project bacter by tgvaughan.
the class DifferenceFromTrueACG method countRecoveredTrueConvs.
public static int countRecoveredTrueConvs(Map<Conversion, Integer> convHist, int nSampledACGs, double minimumSupport) {
int count = 0;
int sampleThresh = (int) Math.round(nSampledACGs * minimumSupport);
for (Conversion conv : convHist.keySet()) {
if (convHist.get(conv) >= sampleThresh)
count += 1;
}
return count;
}
use of bacter.Conversion 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;
}
use of bacter.Conversion in project bacter by tgvaughan.
the class MergeSplitConversion method splitProposal.
/**
* Perform split component of merge/split move.
*
* @param locus locus on which to apply move
* @return log of move HR
*/
private double splitProposal(Locus locus) {
double logHGF = 0.0;
if (acg.getConvCount(locus) == 0)
return Double.NEGATIVE_INFINITY;
Conversion conv1 = acg.getConversions(locus).get(Randomizer.nextInt(acg.getConvCount(locus)));
logHGF -= Math.log(1.0 / acg.getConvCount(locus));
Conversion conv2 = new Conversion();
conv2.setLocus(locus);
conv2.setNode1(conv1.getNode1());
conv2.setNode2(conv1.getNode2());
int s1, s2, e1, e2;
int m1 = conv1.getStartSite() + Randomizer.nextInt(conv1.getSiteCount());
int m2 = conv1.getStartSite() + Randomizer.nextInt(conv1.getSiteCount());
if (Randomizer.nextBoolean()) {
s1 = conv1.getStartSite();
s2 = m1;
} else {
s1 = m1;
s2 = conv1.getStartSite();
}
if (Randomizer.nextBoolean()) {
e1 = conv1.getEndSite();
e2 = m2;
} else {
e1 = m2;
e2 = conv1.getEndSite();
}
if (e1 < s1 || e2 < s2)
return Double.NEGATIVE_INFINITY;
logHGF -= 2.0 * Math.log(0.5 / (conv1.getSiteCount()));
conv1.setStartSite(s1);
conv1.setEndSite(e1);
conv2.setStartSite(s2);
conv2.setEndSite(e2);
conv2.setHeight1(conv1.getNode1().getHeight() + Randomizer.nextDouble() * conv1.getNode1().getLength());
logHGF -= Math.log(1.0 / conv1.getNode1().getLength());
if (conv1.getNode2().isRoot()) {
double lambda = 1.0 / (conv1.getHeight2() - conv1.getNode2().getHeight());
conv2.setHeight1(conv1.getNode2().getHeight() + Randomizer.nextExponential(lambda));
logHGF -= -lambda * (conv2.getHeight2() - conv2.getNode2().getHeight()) + Math.log(lambda);
} else {
conv2.setHeight2(conv1.getNode2().getHeight() + Randomizer.nextDouble() * conv1.getNode2().getLength());
logHGF -= Math.log(1.0 / conv1.getNode2().getLength());
}
if (conv2.getHeight2() < conv2.getHeight1())
return Double.NEGATIVE_INFINITY;
acg.addConversion(conv2);
logHGF += Math.log(1.0 / (acg.getConvCount(locus) * (acg.getConvCount(locus) - 1)));
return logHGF;
}
Aggregations