Examples with TreeIntervals dr.evomodel.coalescent.TreeIntervals used on opensource projects

Example 1 with TreeIntervals

use of dr.evomodel.coalescent.TreeIntervals in project beast-mcmc by beast-dev.

the class RandomWalkGenerator method main.

public static void main(String[] args) throws Exception {
    // Parameter data = new Parameter.Default(100);
    // Parameter firstElemPrec = new Parameter.Default(0.25);
    // Parameter prec = new Parameter.Default(4.0);
    // RandomWalkGenerator gen = new RandomWalkGenerator(data, firstElemPrec, prec);
    // 
    // System.out.println("Data = " + Arrays.toString(data.getParameterValues()));
    // 
    // double[] newData = gen.nextRandom();
    // 
    // System.out.println("New Data = " + Arrays.toString(newData));
    // 
    // System.out.print("dat = c(");
    // for (int i = 0; i < newData.length; i++) {
    // data.setParameterValue(i, newData[i]);
    // System.out.printf("%.3f,", newData[i]);
    // }
    // System.out.println();
    // 
    // System.out.println("New Data too = " + Arrays.toString(data.getParameterValues()));
    // 
    // System.out.println("Log-likelihood = " + gen.getLogLikelihood());
    // 
    // double[] values = new double[5];
    // System.out.println("New Data = " + Arrays.toString(values));
    NewickImporter importer = new NewickImporter("((((t75_0:8.196509369,(t76_0:5.19700607,t77_0:5.117706795):2.208323646):3.611573553,(t78_0:10.09414276,(t79_0:1.555758051,t80_0:0.6394209368):8.440019033):0.7958799608):36.90916229,(t81_0:6.965230429,t82_0:6.660700811):39.60826152):84.84805599,(((((t83_0:8.461009149,(((t84_0:6.681657444,t85_0:6.329850896):0.4970932551,(t86_0:5.709749921,t87_0:5.207473638):1.032134073):0.2907208781,t88_0:5.860923226):0.6714693579):0.9286485984,(t89_0:4.35176728,(t90_0:3.383962439,t91_0:2.365432603):0.5895519359):2.156672574):0.8045815339,(((t92_0:2.347375928,(t93_0:1.806692774,t95_0:0.2737786494):0.06086363256):1.778880139,t96_0:1.745952006):0.1290466242,t97_0:1.831227907):1.231195333):2.445021787,((t98_0:0.8797524272,t99_0:0.8146248384):2.02910329,t100_0:2.379928659):2.546864865):8.217636348,(t94_0:0.1096345611,(((((t72_0:3.137047799,(t65_0:0.7428752174,(((t54_0:2.298056151,t55_0:1.386634319):0.07055842264,(t56_0:0.6211411999,t57_0:0.5106794843):0.7290213103):1.302958439,(((t32_0:0.01892218578,(((t28_0:0.8330661569,((((t25_0:0.1067751399,t1_0:9.307914765):0.2959298825,t2_0:9.425199804):0.07190576719,t3_0:8.826132669):0.2342818567,(t4_0:8.321528621,t5_0:8.0599082):0.2509419834):0.9408918394):0.7425197124,t6_0:9.694992029):2.36508096,t7_0:11.2758228):0.1994631268):8.443974525,((t8_0:10.76752205,((((t9_0:1.493044826,(t10_0:0.6698090384,t11_0:0.2989691445):0.7176777198):2.286111176,t12_0:2.654871592):1.588547853,t13_0:4.203655954):3.489759422,t14_0:7.542791391):1.88216266):8.014313382,((((t15_0:4.455846169,(t16_0:0.9348779818,t17_0:0.7092212787):2.927543424):0.8337954126,t18_0:4.371535283):2.561638646,(t19_0:2.317319668,(t20_0:0.2652308407,t21_0:0.2469742062):1.962177717):3.646397641):4.201926648,(((t22_0:3.626205298,(t23_0:1.814394049,t24_0:1.722932275):1.576199638):0.6832367053,(t26_0:0.143315467,t27_0:0.135478801):1.855008):4.430521602,(t29_0:0.277418938,t30_0:0.2751244534):2.656950561):1.277291788):4.44558088):0.4877611393):1.124425501,(((t31_0:2.199804416,t33_0:1.72439591):5.99779354,t34_0:7.697646612):0.2314130313,t35_0:7.7292073):1.578317139):0.6787721178):5.713754691):4.899357656):2.690189889,((((((t36_0:8.464440253,((t37_0:1.516370472,(t38_0:0.4518582713,t39_0:0.3776256558):1.039775741):5.008595559,(t40_0:1.662591706,t41_0:1.311689096):4.010173892):1.134735834):0.9467574058,t42_0:7.365820535):0.4308090259,((((t43_0:4.977793922,t44_0:4.1359021):0.04952307756,t45_0:4.176477429):0.004420431545,t46_0:3.828337992):0.1806606434,((t47_0:2.488665975,t48_0:2.460902709):0.2186504359,(t49_0:1.415777402,t50_0:1.358357241):0.6077191637):0.7478779126):1.81853207):7.625751786,(t51_0:1.460127624,t52_0:1.344261146):10.22672898):4.701163396,((t53_0:9.827071398,(t58_0:1.163158974,t59_0:0.8080003353):4.25790361):5.797987449,(t60_0:8.676431981,(t61_0:0.7238168631,t62_0:0.5772959095):6.677580493):1.790578221):0.5861914293):0.02660969378,t63_0:9.157389255):1.055890624):1.298920987,((t64_0:1.403696795,t66_0:0.17472493):4.761586473,(t67_0:0.8682056765,t68_0:0.7870286654):3.996084715):3.608662392):1.291247032,(t69_0:3.938693436,(t70_0:1.363501591,t71_0:1.047508982):1.675768049):5.725301901):1.748629464,(t73_0:0.4162471532,t74_0:0.3326790869):5.542586246):2.615577619):15.02704211):110.212517)");
    Tree tree = importer.importNextTree();
    Parameter N0 = new Parameter.Default(new double[] { 10831.41518, 18564.20864, 4012.90013, 4844.7051, 2191.07171, 864.58823, 1954.18284, 7233.75831, 1624.29542, 1204.22108, 2424.79383, 3236.60901, 741.81022, 309.39496, 739.2772, 436.22362, 793.08349, 708.42725, 470.52578, 731.14609, 609.88306, 380.00053, 333.14087, 445.31313, 385.08516, 310.62223, 117.55287, 15.35954, 5.30181, 9.55051, 31.10367, 59.87409, 53.83162, 28.2465, 57.93228, 26.31053, 61.78336, 98.70132, 104.77869, 255.49281, 715.66912, 324.43738, 163.67341, 89.23281, 91.41586, 133.36657, 71.08849, 93.51971, 329.73119, 595.42027, 274.05989, 192.55929, 200.60026, 180.83781, 225.98434, 96.75023, 108.98978, 169.2301, 174.38823, 206.05087, 24.43085, 32.28332, 29.62218, 22.61263, 67.1912, 12.66908, 14.31542, 17.22216, 40.34572, 29.89561, 269.60226, 792.19973, 1639.95258, 495.38584, 438.89219, 266.81732, 943.54934, 1430.35456, 454.51251, 310.67618, 500.64104, 308.28057, 289.18023, 1065.39065, 633.67468, 2127.13342, 4534.01476, 1561.61864, 1958.98386, 2209.80371, 655.44965, 352.61598, 107.18735, 53.62083, 124.80736, 189.3125, 67.15406, 7.91813, 15.98184, 26.89309 });
    double[] epochLengths = new double[] { 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024, 1.699024 };
    DemographicModel Ne = new PiecewisePopulationModel("Ne(t)", N0, epochLengths, false, Units.Type.DAYS);
    TreeIntervals intervalList = new TreeIntervals(tree, null, null);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(intervalList, Ne);
    double logLik = coalescent.getLogLikelihood();
    System.out.printf("Loglik = %f\n", logLik);
    Parameter prec = new Parameter.Default(2.0);
    Parameter logPop = new TransformedParameter(N0, new Transform.LogTransform());
    RandomWalkGenerator rwg = new RandomWalkGenerator(logPop, new Parameter.Default(0.01), prec);
    double logFieldPri = rwg.getLogLikelihood();
    System.out.printf("LogFieldPri = %f\n", logFieldPri);
// ArrayList<Double> al = new ArrayList<Double>();
// List<Double> lst = new Collections.synchronizedList(al);
}

