Search in sources :

Example 1 with SimpleSiteList

use of dr.evolution.alignment.SimpleSiteList in project beast-mcmc by beast-dev.

the class BASTADensityTester method main.

public static void main(String[] args) {
    MathUtils.setSeed(666);
    // turn off logging to avoid screen noise...
    Logger logger = Logger.getLogger("dr");
    logger.setUseParentHandlers(false);
    System.out.println("EXAMPLE 1: 4 taxa with 2 demes");
    SimpleAlignment alignment = createAlignment(sequences, Nucleotides.INSTANCE);
    TreeModel treeModel;
    try {
        treeModel = createSpecifiedTree("(((AB192965Japan2004:0.42343888910376215,AF189018Indonesia2005:1.4234388891037622)5:9.725099517053918,AF071228Spain1997:3.14853840615768)6:1.7275747971782511,AF105975Portugal1995:2.8761132033359313)");
    } catch (Exception e) {
        throw new RuntimeException("Unable to parse Newick tree");
    }
    List<String> states = new ArrayList<String>();
    states.add("Asia");
    states.add("West_Medit");
    GeneralDataType dataType = new GeneralDataType(states);
    Taxa taxonList = new Taxa();
    Taxon taxonOne = new Taxon(sequences[0][0]);
    taxonOne.setAttribute(DEME, states.get(0));
    taxonList.addTaxon(taxonOne);
    Taxon taxonTwo = new Taxon(sequences[0][1]);
    taxonTwo.setAttribute(DEME, states.get(1));
    taxonList.addTaxon(taxonTwo);
    Taxon taxonThree = new Taxon(sequences[0][2]);
    taxonThree.setAttribute(DEME, states.get(1));
    taxonList.addTaxon(taxonThree);
    Taxon taxonFour = new Taxon(sequences[0][3]);
    taxonFour.setAttribute(DEME, states.get(0));
    taxonList.addTaxon(taxonFour);
    int[] pattern = new int[sequences[0].length];
    pattern[0] = 0;
    pattern[1] = 1;
    pattern[2] = 1;
    pattern[3] = 0;
    SimpleSiteList patterns = new SimpleSiteList(dataType, taxonList);
    patterns.addPattern(pattern);
    double[] freqValues = new double[2];
    freqValues[0] = 0.5;
    freqValues[1] = 0.5;
    FrequencyModel freqModel = new FrequencyModel(dataType, freqValues);
    Parameter ratesParameter = new Parameter.Default("migration rates", 2);
    ratesParameter.setParameterValue(0, 1.0);
    ratesParameter.setParameterValue(1, 1.0);
    Parameter popSizesParameter = new Parameter.Default("population sizes", 2);
    popSizesParameter.setParameterValue(0, 1.0);
    popSizesParameter.setParameterValue(1, 1.0);
    // GeneralSubstitutionModel migrationModel = new GeneralSubstitutionModel("migrationModel", dataType, null, ratesParameter, 0);
    SVSComplexSubstitutionModel migrationModel = new SVSComplexSubstitutionModel("migrationModel", dataType, freqModel, ratesParameter, null);
    TaxonList includeSubtree = null;
    List<TaxonList> excludeSubtrees = new ArrayList<TaxonList>();
    int subIntervals = 2;
    BranchRateModel branchRateModel = new DefaultBranchRateModel();
    StructuredCoalescentLikelihood structured = null;
    try {
        structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
    } catch (TreeUtils.MissingTaxonException missing) {
        System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
    }
    System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
    System.out.println("EXAMPLE 2: 5 taxa (2 identical sampling dates) with 2 demes");
    try {
        treeModel = createSpecifiedTree("((((AB192965Japan2004:0.48972300366304955,AJ842306France2004:0.48972300366304955):1.2052625264106087,AF189018Indonesia2005:2.6949855300736583):6.86618775860274,AF071228Spain1997:1.5611732886763985):1.4313107562768952,AF105975Portugal1995:0.9924840449532937)");
    } catch (Exception e) {
        throw new RuntimeException("Unable to parse Newick tree");
    }
    Taxon taxonFive = new Taxon(sequencesTwo[0][4]);
    taxonFive.setAttribute(DEME, states.get(1));
    taxonList.addTaxon(taxonFive);
    pattern = new int[sequencesTwo[0].length];
    pattern[0] = 0;
    pattern[1] = 1;
    pattern[2] = 1;
    pattern[3] = 0;
    pattern[4] = 1;
    patterns = new SimpleSiteList(dataType, taxonList);
    patterns.addPattern(pattern);
    try {
        structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
    } catch (TreeUtils.MissingTaxonException missing) {
        System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
    }
    System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
    System.out.println("EXAMPLE 3: 6 taxa (2 identical sampling dates) with 3 demes");
    try {
        treeModel = createSpecifiedTree("(((AB192965Japan2004:4.509764060335954,((AF189018Indonesia2005:1.496797606742874,AJ842306France2004:0.4967976067428741):0.3617175987505199,FM212660Cameroon2007:3.858515205493394):3.6512488548425597):4.957323993357177,AF071228Spain1997:2.4670880536931303):0.9760399779205056,AF105975Portugal1995:1.443128031613636)");
    } catch (Exception e) {
        throw new RuntimeException("Unable to parse Newick tree");
    }
    states = new ArrayList<String>();
    states.add("Asia");
    states.add("West_Medit");
    states.add("Africa");
    dataType = new GeneralDataType(states);
    taxonOne.setAttribute(DEME, states.get(0));
    taxonTwo.setAttribute(DEME, states.get(1));
    taxonThree.setAttribute(DEME, states.get(1));
    taxonFour.setAttribute(DEME, states.get(0));
    taxonFive.setAttribute(DEME, states.get(1));
    Taxon taxonSix = new Taxon(sequencesThree[0][5]);
    taxonSix.setAttribute(DEME, states.get(2));
    taxonList.addTaxon(taxonSix);
    pattern = new int[sequencesThree[0].length];
    pattern[0] = 0;
    pattern[1] = 1;
    pattern[2] = 1;
    pattern[3] = 0;
    pattern[4] = 1;
    pattern[5] = 2;
    patterns = new SimpleSiteList(dataType, taxonList);
    patterns.addPattern(pattern);
    freqValues = new double[3];
    freqValues[0] = 1.0 / 3.0;
    freqValues[1] = 1.0 / 3.0;
    freqValues[2] = 1 - freqValues[0] - freqValues[1];
    freqModel = new FrequencyModel(dataType, freqValues);
    ratesParameter = new Parameter.Default("migration rates", 6);
    ratesParameter.setParameterValue(0, 1.0);
    ratesParameter.setParameterValue(1, 1.0);
    ratesParameter.setParameterValue(2, 1.0);
    ratesParameter.setParameterValue(3, 1.0);
    ratesParameter.setParameterValue(4, 1.0);
    ratesParameter.setParameterValue(5, 1.0);
    popSizesParameter = new Parameter.Default("population sizes", 3);
    popSizesParameter.setParameterValue(0, 1.0);
    popSizesParameter.setParameterValue(1, 1.0);
    popSizesParameter.setParameterValue(2, 1.0);
    // GeneralSubstitutionModel migrationModel = new GeneralSubstitutionModel("migrationModel", dataType, null, ratesParameter, 0);
    migrationModel = new SVSComplexSubstitutionModel("migrationModel", dataType, freqModel, ratesParameter, null);
    try {
        structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
    } catch (TreeUtils.MissingTaxonException missing) {
        System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
    }
    System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
    System.out.println("EXAMPLE 4: 6 taxa (2 identical sampling dates) with 3 demes and simple branch lengths");
    try {
        treeModel = createSpecifiedTree("(((AB192965Japan2004:4.5,((AF189018Indonesia2005:1.5,AJ842306France2004:0.5):0.25,FM212660Cameroon2007:3.75):3.75):5.0,AF071228Spain1997:2.5):1.0,AF105975Portugal1995:1.5)");
    } catch (Exception e) {
        throw new RuntimeException("Unable to parse Newick tree");
    }
    try {
        structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
    } catch (TreeUtils.MissingTaxonException missing) {
        System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
    }
    System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
    System.out.println("EXAMPLE 5: 6 taxa (2 identical sampling dates) with 3 demes and simple (but different than in EXAMPLE 4) branch lengths");
    try {
        treeModel = createSpecifiedTree("(((((AF189018Indonesia2005:1.5,AJ842306France2004:0.5):0.75,AB192965Japan2004:1.25):3.25,FM212660Cameroon2007:7.5):5.0,AF071228Spain1997:2.5):1.0,AF105975Portugal1995:1.5)");
    } catch (Exception e) {
        throw new RuntimeException("Unable to parse Newick tree");
    }
    try {
        structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
    } catch (TreeUtils.MissingTaxonException missing) {
        System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
    }
    System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) ArrayList(java.util.ArrayList) Logger(java.util.logging.Logger) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) TreeModel(dr.evomodel.tree.TreeModel) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) SimpleAlignment(dr.evolution.alignment.SimpleAlignment) TreeUtils(dr.evolution.tree.TreeUtils) SimpleSiteList(dr.evolution.alignment.SimpleSiteList) GeneralDataType(dr.evolution.datatype.GeneralDataType) SVSComplexSubstitutionModel(dr.evomodel.substmodel.SVSComplexSubstitutionModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) Parameter(dr.inference.model.Parameter)

