Search in sources :

Example 1 with Patterns

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

the class CheckPointTreeModifier method incorporateAdditionalTaxa.

/**
     * Add the remaining taxa, which can be identified through the TreeDataLikelihood XML elements.
     */
public ArrayList<NodeRef> incorporateAdditionalTaxa(CheckPointUpdaterApp.UpdateChoice choice, BranchRates rateModel) {
    System.out.println("Tree before adding taxa:\n" + treeModel.toString() + "\n");
    ArrayList<NodeRef> newTaxaNodes = new ArrayList<NodeRef>();
    for (String str : newTaxaNames) {
        for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
            if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId().equals(str)) {
                newTaxaNodes.add(treeModel.getExternalNode(i));
                //always take into account Taxon dates vs. dates set through a TreeModel
                System.out.println(treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " with height " + treeModel.getNodeHeight(treeModel.getExternalNode(i)) + " or " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getHeight());
            }
        }
    }
    System.out.println("newTaxaNodes length = " + newTaxaNodes.size());
    ArrayList<Taxon> currentTaxa = new ArrayList<Taxon>();
    for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
        boolean taxonFound = false;
        for (String str : newTaxaNames) {
            if (str.equals((treeModel.getNodeTaxon(treeModel.getExternalNode(i))).getId())) {
                taxonFound = true;
            }
        }
        if (!taxonFound) {
            System.out.println("Adding " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " to list of current taxa");
            currentTaxa.add(treeModel.getNodeTaxon(treeModel.getExternalNode(i)));
        }
    }
    System.out.println("Current taxa count = " + currentTaxa.size());
    //iterate over both current taxa and to be added taxa
    boolean originTaxon = true;
    for (Taxon taxon : currentTaxa) {
        if (taxon.getHeight() == 0.0) {
            originTaxon = false;
            System.out.println("Current taxon " + taxon.getId() + " has node height 0.0");
        }
    }
    for (NodeRef newTaxon : newTaxaNodes) {
        if (treeModel.getNodeTaxon(newTaxon).getHeight() == 0.0) {
            originTaxon = false;
            System.out.println("New taxon " + treeModel.getNodeTaxon(newTaxon).getId() + " has node height 0.0");
        }
    }
    //check the Tree(Data)Likelihoods in the connected set of likelihoods
    //focus on TreeDataLikelihood, which has getTree() to get the tree for each likelihood
    //also get the DataLikelihoodDelegate from TreeDataLikelihood
    ArrayList<TreeDataLikelihood> likelihoods = new ArrayList<TreeDataLikelihood>();
    ArrayList<Tree> trees = new ArrayList<Tree>();
    ArrayList<DataLikelihoodDelegate> delegates = new ArrayList<DataLikelihoodDelegate>();
    for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
        if (likelihood instanceof TreeDataLikelihood) {
            likelihoods.add((TreeDataLikelihood) likelihood);
            trees.add(((TreeDataLikelihood) likelihood).getTree());
            delegates.add(((TreeDataLikelihood) likelihood).getDataLikelihoodDelegate());
        }
    }
    //suggested to go through TreeDataLikelihoodParser and give it an extra option to create a HashMap
    //keyed by the tree; am currently not overly fond of this approach
    ArrayList<PatternList> patternLists = new ArrayList<PatternList>();
    for (DataLikelihoodDelegate del : delegates) {
        if (del instanceof BeagleDataLikelihoodDelegate) {
            patternLists.add(((BeagleDataLikelihoodDelegate) del).getPatternList());
        } else if (del instanceof MultiPartitionDataLikelihoodDelegate) {
            MultiPartitionDataLikelihoodDelegate mpdld = (MultiPartitionDataLikelihoodDelegate) del;
            List<PatternList> list = mpdld.getPatternLists();
            for (PatternList pList : list) {
                patternLists.add(pList);
            }
        }
    }
    if (patternLists.size() == 0) {
        throw new RuntimeException("No patterns detected. Please make sure the XML file is BEAST 1.9 compatible.");
    }
    //aggregate all patterns to create distance matrix
    //TODO What about different trees for different partitions?
    Patterns patterns = new Patterns(patternLists.get(0));
    if (patternLists.size() > 1) {
        for (int i = 1; i < patternLists.size(); i++) {
            patterns.addPatterns(patternLists.get(i));
        }
    }
    //set the patterns for the distance matrix computations
    choice.setPatterns(patterns);
    //add new taxa one at a time
    System.out.println("Adding " + newTaxaNodes.size() + " taxa ...");
    for (NodeRef newTaxon : newTaxaNodes) {
        treeModel.setNodeHeight(newTaxon, treeModel.getNodeTaxon(newTaxon).getHeight());
        System.out.println("\nadding Taxon: " + newTaxon + " (height = " + treeModel.getNodeHeight(newTaxon) + ")");
        //check if this taxon has a more recent sampling date than all other nodes in the current TreeModel
        double offset = checkCurrentTreeNodes(newTaxon, treeModel.getRoot());
        System.out.println("Sampling date offset when adding " + newTaxon + " = " + offset);
        //AND set its current node height to 0.0 IF no originTaxon has been found
        if (offset < 0.0) {
            if (!originTaxon) {
                System.out.println("Updating all node heights with offset " + Math.abs(offset));
                updateAllTreeNodes(Math.abs(offset), treeModel.getRoot());
                treeModel.setNodeHeight(newTaxon, 0.0);
            }
        } else if (offset == 0.0) {
            if (!originTaxon) {
                treeModel.setNodeHeight(newTaxon, 0.0);
            }
        }
        //get the closest Taxon to the Taxon that needs to be added
        //take into account which taxa can currently be chosen
        Taxon closest = choice.getClosestTaxon(treeModel.getNodeTaxon(newTaxon), currentTaxa);
        System.out.println("\nclosest Taxon: " + closest + " with original height: " + closest.getHeight());
        //get the distance between these two taxa
        double distance = choice.getDistance(treeModel.getNodeTaxon(newTaxon), closest);
        System.out.println("at distance: " + distance);
        //TODO what if distance == 0.0 ??? how to choose closest taxon then (in absence of geo info)?
        //find the NodeRef for the closest Taxon (do not rely on node numbering)
        NodeRef closestRef = null;
        //careful with node numbering and subtract number of new taxa
        for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
            if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)) == closest) {
                closestRef = treeModel.getExternalNode(i);
            }
        }
        System.out.println(closestRef + " with height " + treeModel.getNodeHeight(closestRef));
        //System.out.println("trying to set node height: " + closestRef + " from " + treeModel.getNodeHeight(closestRef) + " to " + closest.getHeight());
        //treeModel.setNodeHeight(closestRef, closest.getHeight());
        double timeForDistance = distance / rateModel.getBranchRate(treeModel, closestRef);
        System.out.println("timeForDistance = " + timeForDistance);
        //get parent node of branch that will be split
        NodeRef parent = treeModel.getParent(closestRef);
        //determine height of new node
        double insertHeight;
        if (treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon)) {
            insertHeight = treeModel.getNodeHeight(closestRef) + timeForDistance / 2.0;
            System.out.println("treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon): " + insertHeight);
            if (insertHeight >= treeModel.getNodeHeight(parent)) {
                insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
            }
        } else {
            double remainder = (timeForDistance - Math.abs(treeModel.getNodeHeight(closestRef) - treeModel.getNodeHeight(newTaxon))) / 2.0;
            if (remainder > 0) {
                insertHeight = Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)) + remainder;
                System.out.println("remainder > 0: " + insertHeight);
                if (insertHeight >= treeModel.getNodeHeight(parent)) {
                    insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
                }
            } else {
                insertHeight = EPSILON * (treeModel.getNodeHeight(parent) - Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)));
                insertHeight += Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon));
                System.out.println("remainder <= 0: " + insertHeight);
            }
        }
        System.out.println("insert at height: " + insertHeight);
        //pass on all the necessary variables to a method that adds the new taxon to the tree
        addTaxonAlongBranch(newTaxon, parent, closestRef, insertHeight);
        //option to print tree after each taxon addition
        System.out.println("\nTree after adding taxon " + newTaxon + ":\n" + treeModel.toString());
        //add newly added Taxon to list of current taxa
        currentTaxa.add(treeModel.getNodeTaxon(newTaxon));
    }
    return newTaxaNodes;
}
Also used : Likelihood(dr.inference.model.Likelihood) TreeDataLikelihood(dr.evomodel.treedatalikelihood.TreeDataLikelihood) Taxon(dr.evolution.util.Taxon) ArrayList(java.util.ArrayList) PatternList(dr.evolution.alignment.PatternList) BeagleDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.BeagleDataLikelihoodDelegate) NodeRef(dr.evolution.tree.NodeRef) TreeDataLikelihood(dr.evomodel.treedatalikelihood.TreeDataLikelihood) BeagleDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.BeagleDataLikelihoodDelegate) MultiPartitionDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.MultiPartitionDataLikelihoodDelegate) DataLikelihoodDelegate(dr.evomodel.treedatalikelihood.DataLikelihoodDelegate) Tree(dr.evolution.tree.Tree) PatternList(dr.evolution.alignment.PatternList) ArrayList(java.util.ArrayList) List(java.util.List) MultiPartitionDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.MultiPartitionDataLikelihoodDelegate) Patterns(dr.evolution.alignment.Patterns)