Example 2 with TreeIntervals

use of dr.evomodel.coalescent.TreeIntervals in project beast-mcmc by beast-dev.

the class BigFastTreeTreeIntervalsTest method testOrderChange.

// 0.5
// 1.0           +------ 0 * [1.5]
// +----------------------- |(7) *  [2]
// 1.0           |                        |                                             +-- 1 * [0]
// +----------------------- |(8) *  [3]              +---------------------------------------------|(6) [0.1]
// |                        |                                                                      +--- 2 * [0]
// |                        |                1.5
// |(10) &  [4]             +-------------------------------- 3 * [1.5]
// |                                                       2.0
// |        1.5                   +----------------------------------------------------- 4 & [0.5]
// +-----------------------------|(9) & [2.5]     1.51
// +-------------------------------------------- 5 & [0.99]
public void testOrderChange() throws TreeUtils.MissingTaxonException, IOException, Importer.ImportException {
    NewickImporter importer = new NewickImporter("(((0:0.5,(1:1.0,2:1.0)n6:1.0)n7:1.0,3:1.5)n8:1.0,(4:2.0,5:1.51)n9:1.5)n10;");
    tree = new BigFastTreeModel(importer.importTree(null));
    IntervalList bigFastIntervals = new BigFastTreeIntervals(tree);
    IntervalList intervals = new TreeIntervals(tree, null, null);
    bigFastIntervals.calculateIntervals();
    // node             1-2,6,  4  , 5  ,  0  ,  3, 7  , 9  , 8  , 10
    boolean pass = true;
    NodeRef parent = tree.getNode(6);
    NodeRef sibling = tree.getNode(2);
    NodeRef grandParent = tree.getParent(parent);
    NodeRef j = tree.getNode(2);
    tree.beginTreeEdit();
    tree.removeChild(grandParent, parent);
    tree.removeChild(parent, sibling);
    tree.addChild(grandParent, sibling);
    NodeRef jParent = tree.getParent(j);
    tree.removeChild(jParent, j);
    tree.addChild(jParent, parent);
    tree.addChild(parent, j);
    tree.setNodeHeight(parent, 0.1);
    tree.endTreeEdit();
    bigFastIntervals.calculateIntervals();
    List<Integer> missed = new ArrayList<>();
    for (int i = 0; i < bigFastIntervals.getIntervalCount(); i++) {
        if (bigFastIntervals.getInterval(i) != intervals.getInterval(i)) {
            System.out.println("expected: " + intervals.getInterval(i) + " got: " + bigFastIntervals.getInterval(i));
            missed.add(i);
        }
    }
    System.out.println(missed);
    assertEquals(0, missed.size());
}