Example 2 with SimpleSiteList

use of dr.evolution.alignment.SimpleSiteList in project beast-mcmc by beast-dev.

the class AttributePatternsParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    String attributeName = xo.getStringAttribute(ATTRIBUTE);
    String secondaryAttributeName = xo.getAttribute(SECONDARY_ATTRIBUTE, (String) null);
    TaxonList taxa = (TaxonList) xo.getChild(TaxonList.class);
    DataType dataType = DataTypeUtils.getDataType(xo);
    if (dataType == null) {
        throw new XMLParseException("dataType expected for attributePatterns element");
    }
    // using a SimpleSiteList rather than Patterns to allow ancestral reconstruction
    SimpleSiteList patterns = new SimpleSiteList(dataType, taxa);
    int[] pattern = new int[taxa.getTaxonCount()];
    boolean attributeFound = true;
    for (int i = 0; i < taxa.getTaxonCount(); i++) {
        Taxon taxon = taxa.getTaxon(i);
        if (secondaryAttributeName == null || secondaryAttributeName.isEmpty()) {
            Object value = taxon.getAttribute(attributeName);
            if (value != null) {
                int state = dataType.getState(value.toString());
                if (state < 0) {
                    throw new XMLParseException("State for attribute, " + attributeName + ", in taxon, " + taxon.getId() + ", is unknown: " + value.toString());
                }
                pattern[i] = state;
            } else {
                attributeFound = false;
            }
        } else {
            Object value1 = taxon.getAttribute(attributeName);
            Object value2 = taxon.getAttribute(secondaryAttributeName);
            if (value1 != null && value2 != null) {
                String code = value1.toString() + CompositeDataTypeParser.COMPOSITE_STATE_SEPARATOR + value2.toString();
                int state = dataType.getState(code);
                if (state < 0) {
                    throw new XMLParseException("State for attributes, " + attributeName + " & " + secondaryAttributeName + ", in taxon, " + taxon.getId() + ", is unknown: " + code);
                }
                pattern[i] = state;
            } else {
                attributeFound = false;
            }
        }
    }
    if (!attributeFound) {
        throw new XMLParseException("The attribute, " + attributeName + " was not found in all taxa. Check the name of the attribute.");
    }
    patterns.addPattern(pattern);
    if (xo.hasAttribute(XMLParser.ID)) {
        Logger.getLogger("dr.evoxml").info("Read attribute patterns, '" + xo.getId() + "' for attribute, " + attributeName);
    } else {
        Logger.getLogger("dr.evoxml").info("Read attribute patterns for attribute, " + attributeName);
    }
    return patterns;
}
Also used : TaxonList(dr.evolution.util.TaxonList) Taxon(dr.evolution.util.Taxon) DataType(dr.evolution.datatype.DataType) SimpleSiteList(dr.evolution.alignment.SimpleSiteList)

Aggregations

SimpleSiteList (dr.evolution.alignment.SimpleSiteList)2 SimpleAlignment (dr.evolution.alignment.SimpleAlignment)1 DataType (dr.evolution.datatype.DataType)1 GeneralDataType (dr.evolution.datatype.GeneralDataType)1 TreeUtils (dr.evolution.tree.TreeUtils)1 Taxon (dr.evolution.util.Taxon)1 TaxonList (dr.evolution.util.TaxonList)1 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)1 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)1 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)1 SVSComplexSubstitutionModel (dr.evomodel.substmodel.SVSComplexSubstitutionModel)1 DefaultTreeModel (dr.evomodel.tree.DefaultTreeModel)1 TreeModel (dr.evomodel.tree.TreeModel)1 Parameter (dr.inference.model.Parameter)1 ArrayList (java.util.ArrayList)1 Logger (java.util.logging.Logger)1