Examples with RealParameter beast.core.parameter.RealParameter used on opensource projects

Example 36 with RealParameter

use of beast.core.parameter.RealParameter in project MultiTypeTree by tgvaughan.

the class MultiTypeTreeScale method proposal.

@Override
public double proposal() {
    // Choose scale factor:
    double u = Randomizer.nextDouble();
    double f = u * scaleFactorInput.get() + (1.0 - u) / scaleFactorInput.get();
    // Keep track of Hastings ratio:
    double logf = Math.log(f);
    double logHR = -2 * logf;
    // Scale colour change times on external branches:
    if (!useOldTreeScalerInput.get()) {
        for (Node leaf : mtTree.getExternalNodes()) {
            double lold = leaf.getParent().getHeight() - leaf.getHeight();
            double lnew = f * leaf.getParent().getHeight() - leaf.getHeight();
            for (int c = 0; c < ((MultiTypeNode) leaf).getChangeCount(); c++) {
                double oldTime = ((MultiTypeNode) leaf).getChangeTime(c);
                double newTime = leaf.getHeight() + (oldTime - leaf.getHeight()) * lnew / lold;
                ((MultiTypeNode) leaf).setChangeTime(c, newTime);
            }
            logHR += ((MultiTypeNode) leaf).getChangeCount() * Math.log(lnew / lold);
        }
    } else {
        for (Node leaf : mtTree.getExternalNodes()) {
            for (int c = 0; c < ((MultiTypeNode) leaf).getChangeCount(); c++) {
                double oldTime = ((MultiTypeNode) leaf).getChangeTime(c);
                ((MultiTypeNode) leaf).setChangeTime(c, f * oldTime);
                logHR += logf;
            }
        }
    }
    // Scale internal node heights and colour change times:
    for (Node node : mtTree.getInternalNodes()) {
        node.setHeight(node.getHeight() * f);
        logHR += logf;
        for (int c = 0; c < ((MultiTypeNode) node).getChangeCount(); c++) {
            double oldTime = ((MultiTypeNode) node).getChangeTime(c);
            ((MultiTypeNode) node).setChangeTime(c, f * oldTime);
            logHR += logf;
        }
    }
    // Reject invalid tree scalings:
    if (f < 1.0) {
        for (Node leaf : mtTree.getExternalNodes()) {
            if (leaf.getParent().getHeight() < leaf.getHeight())
                return Double.NEGATIVE_INFINITY;
            if (((MultiTypeNode) leaf).getChangeCount() > 0 && ((MultiTypeNode) leaf).getChangeTime(0) < leaf.getHeight())
                if (!useOldTreeScalerInput.get())
                    throw new IllegalStateException("Scaled colour change time " + "has dipped below age of leaf - this should never " + "happen!");
                else
                    return Double.NEGATIVE_INFINITY;
        }
    }
    // Scale parameters:
    for (int pidx = 0; pidx < parametersInput.get().size(); pidx++) {
        RealParameter param = parametersInput.get().get(pidx);
        for (int i = 0; i < param.getDimension(); i++) {
            if (!indicatorsUsed || indicatorsInput.get().get(pidx).getValue(i)) {
                double oldValue = param.getValue(i);
                double newValue = oldValue * f;
                if (newValue < param.getLower() || newValue > param.getUpper())
                    return Double.NEGATIVE_INFINITY;
                param.setValue(i, newValue);
                logHR += logf;
            }
        }
    }
    // Scale parameters inversely:
    for (int pidx = 0; pidx < parametersInverseInput.get().size(); pidx++) {
        RealParameter param = parametersInverseInput.get().get(pidx);
        for (int i = 0; i < param.getDimension(); i++) {
            if (!indicatorsInverseUsed || indicatorsInverseInput.get().get(pidx).getValue(i)) {
                double oldValue = param.getValue(i);
                double newValue = oldValue / f;
                if (newValue < param.getLower() || newValue > param.getUpper())
                    return Double.NEGATIVE_INFINITY;
                param.setValue(i, newValue);
                logHR -= logf;
            }
        }
    }
    // Return Hastings ratio:
    return logHR;
}

Example 37 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodCaching.

@Test
public void testLikelihoodCaching() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 5.0 / locus.getSiteCount(), "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
    State state = new State();
    state.initByName("stateNode", acg);
    state.initialise();
    System.out.println(acg);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", new RealParameter("1"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    // Calculate likelihood 1:
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logP1 = argLikelihood.calculateLogP();
    double logP1prime = argLikelihoodSlow.calculateLogP();
    double relError = 2.0 * Math.abs(logP1 - logP1prime) / Math.abs(logP1 + logP1prime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP1, logP1prime, relError);
    assertTrue(relError < 1e-13);
    // Add a single recombination event
    Node node1 = acg.getExternalNodes().get(0);
    Node node2 = node1.getParent();
    double height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
    double height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
    int startLocus = 500;
    int endLocus = 600;
    Conversion recomb1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb1);
    double logP2 = argLikelihood.calculateLogP();
    double logP2prime = argLikelihoodSlow.calculateLogP();
    relError = 2.0 * Math.abs(logP2 - logP2prime) / Math.abs(logP2 + logP2prime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP2, logP2prime, relError);
    assertTrue(relError < 1e-13);
    state.restore();
    double logP3 = argLikelihood.calculateLogP();
    double logP3prime = argLikelihoodSlow.calculateLogP();
    relError = 2.0 * Math.abs(logP3 - logP3prime) / Math.abs(logP3 + logP3prime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP3, logP3prime, relError);
    assertTrue(relError < 1e-13);
}

Example 38 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class SimulatedAlignmentTest method test.

@Test
public void test() throws Exception {
    // Randomizer.setSeed(26);
    Randomizer.setSeed(7);
    Locus locus = new Locus("locus", 100000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 1.0 / 100000, "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
    System.out.println(acg);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", new RealParameter("1.0"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    for (Region region : acg.getRegions(locus)) System.out.println(new MarginalTree(acg, region));
    // (Should be enough info here for precise agreement)
    for (Region region : acg.getRegions(locus)) {
        Alignment margAlign = createMarginalAlignment(alignment, region);
        ClusterTree upgmaTree = new ClusterTree();
        upgmaTree.initByName("clusterType", "upgma", "taxa", margAlign);
        MarginalTree marginalTree = new MarginalTree(acg, region);
        // outfMarg.println(marginalTree.getRoot() + ";");
        // outfMarg.flush();
        // 
        // outfUPGMA.println(upgmaTree.getRoot() + ";");
        // outfUPGMA.flush();
        assertTrue(topologiesEquivalent(marginalTree.getRoot(), upgmaTree.getRoot()));
    }
// outfMarg.close();
// outfUPGMA.close();
}

Example 39 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class SkylinePopulationFunctionTest method test3.

@Test
public void test3() throws Exception {
    SkylinePopulationFunction skyline = new SkylinePopulationFunction();
    skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0"), "piecewiseLinear", true);
    for (double t = 0.0; t < 10; t += 0.01) assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
}

Example 40 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class SkylinePopulationFunctionTest method test4.

@Test
public void test4() throws Exception {
    SkylinePopulationFunction skyline = new SkylinePopulationFunction();
    skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0"), "piecewiseLinear", true);
    for (CFEventList.Event cfEvent : acg.getCFEvents()) {
        double t = cfEvent.getHeight();
        assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
    }
}