Example 2 with Patterns

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

the class OldTreesPanel method createTree.

private void createTree() {
    if (generateTreeDialog == null) {
        generateTreeDialog = new GenerateTreeDialog(frame);
    }
    int result = generateTreeDialog.showDialog(options);
    if (result != JOptionPane.CANCEL_OPTION) {
        GenerateTreeDialog.MethodTypes methodType = generateTreeDialog.getMethodType();
        PartitionData partition = generateTreeDialog.getDataPartition();
        Patterns patterns = new Patterns(partition.getAlignment());
        DistanceMatrix distances = new F84DistanceMatrix(patterns);
        Tree tree;
        TemporalRooting temporalRooting;
        switch(methodType) {
            case NJ:
                tree = new NeighborJoiningTree(distances);
                temporalRooting = new TemporalRooting(tree);
                tree = temporalRooting.findRoot(tree, TemporalRooting.RootingFunction.CORRELATION);
                break;
            case UPGMA:
                tree = new UPGMATree(distances);
                temporalRooting = new TemporalRooting(tree);
                break;
            default:
                throw new IllegalArgumentException("unknown method type");
        }
        tree.setId(generateTreeDialog.getName());
        options.userTrees.add(tree);
        treesTableModel.fireTableDataChanged();
        int row = options.userTrees.size() - 1;
        treesTable.getSelectionModel().setSelectionInterval(row, row);
    }
    fireTreePriorsChanged();
}
Also used : F84DistanceMatrix(dr.evolution.distance.F84DistanceMatrix) NeighborJoiningTree(dr.evolution.tree.NeighborJoiningTree) PartitionData(dr.app.beauti.options.PartitionData) UPGMATree(dr.evolution.tree.UPGMATree) NeighborJoiningTree(dr.evolution.tree.NeighborJoiningTree) Tree(dr.evolution.tree.Tree) TemporalRooting(dr.app.pathogen.TemporalRooting) F84DistanceMatrix(dr.evolution.distance.F84DistanceMatrix) DistanceMatrix(dr.evolution.distance.DistanceMatrix) UPGMATree(dr.evolution.tree.UPGMATree) Patterns(dr.evolution.alignment.Patterns)

