use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.
the class CheckPointTreeModifier method incorporateAdditionalTaxa.
/**
* Add the remaining taxa, which can be identified through the TreeDataLikelihood XML elements.
*/
public ArrayList<NodeRef> incorporateAdditionalTaxa(CheckPointUpdaterApp.UpdateChoice choice, BranchRates rateModel) {
System.out.println("Tree before adding taxa:\n" + treeModel.toString() + "\n");
ArrayList<NodeRef> newTaxaNodes = new ArrayList<NodeRef>();
for (String str : newTaxaNames) {
for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId().equals(str)) {
newTaxaNodes.add(treeModel.getExternalNode(i));
//always take into account Taxon dates vs. dates set through a TreeModel
System.out.println(treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " with height " + treeModel.getNodeHeight(treeModel.getExternalNode(i)) + " or " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getHeight());
}
}
}
System.out.println("newTaxaNodes length = " + newTaxaNodes.size());
ArrayList<Taxon> currentTaxa = new ArrayList<Taxon>();
for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
boolean taxonFound = false;
for (String str : newTaxaNames) {
if (str.equals((treeModel.getNodeTaxon(treeModel.getExternalNode(i))).getId())) {
taxonFound = true;
}
}
if (!taxonFound) {
System.out.println("Adding " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " to list of current taxa");
currentTaxa.add(treeModel.getNodeTaxon(treeModel.getExternalNode(i)));
}
}
System.out.println("Current taxa count = " + currentTaxa.size());
//iterate over both current taxa and to be added taxa
boolean originTaxon = true;
for (Taxon taxon : currentTaxa) {
if (taxon.getHeight() == 0.0) {
originTaxon = false;
System.out.println("Current taxon " + taxon.getId() + " has node height 0.0");
}
}
for (NodeRef newTaxon : newTaxaNodes) {
if (treeModel.getNodeTaxon(newTaxon).getHeight() == 0.0) {
originTaxon = false;
System.out.println("New taxon " + treeModel.getNodeTaxon(newTaxon).getId() + " has node height 0.0");
}
}
//check the Tree(Data)Likelihoods in the connected set of likelihoods
//focus on TreeDataLikelihood, which has getTree() to get the tree for each likelihood
//also get the DataLikelihoodDelegate from TreeDataLikelihood
ArrayList<TreeDataLikelihood> likelihoods = new ArrayList<TreeDataLikelihood>();
ArrayList<Tree> trees = new ArrayList<Tree>();
ArrayList<DataLikelihoodDelegate> delegates = new ArrayList<DataLikelihoodDelegate>();
for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
if (likelihood instanceof TreeDataLikelihood) {
likelihoods.add((TreeDataLikelihood) likelihood);
trees.add(((TreeDataLikelihood) likelihood).getTree());
delegates.add(((TreeDataLikelihood) likelihood).getDataLikelihoodDelegate());
}
}
//suggested to go through TreeDataLikelihoodParser and give it an extra option to create a HashMap
//keyed by the tree; am currently not overly fond of this approach
ArrayList<PatternList> patternLists = new ArrayList<PatternList>();
for (DataLikelihoodDelegate del : delegates) {
if (del instanceof BeagleDataLikelihoodDelegate) {
patternLists.add(((BeagleDataLikelihoodDelegate) del).getPatternList());
} else if (del instanceof MultiPartitionDataLikelihoodDelegate) {
MultiPartitionDataLikelihoodDelegate mpdld = (MultiPartitionDataLikelihoodDelegate) del;
List<PatternList> list = mpdld.getPatternLists();
for (PatternList pList : list) {
patternLists.add(pList);
}
}
}
if (patternLists.size() == 0) {
throw new RuntimeException("No patterns detected. Please make sure the XML file is BEAST 1.9 compatible.");
}
//aggregate all patterns to create distance matrix
//TODO What about different trees for different partitions?
Patterns patterns = new Patterns(patternLists.get(0));
if (patternLists.size() > 1) {
for (int i = 1; i < patternLists.size(); i++) {
patterns.addPatterns(patternLists.get(i));
}
}
//set the patterns for the distance matrix computations
choice.setPatterns(patterns);
//add new taxa one at a time
System.out.println("Adding " + newTaxaNodes.size() + " taxa ...");
for (NodeRef newTaxon : newTaxaNodes) {
treeModel.setNodeHeight(newTaxon, treeModel.getNodeTaxon(newTaxon).getHeight());
System.out.println("\nadding Taxon: " + newTaxon + " (height = " + treeModel.getNodeHeight(newTaxon) + ")");
//check if this taxon has a more recent sampling date than all other nodes in the current TreeModel
double offset = checkCurrentTreeNodes(newTaxon, treeModel.getRoot());
System.out.println("Sampling date offset when adding " + newTaxon + " = " + offset);
//AND set its current node height to 0.0 IF no originTaxon has been found
if (offset < 0.0) {
if (!originTaxon) {
System.out.println("Updating all node heights with offset " + Math.abs(offset));
updateAllTreeNodes(Math.abs(offset), treeModel.getRoot());
treeModel.setNodeHeight(newTaxon, 0.0);
}
} else if (offset == 0.0) {
if (!originTaxon) {
treeModel.setNodeHeight(newTaxon, 0.0);
}
}
//get the closest Taxon to the Taxon that needs to be added
//take into account which taxa can currently be chosen
Taxon closest = choice.getClosestTaxon(treeModel.getNodeTaxon(newTaxon), currentTaxa);
System.out.println("\nclosest Taxon: " + closest + " with original height: " + closest.getHeight());
//get the distance between these two taxa
double distance = choice.getDistance(treeModel.getNodeTaxon(newTaxon), closest);
System.out.println("at distance: " + distance);
//TODO what if distance == 0.0 ??? how to choose closest taxon then (in absence of geo info)?
//find the NodeRef for the closest Taxon (do not rely on node numbering)
NodeRef closestRef = null;
//careful with node numbering and subtract number of new taxa
for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)) == closest) {
closestRef = treeModel.getExternalNode(i);
}
}
System.out.println(closestRef + " with height " + treeModel.getNodeHeight(closestRef));
//System.out.println("trying to set node height: " + closestRef + " from " + treeModel.getNodeHeight(closestRef) + " to " + closest.getHeight());
//treeModel.setNodeHeight(closestRef, closest.getHeight());
double timeForDistance = distance / rateModel.getBranchRate(treeModel, closestRef);
System.out.println("timeForDistance = " + timeForDistance);
//get parent node of branch that will be split
NodeRef parent = treeModel.getParent(closestRef);
//determine height of new node
double insertHeight;
if (treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon)) {
insertHeight = treeModel.getNodeHeight(closestRef) + timeForDistance / 2.0;
System.out.println("treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon): " + insertHeight);
if (insertHeight >= treeModel.getNodeHeight(parent)) {
insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
}
} else {
double remainder = (timeForDistance - Math.abs(treeModel.getNodeHeight(closestRef) - treeModel.getNodeHeight(newTaxon))) / 2.0;
if (remainder > 0) {
insertHeight = Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)) + remainder;
System.out.println("remainder > 0: " + insertHeight);
if (insertHeight >= treeModel.getNodeHeight(parent)) {
insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
}
} else {
insertHeight = EPSILON * (treeModel.getNodeHeight(parent) - Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)));
insertHeight += Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon));
System.out.println("remainder <= 0: " + insertHeight);
}
}
System.out.println("insert at height: " + insertHeight);
//pass on all the necessary variables to a method that adds the new taxon to the tree
addTaxonAlongBranch(newTaxon, parent, closestRef, insertHeight);
//option to print tree after each taxon addition
System.out.println("\nTree after adding taxon " + newTaxon + ":\n" + treeModel.toString());
//add newly added Taxon to list of current taxa
currentTaxa.add(treeModel.getNodeTaxon(newTaxon));
}
return newTaxaNodes;
}
use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.
the class BeagleTreeLikelihood method main.
public static void main(String[] args) {
try {
MathUtils.setSeed(666);
System.out.println("Test case 1: simulateOnePartition");
int sequenceLength = 1000;
ArrayList<Partition> partitionsList = new ArrayList<Partition>();
// create tree
NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
Tree tree = importer.importTree(null);
TreeModel treeModel = new TreeModel(tree);
// create Frequency Model
Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
// create branch model
Parameter kappa1 = new Parameter.Default(1, 1);
Parameter kappa2 = new Parameter.Default(1, 1);
HKY hky1 = new HKY(kappa1, freqModel);
HKY hky2 = new HKY(kappa2, freqModel);
HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
substitutionModels.add(hky1);
substitutionModels.add(hky2);
List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
freqModels.add(freqModel);
Parameter epochTimes = new Parameter.Default(1, 20);
// create branch rate model
Parameter rate = new Parameter.Default(1, 0.001);
BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
// create site model
GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
// create partition
Partition partition1 = new //
Partition(//
treeModel, //
homogenousBranchSubstitutionModel, //
siteRateModel, //
branchRateModel, //
freqModel, // from
0, // to
sequenceLength - 1, // every
1);
partitionsList.add(partition1);
// feed to sequence simulator and generate data
BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
Alignment alignment = simulator.simulate(false, false);
BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
} catch (Exception e) {
e.printStackTrace();
System.exit(-1);
}
// END: try-catch block
}
use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.
the class TiterImporter method importTrees.
/**
* importTrees.
*/
public Tree[] importTrees(TaxonList taxonList) throws IOException, ImportException {
boolean done = false;
ArrayList<FlexibleTree> array = new ArrayList<FlexibleTree>();
do {
try {
skipUntil("(");
unreadCharacter('(');
FlexibleNode root = readInternalNode(taxonList);
FlexibleTree tree = new FlexibleTree(root, false, true);
array.add(tree);
if (taxonList == null) {
taxonList = tree;
}
if (readCharacter() != ';') {
throw new BadFormatException("Expecting ';' after tree");
}
} catch (EOFException e) {
done = true;
}
} while (!done);
Tree[] trees = new Tree[array.size()];
array.toArray(trees);
return trees;
}
use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.
the class OldTreesPanel method createTree.
private void createTree() {
if (generateTreeDialog == null) {
generateTreeDialog = new GenerateTreeDialog(frame);
}
int result = generateTreeDialog.showDialog(options);
if (result != JOptionPane.CANCEL_OPTION) {
GenerateTreeDialog.MethodTypes methodType = generateTreeDialog.getMethodType();
PartitionData partition = generateTreeDialog.getDataPartition();
Patterns patterns = new Patterns(partition.getAlignment());
DistanceMatrix distances = new F84DistanceMatrix(patterns);
Tree tree;
TemporalRooting temporalRooting;
switch(methodType) {
case NJ:
tree = new NeighborJoiningTree(distances);
temporalRooting = new TemporalRooting(tree);
tree = temporalRooting.findRoot(tree, TemporalRooting.RootingFunction.CORRELATION);
break;
case UPGMA:
tree = new UPGMATree(distances);
temporalRooting = new TemporalRooting(tree);
break;
default:
throw new IllegalArgumentException("unknown method type");
}
tree.setId(generateTreeDialog.getName());
options.userTrees.add(tree);
treesTableModel.fireTableDataChanged();
int row = options.userTrees.size() - 1;
treesTable.getSelectionModel().setSelectionInterval(row, row);
}
fireTreePriorsChanged();
}
use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.
the class TimeSlicer method readAndAnalyzeTrees.
// private List<Tree> importTrees(String treeFileName, int burnin) throws IOException, Importer.ImportException {
//
// int totalTrees = 10000;
//
// progressStream.println("Reading trees (bar assumes 10,000 trees)...");
// progressStream.println("0 25 50 75 100");
// progressStream.println("|--------------|--------------|--------------|--------------|");
//
// int stepSize = totalTrees / 60;
// if (stepSize < 1) stepSize = 1;
//
// List<Tree> treeList = new ArrayList<Tree>();
// BufferedReader reader1 = new BufferedReader(new FileReader(treeFileName));
//
// String line1 = reader1.readLine();
// TreeImporter importer1;
// if (line1.toUpperCase().startsWith("#NEXUS")) {
// importer1 = new NexusImporter(new FileReader(treeFileName));
// } else {
// importer1 = new NewickImporter(new FileReader(treeFileName));
// }
// totalTrees = 0;
// while (importer1.hasTree()) {
// Tree treeTime = importer1.importNextTree();
//
// if (totalTrees > burnin)
// treeList.add(treeTime);
//
// if (totalTrees > 0 && totalTrees % stepSize == 0) {
// progressStream.print("*");
// progressStream.flush();
// }
// totalTrees++;
// }
// return treeList;
// }
private void readAndAnalyzeTrees(String treeFileName, int burnin, int skipEvery, String[] traits, double[] slices, boolean impute, boolean trueNoise, Normalization normalize, boolean divideByBranchLength, BranchSet branchset, Set taxaSet) throws IOException, Importer.ImportException {
int totalTrees = 10000;
int totalStars = 0;
progressStream.println("Reading and analyzing trees (bar assumes 10,000 trees)...");
progressStream.println("0 25 50 75 100");
progressStream.println("|--------------|--------------|--------------|--------------|");
int stepSize = totalTrees / 60;
if (stepSize < 1)
stepSize = 1;
BufferedReader reader1 = new BufferedReader(new FileReader(treeFileName));
String line1 = reader1.readLine();
TreeImporter importer1;
if (line1.toUpperCase().startsWith("#NEXUS")) {
importer1 = new NexusImporter(new FileReader(treeFileName));
} else {
importer1 = new NewickImporter(new FileReader(treeFileName));
}
totalTrees = 0;
while (importer1.hasTree()) {
Tree treeTime = importer1.importNextTree();
if (totalTrees % skipEvery == 0) {
treesRead++;
if (totalTrees >= burnin) {
analyzeTree(treeTime, traits, slices, impute, trueNoise, normalize, divideByBranchLength, branchset, taxaSet);
}
}
if (totalTrees > 0 && totalTrees % stepSize == 0) {
progressStream.print("*");
totalStars++;
if (totalStars % 61 == 0)
progressStream.print("\n");
progressStream.flush();
}
totalTrees++;
}
progressStream.print("\n");
}
Aggregations