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

Example 41 with BranchRateModel

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

the class SequenceSimulator method main.

// getDefaultSiteModel
public static void main(String[] args) {
    try {
        int nReplications = 10;
        // create tree
        NewickImporter importer = new NewickImporter("((A:1.0,B:1.0)AB:1.0,(C:1.0,D:1.0)CD:1.0)ABCD;");
        Tree tree = importer.importTree(null);
        // create site model
        SiteModel siteModel = getDefaultSiteModel();
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // feed to sequence simulator and generate leaves
        SequenceSimulator treeSimulator = new SequenceSimulator(tree, siteModel, branchRateModel, nReplications);
        Sequence ancestralSequence = new Sequence();
        ancestralSequence.appendSequenceString("TCAGGTCAAG");
        treeSimulator.setAncestralSequence(ancestralSequence);
        System.out.println(treeSimulator.simulate().toString());
    } catch (Exception e) {
        e.printStackTrace();
    }
//END: try-catch block
}

Example 42 with BranchRateModel

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

the class ContinuousDiffusionStatistic method getStatisticValue.

public double getStatisticValue(int dim) {
    double treeLength = 0;
    double treeDistance = 0;
    double totalMaxDistanceFromRoot = 0;
    // can only be used when cumulative and not associated with discrete state (not based on the distances on the branches from the root up that point)
    double maxDistanceFromRootCumulative = 0;
    double maxBranchDistanceFromRoot = 0;
    // can only be used when cumulative and not associated with discrete state (not based on the distances on the branches from the root up that point)
    double maxDistanceOverTimeFromRootWA = 0;
    double maxBranchDistanceOverTimeFromRootWA = 0;
    List<Double> rates = new ArrayList<Double>();
    List<Double> distances = new ArrayList<Double>();
    List<Double> times = new ArrayList<Double>();
    List<Double> traits = new ArrayList<Double>();
    List<double[]> traits2D = 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();
        String traitName = traitLikelihood.getTraitName();
        for (int i = 0; i < tree.getNodeCount(); i++) {
            NodeRef node = tree.getNode(i);
            if (node != tree.getRoot()) {
                NodeRef parentNode = tree.getParent(node);
                boolean testNode = true;
                if (branchset.equals(BranchSet.CLADE)) {
                    try {
                        testNode = inClade(tree, node, taxonList);
                    } catch (TreeUtils.MissingTaxonException mte) {
                        throw new RuntimeException(mte.toString());
                    }
                } else if (branchset.equals(BranchSet.BACKBONE)) {
                    if (backboneTime > 0) {
                        testNode = onAncestralPathTime(tree, node, backboneTime);
                    } else {
                        try {
                            testNode = onAncestralPathTaxa(tree, node, taxonList);
                        } catch (TreeUtils.MissingTaxonException mte) {
                            throw new RuntimeException(mte.toString());
                        }
                    }
                }
                if (testNode) {
                    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);
                        //                        System.out.println(rate);
                        MultivariateDiffusionModel diffModel = traitLikelihood.diffusionModel;
                        double[] precision = diffModel.getPrecisionParameter().getParameterValues();
                        History history = null;
                        if (stateString != null) {
                            history = setUpHistory(markovJumpLikelihood.getHistoryForNode(tree, node, SITE), markovJumpLikelihood.getStatesForNode(tree, node)[SITE], markovJumpLikelihood.getStatesForNode(tree, parentNode)[SITE], timeLow, timeUp);
                        }
                        if (tree.getNodeHeight(parentNode) > upperHeight) {
                            timeUp = upperHeight;
                            traitUp = imputeValue(trait, parentTrait, upperHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, trueNoise);
                            if (stateString != null) {
                                history.truncateUpper(timeUp);
                            }
                        }
                        if (tree.getNodeHeight(node) < lowerHeight) {
                            timeLow = lowerHeight;
                            traitLow = imputeValue(trait, parentTrait, lowerHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, trueNoise);
                            if (stateString != null) {
                                history.truncateLower(timeLow);
                            }
                        }
                        if (dimension > traitLow.length) {
                            System.err.println("specified trait dimension for continuous trait summary, " + dimension + ", is > dimensionality of trait, " + traitLow.length + ". No trait summarized.");
                        } else {
                            traits.add(traitLow[(dimension - 1)]);
                        }
                        if (traitLow.length == 2) {
                            traits2D.add(traitLow);
                        }
                        double time;
                        if (stateString != null) {
                            time = history.getStateTime(stateString);
                        //                            System.out.println("tine before = "+(timeUp - timeLow)+", time after= "+time);
                        } else {
                            time = timeUp - timeLow;
                        }
                        treeLength += time;
                        times.add(time);
                        //setting up continuous trait values for heights in discrete trait history
                        if (stateString != null) {
                            history.setTraitsforHeights(traitUp, traitLow, precision, rate, trueNoise);
                        }
                        double[] rootTrait = traitLikelihood.getTraitForNode(tree, tree.getRoot(), traitName);
                        double timeFromRoot = (tree.getNodeHeight(tree.getRoot()) - timeLow);
                        if (useGreatCircleDistances && (trait.length == 2)) {
                            // Great Circle distance
                            double distance;
                            if (stateString != null) {
                                distance = history.getStateGreatCircleDistance(stateString);
                            } else {
                                distance = getGreatCircleDistance(traitLow, traitUp);
                            }
                            distances.add(distance);
                            if (time > 0) {
                                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 tempDistanceFromRootLow = rootCoord.distance(new SphericalPolarCoordinates(traitUp[0], traitUp[1]));
                            if (tempDistanceFromRootLow > totalMaxDistanceFromRoot) {
                                totalMaxDistanceFromRoot = tempDistanceFromRootLow;
                                if (stateString != null) {
                                    double[] stateTimeDistance = getStateTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, precision, branchRates, true);
                                    if (stateTimeDistance[0] > 0) {
                                        maxDistanceFromRootCumulative = tempDistanceFromRootLow * (stateTimeDistance[0] / timeFromRoot);
                                        maxDistanceOverTimeFromRootWA = maxDistanceFromRootCumulative / stateTimeDistance[0];
                                        maxBranchDistanceFromRoot = stateTimeDistance[1];
                                        maxBranchDistanceOverTimeFromRootWA = stateTimeDistance[1] / stateTimeDistance[0];
                                    }
                                } else {
                                    maxDistanceFromRootCumulative = tempDistanceFromRootLow;
                                    maxDistanceOverTimeFromRootWA = tempDistanceFromRootLow / timeFromRoot;
                                    double[] timeDistance = getTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, true);
                                    maxBranchDistanceFromRoot = timeDistance[1];
                                    maxBranchDistanceOverTimeFromRootWA = timeDistance[1] / timeDistance[0];
                                }
                                //distance between traitLow and traitUp for maxDistanceFromRootCumulative
                                if (timeUp == upperHeight) {
                                    if (time > 0) {
                                        maxDistanceFromRootCumulative = distance;
                                        maxDistanceOverTimeFromRootWA = distance / time;
                                        maxBranchDistanceFromRoot = distance;
                                        maxBranchDistanceOverTimeFromRootWA = distance / time;
                                    }
                                }
                            }
                        } else {
                            double distance;
                            if (stateString != null) {
                                distance = history.getStateNativeDistance(stateString);
                            } else {
                                distance = getNativeDistance(traitLow, traitUp);
                            }
                            distances.add(distance);
                            if (time > 0) {
                                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 > totalMaxDistanceFromRoot) {
                                totalMaxDistanceFromRoot = tempDistanceFromRoot;
                                if (stateString != null) {
                                    double[] stateTimeDistance = getStateTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, precision, branchRates, false);
                                    if (stateTimeDistance[0] > 0) {
                                        maxDistanceFromRootCumulative = tempDistanceFromRoot * (stateTimeDistance[0] / timeFromRoot);
                                        maxDistanceOverTimeFromRootWA = maxDistanceFromRootCumulative / stateTimeDistance[0];
                                        maxBranchDistanceFromRoot = stateTimeDistance[1];
                                        maxBranchDistanceOverTimeFromRootWA = stateTimeDistance[1] / stateTimeDistance[0];
                                    }
                                } else {
                                    maxDistanceFromRootCumulative = tempDistanceFromRoot;
                                    maxDistanceOverTimeFromRootWA = tempDistanceFromRoot / timeFromRoot;
                                    double[] timeDistance = getTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, false);
                                    maxBranchDistanceFromRoot = timeDistance[1];
                                    maxBranchDistanceOverTimeFromRootWA = timeDistance[1] / timeDistance[0];
                                }
                                //distance between traitLow and traitUp for maxDistanceFromRootCumulative
                                if (timeUp == upperHeight) {
                                    if (time > 0) {
                                        maxDistanceFromRootCumulative = distance;
                                        maxDistanceOverTimeFromRootWA = distance / time;
                                        maxBranchDistanceFromRoot = distance;
                                        maxBranchDistanceOverTimeFromRootWA = 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) {
            final double mean = DiscreteStatistics.mean(toArray(rates));
            return Math.sqrt(DiscreteStatistics.variance(toArray(rates), mean)) / mean;
        //weighted average
        } else {
            return treeDistance / treeLength;
        }
    } else if (summaryStat == summaryStatistic.TRAIT) {
        if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(traits));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(traits));
            return Math.sqrt(DiscreteStatistics.variance(toArray(traits), mean)) / mean;
        // default is average. A warning is thrown by the parser when trying to use WEIGHTED_AVERAGE
        } else {
            return DiscreteStatistics.mean(toArray(traits));
        }
    } else if (summaryStat == summaryStatistic.TRAIT2DAREA) {
        double area = getAreaFrom2Dtraits(traits2D, 0.99);
        return area;
    } 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;
        }
    //wavefront distance
    //TODO: restrict to non state-specific wavefrontDistance/rate
    } else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE) {
        return maxDistanceFromRootCumulative;
    //            return maxBranchDistanceFromRoot;
    } else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE_PHYLO) {
        return maxBranchDistanceFromRoot;
    //wavefront rate, only weighted average TODO: extend for average, median, COEFFICIENT_OF_VARIATION?
    } else if (summaryStat == summaryStatistic.WAVEFRONT_RATE) {
        return maxDistanceOverTimeFromRootWA;
    //            return maxBranchDistanceOverTimeFromRootWA;
    } else if (summaryStat == summaryStatistic.DIFFUSION_DISTANCE) {
        return treeDistance;
    //DIFFUSION_TIME
    } else if (summaryStat == summaryStatistic.DISTANCE_TIME_CORRELATION) {
        if (summaryMode == Mode.SPEARMAN) {
            return getSpearmanRho(convertDoubles(times), convertDoubles(distances));
        } else if (summaryMode == Mode.R_SQUARED) {
            Regression r = new Regression(convertDoubles(times), convertDoubles(distances));
            return r.getRSquared();
        } else {
            Regression r = new Regression(convertDoubles(times), convertDoubles(distances));
            return r.getCorrelationCoefficient();
        }
    } else {
        return treeLength;
    }
}

