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

Example 26 with SiteModel

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

the class UnorderedAlignmentsTest method testUnorderedAlignment.

@Test
public void testUnorderedAlignment() throws Exception {
    TaxonSet taxa = getTaxa();
    Tree tree = getTree(taxa);
    SiteModel siteModel = getSiteModel();
    double logP = 0.0;
    double shuffledLogP = 0.0;
    for (int i = 0; i < 3; ++i) {
        Alignment data = getAlignment(taxa, tree, siteModel);
        // First calculate in order
        TreeLikelihood likelihood = new TreeLikelihood();
        likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
        logP += likelihood.calculateLogP();
        // Now calculate again, with shuffled taxon order
        Collections.shuffle(data.sequenceInput.get());
        likelihood = new TreeLikelihood();
        likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
        shuffledLogP += likelihood.calculateLogP();
    }
    assertEquals(logP, shuffledLogP, 1E-9);
}

Example 27 with SiteModel

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

the class TreeLikelihood method initAndValidate.

@Override
public void initAndValidate() {
    // sanity check: alignment should have same #taxa as tree
    if (dataInput.get().getTaxonCount() != treeInput.get().getLeafNodeCount()) {
        throw new IllegalArgumentException("The number of nodes in the tree does not match the number of sequences");
    }
    beagle = null;
    beagle = new BeagleTreeLikelihood();
    try {
        beagle.initByName("data", dataInput.get(), "tree", treeInput.get(), "siteModel", siteModelInput.get(), "branchRateModel", branchRateModelInput.get(), "useAmbiguities", m_useAmbiguities.get(), "useTipLikelihoods", m_useTipLikelihoods.get(), "scaling", scaling.get().toString());
        if (beagle.beagle != null) {
            // a Beagle instance was found, so we use it
            return;
        }
    } catch (Exception e) {
    // ignore
    }
    // No Beagle instance was found, so we use the good old java likelihood core
    beagle = null;
    int nodeCount = treeInput.get().getNodeCount();
    if (!(siteModelInput.get() instanceof SiteModel.Base)) {
        throw new IllegalArgumentException("siteModel input should be of type SiteModel.Base");
    }
    m_siteModel = (SiteModel.Base) siteModelInput.get();
    m_siteModel.setDataType(dataInput.get().getDataType());
    substitutionModel = m_siteModel.substModelInput.get();
    if (branchRateModelInput.get() != null) {
        branchRateModel = branchRateModelInput.get();
    } else {
        branchRateModel = new StrictClockModel();
    }
    m_branchLengths = new double[nodeCount];
    storedBranchLengths = new double[nodeCount];
    int stateCount = dataInput.get().getMaxStateCount();
    int patterns = dataInput.get().getPatternCount();
    if (stateCount == 4) {
        likelihoodCore = new BeerLikelihoodCore4();
    } else {
        likelihoodCore = new BeerLikelihoodCore(stateCount);
    }
    String className = getClass().getSimpleName();
    Alignment alignment = dataInput.get();
    Log.info.println(className + "(" + getID() + ") uses " + likelihoodCore.getClass().getSimpleName());
    Log.info.println("  " + alignment.toString(true));
    // print startup messages via Log.print*
    proportionInvariant = m_siteModel.getProportionInvariant();
    m_siteModel.setPropInvariantIsCategory(false);
    if (proportionInvariant > 0) {
        calcConstantPatternIndices(patterns, stateCount);
    }
    initCore();
    patternLogLikelihoods = new double[patterns];
    m_fRootPartials = new double[patterns * stateCount];
    matrixSize = (stateCount + 1) * (stateCount + 1);
    probabilities = new double[(stateCount + 1) * (stateCount + 1)];
    Arrays.fill(probabilities, 1.0);
    if (dataInput.get().isAscertained) {
        useAscertainedSitePatterns = true;
    }
}

Example 28 with SiteModel

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

the class SequenceSimulator method main.