Example 3 with Patterns

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

the class BalancedBeagleTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(BeagleTreeLikelihoodParser.USE_AMBIGUITIES, false);
    /*int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
        if (instanceCount < 1) {
            instanceCount = 1;
        }

        String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
        if (ic != null && ic.length() > 0) {
            instanceCount = Integer.parseInt(ic);
        }*/
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    GammaSiteRateModel siteRateModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        if (substitutionModel == null) {
            substitutionModel = siteRateModel.getSubstitutionModel();
        }
        if (substitutionModel == null) {
            throw new XMLParseException("No substitution model available for TreeLikelihood: " + xo.getId());
        }
        branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    // if (xo.getChild(TipStatesModel.class) != null) {
    // throw new XMLParseException("Sequence Error Models are not supported under BEAGLE yet. Please use Native BEAST Likelihood.");
    // }
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    if (xo.hasAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME)) {
        // scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    boolean delayScaling = true;
    Map<Set<String>, Parameter> partialsRestrictions = null;
    if (xo.hasChildNamed(PARTIALS_RESTRICTION)) {
        XMLObject cxo = xo.getChild(PARTIALS_RESTRICTION);
        TaxonList taxonList = (TaxonList) cxo.getChild(TaxonList.class);
        // Parameter parameter = (Parameter) cxo.getChild(Parameter.class);
        try {
            TreeUtils.getLeavesForTaxa(treeModel, taxonList);
        } catch (TreeUtils.MissingTaxonException e) {
            throw new XMLParseException("Unable to parse taxon list: " + e.getMessage());
        }
        throw new XMLParseException("Restricting internal nodes is not yet implemented.  Contact Marc");
    }
    /*if (instanceCount == 1 || patternList.getPatternCount() < instanceCount) {
            return createTreeLikelihood(
                    patternList,
                    treeModel,
                    branchModel,
                    siteRateModel,
                    branchRateModel,
                    tipStatesModel,
                    useAmbiguities,
                    scalingScheme,
                    partialsRestrictions,
                    xo
            );
        }*/
    // first run a test for instanceCount == 1
    System.err.println("\nTesting instanceCount == 1");
    Likelihood baseLikelihood = createTreeLikelihood(patternList, treeModel, branchModel, siteRateModel, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
    double start = System.nanoTime();
    for (int i = 0; i < TEST_RUNS; i++) {
        baseLikelihood.makeDirty();
        baseLikelihood.getLogLikelihood();
    }
    double end = System.nanoTime();
    double baseResult = end - start;
    System.err.println("Evaluation took: " + baseResult);
    if (!(patternList instanceof SitePatterns)) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with BEAUti-selected codon partitioning.");
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    // List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    List<Likelihood> likelihoods = null;
    CompoundLikelihood compound = null;
    int instanceCount = 2;
    boolean optimal = false;
    while (optimal == false) {
        System.err.println("\nCreating instanceCount == " + instanceCount);
        likelihoods = new ArrayList<Likelihood>();
        for (int i = 0; i < instanceCount; i++) {
            Patterns subPatterns = new Patterns((SitePatterns) patternList, 0, 0, 1, i, instanceCount);
            AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
            treeLikelihood.setId(xo.getId() + "_" + instanceCount);
            likelihoods.add(treeLikelihood);
        }
        // construct compoundLikelihood
        compound = new CompoundLikelihood(instanceCount, likelihoods);
        // test timings
        System.err.println("\nTesting instanceCount == " + instanceCount);
        start = System.nanoTime();
        for (int i = 0; i < TEST_RUNS; i++) {
            compound.makeDirty();
            compound.getLogLikelihood();
        }
        end = System.nanoTime();
        double newResult = end - start;
        System.err.println("Evaluation took: " + newResult);
        if (baseResult / newResult > TEST_CUTOFF) {
            instanceCount++;
            baseResult = newResult;
        } else {
            optimal = true;
            instanceCount--;
            System.err.println("\nCreating final BeagleTreeLikelihood with instanceCount: " + instanceCount);
            likelihoods = new ArrayList<Likelihood>();
            for (int i = 0; i < instanceCount; i++) {
                Patterns subPatterns = new Patterns((SitePatterns) patternList, 0, 0, 1, i, instanceCount);
                AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
                treeLikelihood.setId(xo.getId() + "_" + instanceCount);
                likelihoods.add(treeLikelihood);
            }
            // construct compoundLikelihood
            compound = new CompoundLikelihood(instanceCount, likelihoods);
        }
    }
    return compound;
/*for (int i = 0; i < instanceCount; i++) {

            Patterns subPatterns = new Patterns((SitePatterns)patternList, 0, 0, 1, i, instanceCount);

            AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(
                    subPatterns,
                    treeModel,
                    branchModel,
                    siteRateModel,
                    branchRateModel,
                    null,
                    useAmbiguities,
                    scalingScheme,
                    partialsRestrictions,
                    xo);
            treeLikelihood.setId(xo.getId() + "_" + instanceCount);
            likelihoods.add(treeLikelihood);
        }

        return new CompoundLikelihood(likelihoods);*/
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Set(java.util.Set) CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) BeagleTreeLikelihood(dr.evomodel.treelikelihood.BeagleTreeLikelihood) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PatternList(dr.evolution.alignment.PatternList) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PartialsRescalingScheme(dr.evomodel.treelikelihood.PartialsRescalingScheme) BranchModel(dr.evomodel.branchmodel.BranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) TreeModel(dr.evomodel.tree.TreeModel) Patterns(dr.evolution.alignment.Patterns) SitePatterns(dr.evolution.alignment.SitePatterns) TreeUtils(dr.evolution.tree.TreeUtils) SitePatterns(dr.evolution.alignment.SitePatterns) TaxonList(dr.evolution.util.TaxonList) CompoundLikelihood(dr.inference.model.CompoundLikelihood) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) Parameter(dr.inference.model.Parameter)

