Examples with TreeModel dr.evomodel.tree.TreeModel used on opensource projects

Example 11 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class ImportancePruneAndRegraftTestProblem method testDoOperation.

/**
	 * Test method for {@link SimpleMCMCOperator#doOperation()}.
	 * @throws ImportException 
	 * @throws IOException 
	 */
public void testDoOperation() throws IOException, ImportException {
    // probability of picking (A,B) node is 1/(2n-3) = 1/7
    // probability of swapping with D is 1/2
    // total = 1/14
    //probability of picking {D} node is 1/(2n-3) = 1/7
    //probability of picking {A,B} is 1/5
    // total = 1/35
    //total = 1/14 + 1/35 = 7/70 = 0.1
    System.out.println("Test 1: Forward");
    String treeMatch = "(((D,C),(A,B)),E);";
    int count = 0;
    int reps = 1000000;
    for (int i = 0; i < reps; i++) {
        TreeModel treeModel = new TreeModel("treeModel", tree5);
        ImportancePruneAndRegraft operator = new ImportancePruneAndRegraft(treeModel, 1.0, 0);
        operator.doOperation();
        String tree = TreeUtils.newickNoLengths(treeModel);
        if (tree.equals(treeMatch)) {
            count += 1;
        }
    }
    double p_1 = (double) count / (double) reps;
    System.out.println("Number of proposals:\t" + count);
    System.out.println("Number of tries:\t" + reps);
    System.out.println("Number of ratio:\t" + p_1);
    System.out.println("Number of expected ratio:\t" + 0.1);
    assertExpectation(0.1, p_1, reps);
    // lets see what the backward probability is for the hastings ratio
    // (((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0) -> ((((A,B),C),D),E)
    // probability of picking (A,B) node is 1/(2n-3) = 1/7
    // probability of swapping with D is 1/3
    // total = 1/21
    //probability of picking {D} node is 1/(2n-2) = 1/7
    //probability of picking {A,B} is 1/4
    // total = 1/28
    //total = 1/21 + 1/28 = 7/84 = 0.08333333
    System.out.println("Test 2: Backward");
    treeMatch = "((((A,B),C),D),E);";
    NewickImporter importer = new NewickImporter("(((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0);");
    FlexibleTree tree5_2 = (FlexibleTree) importer.importTree(null);
    count = 0;
    for (int i = 0; i < reps; i++) {
        TreeModel treeModel = new TreeModel("treeModel", tree5_2);
        ImportancePruneAndRegraft operator = new ImportancePruneAndRegraft(treeModel, 1.0, 1);
        operator.doOperation();
        String tree = TreeUtils.newickNoLengths(treeModel);
        if (tree.equals(treeMatch)) {
            count += 1;
        }
    }
    double p_2 = (double) count / (double) reps;
    System.out.println("Number of proposals:\t" + count);
    System.out.println("Number of tries:\t" + reps);
    System.out.println("Number of ratio:\t" + p_2);
    System.out.println("Number of expected ratio:\t" + 0.0833333);
    assertExpectation(0.0833333, p_2, reps);
}

Example 12 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class ImportanceSubtreeSwapTestProblem method testDoOperation.

/**
	 * Test method for {@link SimpleMCMCOperator#doOperation()}.
	 * @throws ImportException 
	 * @throws IOException 
	 */
