use of bacter.Conversion in project bacter by tgvaughan.
the class ACGScaler method proposal.
@Override
public double proposal() {
// Keep track of number of positively scaled elements minus
// negatively scaled elements.
int count = 0;
// Choose scaling factor:
double f = scaleParam + Randomizer.nextDouble() * (1.0 / scaleParam - scaleParam);
// Scale clonal frame:
if (rootOnly) {
acg.getRoot().setHeight(acg.getRoot().getHeight() * f);
count += 1;
} else {
for (Node node : acg.getInternalNodes()) {
node.setHeight(node.getHeight() * f);
count += 1;
}
}
// Scale recombinant edges:
for (Locus locus : acg.getLoci()) {
for (Conversion recomb : acg.getConversions(locus)) {
if (!rootOnly || recomb.getNode1().getParent().isRoot()) {
recomb.setHeight1(recomb.getHeight1() * f);
count += 1;
}
if (!rootOnly || recomb.getNode2().isRoot() || recomb.getNode2().getParent().isRoot()) {
recomb.setHeight2(recomb.getHeight2() * f);
count += 1;
}
if (recomb.getHeight1() < recomb.getNode1().getHeight() || recomb.getHeight2() < recomb.getNode2().getHeight() || recomb.getHeight1() > recomb.getHeight2())
return Double.NEGATIVE_INFINITY;
}
}
// Check for illegal node heights:
if (rootOnly) {
for (Node node : acg.getRoot().getChildren()) {
if (node.getHeight() > node.getParent().getHeight())
return Double.NEGATIVE_INFINITY;
}
} else {
for (Node node : acg.getExternalNodes()) {
if (node.getHeight() > node.getParent().getHeight())
return Double.NEGATIVE_INFINITY;
}
}
for (RealParameter param : parametersInput.get()) {
try {
param.startEditing(null);
param.scale(f);
} catch (Exception e) {
// bounds. Needs to change!
return Double.NEGATIVE_INFINITY;
}
count += param.getDimension();
}
for (RealParameter paramInv : parametersInverseInput.get()) {
try {
paramInv.startEditing(null);
paramInv.scale(1.0 / f);
} catch (Exception e) {
// bounds. Needs to change!
return Double.NEGATIVE_INFINITY;
}
count -= paramInv.getDimension();
}
// Return log of Hastings ratio:
return (count - 2) * Math.log(f);
}
use of bacter.Conversion in project bacter by tgvaughan.
the class AddRemoveConversion method drawNewConversion.
/**
* Add new conversion to ACG, returning the probability density of the
* new edge and converted region location.
*
* @return log of proposal density
*/
public double drawNewConversion() {
Conversion newConversion = new Conversion();
double logP = attachEdge(newConversion) + drawAffectedRegion(newConversion);
acg.addConversion(newConversion);
return logP;
}
use of bacter.Conversion in project bacter by tgvaughan.
the class AddRemoveDetour method removeDetour.
/**
* Detour deletion variant of move.
*
* @return log HGF
*/
private double removeDetour() {
double logHGF = 0.0;
// Choose non-root detour edge at random
Node detour = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
logHGF -= Math.log(1.0 / (acg.getNodeCount() - 1));
List<Conversion> convApotentials = new ArrayList<>();
List<Conversion> convBpotentials = new ArrayList<>();
for (Locus locus : acg.getLoci()) {
for (Conversion conv : acg.getConversions(locus)) {
if (conv.getNode2() == detour && conv.getNode1() != detour)
convApotentials.add(conv);
if (conv.getNode1() == detour && conv.getNode2() != detour)
convBpotentials.add(conv);
}
}
if (convApotentials.isEmpty() || convBpotentials.isEmpty())
return Double.NEGATIVE_INFINITY;
Conversion convA = convApotentials.get(Randomizer.nextInt(convApotentials.size()));
Conversion convB = convBpotentials.get(Randomizer.nextInt(convBpotentials.size()));
if (convA.getHeight2() > convB.getHeight1())
return Double.NEGATIVE_INFINITY;
logHGF -= Math.log(1.0 / (convApotentials.size() * convBpotentials.size()));
double tLowerBound = convA.getHeight1();
double tUpperBound = convB.getHeight2();
Conversion conv = new Conversion();
conv.setNode1(convA.getNode1());
conv.setHeight1(convA.getHeight1());
conv.setNode2(convB.getNode2());
conv.setHeight2(convB.getHeight2());
conv.setStartSite(convA.getStartSite());
conv.setEndSite(convA.getEndSite());
conv.setLocus(convA.getLocus());
acg.deleteConversion(convA);
acg.deleteConversion(convB);
acg.addConversion(conv);
logHGF += Math.log(1.0 / acg.getTotalConvCount()) + Math.log(1.0 / (acg.getNodeCount() - 1)) + Math.log(2.0) - 2.0 * Math.log(tUpperBound - tLowerBound) + getAffectedRegionProb(convB);
return logHGF;
}
use of bacter.Conversion in project bacter by tgvaughan.
the class AddRemoveRedundantConversion method getRedundantConversions.
/**
* Obtain list of redundant conversions.
*
* @param cfNode node identifying edge on CF
* @return conversion list
*/
private List<Conversion> getRedundantConversions(Node cfNode) {
Node cfParent = cfNode.getParent();
List<Conversion> redundantConversions = new ArrayList<>();
double maxL = acg.getRoot().getHeight() * apertureInput.get();
for (Locus locus : acg.getLoci()) {
for (Conversion conv : acg.getConversions(locus)) {
if (((conv.getNode1() == cfNode || conv.getNode1().getParent() == cfNode) && Math.abs(conv.getHeight1() - cfNode.getHeight()) < maxL) && (conv.getNode2() == cfParent || conv.getNode2().getParent() == cfParent) && Math.abs(conv.getHeight2() - cfParent.getHeight()) < maxL) {
redundantConversions.add(conv);
}
}
}
return redundantConversions;
}
use of bacter.Conversion in project bacter by tgvaughan.
the class AddRemoveRedundantConversion method proposal.
@Override
public double proposal() {
double logHGF = 0;
// Select non-root CF node
Node cfNode = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
Node cfParent = cfNode.getParent();
double maxL = apertureInput.get() * acg.getRoot().getHeight();
if (Randomizer.nextBoolean()) {
// Add redundant conversion
Conversion newConv = new Conversion();
double L = Math.min(getEdgeLength(cfNode), maxL);
for (Node child : cfNode.getChildren()) L += Math.min(getEdgeLength(child), maxL);
logHGF -= Math.log(1.0 / L);
double fromPoint = L * Randomizer.nextDouble();
if (fromPoint < Math.min(getEdgeLength(cfNode), maxL)) {
newConv.setNode1(cfNode);
newConv.setHeight1(cfNode.getHeight() + fromPoint);
} else {
fromPoint -= Math.min(getEdgeLength(cfNode), maxL);
for (Node child : cfNode.getChildren()) {
if (fromPoint < Math.min(getEdgeLength(child), maxL)) {
newConv.setNode1(child);
newConv.setHeight1(cfNode.getHeight() - fromPoint);
break;
}
fromPoint -= Math.min(getEdgeLength(child), maxL);
}
}
L = Math.min(getEdgeLength(cfParent), maxL);
for (Node child : cfParent.getChildren()) L += Math.min(getEdgeLength(child), maxL);
logHGF -= Math.log(1.0 / L);
double toPoint = L * Randomizer.nextDouble();
if (toPoint < Math.min(getEdgeLength(cfParent), maxL)) {
newConv.setNode2(cfParent);
newConv.setHeight2(cfParent.getHeight() + toPoint);
} else {
toPoint -= Math.min(getEdgeLength(cfParent), maxL);
for (Node child : cfParent.getChildren()) {
if (toPoint < Math.min(getEdgeLength(child), maxL)) {
newConv.setNode2(child);
newConv.setHeight2(cfParent.getHeight() - toPoint);
break;
}
toPoint -= Math.min(getEdgeLength(child), maxL);
}
}
if (newConv.getHeight1() > newConv.getHeight2())
return Double.NEGATIVE_INFINITY;
logHGF -= drawAffectedRegion(newConv);
// Add conversion
acg.addConversion(newConv);
// Add probability of reverse move deleting this conversion
// to HGF:
logHGF += Math.log(1.0 / getRedundantConversions(cfNode).size());
} else {
// Remove
// Identify redundant conversions
List<Conversion> redundantConversions = getRedundantConversions(cfNode);
if (redundantConversions.size() == 0)
return Double.NEGATIVE_INFINITY;
// Choose conversion to remove
Conversion conv = redundantConversions.get(Randomizer.nextInt(redundantConversions.size()));
logHGF -= Math.log(1.0 / redundantConversions.size());
// Add probability of reverse move generating this conversion
// to HGF:
double L = Math.min(getEdgeLength(cfNode), maxL);
for (Node child : cfNode.getChildren()) L += Math.min(getEdgeLength(child), maxL);
logHGF += Math.log(1.0 / L);
L = Math.min(getEdgeLength(cfParent), maxL);
for (Node child : cfParent.getChildren()) L += Math.min(getEdgeLength(child), maxL);
logHGF += Math.log(1.0 / L);
logHGF += getAffectedRegionProb(conv);
// Remove conversion
acg.deleteConversion(conv);
}
return logHGF;
}
Aggregations