Example 43 with BranchRateModel

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

the class MicrosatelliteSimulatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Microsatellite msatDataType = (Microsatellite) xo.getChild(Microsatellite.class);
    Taxa taxa = (Taxa) xo.getChild(Taxa.class);
    Tree tree = (Tree) xo.getChild(Tree.class);
    MicrosatelliteModel msatModel = (MicrosatelliteModel) xo.getChild(MicrosatelliteModel.class);
    BranchRateModel brModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (brModel == null) {
        brModel = new DefaultBranchRateModel();
    }
    MicrosatelliteSimulator msatSim = new MicrosatelliteSimulator(msatDataType, taxa, tree, new GammaSiteModel(msatModel), brModel);
    Patterns patterns = msatSim.simulateMsatPattern();
    String msatPatId = xo.getAttribute("id", "simMsatPat");
    patterns.setId(msatPatId);
    MicrosatellitePatternParser.printDetails(patterns);
    MicrosatellitePatternParser.printMicrosatContent(patterns);
    return patterns;
}

Example 44 with BranchRateModel

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

the class VariableBranchCompleteHistorySimulatorTest method testHKYVariableSimulation.

public void testHKYVariableSimulation() {
    System.out.println("Starting HKY variable branch simulation");
    Parameter kappa = new Parameter.Default(1, 2.0);
    double[] pi = { 0.45, 0.05, 0.25, 0.25 };
    Parameter freqs = new Parameter.Default(pi);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    int stateCount = hky.getDataType().getStateCount();
    Parameter mu = new Parameter.Default(1, 0.5);
    Parameter alpha = new Parameter.Default(1, 0.5);
    GammaSiteRateModel siteModel = new GammaSiteRateModel("gammaModel", mu, alpha, 4, null);
    siteModel.setSubstitutionModel(hky);
    BranchRateModel branchRateModel = new DefaultBranchRateModel();
    double analyticResult = TreeUtils.getTreeLength(tree, tree.getRoot()) * mu.getParameterValue(0);
    int nSites = 200;
    double[] register1 = new double[stateCount * stateCount];
    double[] register2 = new double[stateCount * stateCount];
    // Count all jumps
    MarkovJumpsCore.fillRegistrationMatrix(register1, stateCount);
    // Move some jumps from 1 to 2
    register1[1 * stateCount + 2] = 0;
    register2[1 * stateCount + 2] = 1;
    register1[1 * stateCount + 3] = 0;
    register2[1 * stateCount + 3] = 1;
    register1[2 * stateCount + 3] = 0;
    register2[2 * stateCount + 3] = 1;
    double[] branchValues = { 10.0, 10.0, 10.0, 10.0, 10.0 };
    Parameter branchValuesParam = new Parameter.Default(branchValues);
    runSimulation(N, tree, siteModel, branchRateModel, nSites, new double[][] { register1, register2 }, analyticResult, kappa, branchValuesParam);
}