use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class SkylinePopulationFunctionTest method test2.
@Test
public void test2() throws Exception {
SkylinePopulationFunction skyline = new SkylinePopulationFunction();
skyline.initByName("acg", acg, "popSizes", new RealParameter("5.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0 0 0 0"));
for (CFEventList.Event cfEvent : acg.getCFEvents()) {
double t = cfEvent.getHeight();
assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
}
}
use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class SkylinePopulationFunctionTest method test1.
@Test
public void test1() throws Exception {
SkylinePopulationFunction skyline = new SkylinePopulationFunction();
skyline.initByName("acg", acg, "popSizes", new RealParameter("5.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0 0 0 0"));
for (double t = 0.0; t < 10; t += 0.01) assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
}
use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class AddRemoveConversionTest method testHR.
/**
* Tests that probability density of forward move calculated
* by drawNewConversion() matches probability density of backward
* move calculated by getConversionProb().
*
* @throws Exception
*/
@Test
public void testHR() throws Exception {
ConstantPopulation popFunc = new ConstantPopulation();
popFunc.initByName("popSize", new RealParameter("1.0"));
Locus locus = new Locus("locus", 10000);
TaxonSet taxonSet = getTaxonSet(10);
SimulatedACG acg = new SimulatedACG();
acg.initByName("rho", 1.0 / locus.getSiteCount(), "delta", 50.0, "locus", locus, "taxonset", taxonSet, "populationModel", popFunc);
AddRemoveConversion operator = new AddRemoveConversion();
// Loop until a valid proposal is made
double logP1;
List<Conversion> oldConversions;
do {
operator.initByName("weight", 1.0, "acg", acg, "delta", new RealParameter("50.0"), "populationModel", popFunc);
oldConversions = Lists.newArrayList(acg.getConversions(locus));
logP1 = operator.drawNewConversion();
} while (Double.isInfinite(logP1));
System.out.println("logP1 = " + logP1);
// Identify new recomination
Conversion newRecomb = null;
for (Conversion recomb : acg.getConversions(locus)) {
if (!oldConversions.contains(recomb))
newRecomb = recomb;
}
assertNotNull(newRecomb);
double logP2 = operator.getConversionProb(newRecomb);
System.out.println("logP2 = " + logP2);
assertTrue(Math.abs(logP1 - logP2) < 1e-10);
}
use of beast.core.parameter.RealParameter 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 beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class AddRemoveConversion method main.
public static void main(String[] args) throws Exception {
ConversionGraph acg = new ConversionGraph();
ConstantPopulation popFunc = new ConstantPopulation();
AddRemoveConversion operator = new AddRemoveConversion();
operator.initByName("weight", 1.0, "acg", acg, "populationModel", popFunc, "rho", new RealParameter(Double.toString(1.0 / 10000.0)), "delta", new RealParameter("50.0"));
popFunc.initByName("popSize", new RealParameter("1.0"));
TaxonSet taxonSet = new TaxonSet();
taxonSet.taxonsetInput.setValue(new Taxon("t1"), taxonSet);
taxonSet.taxonsetInput.setValue(new Taxon("t2"), taxonSet);
Locus locus = new Locus("locus", 10000);
try (PrintStream ps = new PrintStream("out.txt")) {
for (int i = 0; i < 100000; i++) {
acg.initByName("locus", locus, "taxonset", taxonSet, "fromString", "(0:1.0,1:1.0)2:0.0;");
operator.drawNewConversion();
ps.println(acg.getConversions(locus).get(0).getStartSite() + " " + acg.getConversions(locus).get(0).getEndSite());
}
}
}
Aggregations