Search in sources :

Example 31 with RealParameter

use of beast.core.parameter.RealParameter 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);
}
Also used : Operator(beast.core.Operator) MultiTypeTreeStatLogger(multitypetree.util.MultiTypeTreeStatLogger) MCMC(beast.core.MCMC) RealParameter(beast.core.parameter.RealParameter) SCMigrationModel(beast.evolution.tree.SCMigrationModel) State(beast.core.State) MultiTypeTreeFromNewick(beast.evolution.tree.MultiTypeTreeFromNewick) TypeSet(beast.evolution.tree.TypeSet) StructuredCoalescentTreeDensity(multitypetree.distributions.StructuredCoalescentTreeDensity) Test(org.junit.Test)

Example 32 with RealParameter

use of beast.core.parameter.RealParameter 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();
        }
    }
}
Also used : PrintStream(java.io.PrintStream) State(beast.core.State) RealParameter(beast.core.parameter.RealParameter)

Example 33 with RealParameter

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

the class InitMigrationModelConnector method customConnector.

public static boolean customConnector(BeautiDoc doc) {
    for (BEASTInterface p : doc.getPartitions("Tree")) {
        TreeLikelihood treeLikelihood = (TreeLikelihood) p;
        StructuredCoalescentMultiTypeTree tree = (StructuredCoalescentMultiTypeTree) treeLikelihood.treeInput.get();
        String pID = BeautiDoc.parsePartition(tree.getID());
        SCMigrationModel migModel = (SCMigrationModel) doc.pluginmap.get("migModel.t:" + pID);
        SCMigrationModel migModelInit = (SCMigrationModel) doc.pluginmap.get("migModelInit.t:" + pID);
        String rateMatrixStr = getParameterString((RealParameter) migModel.rateMatrixInput.get());
        String popSizesStr = getParameterString((RealParameter) migModel.popSizesInput.get());
        // Ensure model has appropriate number of demes
        int uniqueTraitCount = uniqueTraitsInData(tree).size();
        StringBuilder rateMatrixStrBuilder = new StringBuilder();
        StringBuilder popSizesStrBuilder = new StringBuilder();
        migModel.getTypeSet().initAndValidate();
        if (migModel.popSizesInput.get().getDimension() != migModel.getNTypes()) {
            for (int i = 0; i < migModel.getNTypes(); i++) {
                popSizesStrBuilder.append(" 1.0");
                for (int j = 0; j < migModel.getNTypes(); j++) {
                    if (j == i)
                        continue;
                    rateMatrixStrBuilder.append(" 1.0");
                }
            }
            popSizesStr = popSizesStrBuilder.toString();
            rateMatrixStr = rateMatrixStrBuilder.toString();
            ((RealParameter) migModel.popSizesInput.get()).setDimension(migModel.getNTypes());
            ((RealParameter) migModel.popSizesInput.get()).valuesInput.setValue(popSizesStr, (RealParameter) migModel.popSizesInput.get());
            ((RealParameter) migModel.rateMatrixInput.get()).setDimension(migModel.getNTypes() * (migModel.getNTypes() - 1));
            ((RealParameter) migModel.rateMatrixInput.get()).valuesInput.setValue(rateMatrixStr, (RealParameter) migModel.rateMatrixInput.get());
            ((RealParameter) migModel.popSizesInput.get()).initAndValidate();
            ((RealParameter) migModel.rateMatrixInput.get()).initAndValidate();
            migModel.initAndValidate();
        }
        ((RealParameter) migModelInit.popSizesInput.get()).setDimension(migModel.getNTypes());
        ((RealParameter) migModelInit.popSizesInput.get()).valuesInput.setValue(popSizesStr, (RealParameter) migModelInit.popSizesInput.get());
        ((RealParameter) migModelInit.rateMatrixInput.get()).setDimension(migModel.getNTypes() * (migModel.getNTypes() - 1));
        ((RealParameter) migModelInit.rateMatrixInput.get()).valuesInput.setValue(rateMatrixStr, (RealParameter) migModelInit.rateMatrixInput.get());
        try {
            ((RealParameter) migModelInit.popSizesInput.get()).initAndValidate();
            ((RealParameter) migModelInit.rateMatrixInput.get()).initAndValidate();
            migModelInit.initAndValidate();
        } catch (Exception ex) {
            System.err.println("Error configuring initial migration model.");
        }
    }
    return false;
}
Also used : TreeLikelihood(beast.evolution.likelihood.TreeLikelihood) RealParameter(beast.core.parameter.RealParameter) BEASTInterface(beast.core.BEASTInterface) SCMigrationModel(beast.evolution.tree.SCMigrationModel) StructuredCoalescentMultiTypeTree(beast.evolution.tree.StructuredCoalescentMultiTypeTree)

Example 34 with RealParameter

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

the class StructuredCoalescentUntypedTree method main.

/**
 * Generates an ensemble of trees from the structured coalescent for testing
 * coloured tree-space samplers.
 *
 * @param argv
 * @throws Exception
 */
