Examples with MCMC beast.core.MCMC used on opensource projects

Example 16 with MCMC

use of beast.core.MCMC in project MultiTypeTree by tgvaughan.

the class TWB_TS_Test method testTWB2.

@Test
public void testTWB2() throws Exception {
    System.out.println("TWB_test 2");
    // Fix seed.
    Randomizer.setSeed(42);
    // Assemble initial MultiTypeTree
    String newickStr = "((1[&deme=1]:1,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
    MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
    mtTree.initByName("value", newickStr, "typeLabel", "deme");
    // Assemble migration model:
    RealParameter rateMatrix = new RealParameter("0.1 0.1");
    RealParameter popSizes = new RealParameter("7.0 7.0");
    SCMigrationModel migModel = new SCMigrationModel();
    migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
    // Assemble distribution:
    StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
    distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree);
    // Set up state:
    State state = new State();
    state.initByName("stateNode", mtTree);
    // Set up operator:
    TypedWilsonBalding operatorTWB = new TypedWilsonBalding();
    operatorTWB.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "alpha", 0.2);
    Operator operatorMTTS = new MultiTypeTreeScale();
    operatorMTTS.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "scaleFactor", 0.8, "useOldTreeScaler", false);
    // Set up stat analysis logger:
    MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
    logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.1, "logEvery", 1000);
    // Set up MCMC:
    MCMC mcmc = new MCMC();
    mcmc.initByName("chainLength", "1000000", "state", state, "distribution", distribution, "operator", operatorTWB, "logger", logger);
    // Run MCMC:
    mcmc.run();
    System.out.format("height mean = %s\n", logger.getHeightMean());
    System.out.format("height var = %s\n", logger.getHeightVar());
    System.out.format("height ESS = %s\n", logger.getHeightESS());
    // Direct simulation:
    double[] heights = UtilMethods.getSimulatedHeights(migModel, new IntegerParameter("1 0 0"));
    double simHeightMean = DiscreteStatistics.mean(heights);
    double simHeightVar = DiscreteStatistics.variance(heights);
    // Compare analysis results with truth:
    boolean withinTol = (logger.getHeightESS() > 400) && (Math.abs(logger.getHeightMean() - simHeightMean) < 1.0) && (Math.abs(logger.getHeightVar() - simHeightVar) < 30);
    Assert.assertTrue(withinTol);
}

Example 17 with MCMC

use of beast.core.MCMC in project beast2 by CompEvol.

the class BeautiPanel method cloneFrom.

/**
 * Clones partition identified by sourceID to targetID and type (Site/Clock/Tree model)
 * as stored in config.
 * @param sourceID
 * @param targetID
 */
