use of bacter.Locus 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);
}
}
use of bacter.Locus in project bacter by tgvaughan.
the class SimulatedAlignment method initAndValidate.
@Override
public void initAndValidate() {
acg = acgInput.get();
siteModel = siteModelInput.get();
Locus locus;
if (acg.getLoci().size() == 1)
locus = acg.getLoci().get(0);
else {
locus = locusInput.get();
if (locus == null)
throw new IllegalArgumentException("Must specify locus" + " when simulating alignment from ACG associated" + " with multiple loci.");
}
// We can't wait for Alignment.initAndValidate() to get the
// data type for us.
grabDataType();
// Simulate alignment
simulate(locus);
super.initAndValidate();
// Write simulated alignment to disk if requested:
if (outputFileNameInput.get() != null) {
try (PrintStream pstream = new PrintStream(outputFileNameInput.get())) {
if (useNexusInput.get()) {
NexusBuilder nb = new NexusBuilder();
nb.append(new TaxaBlock(acg.getTaxonset()));
nb.append(new CharactersBlock(this));
nb.write(pstream);
} else {
for (String taxonName : acg.getTaxaNames()) {
pstream.print(">" + taxonName + "\n");
pstream.print(getSequenceAsString(taxonName) + "\n");
}
}
} catch (FileNotFoundException e) {
e.printStackTrace();
}
}
}
use of bacter.Locus 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);
}
}
}
use of bacter.Locus in project bacter by tgvaughan.
the class AddRemoveDetour method addDetour.
/**
* Detour addition variant of move.
*
* @return log HGF
*/
private double addDetour() {
double logHGF = 0.0;
if (acg.getTotalConvCount() == 0)
return Double.NEGATIVE_INFINITY;
// Select conversion at random
Conversion conv = chooseConversion();
logHGF -= Math.log(1.0 / acg.getTotalConvCount());
// Select detour times:
double t1 = conv.getHeight1() + Randomizer.nextDouble() * (conv.getHeight2() - conv.getHeight1());
double t2 = conv.getHeight1() + Randomizer.nextDouble() * (conv.getHeight2() - conv.getHeight1());
double tLower, tUpper;
if (t1 < t2) {
tLower = t1;
tUpper = t2;
} else {
tLower = t2;
tUpper = t1;
}
logHGF -= Math.log(2.0) - 2.0 * Math.log(conv.getHeight2() - conv.getHeight1());
// Select non-root node below detour edge at random
Node detour = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
logHGF -= Math.log(1.0 / (acg.getNodeCount() - 1));
// Abort if selected detour edge does not contain tLower and tUpper
if (detour.getHeight() > tLower || detour.getParent().getHeight() < tUpper)
return Double.NEGATIVE_INFINITY;
// Abort if conv is already attached to selected detour edge:
if (detour == conv.getNode1() || detour == conv.getNode2())
return Double.NEGATIVE_INFINITY;
Conversion convA = new Conversion();
convA.setLocus(conv.getLocus());
convA.setNode1(conv.getNode1());
convA.setHeight1(conv.getHeight1());
convA.setNode2(detour);
convA.setHeight2(tLower);
convA.setStartSite(conv.getStartSite());
convA.setEndSite(conv.getEndSite());
Conversion convB = new Conversion();
convB.setNode1(detour);
convB.setHeight1(tUpper);
convB.setNode2(conv.getNode2());
convB.setHeight2(conv.getHeight2());
logHGF -= drawAffectedRegion(convB);
acg.deleteConversion(conv);
acg.addConversion(convA);
acg.addConversion(convB);
// Count number of node1s and node2s attached to detour edge
int node1Count = 0;
int node2Count = 0;
for (Locus locus : acg.getLoci()) {
for (Conversion thisConv : acg.getConversions(locus)) {
if (thisConv.getNode1() == detour && thisConv.getNode2() != detour)
node1Count += 1;
if (thisConv.getNode2() == detour && thisConv.getNode1() != detour)
node2Count += 1;
}
}
// Incorporate probability of reverse move:
logHGF += Math.log(1.0 / ((acg.getNodeCount() - 1) * node1Count * node2Count));
return logHGF;
}
use of bacter.Locus in project bacter by tgvaughan.
the class AddRemoveRedundantConversion method drawAffectedRegion.
private double drawAffectedRegion(Conversion conv) {
double logP = 0.0;
Locus locus = acg.getLoci().get(Randomizer.nextInt(acg.getLoci().size()));
conv.setLocus(locus);
logP += Math.log(1.0 / acg.getLoci().size());
if (!acg.wholeLocusModeOn()) {
int site1 = Randomizer.nextInt(locus.getSiteCount());
int site2 = Randomizer.nextInt(locus.getSiteCount());
if (site1 < site2) {
conv.setStartSite(site1);
conv.setEndSite(site2);
} else {
conv.setStartSite(site2);
conv.setEndSite(site1);
}
logP += 2.0 * Math.log(1.0 / locus.getSiteCount());
if (site1 != site2)
logP += Math.log(2.0);
} else {
conv.setStartSite(0);
conv.setEndSite(locus.getSiteCount() - 1);
}
return logP;
}
Aggregations