Example 4 with Patterns

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

the class BeagleTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(USE_AMBIGUITIES, false);
    int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
    if (instanceCount < 1) {
        instanceCount = 1;
    }
    String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
    if (ic != null && ic.length() > 0) {
        instanceCount = Integer.parseInt(ic);
    }
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    MutableTreeModel treeModel = (MutableTreeModel) xo.getChild(MutableTreeModel.class);
    GammaSiteRateModel siteRateModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        if (substitutionModel == null) {
            substitutionModel = siteRateModel.getSubstitutionModel();
        }
        if (substitutionModel == null) {
            throw new XMLParseException("No substitution model available for TreeLikelihood: " + xo.getId());
        }
        branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (branchModel instanceof EpochBranchModel && rootFreqModel != null) {
        EpochBranchModel epochBranchModel = (EpochBranchModel) branchModel;
        epochBranchModel.setRootFrequencyModel(rootFreqModel);
    }
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    // if (xo.getChild(TipStatesModel.class) != null) {
    // throw new XMLParseException("Sequence Error Models are not supported under BEAGLE yet. Please use Native BEAST Likelihood.");
    // }
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    boolean delayScaling = true;
    if (xo.hasAttribute(SCALING_SCHEME)) {
        scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    if (xo.hasAttribute(DELAY_SCALING)) {
        delayScaling = xo.getBooleanAttribute(DELAY_SCALING);
    }
    Map<Set<String>, Parameter> partialsRestrictions = null;
    if (xo.hasChildNamed(PARTIALS_RESTRICTION)) {
        XMLObject cxo = xo.getChild(PARTIALS_RESTRICTION);
        TaxonList taxonList = (TaxonList) cxo.getChild(TaxonList.class);
        // Parameter parameter = (Parameter) cxo.getChild(Parameter.class);
        try {
            TreeUtils.getLeavesForTaxa(treeModel, taxonList);
        } catch (TreeUtils.MissingTaxonException e) {
            throw new XMLParseException("Unable to parse taxon list: " + e.getMessage());
        }
        throw new XMLParseException("Restricting internal nodes is not yet implemented.  Contact Marc");
    }
    int beagleThreadCount = -1;
    if (System.getProperty(BEAGLE_THREAD_COUNT) != null) {
        beagleThreadCount = Integer.parseInt(System.getProperty(BEAGLE_THREAD_COUNT));
    }
    if (beagleThreadCount == -1) {
        // the default is -1 threads (automatic thread pool size) but an XML attribute can override it
        int threadCount = xo.getAttribute(THREADS, -1);
        if (System.getProperty(THREAD_COUNT) != null) {
            threadCount = Integer.parseInt(System.getProperty(THREAD_COUNT));
        }
        // Todo: allow for different number of threads per beagle instance according to pattern counts
        if (threadCount >= 0) {
            System.setProperty(BEAGLE_THREAD_COUNT, Integer.toString(threadCount / instanceCount));
        }
    }
    if (instanceCount == 1 || patternList.getPatternCount() < instanceCount) {
        return createTreeLikelihood(patternList, treeModel, branchModel, siteRateModel, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    for (int i = 0; i < instanceCount; i++) {
        Patterns subPatterns = new Patterns(patternList, i, instanceCount);
        AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
        treeLikelihood.setId(xo.getId() + "_" + instanceCount);
        likelihoods.add(treeLikelihood);
    }
    return new CompoundLikelihood(likelihoods);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Set(java.util.Set) CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) BeagleTreeLikelihood(dr.evomodel.treelikelihood.BeagleTreeLikelihood) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PatternList(dr.evolution.alignment.PatternList) ArrayList(java.util.ArrayList) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PartialsRescalingScheme(dr.evomodel.treelikelihood.PartialsRescalingScheme) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) BranchModel(dr.evomodel.branchmodel.BranchModel) EpochBranchModel(dr.evomodel.branchmodel.EpochBranchModel) MutableTreeModel(dr.evolution.tree.MutableTreeModel) Patterns(dr.evolution.alignment.Patterns) TreeUtils(dr.evolution.tree.TreeUtils) EpochBranchModel(dr.evomodel.branchmodel.EpochBranchModel) TaxonList(dr.evolution.util.TaxonList) CompoundLikelihood(dr.inference.model.CompoundLikelihood) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) Parameter(dr.inference.model.Parameter)