Example 3 with TreeIntervals

use of dr.evomodel.coalescent.TreeIntervals in project beast-mcmc by beast-dev.

the class StrictClockTest method testStrictClock.

public void testStrictClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 380.0, 0, 38000.0);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    TreeIntervals intervalList = new TreeIntervals(treeModel, null, null);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(intervalList, constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter rateParameter = new Parameter.Default(StrictClockBranchRates.RATE, 2.3E-5, 0, 100.0);
    StrictClockBranchRates branchRateModel = new StrictClockBranchRates(rateParameter);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 100.0);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    // siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    // treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(rateParameter, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter allInternalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(true, true, false);
    operator = new UpDownOperator(new Scalable[] { new Scalable.Default(rateParameter) }, new Scalable[] { new Scalable.Default(allInternalHeights) }, 0.75, 3.0, AdaptationMode.ADAPTATION_ON);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = ((DefaultTreeModel) treeModel).getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(((DefaultTreeModel) treeModel), 15.0, 1.0, true, false, false, false, AdaptationMode.ADAPTATION_ON, AdaptableMCMCOperator.DEFAULT_ADAPTATION_TARGET);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 15.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 3.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 3.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    // CompoundLikelihood
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(coalescent);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    likelihoods.clear();
    likelihoods.add(treeLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    posterior.setId(CompoundLikelihoodParser.POSTERIOR);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 500, false);
    loggers[0].add(posterior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(rateParameter);
    loggers[0].add(popSize);
    loggers[0].add(kappa);
    loggers[0].add(coalescent);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[1].add(posterior);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(rateParameter);
    loggers[1].add(coalescent);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("RandomLocalClockTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    // <expectation name="posterior" value="-3928.71"/>
    // <expectation name="clock.rate" value="8.04835E-4"/>
    // <expectation name="constant.popSize" value="37.3762"/>
    // <expectation name="hky.kappa" value="18.2782"/>
    // <expectation name="treeModel.rootHeight" value="69.0580"/>
    // <expectation name="treeLikelihood" value="-3856.59"/>
    // <expectation name="coalescent" value="-72.1285"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -3928.71);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -3856.59);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 69.0580);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 18.2782);
    TraceCorrelation rateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(StrictClockBranchRates.RATE));
    assertExpectation(StrictClockBranchRates.RATE, rateStats, 8.04835E-4);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 37.3762);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, -72.1285);
}

