Examples with State beast.core.State used on opensource projects

Example 11 with State

use of beast.core.State in project MultiTypeTree by tgvaughan.

the class TWBR_TS_Test method test1.

@Test
public void test1() throws Exception {
    System.out.println("TWBR_test 1");
    // Fix seed.
    Randomizer.setSeed(42);
    // Assemble initial MultiTypeTree
    String newickStr = "((1[&deme=0]:1,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
    MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
    mtTree.initByName("value", newickStr, "typeLabel", "deme");
    // Assemble migration model:
    RealParameter rateMatrix = new RealParameter("0.1 0.1");
    RealParameter popSizes = new RealParameter("7.0 7.0");
    SCMigrationModel migModel = new SCMigrationModel();
    migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
    // Assemble distribution:
    StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
    distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree);
    // Set up state:
    State state = new State();
    state.initByName("stateNode", mtTree);
    // Set up operator:
    TypedWilsonBaldingRandom operatorTWBR = new TypedWilsonBaldingRandom();
    operatorTWBR.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "mu", 0.2, "alpha", 0.2);
    Operator operatorMTTS = new MultiTypeTreeScale();
    operatorMTTS.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "scaleFactor", 0.8, "useOldTreeScaler", false);
    // Set up stat analysis logger:
    MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
    logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.2, "logEvery", 1000);
    // Set up MCMC:
    MCMC mcmc = new MCMC();
    mcmc.initByName("chainLength", "1000000", "state", state, "distribution", distribution, "operator", operatorTWBR, "operator", operatorMTTS, "logger", logger);
    // Run MCMC:
    mcmc.run();
    System.out.format("height mean = %s\n", logger.getHeightMean());
    System.out.format("height var = %s\n", logger.getHeightVar());
    System.out.format("height ESS = %s\n", logger.getHeightESS());
    // Compare analysis results with truth:
    boolean withinTol = (logger.getHeightESS() > 500) && (Math.abs(logger.getHeightMean() - 19) < 0.5) && (Math.abs(logger.getHeightVar() - 300) < 50);
    Assert.assertTrue(withinTol);
}

Example 12 with State

use of beast.core.State in project MultiTypeTree by tgvaughan.

the class UniformizationRetypeOperator method main.

/**
 * Main method for debugging.
 *
 * @param args
 * @throws java.lang.Exception
 */
public static void main(String[] args) throws Exception {
    // Generate an ensemble of paths along a branch of a tree.
    // Assemble initial MultiTypeTree
    String newickStr = "((1[deme='3']:50)[deme='1']:150,2[deme='1']:200)[deme='1']:0;";
    MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
    mtTree.initByName("newick", newickStr, "typeLabel", "deme", "nTypes", 4);
    // Assemble migration model:
    RealParameter rateMatrix = new RealParameter("0.20 0.02 0.03 0.04 " + "0.05 0.06 0.07 0.08 " + "0.09 0.10 0.11 0.12");
    RealParameter popSizes = new RealParameter("7.0 7.0 7.0 7.0");
    SCMigrationModel migModel = new SCMigrationModel();
    migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes);
    // Set up state:
    State state = new State();
    state.initByName("stateNode", mtTree);
    UniformizationRetypeOperator op = new UniformizationRetypeOperator() {

        @Override
        public void initAndValidate() {
        }

        @Override
        public double proposal() {
            // To change body of generated methods, choose Tools | Templates.
            throw new UnsupportedOperationException("Not supported yet.");
        }
    };
    op.initByName("multiTypeTree", mtTree, "migrationModel", migModel);
    op.mtTree = mtTree;
    migModel.getQ(false).print();
    try (PrintStream outfile = new PrintStream("counts.txt")) {
        outfile.print("totalCounts");
        for (int c = 0; c < migModel.getNTypes(); c++) outfile.print(" counts" + c);
        outfile.println();
        for (int i = 0; i < 10000; i++) {
            MultiTypeNode srcNode = (MultiTypeNode) mtTree.getRoot().getLeft();
            op.retypeBranch(srcNode);
            outfile.print(srcNode.getChangeCount());
            int[] counts = new int[4];
            for (int j = 0; j < srcNode.getChangeCount(); j++) {
                counts[srcNode.getChangeType(j)] += 1;
            }
            for (int c = 0; c < migModel.getNTypes(); c++) outfile.print(" " + counts[c]);
            outfile.println();
        }
    }
}

Example 13 with State

use of beast.core.State 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 14 with State

use of beast.core.State in project MultiTypeTree by tgvaughan.

the class Ewing_Test method test2.