// printUsageAndExit
@SuppressWarnings("unchecked")
public static void main(String[] args) {
    try {
        // parse arguments
        if (args.length < 2) {
            printUsageAndExit();
        }
        String fileName = args[0];
        int replications = Integer.parseInt(args[1]);
        PrintStream out = System.out;
        if (args.length == 3) {
            File file = new File(args[2]);
            out = new PrintStream(file);
        }
        // grab the file
        String xml = "";
        BufferedReader fin = new BufferedReader(new FileReader(fileName));
        while (fin.ready()) {
            xml += fin.readLine();
        }
        fin.close();
        // parse the xml
        XMLParser parser = new XMLParser();
        BEASTInterface beastObject = parser.parseFragment(xml, true);
        // find relevant objects from the model
        TreeLikelihood treeLikelihood = getTreeLikelihood(beastObject);
        if (treeLikelihood == null) {
            throw new IllegalArgumentException("No treelikelihood found in file. Giving up now.");
        }
        Alignment data = ((Input<Alignment>) treeLikelihood.getInput("data")).get();
        Tree tree = ((Input<Tree>) treeLikelihood.getInput("tree")).get();
        SiteModel pSiteModel = ((Input<SiteModel>) treeLikelihood.getInput("siteModel")).get();
        BranchRateModel pBranchRateModel = ((Input<BranchRateModel>) treeLikelihood.getInput("branchRateModel")).get();
        // feed to sequence simulator and generate leaves
        SequenceSimulator treeSimulator = new SequenceSimulator();
        treeSimulator.init(data, tree, pSiteModel, pBranchRateModel, replications);
        XMLProducer producer = new XMLProducer();
        Alignment alignment = treeSimulator.simulate();
        xml = producer.toRawXML(alignment);
        out.println("<beast version='2.0'>");
        out.println(xml);
        out.println("</beast>");
    } catch (Exception e) {
        e.printStackTrace();
    }
}

Example 29 with SiteModel

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

the class SiteModelInputEditor method createMutationRateEditor.

public InputEditor createMutationRateEditor() {
    SiteModel sitemodel = ((SiteModel) m_input.get());
    final Input<?> input = sitemodel.muParameterInput;
    ParameterInputEditor mutationRateEditor = new ParameterInputEditor(doc);
    mutationRateEditor.init(input, sitemodel, -1, ExpandOption.FALSE, true);
    mutationRateEditor.getEntry().setEnabled(!doc.autoUpdateFixMeanSubstRate);
    return mutationRateEditor;
}

Example 30 with SiteModel

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

the class SiteModelInputEditor method setUpOperator.

/**
 * set up relative weights and parameter input *
 */
public void setUpOperator() {
    boolean isAllClocksAreEqual = true;
    try {
        boolean hasOneEstimatedRate = customConnector(doc);
        if (doc.autoUpdateFixMeanSubstRate) {
            fixMeanRatesCheckBox.setSelected(hasOneEstimatedRate);
            doFixMeanRates(hasOneEstimatedRate);
        }
        try {
            double commonClockRate = -1;
            CompoundDistribution likelihood = (CompoundDistribution) doc.pluginmap.get("likelihood");
            for (Distribution d : likelihood.pDistributions.get()) {
                GenericTreeLikelihood treelikelihood = (GenericTreeLikelihood) d;
                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()) {
                        if (commonClockRate < 0) {
                            commonClockRate = mutationRate.valuesInput.get().get(0);
                        } else {
                            if (Math.abs(commonClockRate - mutationRate.valuesInput.get().get(0)) > 1e-10) {
                                isAllClocksAreEqual = false;
                            }
                        }
                    }
                }
            }
        } catch (Exception e) {
        }
        List<RealParameter> parameters = operator.parameterInput.get();
        if (!fixMeanRatesCheckBox.isSelected()) {
            fixMeanRatesValidateLabel.setVisible(false);
            repaint();
            return;
        }
        if (parameters.size() == 0) {
            fixMeanRatesValidateLabel.setVisible(true);
            fixMeanRatesValidateLabel.m_circleColor = Color.red;
            fixMeanRatesValidateLabel.setToolTipText("The model is invalid: At least one substitution rate should be estimated.");
            repaint();
            return;
        }
        if (!isAllClocksAreEqual) {
            fixMeanRatesValidateLabel.setVisible(true);
            fixMeanRatesValidateLabel.m_circleColor = Color.orange;
            fixMeanRatesValidateLabel.setToolTipText("Not all substitution rates are equal. Are you sure this is what you want?");
        } else if (parameters.size() == 1) {
            fixMeanRatesValidateLabel.setVisible(true);
            fixMeanRatesValidateLabel.m_circleColor = Color.orange;
            fixMeanRatesValidateLabel.setToolTipText("At least 2 clock models should have their rate estimated");
        } else if (parameters.size() < doc.getPartitions("SiteModel").size()) {
            fixMeanRatesValidateLabel.setVisible(true);
            fixMeanRatesValidateLabel.m_circleColor = Color.orange;
            fixMeanRatesValidateLabel.setToolTipText("Not all partitions have their rate estimated");
        } else {
            fixMeanRatesValidateLabel.setVisible(false);
        }
        repaint();
    } catch (Exception e) {
        e.printStackTrace();
    }
}