public void cloneFrom(String sourceID, String targetID) {
    if (sourceID.equals(targetID)) {
        return;
    }
    String type = config.hasPartitionsInput.get().toString();
    java.util.List<BEASTInterface> list = doc.getPartitions(type);
    int source = -1, target = -1;
    for (int i = 0; i < list.size(); i++) {
        BEASTInterface partition = list.get(i);
        if (type.equals("SiteModel")) {
            partition = (BEASTInterface) ((GenericTreeLikelihood) partition).siteModelInput.get();
        } else if (type.equals("ClockModel")) {
            partition = ((GenericTreeLikelihood) partition).branchRateModelInput.get();
        } else if (type.equals("Tree")) {
            partition = (BEASTInterface) ((GenericTreeLikelihood) partition).treeInput.get();
        }
        String partitionID = partition.getID();
        partitionID = partitionID.substring(partitionID.lastIndexOf('.') + 1);
        if (partitionID.length() > 1 && partitionID.charAt(1) == ':') {
            partitionID = partitionID.substring(2);
        }
        if (partitionID.equals(sourceID)) {
            source = i;
        }
        if (partitionID.equals(targetID)) {
            target = i;
        }
    }
    if (target == -1) {
        throw new RuntimeException("Programmer error: sourceID and targetID should be in list");
    }
    CompoundDistribution likelihoods = (CompoundDistribution) doc.pluginmap.get("likelihood");
    GenericTreeLikelihood likelihoodSource = (GenericTreeLikelihood) likelihoods.pDistributions.get().get(source);
    GenericTreeLikelihood likelihood = (GenericTreeLikelihood) likelihoods.pDistributions.get().get(target);
    PartitionContext oldContext = doc.getContextFor(likelihoodSource);
    PartitionContext newContext = doc.getContextFor(likelihood);
    // this ensures the config.sync does not set any input value
    config._input.setValue(null, config);
    if (type.equals("SiteModel")) {
        SiteModelInterface siteModelSource = likelihoodSource.siteModelInput.get();
        SiteModelInterface siteModel = null;
        try {
            siteModel = (SiteModel.Base) BeautiDoc.deepCopyPlugin((BEASTInterface) siteModelSource, likelihood, (MCMC) doc.mcmc.get(), oldContext, newContext, doc, null);
        } catch (RuntimeException e) {
            JOptionPane.showMessageDialog(this, "Could not clone " + sourceID + " to " + targetID + " " + e.getMessage());
            return;
        }
        likelihood.siteModelInput.setValue(siteModel, likelihood);
        return;
    } else if (type.equals("ClockModel")) {
        BranchRateModel clockModelSource = likelihoodSource.branchRateModelInput.get();
        BranchRateModel clockModel = null;
        try {
            clockModel = (BranchRateModel) BeautiDoc.deepCopyPlugin((BEASTInterface) clockModelSource, likelihood, (MCMC) doc.mcmc.get(), oldContext, newContext, doc, null);
        } catch (Exception e) {
            JOptionPane.showMessageDialog(this, "Could not clone " + sourceID + " to " + targetID + " " + e.getMessage());
            return;
        }
        // make sure that *if* the clock model has a tree as input, it is
        // the same as for the likelihood
        TreeInterface tree = null;
        try {
            for (Input<?> input : ((BEASTInterface) clockModel).listInputs()) {
                if (input.getName().equals("tree")) {
                    tree = (TreeInterface) input.get();
                }
            }
        } catch (IllegalArgumentException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
        if (tree != null && tree != likelihood.treeInput.get()) {
            // likelihood.treeInput.setValue(tree, likelihood);
            JOptionPane.showMessageDialog(null, "Cannot clone clock model with different trees");
            return;
        }
        likelihood.branchRateModelInput.setValue(clockModel, likelihood);
        return;
    } else if (type.equals("Tree")) {
        TreeInterface tree = null;
        TreeInterface treeSource = likelihoodSource.treeInput.get();
        try {
            tree = (TreeInterface) BeautiDoc.deepCopyPlugin((BEASTInterface) treeSource, likelihood, (MCMC) doc.mcmc.get(), oldContext, newContext, doc, null);
        } catch (Exception e) {
            JOptionPane.showMessageDialog(this, "Could not clone " + sourceID + " to " + targetID + " " + e.getMessage());
            return;
        }
        // sanity check: make sure taxon sets are compatible
        Taxon.assertSameTaxa(tree.getID(), tree.getTaxonset().getTaxaNames(), likelihood.dataInput.get().getID(), likelihood.dataInput.get().getTaxaNames());
        likelihood.treeInput.setValue(tree, likelihood);
        return;
    } else {
        throw new RuntimeException("Programmer error calling cloneFrom: Should only clone Site/Clock/Tree model");
    }
}

Example 18 with MCMC

use of beast.core.MCMC in project beast2 by CompEvol.

the class BeautiBase method assertOperatorsEqual.

void assertOperatorsEqual(String... ids) {
    System.err.println("assertOperatorsEqual");
    MCMC mcmc = (MCMC) doc.mcmc.get();
    List<Operator> operators = mcmc.operatorsInput.get();
    asserListsEqual(operators, ids);
}

Example 19 with MCMC

use of beast.core.MCMC in project beast2 by CompEvol.

the class BeautiDoc method addMRCAPrior.

public void addMRCAPrior(MRCAPrior mrcaPrior) {
    Tree tree = (Tree) pluginmap.get("Tree.t:" + alignments.get(0).getID());
    if (tree == null) {
        for (String key : pluginmap.keySet()) {
            if (key.startsWith("Tree.t:")) {
                tree = (Tree) pluginmap.get(key);
                break;
            }
        }
    }
    // make sure we have the appropriate tree:
    if (alignments.size() > 1) {
        String testTaxon = mrcaPrior.taxonsetInput.get().toString().split("\n")[1].trim();
        // = tree.getTaxaNames();
        String[] taxaNames;
        int index = -1;
        int j = 0;
        while (index < 0 && j++ < alignments.size()) {
            tree = (Tree) pluginmap.get("Tree.t:" + alignments.get(j - 1).getID());
            taxaNames = tree.getTaxaNames();
            for (int i = 0; i < taxaNames.length && index < 0; i++) if (testTaxon.equals(taxaNames[i]))
                index = i;
        }
    }
    CompoundDistribution prior = (CompoundDistribution) pluginmap.get("prior");
    mrcaPrior.treeInput.setValue(tree, mrcaPrior);
    ParametricDistribution distr = mrcaPrior.distInput.get();
    TaxonSet t = mrcaPrior.taxonsetInput.get();
    if (taxaset.keySet().contains(t.getID())) {
        Log.warning.println("taxonset " + t.getID() + " already exists: MRCAPrior " + mrcaPrior.getID() + " can not be added");
    } else {
        taxaset.put(t.getID(), t);
        // ensure TaxonSets are not duplicated
        List<Taxon> taxa = t.taxonsetInput.get();
        for (int i = 0; i < taxa.size(); i++) {
            if (taxaset.containsKey(taxa.get(i).getID())) {
                taxa.set(i, taxaset.get(taxa.get(i).getID()));
            } else {
                taxaset.put(taxa.get(i).getID(), taxa.get(i));
            }
        }
        if (distr instanceof Normal && (Double.isInfinite(((Normal) distr).sigmaInput.get().getValue()))) {
        // it is a 'fixed' calibration, no need to add a distribution
        } else {
            prior.pDistributions.setValue(mrcaPrior, prior);
            connect(mrcaPrior, "tracelog", "log");
        }
    }
    if (t.taxonsetInput.get().size() == 1 && distr != null) {
        // only add operators if it is NOT a 'fixed' calibration
        if (!(distr instanceof Normal && (Double.isInfinite(((Normal) distr).sigmaInput.get().getValue())))) {
            TipDatesRandomWalker operator = new TipDatesRandomWalker();
            t.initAndValidate();
            operator.initByName("taxonset", t, "weight", 1.0, "tree", tree, "windowSize", 1.0);
            operator.setID("TipDatesRandomWalker." + t.getID());
            MCMC mcmc = (MCMC) this.mcmc.get();
            mcmc.operatorsInput.setValue(operator, mcmc);
        }
        // set up date trait
        double date = distr.getMean();
        TraitSet dateTrait = null;
        for (TraitSet ts : tree.m_traitList.get()) {
            if (ts.isDateTrait()) {
                dateTrait = ts;
            }
        }
        if (dateTrait == null) {
            dateTrait = new TraitSet();
            dateTrait.setID("traitsetDate");
            registerPlugin(dateTrait);
            dateTrait.initByName("traitname", TraitSet.DATE_BACKWARD_TRAIT, "taxa", tree.getTaxonset(), "value", t.taxonsetInput.get().get(0).getID() + "=" + date);
            tree.m_traitList.setValue(dateTrait, tree);
            tree.initAndValidate();
        } else {
            dateTrait.traitsInput.setValue(dateTrait.traitsInput.get() + ",\n" + t.taxonsetInput.get().get(0).getID() + "=" + date, dateTrait);
        }
        mrcaPrior.onlyUseTipsInput.setValue(true, mrcaPrior);
    }
}

Example 20 with MCMC

use of beast.core.MCMC in project beast2 by CompEvol.

the class BeautiDoc method setUpActivePlugins.

// scrubAll
protected void setUpActivePlugins() {
    posteriorPredecessors = new ArrayList<>();
    collectPredecessors(((MCMC) mcmc.get()).posteriorInput.get(), posteriorPredecessors);
    likelihoodPredecessors = new ArrayList<>();
    if (pluginmap.containsKey("likelihood")) {
        collectPredecessors(pluginmap.get("likelihood"), likelihoodPredecessors);
    }
    Log.trace.print("InPosterior=");
    for (BEASTInterface o : posteriorPredecessors) {
        pluginmap.put(o.getID(), o);
        Log.trace.print(o.getID() + " ");
    // if (!pluginmap.containsKey(o)) {
    // System.err.println("MISSING: " + o.getID());
    // }
    }
    Log.trace.println();
}