@Test
public void test2() throws Exception {
    System.out.println("Ewing_test 2");
    // Fix seed.
    Randomizer.setSeed(42);
    // Assemble initial MultiTypeTree
    String newickStr = "(((1[&deme=1]:0.5)[&deme=0]:0.5,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
    MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
    mtTree.initByName("value", newickStr, "typeLabel", "deme");
    // Assemble migration model:
    RealParameter rateMatrix = new RealParameter("0.1 0.1");
    RealParameter popSizes = new RealParameter("7.0 7.0");
    SCMigrationModel migModel = new SCMigrationModel();
    migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
    // Assemble distribution:
    StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
    distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree, "checkValidity", true);
    // Set up state:
    State state = new State();
    state.initByName("stateNode", mtTree);
    // Set up operators:
    TypedWilsonBaldingEasy twbOperator = new TypedWilsonBaldingEasy();
    twbOperator.initByName("weight", 1.0, "migrationModel", migModel, "multiTypeTree", mtTree);
    TypedSubtreeExchangeEasy tsxOperator = new TypedSubtreeExchangeEasy();
    tsxOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "isNarrow", true);
    MultiTypeUniform mtuOperator = new MultiTypeUniform();
    mtuOperator.initByName("weight", 1.0, "migrationModel", migModel, "multiTypeTree", mtTree);
    MultiTypeTreeScale mttsOperator = new MultiTypeTreeScale();
    mttsOperator.initByName("weight", 1.0, "scaleFactor", 0.8, "useOldTreeScaler", true, "migrationModel", migModel, "multiTypeTree", mtTree);
    TypePairBirthDeath tpbdOperator = new TypePairBirthDeath();
    tpbdOperator.initByName("weight", 1.0, "migrationModel", migModel, "multiTypeTree", mtTree);
    TypeMergeSplit tmsOperator = new TypeMergeSplit();
    tmsOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "includeRoot", true);
    // Set up stat analysis logger:
    MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
    logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.1, "logEvery", 1000);
    // Set up MCMC:
    MCMC mcmc = new MCMC();
    mcmc.initByName("chainLength", "10000000", "state", state, "distribution", distribution, "operator", twbOperator, "operator", tsxOperator, "operator", mtuOperator, "operator", mttsOperator, "operator", tpbdOperator, "operator", tmsOperator, "logger", logger);
    // Run MCMC:
    mcmc.run();
    System.out.format("height mean = %s\n", logger.getHeightMean());
    System.out.format("height var = %s\n", logger.getHeightVar());
    System.out.format("height ESS = %s\n", logger.getHeightESS());
    // Compare analysis results with truth:
    boolean withinTol = (logger.getHeightESS() > 1000) && (Math.abs(logger.getHeightMean() - 23) < 0.5) && (Math.abs(logger.getHeightVar() - 300) < 30.0);
    Assert.assertTrue(withinTol);
}

Example 15 with State

use of beast.core.State in project MultiTypeTree by tgvaughan.

the class Ewing_Test method test1.

@Test
public void test1() throws Exception {
    System.out.println("Ewing_test1");
    // Fix seed.
    Randomizer.setSeed(1);
    // Assemble initial MultiTypeTree
    String newickStr = "((1[&deme=0]:1,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
    MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
    mtTree.initByName("value", newickStr, "typeLabel", "deme");
    // Assemble migration model:
    RealParameter rateMatrix = new RealParameter("0.1 0.1");
    RealParameter popSizes = new RealParameter("7.0 7.0");
    SCMigrationModel migModel = new SCMigrationModel();
    migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
    // Assemble distribution:
    StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
    distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree);
    // Set up state:
    State state = new State();
    state.initByName("stateNode", mtTree);
    // Set up operators:
    TypedWilsonBaldingEasy twbOperator = new TypedWilsonBaldingEasy();
    twbOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
    TypedSubtreeExchangeEasy tsxOperator = new TypedSubtreeExchangeEasy();
    tsxOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "isNarrow", true);
    MultiTypeUniform mtuOperator = new MultiTypeUniform();
    mtuOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
    MultiTypeTreeScale mttsOperator = new MultiTypeTreeScale();
    mttsOperator.initByName("weight", 1.0, "scaleFactor", 0.8, "useOldTreeScaler", true, "multiTypeTree", mtTree, "migrationModel", migModel);
    TypePairBirthDeath tpbdOperator = new TypePairBirthDeath();
    tpbdOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel);
    TypeMergeSplit tmsOperator = new TypeMergeSplit();
    tmsOperator.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "includeRoot", true);
    // Set up stat analysis logger:
    MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
    logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.1, "logEvery", 1000);
    // Logger fileLogger = new Logger();
    // TreeHeightLogger thLogger = new TreeHeightLogger();
    // thLogger.initByName("tree", mtTree);
    // fileLogger.initByName(
    // "fileName", "test.log",
    // "logEvery", 1000,
    // "model", distribution,
    // "log", thLogger);
    // Set up MCMC:
    MCMC mcmc = new MCMC();
    mcmc.initByName("chainLength", "10000000", "state", state, "distribution", distribution, "operator", twbOperator, "operator", tsxOperator, "operator", mtuOperator, "operator", mttsOperator, "operator", tpbdOperator, "operator", tmsOperator, "logger", logger);
    // Run MCMC:
    mcmc.run();
    System.out.format("height mean = %s\n", logger.getHeightMean());
    System.out.format("height var = %s\n", logger.getHeightVar());
    System.out.format("height ESS = %s\n", logger.getHeightESS());
    // Compare analysis results with truth:
    boolean withinTol = (logger.getHeightESS() > 3000) && (Math.abs(logger.getHeightMean() - 19.15) < 0.5) && (Math.abs(logger.getHeightVar() - 310) < 20);
    Assert.assertTrue(withinTol);
}