Example 4 with TreeIntervals

use of dr.evomodel.coalescent.TreeIntervals in project beast-mcmc by beast-dev.

the class RandomLocalClockTestProblem method testRandomLocalClock.

public void testRandomLocalClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 0.077, 0, Double.POSITIVE_INFINITY);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    TreeIntervals intervalList = new TreeIntervals(treeModel, null, null);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(intervalList, constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter ratesParameter = new Parameter.Default(RandomLocalClockModelParser.RATES, 10, 1);
    Parameter rateIndicatorParameter = new Parameter.Default(RandomLocalClockModelParser.RATE_INDICATORS, 10, 1);
    Parameter meanRateParameter = new Parameter.Default(RandomLocalClockModelParser.CLOCK_RATE, 1, 1.0);
    RandomLocalClockModel branchRateModel = new RandomLocalClockModel(treeModel, meanRateParameter, rateIndicatorParameter, ratesParameter, false, 0.5);
    SumStatistic rateChanges = new SumStatistic("rateChangeCount", true, null);
    rateChanges.addStatistic(rateIndicatorParameter);
    RateStatistic meanRate = new RateStatistic("meanRate", treeModel, branchRateModel, true, true, RateStatisticParser.MEAN);
    RateStatistic coefficientOfVariation = new RateStatistic(RateStatisticParser.COEFFICIENT_OF_VARIATION, treeModel, branchRateModel, true, true, RateStatisticParser.COEFFICIENT_OF_VARIATION);
    RateCovarianceStatistic covariance = new RateCovarianceStatistic("covariance", treeModel, branchRateModel);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, Double.POSITIVE_INFINITY);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    // siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    // treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(ratesParameter, 0.75);
    operator.setWeight(10.0);
    schedule.addOperator(operator);
    operator = new BitFlipOperator(rateIndicatorParameter, 15.0, true);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = ((DefaultTreeModel) treeModel).getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(((DefaultTreeModel) treeModel), 15.0, 0.0077, true, false, false, false, AdaptationMode.ADAPTATION_ON, AdaptableMCMCOperator.DEFAULT_ADAPTATION_TARGET);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 15.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 3.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 3.0);
    // operator.doOperation();
    schedule.addOperator(operator);
    // CompoundLikelihood
    OneOnXPrior likelihood1 = new OneOnXPrior();
    likelihood1.addData(popSize);
    OneOnXPrior likelihood2 = new OneOnXPrior();
    likelihood2.addData(kappa);
    DistributionLikelihood likelihood3 = new DistributionLikelihood(new GammaDistribution(0.5, 2.0), 0.0);
    likelihood3.addData(ratesParameter);
    DistributionLikelihood likelihood4 = new DistributionLikelihood(new PoissonDistribution(1.0), 0.0);
    likelihood4.addData(rateChanges);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(likelihood1);
    likelihoods.add(likelihood2);
    likelihoods.add(likelihood3);
    likelihoods.add(likelihood4);
    likelihoods.add(coalescent);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    likelihoods.clear();
    likelihoods.add(treeLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    posterior.setId(CompoundLikelihoodParser.POSTERIOR);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    loggers[0].add(posterior);
    loggers[0].add(prior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(kappa);
    // loggers[0].add(meanRate);
    loggers[0].add(rateChanges);
    loggers[0].add(coefficientOfVariation);
    loggers[0].add(covariance);
    loggers[0].add(popSize);
    loggers[0].add(coalescent);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[1].add(posterior);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(meanRate);
    loggers[1].add(rateChanges);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("RandomLocalClockTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    // <expectation name="posterior" value="-1818.26"/>
    // <expectation name="prior" value="-2.70143"/>
    // <expectation name="likelihood" value="-1815.56"/>
    // <expectation name="treeModel.rootHeight" value="6.363E-2"/>
    // <expectation name="constant.popSize" value="9.67405E-2"/>
    // <expectation name="hky.kappa" value="30.0394"/>
    // <expectation name="coefficientOfVariation" value="7.02408E-2"/>
    // covariance 0.47952
    // <expectation name="rateChangeCount" value="0.40786"/>
    // <expectation name="coalescent" value="7.29521"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -1818.26);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.PRIOR));
    assertExpectation(CompoundLikelihoodParser.PRIOR, likelihoodStats, -2.70143);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.56);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 6.363E-2);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 30.0394);
    TraceCorrelation rateChangeStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("rateChangeCount"));
    assertExpectation("rateChangeCount", rateChangeStats, 0.40786);
    TraceCorrelation coefficientOfVariationStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(RateStatisticParser.COEFFICIENT_OF_VARIATION));
    assertExpectation(RateStatisticParser.COEFFICIENT_OF_VARIATION, coefficientOfVariationStats, 7.02408E-2);
    TraceCorrelation covarianceStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("covariance"));
    assertExpectation("covariance", covarianceStats, 0.47952);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 9.67405E-2);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, 7.29521);
}