Example 5 with Patterns

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

the class MultiPartitionDataLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(USE_AMBIGUITIES, false);
    int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
    if (instanceCount < 1) {
        instanceCount = 1;
    }
    String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
    if (ic != null && ic.length() > 0) {
        instanceCount = Integer.parseInt(ic);
    }
    if (DEBUG) {
        System.out.println("instanceCount: " + instanceCount);
    }
    List<PatternList> patternLists = xo.getAllChildren(PatternList.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    List<SiteRateModel> siteRateModels = xo.getAllChildren(SiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    List<BranchModel> branchModels = xo.getAllChildren(BranchModel.class);
    if (branchModels.size() == 0) {
        if (DEBUG) {
            System.out.println("branchModels == null");
        }
        branchModels = new ArrayList<BranchModel>();
        List<SubstitutionModel> substitutionModels = xo.getAllChildren(SubstitutionModel.class);
        if (substitutionModels.size() == 0) {
            if (DEBUG) {
                System.out.println("substitutionModels == null");
            }
            for (SiteRateModel siteRateModel : siteRateModels) {
                SubstitutionModel substitutionModel = ((GammaSiteRateModel) siteRateModel).getSubstitutionModel();
                if (substitutionModel == null) {
                    throw new XMLParseException("No substitution model available for TreeDataLikelihood: " + xo.getId());
                }
                branchModels.add(new HomogeneousBranchModel(substitutionModel, rootFreqModel));
            }
        }
        if (DEBUG) {
            System.out.println("branchModels size: " + branchModels.size());
        }
        for (BranchModel branchModel : branchModels) {
            System.out.println("  " + branchModel.getId() + "  " + branchModel.getModelName());
        }
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (branchRateModel == null) {
        branchRateModel = new DefaultBranchRateModel();
    }
    if (DEBUG) {
        System.out.println("BranchRateModel: " + branchRateModel.getId());
    }
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    boolean delayScaling = true;
    if (xo.hasAttribute(SCALING_SCHEME)) {
        scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    if (xo.hasAttribute(DELAY_SCALING)) {
        delayScaling = xo.getBooleanAttribute(DELAY_SCALING);
    }
    if (instanceCount == 1) {
        if (DEBUG) {
            System.out.println("instanceCount == 1");
        }
        return createTreeDataLikelihood(patternLists, treeModel, branchModels, siteRateModels, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, xo);
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    List<PatternList> patternInstanceLists = new ArrayList<PatternList>();
    for (int j = 0; j < patternLists.size(); j++) {
        for (int i = 0; i < instanceCount; i++) {
            patternInstanceLists.add(new Patterns(patternLists.get(j), i, instanceCount));
        }
    }
    return createTreeDataLikelihood(patternLists, treeModel, branchModels, siteRateModels, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, xo);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) PatternList(dr.evolution.alignment.PatternList) ArrayList(java.util.ArrayList) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) PartialsRescalingScheme(dr.evomodel.treelikelihood.PartialsRescalingScheme) BranchModel(dr.evomodel.branchmodel.BranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) SiteRateModel(dr.evomodel.siteratemodel.SiteRateModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) TreeModel(dr.evomodel.tree.TreeModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) Patterns(dr.evolution.alignment.Patterns)

Aggregations

Patterns (dr.evolution.alignment.Patterns)22 Taxon (dr.evolution.util.Taxon)10 ArrayList (java.util.ArrayList)10 Microsatellite (dr.evolution.datatype.Microsatellite)9 Taxa (dr.evolution.util.Taxa)9 Tree (dr.evolution.tree.Tree)8 TreeModel (dr.evomodel.tree.TreeModel)8 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)6 PatternList (dr.evolution.alignment.PatternList)5 NewickImporter (dr.evolution.io.NewickImporter)4 BranchModel (dr.evomodel.branchmodel.BranchModel)4 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)4 DefaultTreeModel (dr.evomodel.tree.DefaultTreeModel)4 PartialsRescalingScheme (dr.evomodel.treelikelihood.PartialsRescalingScheme)4 Likelihood (dr.inference.model.Likelihood)4 Parameter (dr.inference.model.Parameter)4 GammaSiteModel (dr.oldevomodel.sitemodel.GammaSiteModel)4 AsymmetricQuadraticModel (dr.oldevomodel.substmodel.AsymmetricQuadraticModel)4 SitePatterns (dr.evolution.alignment.SitePatterns)3 TaxonList (dr.evolution.util.TaxonList)3