Examples with SiteModel beast.evolution.sitemodel.SiteModel used on opensource projects

Example 16 with SiteModel

use of beast.evolution.sitemodel.SiteModel in project beast2 by CompEvol.

the class TreeLikelihoodTest method testGTRGLikelihood.

@Test
public void testGTRGLikelihood() throws Exception {
    // Set up GTR model: 4 gamma categories, gamma shape = 0.5
    Alignment data = BEASTTestCase.getAlignment();
    Tree tree = BEASTTestCase.getTree(data);
    Frequencies freqs = new Frequencies();
    freqs.initByName("data", data);
    GeneralSubstitutionModel gsm = new GeneralSubstitutionModel();
    gsm.initByName("rates", "1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0", "frequencies", freqs);
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", "1.0", "gammaCategoryCount", 4, "shape", "0.5", "substModel", gsm);
    TreeLikelihood likelihood = newTreeLikelihood();
    likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
    double logP = 0;
    logP = likelihood.calculateLogP();
    assertEquals(logP, -1949.0360143622, BEASTTestCase.PRECISION);
    likelihood.initByName("useAmbiguities", false, "data", data, "tree", tree, "siteModel", siteModel);
    logP = likelihood.calculateLogP();
    assertEquals(logP, -1949.0360143622, BEASTTestCase.PRECISION);
}

Example 17 with SiteModel

use of beast.evolution.sitemodel.SiteModel in project beast2 by CompEvol.

the class ClusterTreeTest method testUPGMA.

public void testUPGMA() throws Exception {
    Alignment alignment = BEASTTestCase.getAlignment();
    JukesCantor JC = new JukesCantor();
    JC.initAndValidate();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("substModel", JC);
    ClusterTree tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    String expectedNewick = "((((human:0.01903085702575253,(chimp:0.008560512208575313,bonobo:0.008560512208575313):0.010470344817177218):0.007962255880644985,gorilla:0.026993112906397516):0.019197419394211015,orangutan:0.04619053230060853):0.007214240662738673,siamang:0.053404772963347204):0.0";
    String actualNewick = tree.getRoot().toNewick();
    assertEquals(expectedNewick, actualNewick);
    // select 3 sequences
    List<Sequence> seqs = new ArrayList<Sequence>();
    seqs.addAll(alignment.sequenceInput.get());
    List<Sequence> newseqs = alignment.sequenceInput.get();
    newseqs.clear();
    newseqs.add(Sequence.getSequenceByTaxon("bonobo", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("chimp", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("human", seqs));
    alignment.initAndValidate();
    tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    expectedNewick = "((bonobo:0.008560512208575313,chimp:0.008560512208575313):0.010470344817177218,human:0.01903085702575253):0.0";
    System.err.println("Seqs:");
    for (Sequence s : seqs) {
        System.err.println(s.taxonInput.get());
    }
    System.err.println("Newseqs:");
    for (Sequence s : newseqs) {
        System.err.println(s.taxonInput.get());
    }
    actualNewick = tree.getRoot().toNewick();
    assertEquals(expectedNewick, actualNewick);
    // same sequences in different order
    newseqs.clear();
    newseqs.add(Sequence.getSequenceByTaxon("human", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("chimp", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("bonobo", seqs));
    alignment.initAndValidate();
    tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    actualNewick = tree.getRoot().toNewick();
    expectedNewick = "(human:0.01903085702575253,(chimp:0.008560512208575313,bonobo:0.008560512208575313):0.010470344817177218):0.0";
    assertEquals(expectedNewick, actualNewick);
}

Example 18 with SiteModel

use of beast.evolution.sitemodel.SiteModel in project beast2 by CompEvol.

the class SiteModelInputEditor method customConnector.

public static boolean customConnector(BeautiDoc doc) {
    try {
        DeltaExchangeOperator operator = (DeltaExchangeOperator) doc.pluginmap.get("FixMeanMutationRatesOperator");
        if (operator == null) {
            return false;
        }
        List<RealParameter> parameters = operator.parameterInput.get();
        parameters.clear();
        // String weights = "";
        CompoundDistribution likelihood = (CompoundDistribution) doc.pluginmap.get("likelihood");
        boolean hasOneEstimatedRate = false;
        List<String> rateIDs = new ArrayList<>();
        List<Integer> weights = new ArrayList<>();
        for (Distribution d : likelihood.pDistributions.get()) {
            GenericTreeLikelihood treelikelihood = (GenericTreeLikelihood) d;
            Alignment data = treelikelihood.dataInput.get();
            int weight = data.getSiteCount();
            if (data.isAscertained) {
                weight -= data.getExcludedPatternCount();
            }
            if (treelikelihood.siteModelInput.get() instanceof SiteModel) {
                SiteModel siteModel = (SiteModel) treelikelihood.siteModelInput.get();
                RealParameter mutationRate = siteModel.muParameterInput.get();
                // clockRate.m_bIsEstimated.setValue(true, clockRate);
                if (mutationRate.isEstimatedInput.get()) {
                    hasOneEstimatedRate = true;
                    if (rateIDs.indexOf(mutationRate.getID()) == -1) {
                        parameters.add(mutationRate);
                        weights.add(weight);
                        rateIDs.add(mutationRate.getID());
                    } else {
                        int k = rateIDs.indexOf(mutationRate.getID());
                        weights.set(k, weights.get(k) + weight);
                    }
                }
            }
        }
        IntegerParameter weightParameter;
        if (weights.size() == 0) {
            weightParameter = new IntegerParameter();
        } else {
            String weightString = "";
            for (int k : weights) {
                weightString += k + " ";
            }
            weightParameter = new IntegerParameter(weightString);
            weightParameter.setID("weightparameter");
        }
        weightParameter.isEstimatedInput.setValue(false, weightParameter);
        operator.parameterWeightsInput.setValue(weightParameter, operator);
        return hasOneEstimatedRate;
    } catch (Exception e) {
    }
    return false;
}

Example 19 with SiteModel

use of beast.evolution.sitemodel.SiteModel 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 20 with SiteModel

use of beast.evolution.sitemodel.SiteModel in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodAmbiguities.

@Test
public void testLikelihoodAmbiguities() throws Exception {
    Locus locus = new Locus("locus", getAlignment());
    // ConversionGraph
    ConversionGraph acg = new ConversionGraph();
    ClusterTree tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", locus.getAlignment());
    acg.assignFrom(tree);
    acg.initByName("locus", locus);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("substModel", jc);
    // Likelihood
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "tree", acg, "siteModel", siteModel, "useAmbiguities", true);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel, "useAmbiguities", true);
    acg.setEverythingDirty(true);
    double logP = argLikelihood.calculateLogP();
    double logPtrue = argLikelihoodSlow.calculateLogP();
    double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    System.out.println("Test 1.  Truth: " + logPtrue + " Value: " + logP);
    assertTrue(relativeDiff < 1e-14);
    // 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 = 100;
    int endLocus = 200;
    Conversion recomb1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb1);
    logP = argLikelihood.calculateLogP();
    logPtrue = argLikelihoodSlow.calculateLogP();
    relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    assertTrue(relativeDiff < 1e-14);
    // Add another recombination event
    node1 = acg.getExternalNodes().get(0);
    node2 = acg.getNode(20);
    height1 = 0.75 * (node1.getHeight() + node1.getParent().getHeight());
    height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
    startLocus = 250;
    endLocus = 300;
    Conversion recomb2 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb2);
    logP = argLikelihood.calculateLogP();
    logPtrue = argLikelihoodSlow.calculateLogP();
    relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    assertTrue(relativeDiff < 1e-14);
}