Examples with BranchRateModel dr.evomodel.branchratemodel.BranchRateModel used on opensource projects

Example 16 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class RateCovarianceStatisticParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    String name = xo.getAttribute(Statistic.NAME, xo.getId());
    Tree tree = (Tree) xo.getChild(Tree.class);
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    return new RateCovarianceStatistic(name, tree, branchRateModel);
}

Example 17 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class RateStatisticParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    final String name = xo.getAttribute(Statistic.NAME, xo.getId());
    final Tree tree = (Tree) xo.getChild(Tree.class);
    final BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    final boolean internal = xo.getBooleanAttribute("internal");
    final boolean external = xo.getBooleanAttribute("external");
    if (!(internal || external)) {
        throw new XMLParseException("At least one of internal and external must be true!");
    }
    final String mode = xo.getStringAttribute(MODE);
    return new RateStatistic(name, tree, branchRateModel, external, internal, mode);
}

Example 18 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class DiffusionRateStatistic method getStatisticValue.

public double getStatisticValue(int dim) {
    String traitName = traitLikelihoods.get(0).getTraitName();
    double treelength = 0;
    double treeDistance = 0;
    double maxDistanceFromRoot = 0;
    double maxDistanceOverTimeFromRoot = 0;
    //double[] rates =  null;
    List<Double> rates = new ArrayList<Double>();
    //double[] diffusionCoefficients =  null;
    List<Double> diffusionCoefficients = new ArrayList<Double>();
    double waDiffusionCoefficient = 0;
    double lowerHeight = heightLowers[dim];
    double upperHeight = Double.MAX_VALUE;
    if (heightLowers.length == 1) {
        upperHeight = heightUpper;
    } else {
        if (dim > 0) {
            if (!cumulative) {
                upperHeight = heightLowers[dim - 1];
            }
        }
    }
    for (AbstractMultivariateTraitLikelihood traitLikelihood : traitLikelihoods) {
        MultivariateTraitTree tree = traitLikelihood.getTreeModel();
        BranchRateModel branchRates = traitLikelihood.getBranchRateModel();
        for (int i = 0; i < tree.getNodeCount(); i++) {
            NodeRef node = tree.getNode(i);
            if (node != tree.getRoot()) {
                NodeRef parentNode = tree.getParent(node);
                if ((tree.getNodeHeight(parentNode) > lowerHeight) && (tree.getNodeHeight(node) < upperHeight)) {
                    double[] trait = traitLikelihood.getTraitForNode(tree, node, traitName);
                    double[] parentTrait = traitLikelihood.getTraitForNode(tree, parentNode, traitName);
                    double[] traitUp = parentTrait;
                    double[] traitLow = trait;
                    double timeUp = tree.getNodeHeight(parentNode);
                    double timeLow = tree.getNodeHeight(node);
                    double rate = (branchRates != null ? branchRates.getBranchRate(tree, node) : 1.0);
                    MultivariateDiffusionModel diffModel = traitLikelihood.diffusionModel;
                    double[] precision = diffModel.getPrecisionParameter().getParameterValues();
                    if (tree.getNodeHeight(parentNode) > upperHeight) {
                        timeUp = upperHeight;
                        //TODO: implement TrueNoise??
                        traitUp = imputeValue(trait, parentTrait, upperHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, false);
                    }
                    if (tree.getNodeHeight(node) < lowerHeight) {
                        timeLow = lowerHeight;
                        traitLow = imputeValue(trait, parentTrait, lowerHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, false);
                    }
                    double time = timeUp - timeLow;
                    treelength += time;
                    double[] rootTrait = traitLikelihood.getTraitForNode(tree, tree.getRoot(), traitName);
                    if (useGreatCircleDistances && (trait.length == 2)) {
                        // Great Circle distance
                        SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(traitLow[0], traitLow[1]);
                        SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(traitUp[0], traitUp[1]);
                        double distance = coord1.distance(coord2);
                        treeDistance += distance;
                        double dc = Math.pow(distance, 2) / (4 * time);
                        diffusionCoefficients.add(dc);
                        waDiffusionCoefficient += dc * time;
                        rates.add(distance / time);
                        SphericalPolarCoordinates rootCoord = new SphericalPolarCoordinates(rootTrait[0], rootTrait[1]);
                        double tempDistanceFromRoot = rootCoord.distance(coord2);
                        if (tempDistanceFromRoot > maxDistanceFromRoot) {
                            maxDistanceFromRoot = tempDistanceFromRoot;
                            maxDistanceOverTimeFromRoot = tempDistanceFromRoot / (tree.getNodeHeight(tree.getRoot()) - timeLow);
                            //distance between traitLow and traitUp for maxDistanceFromRoot
                            if (timeUp == upperHeight) {
                                maxDistanceFromRoot = distance;
                                maxDistanceOverTimeFromRoot = distance / time;
                            }
                        }
                    } else {
                        double distance = getNativeDistance(traitLow, traitUp);
                        treeDistance += distance;
                        double dc = Math.pow(distance, 2) / (4 * time);
                        diffusionCoefficients.add(dc);
                        waDiffusionCoefficient += dc * time;
                        rates.add(distance / time);
                        double tempDistanceFromRoot = getNativeDistance(traitLow, rootTrait);
                        if (tempDistanceFromRoot > maxDistanceFromRoot) {
                            maxDistanceFromRoot = tempDistanceFromRoot;
                            maxDistanceOverTimeFromRoot = tempDistanceFromRoot / (tree.getNodeHeight(tree.getRoot()) - timeLow);
                            //distance between traitLow and traitUp for maxDistanceFromRoot
                            if (timeUp == upperHeight) {
                                maxDistanceFromRoot = distance;
                                maxDistanceOverTimeFromRoot = distance / time;
                            }
                        }
                    }
                }
            }
        }
    }
    if (summaryStat == summaryStatistic.DIFFUSION_RATE) {
        if (summaryMode == Mode.AVERAGE) {
            return DiscreteStatistics.mean(toArray(rates));
        } else if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(rates));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(rates));
            return Math.sqrt(DiscreteStatistics.variance(toArray(rates), mean)) / mean;
        } else {
            return treeDistance / treelength;
        }
    } else if (summaryStat == summaryStatistic.DIFFUSION_COEFFICIENT) {
        if (summaryMode == Mode.AVERAGE) {
            return DiscreteStatistics.mean(toArray(diffusionCoefficients));
        } else if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(diffusionCoefficients));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(diffusionCoefficients));
            return Math.sqrt(DiscreteStatistics.variance(toArray(diffusionCoefficients), mean)) / mean;
        } else {
            return waDiffusionCoefficient / treelength;
        }
    } else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE) {
        return maxDistanceFromRoot;
    } else {
        return maxDistanceOverTimeFromRoot;
    }
}

