use of beast.evolution.tree.Tree in project beast2 by CompEvol.
the class BeautiDoc method deepCopyPlugin.
/**
* Create a deep copy of a beastObject, but in a different partition context
* First, find all beastObjects that are predecessors of the beastObject to be copied
* that are ancestors of StateNodes
*
* @param beastObject
* @param parent
* @return
*/
public static BEASTInterface deepCopyPlugin(BEASTInterface beastObject, BEASTInterface parent, MCMC mcmc, PartitionContext oldContext, PartitionContext newContext, BeautiDoc doc, List<BEASTInterface> tabooList) {
/**
* taboo = list of beastObjects that should not be copied *
*/
Set<BEASTInterface> taboo = new HashSet<>();
taboo.add(parent);
// add state
taboo.add(mcmc.startStateInput.get());
// add likelihood and prior
if (mcmc.posteriorInput.get() instanceof CompoundDistribution) {
for (Distribution distr : ((CompoundDistribution) mcmc.posteriorInput.get()).pDistributions.get()) {
if (distr instanceof CompoundDistribution) {
taboo.add(distr);
}
}
}
// add posterior
taboo.add(mcmc.posteriorInput.get());
// parent of operators
taboo.add(mcmc);
// add loggers
taboo.addAll(mcmc.loggersInput.get());
// add exception for *BEAST logger (perhaps need to be generalised?)
if (doc.pluginmap.containsKey("SpeciesTreeLoggerX")) {
taboo.add(doc.pluginmap.get("SpeciesTreeLoggerX"));
}
// add trees
for (StateNode node : mcmc.startStateInput.get().stateNodeInput.get()) {
if (node instanceof Tree) {
taboo.add(node);
}
}
// add MRCAPriors
for (String id : doc.pluginmap.keySet()) {
BEASTInterface o = doc.pluginmap.get(id);
if (o instanceof MRCAPrior) {
taboo.add(o);
}
}
if (tabooList != null) {
taboo.addAll(tabooList);
}
// find predecessors of beastObject to be copied
List<BEASTInterface> predecessors = new ArrayList<>();
collectPredecessors(beastObject, predecessors);
// find ancestors of StateNodes that are predecessors + the beastObject
// itself
Set<BEASTInterface> ancestors = new HashSet<>();
collectAncestors(beastObject, ancestors, taboo);
Log.info.print(Arrays.toString(ancestors.toArray()));
for (BEASTInterface beastObject2 : predecessors) {
if (beastObject2 instanceof StateNode) {
Set<BEASTInterface> ancestors2 = new HashSet<>();
collectAncestors(beastObject2, ancestors2, taboo);
ancestors.addAll(ancestors2);
} else if (beastObject2 instanceof Alignment || beastObject2 instanceof FilteredAlignment) {
for (Object output : beastObject2.getOutputs()) {
if (!taboo.contains(output)) {
Set<BEASTInterface> ancestors2 = new HashSet<>();
collectAncestors((BEASTInterface) output, ancestors2, taboo);
ancestors.addAll(ancestors2);
}
}
}
}
// collect priors
predecessors.addAll(ancestors);
for (BEASTInterface o : predecessors) {
if (o instanceof Prior) {
List<BEASTInterface> priorPredecessors = new ArrayList<>();
collectPredecessors(o, priorPredecessors);
ancestors.addAll(priorPredecessors);
}
}
Log.info.print(Arrays.toString(predecessors.toArray()));
for (BEASTInterface p : ancestors) {
Log.info.print("(");
try {
for (BEASTInterface p2 : p.listActiveBEASTObjects()) {
if (ancestors.contains(p2)) {
Log.info.print(p2.getID() + " ");
}
}
} catch (IllegalArgumentException e) {
e.printStackTrace();
}
Log.info.print(") ");
Log.info.println(p.getID());
}
// now the ancestors contain all beastObjects to be copied
// make a copy of all individual BEASTObjects, before connecting them up
Map<String, BEASTInterface> copySet = new HashMap<>();
for (BEASTInterface beastObject2 : ancestors) {
String id = beastObject2.getID();
if (id == null) {
id = beastObject.getClass().getName().replaceAll(".*\\.", "");
int i = 0;
while (doc.pluginmap.containsKey(id + "." + i)) {
i++;
}
id = id + "." + i;
beastObject2.setID(id);
}
String copyID = renameId(id, oldContext, newContext);
if (!id.equals(copyID)) {
if (doc.pluginmap.containsKey(copyID)) {
BEASTInterface org = doc.pluginmap.get(copyID);
copySet.put(id, org);
} else {
BEASTInterface copy;
try {
copy = beastObject2.getClass().newInstance();
copy.setID(copyID);
copySet.put(id, copy);
} catch (InstantiationException | IllegalAccessException e) {
e.printStackTrace();
throw new RuntimeException("Programmer error: every object in the model should have a default constructor that is publicly accessible");
}
}
}
Log.warning.println("Copy: " + id + " -> " + copyID);
}
// set all inputs of copied beastObjects + outputs to taboo
for (BEASTInterface beastObject2 : ancestors) {
String id = beastObject2.getID();
BEASTInterface copy = copySet.get(id);
if (copy != null) {
Log.warning.println("Processing: " + id + " -> " + copy.getID());
// set inputs
for (Input<?> input : beastObject2.listInputs()) {
if (input.get() != null) {
if (input.get() instanceof List) {
// ((List)copy.getInput(input.getName())).clear();
for (Object o : (List<?>) input.get()) {
if (o instanceof BEASTInterface) {
BEASTInterface value = getCopyValue((BEASTInterface) o, copySet, oldContext, newContext, doc);
// make sure it is not already in the list
Object o2 = copy.getInput(input.getName()).get();
boolean alreadyInList = false;
if (o2 instanceof List) {
List<?> currentList = (List<?>) o2;
for (Object v : currentList) {
if (v == value) {
alreadyInList = true;
break;
}
}
}
if (!alreadyInList) {
// add to the list
copy.setInputValue(input.getName(), value);
}
} else {
// it is a primitive value
if (copy instanceof Parameter.Base && input.getName().equals("value")) {
// // prevent appending to parameter values
Parameter.Base<?> p = ((Parameter.Base<?>) copy);
((List<?>) p.valuesInput.get()).clear();
}
copy.setInputValue(input.getName(), input.get());
}
}
} else if (input.get() instanceof BEASTInterface) {
// handle BEASTObject
BEASTInterface value = getCopyValue((BEASTInterface) input.get(), copySet, oldContext, newContext, doc);
copy.setInputValue(input.getName(), value);
} else if (input.get() instanceof String) {
// may need to replace partition info
String s = (String) input.get();
s = s.replaceAll("\\.c:[a-zA-Z0-9_]*", ".c:" + newContext.clockModel);
s = s.replaceAll("\\.s:[a-zA-Z0-9_]*", ".s:" + newContext.siteModel);
s = s.replaceAll("\\.t:[a-zA-Z0-9_]*", ".t:" + newContext.tree);
copy.setInputValue(input.getName(), s);
} else {
// it is a primitive value
copy.setInputValue(input.getName(), input.get());
}
}
}
// set outputs
for (Object output : beastObject2.getOutputs()) {
if (taboo.contains(output) && output != parent) {
BEASTInterface output2 = getCopyValue((BEASTInterface) output, copySet, oldContext, newContext, doc);
for (Input<?> input : ((BEASTInterface) output).listInputs()) {
// do not add state node initialisers automatically
if (input.get() instanceof List && // do not update state node initialisers
!(taboo.contains(output2) && input.getName().equals("init"))) {
List<?> list = (List<?>) input.get();
if (list.contains(beastObject2)) {
List<?> list2 = (List<?>) output2.getInput(input.getName()).get();
if (!list2.contains(copy)) {
output2.setInputValue(input.getName(), copy);
}
}
}
}
}
}
copySet.put(id, copy);
// Log.warning.println(base.operatorsAsString());
}
}
// deep copy must be obtained from copyset, before sorting
// since the sorting changes (deletes items) from the copySet map
BEASTInterface deepCopy = copySet.get(beastObject.getID());
// first need to sort copySet by topology, before we can initAndValidate
// them
List<BEASTInterface> sorted = new ArrayList<>();
Collection<BEASTInterface> values = copySet.values();
while (values.size() > 0) {
for (BEASTInterface copy : values) {
boolean found = false;
for (BEASTInterface beastObject2 : copy.listActiveBEASTObjects()) {
if (values.contains(beastObject2)) {
found = true;
break;
}
}
if (!found) {
sorted.add(copy);
}
}
values.remove(sorted.get(sorted.size() - 1));
}
// initialise copied beastObjects
Set<BEASTInterface> done = new HashSet<>();
for (BEASTInterface copy : sorted) {
try {
if (!done.contains(copy)) {
copy.initAndValidate();
done.add(copy);
}
} catch (Exception e) {
// ignore
Log.warning.print(e.getMessage());
}
if (doc != null) {
doc.addPlugin(copy);
}
}
doc.scrubAll(true, false);
return deepCopy;
}
use of beast.evolution.tree.Tree in project beast2 by CompEvol.
the class BeautiDoc method addTraitSet.
/**
* assigns trait to first available tree *
*/
void addTraitSet(TraitSet trait) {
if (trait != null) {
CompoundDistribution likelihood = (CompoundDistribution) pluginmap.get("likelihood");
for (Distribution d : likelihood.pDistributions.get()) {
if (d instanceof GenericTreeLikelihood) {
try {
// TODO: this might not be a valid type conversion from TreeInterface to Tree
Tree tree = (Tree) ((GenericTreeLikelihood) d).treeInput.get();
tree.m_traitList.setValue(trait, tree);
} catch (Exception e) {
e.printStackTrace();
}
scrubAll(true, false);
return;
}
}
}
}
use of beast.evolution.tree.Tree in project beast2 by CompEvol.
the class BeautiDoc method setClockRate.
void setClockRate() {
boolean needsEstimationBySPTree = false;
if (pluginmap.containsKey("Tree.t:Species")) {
Tree sptree = (Tree) pluginmap.get("Tree.t:Species");
// check whether there is a calibration
for (Object beastObject : sptree.getOutputs()) {
if (beastObject instanceof MRCAPrior) {
MRCAPrior prior = (MRCAPrior) beastObject;
if (prior.distInput.get() != null) {
needsEstimationBySPTree = true;
}
}
}
}
BEASTInterface likelihood = pluginmap.get("likelihood");
if (likelihood instanceof CompoundDistribution) {
int i = 0;
RealParameter firstClock = null;
for (Distribution distr : ((CompoundDistribution) likelihood).pDistributions.get()) {
if (distr instanceof GenericTreeLikelihood) {
GenericTreeLikelihood treeLikelihood = (GenericTreeLikelihood) distr;
boolean needsEstimation = needsEstimationBySPTree;
if (i > 0) {
BranchRateModel.Base model = treeLikelihood.branchRateModelInput.get();
needsEstimation = (model.meanRateInput.get() != firstClock) || firstClock.isEstimatedInput.get();
} else {
// TODO: this might not be a valid type conversion from TreeInterface to Tree
Tree tree = (Tree) treeLikelihood.treeInput.get();
// check whether there are tip dates
if (tree.hasDateTrait()) {
needsEstimation = true;
}
// check whether there is a calibration
for (Object beastObject : tree.getOutputs()) {
if (beastObject instanceof MRCAPrior) {
MRCAPrior prior = (MRCAPrior) beastObject;
if (prior.distInput.get() != null) {
needsEstimation = true;
}
}
}
}
BranchRateModel.Base model = treeLikelihood.branchRateModelInput.get();
if (model != null) {
RealParameter clockRate = model.meanRateInput.get();
clockRate.isEstimatedInput.setValue(needsEstimation, clockRate);
if (firstClock == null) {
firstClock = clockRate;
}
}
i++;
}
}
}
}
use of beast.evolution.tree.Tree in project beast2 by CompEvol.
the class BeautiDoc method scrubAll.
// TreeDistribution getTreePrior(String partition) {
// int k = 0;
// for (Alignment data : alignments) {
// if (data.getID().equals(partition)) {
// return treePriors.get(k);
// }
// k++;
// }
// return null;
// }
public synchronized void scrubAll(boolean useNotEstimatedStateNodes, boolean isInitial) {
try {
if (autoSetClockRate) {
setClockRate();
}
if (autoUpdateFixMeanSubstRate) {
SiteModelInputEditor.customConnector(this);
}
// }
if (pluginmap.containsKey("Tree.t:Species")) {
Tree tree = (Tree) pluginmap.get("Tree.t:Species");
tree.isEstimatedInput.setValue(true, tree);
}
// go through all templates, and process connectors in relevant ones
boolean progress = true;
while (progress) {
warning("============================ start scrubbing ===========================");
progress = false;
setUpActivePlugins();
// process MRCA priors
for (String id : pluginmap.keySet()) {
if (id != null && id.endsWith(".prior")) {
BEASTInterface beastObject = pluginmap.get(id);
if (beastObject instanceof MRCAPrior) {
MRCAPrior prior = (MRCAPrior) beastObject;
if (prior.treeInput.get().isEstimatedInput.get() == false) {
// disconnect
disconnect(beastObject, "prior", "distribution");
} else {
// connect
connect(beastObject, "prior", "distribution");
}
}
}
}
List<BeautiSubTemplate> templates = new ArrayList<>();
templates.add(beautiConfig.partitionTemplate.get());
templates.addAll(beautiConfig.subTemplates);
for (PartitionContext context : possibleContexts) {
applyBeautiRules(templates, isInitial, context);
}
// add 'Species' as special partition name
applyBeautiRules(templates, isInitial, new PartitionContext("Species"));
// if the model changed, some rules that use inposterior() may
// not have been triggered properly
// so we need to check that the model changed, and if so,
// revisit the BeautiConnectors
List<BEASTInterface> posteriorPredecessors2 = new ArrayList<>();
collectPredecessors(((MCMC) mcmc.get()).posteriorInput.get(), posteriorPredecessors2);
if (posteriorPredecessors.size() != posteriorPredecessors2.size()) {
progress = true;
} else {
for (BEASTInterface beastObject : posteriorPredecessors2) {
if (!posteriorPredecessors.contains(beastObject)) {
progress = true;
break;
}
}
}
}
List<BeautiSubTemplate> templates = new ArrayList<>();
templates.add(beautiConfig.hyperPriorTemplate);
for (BEASTInterface beastObject : pluginmap.values()) {
if (beastObject instanceof RealParameter) {
if (beastObject.getID() != null && beastObject.getID().startsWith("parameter.")) {
PartitionContext context = new PartitionContext(beastObject.getID().substring("parameter.".length()));
applyBeautiRules(templates, isInitial, context);
}
}
}
collectClockModels();
// collectTreePriors();
Log.warning.println("PARTITIONS:\n");
Log.warning.println(Arrays.toString(currentPartitions));
determineLinks();
} catch (Exception e) {
Log.err.println(e.getMessage());
}
}
use of beast.evolution.tree.Tree in project beast2 by CompEvol.
the class WilsonBaldingTest method topologyDistribution.
/**
* Test topology distribution.
* @throws Exception
*/
@Test
public void topologyDistribution() throws Exception {
// Fix seed: will hopefully ensure success of test unless something
// goes terribly wrong.
Randomizer.setSeed(42);
// Assemble model:
ConstantPopulation constantPop = new ConstantPopulation();
constantPop.initByName("popSize", new RealParameter("10000.0"));
List<Object> alignmentInitArgs = new ArrayList<Object>();
for (int i = 0; i < 4; i++) {
Sequence thisSeq = new Sequence();
thisSeq.initByName("taxon", String.valueOf(i), "value", "?");
alignmentInitArgs.add("sequence");
alignmentInitArgs.add(thisSeq);
}
Alignment alignment = new Alignment();
alignment.initByName(alignmentInitArgs.toArray());
Tree tree = new RandomTree();
tree.initByName("taxa", alignment, "populationModel", constantPop);
TreeIntervals treeIntervals = new TreeIntervals();
treeIntervals.initByName("tree", tree);
Coalescent coalescentDistrib = new Coalescent();
coalescentDistrib.initByName("treeIntervals", treeIntervals, "populationModel", constantPop);
// Set up state:
State state = new State();
state.initByName("stateNode", tree);
// Set up operator:
WilsonBalding wilsonBalding = new WilsonBalding();
wilsonBalding.initByName("weight", "1", "tree", tree);
// Set up logger:
TreeReport treeReport = new TreeReport();
treeReport.initByName("logEvery", "100", "burnin", "200000", "credibleSetPercentage", "95.0", "log", tree, "silent", true);
// Set up MCMC:
MCMC mcmc = new MCMC();
mcmc.initByName("chainLength", "2000000", "state", state, "distribution", coalescentDistrib, "operator", wilsonBalding, "logger", treeReport);
// Run MCMC:
mcmc.run();
// Obtain analysis results:
TreeTraceAnalysis analysis = treeReport.getAnalysis();
Map<String, Integer> topologyCounts = analysis.getTopologyCounts();
int totalTreesUsed = analysis.getNTrees();
// Test topology distribution against ideal:
double tol = 0.005;
for (int i = 0; i < topologies.length; i++) {
double thisProb = topologyCounts.get(topologies[i]) / (double) totalTreesUsed;
boolean withinTol = (thisProb > probs[i] - tol && thisProb < probs[i] + tol);
Assert.assertTrue(withinTol);
System.err.format("Topology %s rel. freq. %.3f", topologies[i], thisProb);
if (withinTol)
System.err.println(" (Within tolerance " + tol + " of " + String.valueOf(probs[i]) + ")");
else
System.err.println(" (FAILURE: outside tolerance " + tol + " of " + String.valueOf(probs[i]) + ")");
}
}
Aggregations