Example 5 with TreeIntervals

use of dr.evomodel.coalescent.TreeIntervals in project beast-mcmc by beast-dev.

the class BigFastTreeTreeIntervalsTest method testCompareIntervals.

public void testCompareIntervals() throws TreeUtils.MissingTaxonException, IOException, Importer.ImportException {
    NewickImporter importer = new NewickImporter("(Lishui/LS557/2020:0,((Netherlands/Utrecht_10015/2020:0.00006795400000000001,USA/IL-NM073/2020:0.00006799599999999999):0.000033976,England/LOND-D604F/2020:0.000101963):0.000033968,Guangdong/2020XN4459-P0041/2020:0.000000005,(Portugal/PT0063/2020:0,(Spain/Zaragoza2486/2020:0.000102605,Scotland/CVR746/2020:0.000000005,Spain/COV000882/2020:0.000067956,Colombia/INS-79253/2020:0.000101944,Uruguay/UY-4/2020:0.000031515):0.000033979,(Spain/CastillaLaMancha201329/2020:0.000000005,Netherlands/NoordHolland_10011/2020:0.000033987):0.00006799,England/LIVE-9CE87/2020:0.00013727299999999998,Spain/Granada-COV002916/2020:0.000033979999999999997):0.000033968,((USA/VI-CDC-3705/2020:0.000000005,Australia/VIC229/2020:0,USA/MA-MGH-00063/2020:0,(USA/WA-S41/2020:0.000068895,USA/WA-UW114/2020:0.000067978,USA/WA-UW17/2020:0.000000005,(USA/WA-S582/2020:0,USA/WA-UW-1682/2020:0.000000005,USA/WA-S994/2020:0.000101934):0.000033955,USA/WA-S121/2020:0.000000005,USA/WA-S154/2020:0.000067982,USA/WA-UW37/2020:0,USA/WA-S321/2020:0,USA/WA-S445/2020:0,USA/WA-S512/2020:0,USA/WA-S33/2020:0.000033979,Canada/BC_6981299/2020:0.000033972,USA/WA-UW-1294/2020:0.000033972,USA/WA-UW-2247/2020:0.000033988,Australia/VIC140/2020:0.000033984,USA/WA-UW61/2020:0.000033972,Canada/BC_8606204/2020:0.000166157,(USA/WA-S734/2020:0,USA/WA-S844/2020:0.000033983):0.000067965,(USA/WA-S1191/2020:0.000067947,USA/WA-S951/2020:0.000101914):0.000095803,Australia/NSW99/2020:0.000101953,(USA/WA-S317/2020:0.000000005,USA/WA-S721/2020:0.00003397):0.00010195700000000001,USA/WA-UW139/2020:0.000135916,USA/WA-S572/2020:0.000033979999999999997,USA/WA-S279/2020:0.000033972,USA/WA-UW28/2020:0.000034002,USA/WA-S114/2020:0.000033969,(USA/WA-S852/2020:0.000203899,(USA/WA-S568/2020:0,USA/WA-S791/2020:0.00006794599999999999):0.000033983):0.000101964,USA/WA-S842/2020:0.000067951):0.000033986,Singapore/302/2020:0.000101947):0.00016677,(((USA/IL-NM0112/2020:0.00003397,USA/IL-NM053/2020:0.000034229,USA/IL-NM059/2020:0.000101967):0.00003397,USA/WI-UW-218/2020:0.000033995):0.000030539,(USA/UT-QDX-63/2020:0,USA/CA-QDX-111/2020:0,USA/TX-HMH0427/2020:0.000203861):0.000101955):0.00023787300000000002):0.000033959,(((Scotland/CVR3203/2020:0.000000005,Scotland/CVR2246/2020:0.000000005,Scotland/GCVR-1714B2/2020:0.000033975999999999995,Scotland/CVR3514/2020:0.000068628):0.000067954,Australia/NT08/2020:0.000034000999999999995):0.00003397,Spain/COV001440/2020:0,Spain/Alcaniz2449/2020:0.000068985,Spain/COV001548/2020:0,USA/WI-WSLH-200057/2020:0.000000005,Spain/Valencia6/2020:0.0000343,Spain/Granada-COV002944/2020:0.000000005,Spain/COV001929/2020:0.000000005,Spain/COV002049/2020:0.000000005,(Spain/Valencia59/2020:0,Spain/Valencia306/2020:0.000000005):0.000033996,(Spain/COV001117/2020:0.00010265,Spain/COV002055/2020:0.000000005,England/20126000104/2020:0.00006758400000000001):0.000033997,Spain/COV001576/2020:0.000000005,Chile/Santiago-1/2020:0.000000005,Spain/COV000721/2020:0.000000005,(Spain/COV001575/2020:0,Spain/COV001505/2020:0):0.000067968,Spain/Madrid_H12_28/2020:0.000067957,Spain/COV001568/2020:0.000033975,England/CAMB-83357/2020:0.000068619,(Spain/Almeria-COV002842/2020:0.000000005,Spain/Malaga-COV002841/2020:0.000000005):0.000067851):0.000169854,Spain/Madrid_LP16_6193/2020:0.00006795299999999999,Singapore/51/2020:0.000044697,(Thailand/Nonthaburi_193/2020:0,Thailand/Bangkok_237/2020:0,Thailand/Bangkok_238/2020:0,((Thailand/Bangkok-0034/2020:0.000000005,Thailand/Bangkok_2295/2020:0,Thailand/Bangkok-0065/2020:0.000047826,Thailand/Bangkok-CONI-0147/2020:0.000033997):0.000033983,Thailand/SI202769-NT/2020:0.000203899):0.000101951):0.00006797500000000001,Shenzhen/SZTH-002/2020:0.000033999);");
    // NewickImporter importer = new NewickImporter("(((0:0.5,(1:1.0,2:1.0)n6:1.0)n7:1.0,3:1.5)n8:1.0,(4:2.0,5:1.51)n9:1.5)n10;");
    // NewickImporter constraintsImporter = new NewickImporter("(((0:0.5,(1:1.0,2:1.0)n6:1.0)n7:1.0,3:1.5)n8:1.0,(4:2.0,5:1.51)n9:1.5)n10;");
    tree = new DefaultTreeModel(importer.importTree(null));
    IntervalList intervals = new TreeIntervals(tree, null, null);
    BigFastTreeIntervals bigFastTreeIntervals = new BigFastTreeIntervals(tree);
    SubtreeLeapOperator op = new SubtreeLeapOperator(tree, 1, 0.0001, SubtreeLeapOperator.DistanceKernelType.NORMAL, AdaptationMode.ADAPTATION_OFF, 0.2);
    NodeHeightOperator nh = new NodeHeightOperator(tree, 1, 1, NodeHeightOperator.OperatorType.UNIFORM, AdaptationMode.ADAPTATION_OFF, 0.25);
    NodeHeightOperator root = new NodeHeightOperator(tree, 1, 0.75, NodeHeightOperator.OperatorType.SCALEROOT, AdaptationMode.ADAPTATION_OFF, 0.25);
    boolean pass = true;
    MathUtils.setSeed(2);
    for (int i = 0; i < 100000; i++) {
        op.doOperation();
        intervals.calculateIntervals();
        // bigFastIntervals.makeDirty();
        bigFastTreeIntervals.calculateIntervals();
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getInterval(j) != bigFastTreeIntervals.getInterval(j)) {
                System.out.println(i);
                System.out.println("interval wrong");
                pass = false;
                break;
            }
        }
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getLineageCount(j) != bigFastTreeIntervals.getLineageCount(j)) {
                System.out.println(i);
                System.out.println("lineage Counts wrong: " + j);
                System.out.println("expected: " + intervals.getLineageCount(j));
                System.out.println("got " + bigFastTreeIntervals.getLineageCount(j));
                pass = false;
                break;
            }
        }
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getIntervalTime(j) != bigFastTreeIntervals.getIntervalTime(j)) {
                System.out.println(i);
                System.out.println("times wrong");
                pass = false;
                break;
            }
        }
        if (!pass) {
            break;
        }
    }
    assertTrue(pass);
}