Examples with MutableTreeModel dr.evolution.tree.MutableTreeModel used on opensource projects

Example 6 with MutableTreeModel

use of dr.evolution.tree.MutableTreeModel in project beast-mcmc by beast-dev.

the class MultivariateTraitUtils method computeTreeTraitMeanOU.

public static double[] computeTreeTraitMeanOU(FullyConjugateMultivariateTraitLikelihood trait, double[] rootValue, boolean conditionOnRoot) {
    double[] root = trait.getPriorMean();
    MutableTreeModel myTreeModel = trait.getTreeModel();
    double[][] linCombMatrix = computeLinCombMatrix(trait);
    double[] rootMultiplierVect = computeRootMultipliers(trait);
    if (conditionOnRoot) {
        root = rootValue;
    }
    final int nTaxa = myTreeModel.getExternalNodeCount();
    final int branchCount = 2 * nTaxa - 2;
    final int traitDim = trait.dimTrait;
    double[] mean = new double[root.length * nTaxa];
    double[] displacementMeans = new double[branchCount * traitDim];
    double[] tempVect = new double[nTaxa * traitDim];
    NodeRef tempNode;
    for (int k = 0; k < branchCount; k++) {
        tempNode = myTreeModel.getNode(k);
        for (int t = 0; t < traitDim; t++) {
            displacementMeans[k * traitDim + t] = (1 - Math.exp(-trait.getTimeScaledSelection(tempNode))) * trait.getOptimalValue(tempNode)[t];
        }
    }
    // multiply linCombMatrix with displacement means
    for (int i = 0; i < nTaxa; i++) {
        for (int j = 0; j < branchCount; j++) {
            for (int k = 0; k < traitDim; k++) {
                tempVect[i * traitDim + k] = tempVect[i * traitDim + k] + linCombMatrix[i][j] * displacementMeans[j * traitDim + k];
            }
        }
    }
    for (int i = 0; i < nTaxa; ++i) {
        System.arraycopy(root, 0, mean, i * root.length, root.length);
        for (int j = 0; j < traitDim; ++j) {
            mean[i * traitDim + j] = mean[i * traitDim + j] * rootMultiplierVect[i] + tempVect[i * traitDim + j];
        }
    }
    return mean;
}

Example 7 with MutableTreeModel

use of dr.evolution.tree.MutableTreeModel in project beast-mcmc by beast-dev.

the class MultivariateTraitUtils method getShiftContributionToMean.

private static double[] getShiftContributionToMean(NodeRef node, FullyConjugateMultivariateTraitLikelihood trait) {
    MutableTreeModel treeModel = trait.getTreeModel();
    double[] shiftContribution = new double[trait.dimTrait];
    if (!treeModel.isRoot(node)) {
        NodeRef parent = treeModel.getParent(node);
        double[] shiftContributionParent = getShiftContributionToMean(parent, trait);
        for (int i = 0; i < shiftContribution.length; ++i) {
            shiftContribution[i] = trait.getShiftForBranchLength(node)[i] + shiftContributionParent[i];
        }
    }
    return shiftContribution;
}

Example 8 with MutableTreeModel

use of dr.evolution.tree.MutableTreeModel in project beast-mcmc by beast-dev.

the class DataFromTreeTipsParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    TreeTraitParserUtilities utilities = new TreeTraitParserUtilities();
    String traitName = (String) xo.getAttribute(TreeTraitParserUtilities.TRAIT_NAME);
    MutableTreeModel treeModel = (MutableTreeModel) xo.getChild(MutableTreeModel.class);
    TreeTraitParserUtilities.TraitsAndMissingIndices returnValue = utilities.parseTraitsFromTaxonAttributes(xo, traitName, treeModel, true);
    MatrixParameter dataParameter = MatrixParameter.recast(returnValue.traitParameter.getId(), returnValue.traitParameter);
    if (xo.hasChildNamed(TreeTraitParserUtilities.MISSING)) {
        Parameter missing = (Parameter) xo.getChild(TreeTraitParserUtilities.MISSING).getChild(Parameter.class);
        missing.setDimension(dataParameter.getDimension());
        for (int i = 0; i < missing.getDimension(); i++) {
            if (returnValue.missingIndices.contains(i)) {
                missing.setParameterValue(i, 1);
            } else {
                missing.setParameterValue(i, 0);
            }
        }
    }
    return dataParameter;
}

Example 9 with MutableTreeModel

use of dr.evolution.tree.MutableTreeModel 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);
}

Example 10 with MutableTreeModel

use of dr.evolution.tree.MutableTreeModel in project beast-mcmc by beast-dev.

the class EpochBranchModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Logger.getLogger("dr.evomodel").info("\nUsing multi-epoch branch model.");
    MutableTreeModel treeModel = (MutableTreeModel) xo.getChild(MutableTreeModel.class);
    SubstitutionModel ancestralSubstitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
    List<Epoch> epochs = new ArrayList<Epoch>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        if (xo.getChild(i) instanceof XMLObject) {
            XMLObject xoc = (XMLObject) xo.getChild(i);
            if (xoc.getName().equals(EPOCH)) {
                Parameter tt = null;
                if (xoc.hasAttribute(TRANSITION_TIME)) {
                    double t = xoc.getAttribute(TRANSITION_TIME, 0.0);
                    tt = new Parameter.Default(1, t);
                }
                SubstitutionModel s = (SubstitutionModel) xoc.getChild(SubstitutionModel.class);
                if (xoc.hasChildNamed(TRANSITION_TIME)) {
                    if (tt != null) {
                        throw new XMLParseException("An epoch cannot have a transitionTime attribute and a parameter");
                    }
                    tt = (Parameter) xoc.getElementFirstChild(TRANSITION_TIME);
                }
                epochs.add(new Epoch(s, tt));
            }
        }
    }
    Collections.sort(epochs);
    List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
    CompoundParameter transitionTimes = new CompoundParameter("epochTimes");
    for (Epoch epoch : epochs) {
        substitutionModels.add(epoch.substitutionModel);
        transitionTimes.addParameter(epoch.timeParameter);
    }
    substitutionModels.add(ancestralSubstitutionModel);
    return new EpochBranchModel(treeModel, substitutionModels, transitionTimes);
}