public void testDoOperation() throws IOException, ImportException {
    // probability of picking (A,B) node is 1/(2n-3) = 1/7
    // probability of swapping with D is 1/2
    // total = 1/14
    //probability of picking {D} node is 1/(2n-3) = 1/7
    //probability of picking {A,B} is 1/5
    // total = 1/35
    //total = 1/14 + 1/35 = 7/70 = 0.1
    System.out.println("Test 1: Forward");
    String treeMatch = "(((D,C),(A,B)),E);";
    int count = 0;
    int reps = 100000;
    for (int i = 0; i < reps; i++) {
        TreeModel treeModel = new TreeModel("treeModel", tree5);
        ImportanceSubtreeSwap operator = new ImportanceSubtreeSwap(treeModel, 1.0, 0);
        operator.doOperation();
        String tree = TreeUtils.newickNoLengths(treeModel);
        if (tree.equals(treeMatch)) {
            count += 1;
        }
    }
    double p_1 = (double) count / (double) reps;
    System.out.println("Number of proposals:\t" + count);
    System.out.println("Number of tries:\t" + reps);
    System.out.println("Number of ratio:\t" + p_1);
    System.out.println("Number of expected ratio:\t" + 0.1);
    assertExpectation(0.1, p_1, reps);
    // lets see what the backward probability is for the hastings ratio
    // (((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0) -> ((((A,B),C),D),E)
    // probability of picking (A,B) node is 1/(2n-3) = 1/7
    // probability of swapping with D is 1/3
    // total = 1/21
    //probability of picking {D} node is 1/(2n-2) = 1/7
    //probability of picking {A,B} is 1/4
    // total = 1/28
    //total = 1/21 + 1/28 = 7/84 = 0.08333333
    System.out.println("Test 2: Backward");
    treeMatch = "((((A,B),C),D),E);";
    NewickImporter importer = new NewickImporter("(((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0);");
    FlexibleTree tree5_2 = (FlexibleTree) importer.importTree(null);
    count = 0;
    for (int i = 0; i < reps; i++) {
        TreeModel treeModel = new TreeModel("treeModel", tree5_2);
        ImportanceSubtreeSwap operator = new ImportanceSubtreeSwap(treeModel, 1.0, 0);
        operator.doOperation();
        String tree = TreeUtils.newickNoLengths(treeModel);
        if (tree.equals(treeMatch)) {
            count += 1;
        }
    }
    double p_2 = (double) count / (double) reps;
    System.out.println("Number of proposals:\t" + count);
    System.out.println("Number of tries:\t" + reps);
    System.out.println("Number of ratio:\t" + p_2);
    System.out.println("Number of expected ratio:\t" + 0.0833333);
    assertExpectation(0.0833333, p_2, reps);
}

Example 13 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateTwoPartitions.

// END: simulateOnePartition
static void simulateTwoPartitions() {
    try {
        System.out.println("Test case 3: simulateTwoPartitions");
        MathUtils.setSeed(666);
        int sequenceLength = 11;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        3, // every
        1);
        // create partition
        Partition Partition = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        4, // to
        sequenceLength - 1, // every
        1);
        Sequence ancestralSequence = new Sequence();
        ancestralSequence.appendSequenceString("TCAAGTG");
        Partition.setRootSequence(ancestralSequence);
        partitionsList.add(partition1);
        partitionsList.add(Partition);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        System.out.println(simulator.simulate(simulateInPar, false).toString());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 14 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateAminoAcid.

// END: simulateThreePartitions
static void simulateAminoAcid() {
    try {
        System.out.println("Test case 4: simulateAminoAcid");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05 });
        FrequencyModel freqModel = new FrequencyModel(AminoAcids.INSTANCE, freqs);
        // create substitution model
        EmpiricalRateMatrix rateMatrix = Blosum62.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        System.out.println(simulator.simulate(simulateInPar, false).toString());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch
}

Example 15 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateRandomBranchAssignment.

// END: main
static void simulateRandomBranchAssignment() {
    try {
        System.out.println("Test case I dunno which: simulate random branch assignments");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        File treeFile = new File("/home/filip/Dropbox/BeagleSequenceSimulator/SimTree/SimTree.figtree");
        Tree tree = Utils.importTreeFromFile(treeFile);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(Utils.UNIFORM_CODON_FREQUENCIES);
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create base subst model
        Parameter omegaParameter = new Parameter.Default("omega", 1, 1.0);
        Parameter kappaParameter = new Parameter.Default("kappa", 1, 1.0);
        GY94CodonModel baseSubModel = new GY94CodonModel(Codons.UNIVERSAL, omegaParameter, kappaParameter, freqModel);
        RandomBranchModel substitutionModel = new RandomBranchModel(treeModel, baseSubModel, 0.25, false, -1);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        //			Sequence ancestralSequence = new Sequence();
        //			ancestralSequence.appendSequenceString("TCAAGTGAGG");
        //			partition1.setRootSequence(ancestralSequence);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        SimpleAlignment alignment = simulator.simulate(simulateInPar, true);
        // alignment.setOutputType(SimpleAlignment.OutputType.NEXUS);
        alignment.setOutputType(SimpleAlignment.OutputType.XML);
        System.out.println(alignment.toString());
    } catch (Exception e) {
        e.printStackTrace();
    }
// END: try-catch
}