public static void main(String[] argv) throws Exception {
    // Set up migration model.
    RealParameter rateMatrix = new RealParameter();
    rateMatrix.initByName("value", "0.05", "dimension", "12");
    RealParameter popSizes = new RealParameter();
    popSizes.initByName("value", "7.0", "dimension", "4");
    SCMigrationModel migrationModel = new SCMigrationModel();
    migrationModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes);
    // Specify leaf types:
    IntegerParameter leafTypes = new IntegerParameter();
    leafTypes.initByName("value", "0 0 0");
    // Generate ensemble:
    int reps = 1000000;
    double[] heights = new double[reps];
    long startTime = System.currentTimeMillis();
    StructuredCoalescentUntypedTree sctree;
    sctree = new StructuredCoalescentUntypedTree();
    for (int i = 0; i < reps; i++) {
        if (i % 1000 == 0)
            System.out.format("%d reps done\n", i);
        sctree.initByName("migrationModel", migrationModel, "leafTypes", leafTypes, "nTypes", 4);
        heights[i] = sctree.getRoot().getHeight();
    }
    long time = System.currentTimeMillis() - startTime;
    System.out.printf("E[T] = %1.4f +/- %1.4f\n", DiscreteStatistics.mean(heights), DiscreteStatistics.stdev(heights) / Math.sqrt(reps));
    System.out.printf("V[T] = %1.4f\n", DiscreteStatistics.variance(heights));
    System.out.printf("Took %1.2f seconds\n", time / 1000.0);
    try (PrintStream outStream = new PrintStream("heights.txt")) {
        outStream.println("h c");
        for (int i = 0; i < reps; i++) outStream.format("%g\n", heights[i]);
    }
}
Also used : IntegerParameter(beast.core.parameter.IntegerParameter) PrintStream(java.io.PrintStream) RealParameter(beast.core.parameter.RealParameter)

Example 35 with RealParameter

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

the class StructuredCoalescentTreeDensity method main.

/**
 * Test likelihood result. Duplicate of JUnit test for debugging purposes.
 *
 * @param argv
 * @throws java.lang.Exception
 */
public static void main(String[] argv) throws Exception {
    // Assemble test MultiTypeTree:
    String newickStr = "(((A[state=1]:0.25)[state=0]:0.25,B[state=0]:0.5)[state=0]:1.5," + "(C[state=0]:1.0,D[state=0]:1.0)[state=0]:1.0)[state=0]:0.0;";
    MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
    mtTree.initByName("newick", newickStr, "typeLabel", "state", "nTypes", 2);
    // Assemble migration model:
    RealParameter rateMatrix = new RealParameter();
    rateMatrix.initByName("value", "2.0 1.0");
    RealParameter popSizes = new RealParameter();
    popSizes.initByName("value", "5.0 10.0");
    SCMigrationModel migrationModel = new SCMigrationModel();
    migrationModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes);
    // Set up likelihood instance:
    StructuredCoalescentTreeDensity likelihood = new StructuredCoalescentTreeDensity();
    likelihood.initByName("migrationModel", migrationModel, "multiTypeTree", mtTree);
    // Calculated by hand
    double expResult = -16.52831;
    double result = likelihood.calculateLogP();
    System.out.println("Expected result: " + expResult);
    System.out.println("Actual result: " + result);
    System.out.println("Difference: " + String.valueOf(result - expResult));
}
Also used : MultiTypeTreeFromNewick(beast.evolution.tree.MultiTypeTreeFromNewick) RealParameter(beast.core.parameter.RealParameter) SCMigrationModel(beast.evolution.tree.SCMigrationModel)

Aggregations

RealParameter (beast.core.parameter.RealParameter)97 Test (org.junit.Test)39 IntegerParameter (beast.core.parameter.IntegerParameter)23 Tree (beast.evolution.tree.Tree)16 State (beast.core.State)14 SCMigrationModel (beast.evolution.tree.SCMigrationModel)13 Alignment (beast.evolution.alignment.Alignment)11 TypeSet (beast.evolution.tree.TypeSet)11 MCMC (beast.core.MCMC)10 SiteModel (beast.evolution.sitemodel.SiteModel)10 Frequencies (beast.evolution.substitutionmodel.Frequencies)10 ConstantPopulation (beast.evolution.tree.coalescent.ConstantPopulation)10 StructuredCoalescentTreeDensity (multitypetree.distributions.StructuredCoalescentTreeDensity)10 TaxonSet (beast.evolution.alignment.TaxonSet)9 MultiTypeTreeStatLogger (multitypetree.util.MultiTypeTreeStatLogger)9 MultiTypeTreeFromNewick (beast.evolution.tree.MultiTypeTreeFromNewick)8 Node (beast.evolution.tree.Node)8 Operator (beast.core.Operator)7 RandomTree (beast.evolution.tree.RandomTree)7 Locus (bacter.Locus)6