Examples with ExchangeOperator dr.evomodel.operators.ExchangeOperator used on opensource projects

Example 1 with ExchangeOperator

use of dr.evomodel.operators.ExchangeOperator in project beast-mcmc by beast-dev.

the class ExchangeOperatorTest method getOperatorSchedule.

// STATIC METHODS
public OperatorSchedule getOperatorSchedule(TreeModel treeModel) {
    Parameter rootParameter = treeModel.createNodeHeightsParameter(true, false, false);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    ExchangeOperator operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
    ScaleOperator scaleOperator = new ScaleOperator(rootParameter, 0.75, CoercionMode.COERCION_ON, 1.0);
    UniformOperator uniformOperator = new UniformOperator(internalHeights, 1.0);
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    schedule.addOperator(operator);
    schedule.addOperator(scaleOperator);
    schedule.addOperator(uniformOperator);
    return schedule;
}

Example 2 with ExchangeOperator

use of dr.evomodel.operators.ExchangeOperator in project beast-mcmc by beast-dev.

the class NarrowExchangeTest method getOperatorSchedule.

public OperatorSchedule getOperatorSchedule(DefaultTreeModel treeModel) {
    Parameter rootParameter = treeModel.createNodeHeightsParameter(true, false, false);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    ExchangeOperator operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1.0);
    ScaleOperator scaleOperator = new ScaleOperator(rootParameter, 0.75, AdaptationMode.ADAPTATION_ON, 1.0);
    UniformOperator uniformOperator = new UniformOperator(internalHeights, 1.0);
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    schedule.addOperator(operator);
    schedule.addOperator(scaleOperator);
    schedule.addOperator(uniformOperator);
    return schedule;
}

Example 3 with ExchangeOperator

use of dr.evomodel.operators.ExchangeOperator in project beast-mcmc by beast-dev.

the class NarrowExchangeTest method testNarrow.

/**
 * Test method for {@link dr.evomodel.operators.ExchangeOperator#narrow()}.
 */
public void testNarrow() throws IOException, ImportException {
    // probability of picking B node is 1/(2n-4) = 1/6
    // probability of swapping it with C is 1/1
    // total = 1/6
    System.out.println("Test 1: Forward");
    String treeMatch = "((((A,C),B),D),E);";
    int count = 0;
    int reps = 100000;
    for (int i = 0; i < reps; i++) {
        DefaultTreeModel treeModel = new DefaultTreeModel("treeModel", tree5);
        ExchangeOperator operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1);
        operator.doOperation();
        String tree = TreeUtils.newickNoLengths(treeModel);
        if (tree.equals(treeMatch)) {
            count += 1;
        }
    }
    double p_1 = (double) count / (double) reps;
    System.out.println("Number of proposals:\t" + count);
    System.out.println("Number of tries:\t" + reps);
    System.out.println("Number of ratio:\t" + p_1);
    System.out.println("Number of expected ratio:\t" + 1.0 / 6.0);
    assertExpectation(1.0 / 6.0, p_1, reps);
}

Example 4 with ExchangeOperator

use of dr.evomodel.operators.ExchangeOperator 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 5 with ExchangeOperator

use of dr.evomodel.operators.ExchangeOperator in project beast-mcmc by beast-dev.

the class ExchangeOperatorTest method testWideExchangeOperator2.

public void testWideExchangeOperator2() throws IOException, ImportException {
    // probability of picking (A,B) node is 1/(2n-2) = 1/8
    // probability of swapping with D is 1/2
    // total = 1/16
    // probability of picking {D} node is 1/(2n-2) = 1/8
    // probability of picking {A,B} is 1/5
    // total = 1/40
    // total = 1/16 + 1/40 = 0.0625 + 0.025 = 0.0875
    // new test:
    // probability of picking (A,B) node is 1/(2n-2) = 1/8
    // probability of swapping with D is 1/(2n-3) = 1/7
    // total = 1/56
    // probability of picking {D} node is 1/(2n-2) = 1/8
    // probability of picking {A,B} is 1/(2n-3) = 1/7
    // total = 1/56
    // total = 1/56 + 1/56 = 1/28
    System.out.println("Test 1: Forward");
    String treeMatch = "(((D,C),(A,B)),E);";
    int count = 0;
    int reps = 1000000;
    for (int i = 0; i < reps; i++) {
        DefaultTreeModel treeModel = new DefaultTreeModel("treeModel", tree5);
        ExchangeOperator operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
        operator.doOperation();
        String tree = TreeUtils.newickNoLengths(treeModel);
        if (tree.equals(treeMatch)) {
            count += 1;
        }
    }
    double p_1 = (double) count / (double) reps;
    System.out.println("Number of proposals:\t" + count);
    System.out.println("Number of tries:\t" + reps);
    System.out.println("Number of ratio:\t" + p_1);
    System.out.println("Number of expected ratio:\t" + 1.0 / 28.0);
    assertExpectation(1.0 / 28.0, p_1, reps);
    // since this operator is supposed to be symmetric it got a hastings ratio of one
    // this means, it should propose the same move just backwards with the same probability
    // BUT:
    // (((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0) -> ((((A,B),C),D),E)
    // probability of picking (A,B) node is 1/(2n-2) = 1/8
    // probability of swapping with D is 1/3
    // total = 1/24
    // probability of picking {D} node is 1/(2n-2) = 1/8
    // probability of picking {A,B} is 1/4
    // total = 1/32
    // total = 1/24 + 1/32 = 7/96 = 0.07291666666
    // new test:
    // probability of picking (A,B) node is 1/(2n-2) = 1/8
    // probability of swapping with D is 1/(2n-3) = 1/7
    // total = 1/56
    // probability of picking {D} node is 1/(2n-2) = 1/8
    // probability of picking {A,B} is 1/(2n-3) = 1/7
    // total = 1/56
    // total = 1/56 + 1/56 = 1/28
    System.out.println("Test 2: Backward");
    treeMatch = "((((A,B),C),D),E);";
    NewickImporter importer = new NewickImporter("(((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0);");
    FlexibleTree tree5_2 = (FlexibleTree) importer.importTree(null);
    count = 0;
    for (int i = 0; i < reps; i++) {
        DefaultTreeModel treeModel = new DefaultTreeModel("treeModel", tree5_2);
        ExchangeOperator operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
        operator.doOperation();
        String tree = TreeUtils.newickNoLengths(treeModel);
        if (tree.equals(treeMatch)) {
            count += 1;
        }
    }
    double p_2 = (double) count / (double) reps;
    System.out.println("Number of proposals:\t" + count);
    System.out.println("Number of tries:\t" + reps);
    System.out.println("Number of ratio:\t" + p_2);
    System.out.println("Number of expected ratio:\t" + 1.0 / 28.0);
    assertExpectation(1.0 / 28.0, p_2, reps);
}