Example 19 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class MicrosatelliteFullAncestryImportanceSamplingOperatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    final double weight = xo.getDoubleAttribute(MCMCOperator.WEIGHT);
    final Parameter parameter = (Parameter) xo.getChild(Parameter.class);
    final MicrosatelliteSamplerTreeModel msatSamplerTreeModel = (MicrosatelliteSamplerTreeModel) xo.getChild(MicrosatelliteSamplerTreeModel.class);
    final MicrosatelliteModel msatModel = (MicrosatelliteModel) xo.getChild(MicrosatelliteModel.class);
    final BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    return new MicrosatelliteFullAncestryImportanceSamplingOperator(parameter, msatSamplerTreeModel, msatModel, branchRateModel, weight);
}

Example 20 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class CodonPartitionedRobustCountingParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    AncestralStateBeagleTreeLikelihood[] partition = new AncestralStateBeagleTreeLikelihood[3];
    String[] labels = new String[] { FIRST, SECOND, THIRD };
    int patternCount = -1;
    BranchRateModel testBranchRateModel = null;
    for (int i = 0; i < 3; i++) {
        partition[i] = (AncestralStateBeagleTreeLikelihood) xo.getChild(labels[i]).getChild(AncestralStateBeagleTreeLikelihood.class);
        if (i == 0) {
            patternCount = partition[i].getPatternCount();
        } else {
            if (partition[i].getPatternCount() != patternCount) {
                throw new XMLParseException("Codon-partitioned robust counting requires all partitions to have the same length." + " Make sure that partitions include all unique sites and do not strip gaps.");
            }
        }
        // Ensure that siteRateModel has one category
        if (partition[i].getSiteRateModel().getCategoryCount() > 1) {
            throw new XMLParseException("Robust counting currently only implemented for single category models");
        }
        // Ensure that branchRateModel is the same across all partitions
        if (testBranchRateModel == null) {
            testBranchRateModel = partition[i].getBranchRateModel();
        } else if (testBranchRateModel != partition[i].getBranchRateModel()) {
            throw new XMLParseException("Robust counting currently requires the same branch rate model for all partitions");
        }
    }
    TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
    Codons codons = Codons.UNIVERSAL;
    if (xo.hasAttribute(GeneticCode.GENETIC_CODE)) {
        String codeStr = xo.getStringAttribute(GeneticCode.GENETIC_CODE);
        codons = Codons.findByName(codeStr);
    }
    String labelingString = (String) xo.getAttribute(LABELING);
    CodonLabeling codonLabeling = CodonLabeling.parseFromString(labelingString);
    if (codonLabeling == null) {
        throw new XMLParseException("Unrecognized codon labeling '" + labelingString + "'");
    }
    String branchFormatString = xo.getAttribute(BRANCH_FORMAT, StratifiedTraitOutputFormat.SUM_OVER_SITES.getText());
    StratifiedTraitOutputFormat branchFormat = StratifiedTraitOutputFormat.parseFromString(branchFormatString);
    if (branchFormat == null) {
        throw new XMLParseException("Unrecognized branch output format '" + branchFormat + "'");
    }
    String logFormatString = xo.getAttribute(LOG_FORMAT, StratifiedTraitOutputFormat.SUM_OVER_SITES.getText());
    StratifiedTraitOutputFormat logFormat = StratifiedTraitOutputFormat.parseFromString(logFormatString);
    if (logFormat == null) {
        throw new XMLParseException("Unrecognized log output format '" + branchFormat + "'");
    }
    boolean useUniformization = xo.getAttribute(USE_UNIFORMIZATION, false);
    boolean includeExternalBranches = xo.getAttribute(INCLUDE_EXTERNAL, true);
    boolean includeInternalBranches = xo.getAttribute(INCLUDE_INTERNAL, true);
    boolean doUnconditionedPerBranch = xo.getAttribute(DO_UNCONDITIONED_PER_BRANCH, false);
    boolean averageRates = xo.getAttribute(AVERAGE_RATES, true);
    boolean saveCompleteHistory = xo.getAttribute(SAVE_HISTORY, false);
    boolean useNewNeutralModel = xo.getAttribute(USE_NEW_NEUTRAL_MODEL, false);
    String prefix = xo.hasAttribute(PREFIX) ? xo.getStringAttribute(PREFIX) : null;
    return new CodonPartitionedRobustCounting(xo.getId(), tree, partition, codons, codonLabeling, useUniformization, includeExternalBranches, includeInternalBranches, doUnconditionedPerBranch, saveCompleteHistory, averageRates, useNewNeutralModel, branchFormat, logFormat, prefix);
}