#! /usr/bin/env python
# -*- coding: utf-8 -*-
##############################################################################
## DendroPy Phylogenetic Computing Library.
##
## Copyright 2010-2015 Jeet Sukumaran and Mark T. Holder.
## All rights reserved.
##
## See "LICENSE.rst" for terms and conditions of usage.
##
## If you use this work or any portion thereof in published work,
## please cite it as:
##
## Sukumaran, J. and M. T. Holder. 2010. DendroPy: a Python library
## for phylogenetic computing. Bioinformatics 26: 1569-1571.
##
##############################################################################
"""
Statistics, metrics, measurements, and values calculated *between* *two* trees.
"""
import math
import collections
import itertools
from dendropy.utility import error
###############################################################################
## Public Functions
[docs]
def symmetric_difference(tree1, tree2, is_bipartitions_updated=False):
"""
Returns *unweighted* Robinson-Foulds distance between two trees.
Trees need to share the same |TaxonNamespace| reference. The
bipartition bitmasks of the trees must be correct for the current tree
structures (by calling :meth:`Tree.encode_bipartitions()` method) or the
``is_bipartitions_updated`` argument must be |False| to force recalculation
of bipartitions.
Parameters
----------
tree1 : |Tree| object
The first tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree2`` and must have
bipartitions encoded.
tree2 : |Tree| object
The second tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree1`` and must have
bipartitions encoded.
is_bipartitions_updated : bool
If |False|, then the bipartitions on *both* trees will be updated
before comparison. If |True| then the bipartitions will only be
calculated for a |Tree| object if they have not been calculated
before, either explicitly or implicitly.
Returns
-------
d : int
The symmetric difference (a.k.a. the unweighted Robinson-Foulds
distance) between ``tree1`` and ``tree2``.
Examples
--------
::
from dendropy import TaxonNamespace, Tree
from dendropy.calculate import treecompare
tns = TaxonNamespace()
tree1 = Tree.get_from_path(
"t1.nex",
"nexus",
taxon_namespace=tns)
tree2 = Tree.get_from_path(
"t2.nex",
"nexus",
taxon_namespace=tns)
tree1.encode_bipartitions()
tree2.encode_bipartitions()
print(treecompare.symmetric_difference(tree1, tree2))
"""
t = false_positives_and_negatives(
tree1,
tree2,
is_bipartitions_updated=is_bipartitions_updated)
return t[0] + t[1]
[docs]
def unweighted_robinson_foulds_distance(tree1, tree2, is_bipartitions_updated=False):
"""
Alias for ``symmetric_difference()``.
"""
return symmetric_difference(tree1, tree2, is_bipartitions_updated)
[docs]
def weighted_robinson_foulds_distance(
tree1,
tree2,
edge_weight_attr="length",
is_bipartitions_updated=False):
"""
Returns *weighted* Robinson-Foulds distance between two trees based on
``edge_weight_attr``.
Trees need to share the same |TaxonNamespace| reference. The
bipartition bitmasks of the trees must be correct for the current tree
structures (by calling :meth:`Tree.encode_bipartitions()` method) or the
``is_bipartitions_updated`` argument must be |False| to force recalculation of
bipartitions.
Parameters
----------
tree1 : |Tree| object
The first tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree2`` and must have
bipartitions encoded.
tree2 : |Tree| object
The second tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree1`` and must have
bipartitions encoded.
edge_weight_attr : string
Name of attribute on edges of trees to be used as the weight.
is_bipartitions_updated : bool
If |True|, then the bipartitions on *both* trees will be updated before
comparison. If |False| (default) then the bipartitions will only be
calculated for a |Tree| object if they have not been calculated
before, either explicitly or implicitly.
Returns
-------
d : float
The edge-weighted Robinson-Foulds distance between ``tree1`` and ``tree2``.
Examples
--------
::
import dendropy
from dendropy.calculate import treecompare
tns = dendropy.TaxonNamespace()
tree1 = tree.get_from_path(
"t1.nex",
"nexus",
taxon_namespace=tns)
tree2 = tree.get_from_path(
"t2.nex",
"nexus",
taxon_namespace=tns)
tree1.encode_bipartitions()
tree2.encode_bipartitions()
print(treecompare.weighted_robinson_foulds_distance(tree1, tree2))
"""
df = lambda length_diffs: sum([abs(i[0] - i[1]) for i in length_diffs])
return _bipartition_difference(tree1,
tree2,
dist_fn=df,
edge_weight_attr=edge_weight_attr,
value_type=float,
is_bipartitions_updated=is_bipartitions_updated)
[docs]
def false_positives_and_negatives(
reference_tree,
comparison_tree,
is_bipartitions_updated=False):
"""
Counts and returns number of false positive bipar (bipartitions found in
``comparison_tree`` but not in ``reference_tree``) and false negative
bipartitions (bipartitions found in ``reference_tree`` but not in
``comparison_tree``).
Trees need to share the same |TaxonNamespace| reference. The
bipartition bitmasks of the trees must be correct for the current tree
structures (by calling :meth:`Tree.encode_bipartitions()` method) or the
``is_bipartitions_updated`` argument must be |False| to force recalculation of
bipartitions.
Parameters
----------
reference_tree : |Tree| object
The first tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree2`` and must have
bipartitions encoded.
comparison_tree : |Tree| object
The second tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree1`` and must have
bipartitions encoded.
is_bipartitions_updated : bool
If |True|, then the bipartitions on *both* trees will be updated
before comparison. If |False| (default) then the bipartitions
will only be calculated for a |Tree| object if they have not been
calculated before, either explicitly or implicitly.
Returns
-------
t : tuple(int)
A pair of integers, with first integer being the number of false
positives and the second being the number of false negatives.
Examples
--------
::
import dendropy
from dendropy.calculate import treecompare
tns = dendropy.TaxonNamespace()
tree1 = tree.get_from_path(
"t1.nex",
"nexus",
taxon_namespace=tns)
tree2 = tree.get_from_path(
"t2.nex",
"nexus",
taxon_namespace=tns)
tree1.encode_bipartitions()
tree2.encode_bipartitions()
print(treecompare.false_positives_and_negatives(tree1, tree2))
"""
if reference_tree.taxon_namespace is not comparison_tree.taxon_namespace:
raise error.TaxonNamespaceIdentityError(reference_tree, comparison_tree)
if not is_bipartitions_updated:
reference_tree.encode_bipartitions()
comparison_tree.encode_bipartitions()
else:
if reference_tree.bipartition_encoding is None:
reference_tree.encode_bipartitions()
if comparison_tree.bipartition_encoding is None:
comparison_tree.encode_bipartitions()
ref_bipartitions = set(reference_tree.bipartition_encoding)
comparison_bipartitions = set(comparison_tree.bipartition_encoding)
false_positives = comparison_bipartitions.difference(ref_bipartitions)
false_negatives = ref_bipartitions.difference(comparison_bipartitions)
return len(false_positives), len(false_negatives)
[docs]
def euclidean_distance(
tree1,
tree2,
edge_weight_attr="length",
value_type=float,
is_bipartitions_updated=False):
"""
Returns the Euclidean distance (a.k.a. Felsenstein's 2004 "branch length
distance") between two trees based on ``edge_weight_attr``.
Trees need to share the same |TaxonNamespace| reference. The
bipartition bitmasks of the trees must be correct for the current tree
structures (by calling :meth:`Tree.encode_bipartitions()` method) or the
``is_bipartitions_updated`` argument must be |False| to force recalculation of
bipartitions.
Parameters
----------
tree1 : |Tree| object
The first tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree2`` and must have
bipartitions encoded.
tree2 : |Tree| object
The second tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree1`` and must have
bipartitions encoded.
edge_weight_attr : string
Name of attribute on edges of trees to be used as the weight.
is_bipartitions_updated : bool
If |True|, then the bipartitions on *both* trees will be updated
before comparison. If |False| (default) then the bipartitions
will only be calculated for a |Tree| object if they have not been
calculated before, either explicitly or implicitly.
Returns
-------
d : int
The Euclidean distance between ``tree1`` and ``tree2``.
Examples
--------
::
import dendropy
from dendropy.calculate import treecompare
tns = dendropy.TaxonNamespace()
tree1 = tree.get_from_path(
"t1.nex",
"nexus",
taxon_namespace=tns)
tree2 = tree.get_from_path(
"t2.nex",
"nexus",
taxon_namespace=tns)
tree1.encode_bipartitions()
tree2.encode_bipartitions()
print(treecompare.euclidean_distance(tree1, tree2))
"""
df = lambda length_diffs: math.sqrt(sum([pow(i[0] - i[1], 2) for i in length_diffs]))
return _bipartition_difference(tree1,
tree2,
dist_fn=df,
edge_weight_attr=edge_weight_attr,
value_type=value_type,
is_bipartitions_updated=is_bipartitions_updated)
[docs]
def find_missing_bipartitions(reference_tree, comparison_tree, is_bipartitions_updated=False):
"""
Returns a list of bipartitions that are in ``reference_tree``, but
not in ``comparison_tree``.
Trees need to share the same |TaxonNamespace| reference. The
bipartition bitmasks of the trees must be correct for the current tree
structures (by calling :meth:`Tree.encode_bipartitions()` method) or the
``is_bipartitions_updated`` argument must be |False| to force recalculation of
bipartitions.
Parameters
----------
reference_tree : |Tree| object
The first tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree2`` and must have
bipartitions encoded.
comparison_tree : |Tree| object
The second tree of the two trees being compared. This must share the
same |TaxonNamespace| reference as ``tree1`` and must have
bipartitions encoded.
is_bipartitions_updated : bool
If |True|, then the bipartitions on *both* trees will be updated
before comparison. If |False| (default) then the bipartitions
will only be calculated for a |Tree| object if they have not been
calculated before, either explicitly or implicitly.
Returns
-------
s : list[|Bipartition|]
A list of bipartitions that are in the first tree but not in the second.
"""
missing = []
if reference_tree.taxon_namespace is not comparison_tree.taxon_namespace:
raise error.TaxonNamespaceIdentityError(reference_tree, comparison_tree)
if not is_bipartitions_updated:
reference_tree.encode_bipartitions()
comparison_tree.encode_bipartitions()
else:
if reference_tree.bipartition_encoding is None:
reference_tree.encode_bipartitions()
if comparison_tree.bipartition_encoding is None:
comparison_tree.encode_bipartitions()
for bipartition in reference_tree.bipartition_encoding:
if bipartition in comparison_tree.bipartition_encoding:
pass
else:
missing.append(bipartition)
return missing
##############################################################################
### TreeshapeKernel
class TreeShapeKernel(object):
_TreeShapeKernelNodeCache = collections.namedtuple("_TreeShapeKernelNodeCache",
["production", "index", "edge_lengths", "sum_of_square_edge_lengths"])
def __init__(self, **kwargs):
"""
Calculator for tree shape kernel tricking.
References
----------
[1] Poon, A. F., Pond, S. L. K., Bennett, P., Richman, D. D., Brown, A.
J. L., & Frost, S. D. (2007). Adaptation to human populations is
revealed by within-host polymorphisms in HIV-1 and hepatitis C virus.
PLoS Pathog, 3(3), e45.
[2] Poon, A. F., Walker, L. W., Murray, H., McCloskey, R. M., Harrigan,
P. R., & Liang, R. H. (2013). Mapping the shapes of phylogenetic trees
from human and zoonotic RNA viruses. PLoS one, 8(11), e78122.
[3] Poon, A. F., Walker, L. W., Murray, H., McCloskey, R. M., Harrigan,
P. R., & Liang, R. H. (2013). Mapping the shapes of phylogenetic trees
from human and zoonotic RNA viruses. PLoS one, 8(11), e78122.
[4] Poon, A. F. (2015). Phylodynamic inference with kernel ABC and its
application to HIV epidemiology. Molecular biology and evolution,
msv123.
"""
# kernel function
# sigma=1,
# gauss_factor=1,
# decay_factor=0.1,
self.sigma = kwargs.pop("sigma", 1)
self.gauss_factor = kwargs.pop("gauss_factor", 1)
self.decay_factor = kwargs.pop("decay_factor", 0.1)
# cache management
self._tree_cache = {}
def remove_from_cache(self, tree):
del self._tree_cache[tree]
def update_cache(self, tree):
"""
Pre-computes values needed for the kernel trick with this tree and
caches them.
"""
current_tree_cache = {}
for nd in tree.leaf_node_iter():
current_tree_cache[nd] = TreeShapeKernel._TreeShapeKernelNodeCache(
production=0,
index=0,
edge_lengths=None,
sum_of_square_edge_lengths=0)
nd.production = 0
for nd_idx, nd in enumerate(tree.postorder_internal_node_iter()):
nterms = 0
edge_lengths = []
for ch in nd.child_node_iter():
if current_tree_cache[ch].production == 0:
nterms += 1
edge_lengths.append(ch.edge.length)
production = nterms + 1
index = nd_idx
sum_of_square_edge_lengths = sum([elen**2 for elen in edge_lengths])
current_tree_cache[nd] = TreeShapeKernel._TreeShapeKernelNodeCache(
production=production,
index=index,
edge_lengths=edge_lengths,
sum_of_square_edge_lengths=sum_of_square_edge_lengths)
self._tree_cache[tree] = current_tree_cache
return current_tree_cache
def __call__(self,
tree1,
tree2,
is_tree1_cache_updated=False,
is_tree2_cache_updated=False,
):
"""
Recursive function for computing tree convolution
kernel.
Parameters
----------
tree1 : |Tree| instance
First tree to be compared. If it has already been seen by self, its
values will have been cached. If the tree has changed since it has been
seen by self, it will need to be recached, either explicitly before
the calculation by calling 'update_cache' or by specifying
'is_tree1_cache_updated=False'
tree2 : |Tree| instance
Second tree to be compared. If it has already been seen by self, its
values will have been cached. If the tree has changed since it has been
seen by self, it will need to be recached, either explicitly before
the calculation by calling 'update_cache' or by specifying
'is_tree2_cache_updated=True'
is_tree1_cache_updated : bool
If ``tree1`` has not been seen before, then this is ignored as the
cache will be updated regardless. If ``tree1`` has been seen, then
the cached values representing it will be used unless
``is_tree1_cache_updated`` is ``False``.
is_tree2_cache_updated : bool
If ``tree1`` has not been seen before, then this is ignored as the
cache will be updated regardless. If ``tree1`` has been seen, then
the cached values representing it will be used unless
``is_tree1_cache_updated`` is ``False``.
Acknowledgements
----------------
Based in part on:
KAMPHIR
https://github.com/ArtPoon/kamphir.git
KAMPHIR is written and maintained by: Art F.Y. Poon.
With major contributions from: Rosemary McCloskey
Copyright (c) 2015, Art Poon. All rights reserved.
See https://github.com/ArtPoon/kamphir/blob/master/LICENSE.md for
more license information.
Original work adapted from Moschitti (2006) Making tree kernels
practical for natural language learning. Proceedings of the
11th Conference of the European Chapter of the Association
for Computational Linguistics.
"""
internal_nodes2 = list(tree2.postorder_internal_node_iter())
k = 0
if not is_tree1_cache_updated:
tree1_cache = self.update_cache(tree1)
else:
try:
tree1_cache = self._tree_cache[tree1]
except KeyError:
tree1_cache = self.update_cache(tree1)
if not is_tree2_cache_updated:
tree2_cache = self.update_cache(tree2)
else:
try:
tree2_cache = self._tree_cache[tree2]
except KeyError:
tree2_cache = self.update_cache(tree2)
dp_matrix = {}
for ni, tree1_node in enumerate(tree1.postorder_internal_node_iter()):
tree1_cache_node = tree1_cache[tree1_node]
for tree2_node in internal_nodes2:
tree2_cache_node = tree2_cache[tree2_node]
if tree1_cache_node.production != tree2_cache_node.production:
continue
res = self.decay_factor * math.exp( -1. / self.gauss_factor
* (tree1_cache_node.sum_of_square_edge_lengths + tree2_cache_node.sum_of_square_edge_lengths - 2*sum([(tree1_cache_node.edge_lengths[i]*tree2_cache_node.edge_lengths[i]) for i in range(len(tree1_cache_node.edge_lengths))])))
## TODO:
## - (check and) handles cases where unequal number of children
## - how to handle rotation mismatch problems? or do we assume
## trees have equal rotations
# for node_idx in range(2):
# c1 = tree1_node.clades[node_idx]
# c2 = tree2_node.clades[node_idx]
for c1, c2 in zip(tree1_node.child_node_iter(), tree2_node.child_node_iter()):
if tree1_cache[c1].production != tree2_cache[c2].production:
continue
if tree1_cache[c1].production == 0:
# branches are terminal
res *= self.sigma + self.decay_factor
else:
try:
res *= self.sigma + dp_matrix[(tree1_cache[c1].index, tree2_cache[c2].index)]
except KeyError:
res *= self.sigma
dp_matrix[(tree1_cache[tree1_node].index, tree2_cache[tree2_node].index)] = res
k += res
return k
##############################################################################
### AssemblageInducedTree
class AssemblageInducedTreeManager(object):
def __init__(self, *args, **kwargs):
self.is_exchangeable_assemblage_classifications = kwargs.pop("is_exchangeable_assemblage_classifications", True)
self._num_assemblage_classifications = kwargs.pop("num_assemblages", None)
self.induced_tree_factory = kwargs.pop("induced_tree_factory", None)
self.induced_tree_node_factory = kwargs.pop("induced_tree_node_factory", None)
self.skip_null_assemblages = kwargs.pop("skip_null_assemblages", False)
self._tree_assemblage_induced_trees_map = {}
def remove_from_cache(self, tree):
del self._tree_assemblage_induced_trees_map[tree]
def generate_induced_trees(self, tree, assemblage_leaf_sets):
if assemblage_leaf_sets is None:
if self._num_assemblage_classifications is not None:
raise ValueError("Expecting {} assemblage leaf set classifications, but none provided".format(self._num_assemblage_classifications))
else:
if self._num_assemblage_classifications is None:
self._num_assemblage_classifications = len(assemblage_leaf_sets)
elif len(assemblage_leaf_sets) != self._num_assemblage_classifications:
raise ValueError("Expecting {} assemblage leaf set classifications, but only {} specified".format(
self._num_assemblage_classifications,
len(assemblage_leaf_sets)))
induced_trees = []
for idx, assemblage_leaf_set in enumerate(assemblage_leaf_sets):
if len(assemblage_leaf_set) == 0:
if self.skip_null_assemblages:
continue
raise error.NullLeafSetException()
node_filter_fn = lambda nd: nd in assemblage_leaf_set
induced_tree = tree.extract_tree(
node_filter_fn=node_filter_fn,
is_apply_filter_to_leaf_nodes=True,
is_apply_filter_to_internal_nodes=False,
tree_factory=self.induced_tree_factory,
node_factory=self.induced_tree_node_factory)
induced_trees.append(induced_tree)
self._tree_assemblage_induced_trees_map[tree] = induced_trees
return induced_trees
##############################################################################
### AssemblageInducedTreeShapeKernel
class AssemblageInducedTreeShapeKernel(TreeShapeKernel, AssemblageInducedTreeManager):
@staticmethod
def _euclidean_distance(v1, v2, is_weight_values_by_comparison_size=True):
v1_size = len(v1)
v2_size = len(v2)
v1_idx = 0
v2_idx = 0
if v1_size > v2_size:
v1_idx = v1_size - v2_size
weight = float(v2_size)
elif v2_size > v1_size:
v2_idx = v2_size - v1_size
weight = float(v1_size)
else:
weight = float(v1_size)
if not is_weight_values_by_comparison_size:
weight = 1.0
ss = 0.0
while v1_idx < v1_size and v2_idx < v2_size:
ss += pow(v1[v1_idx]/weight - v2[v2_idx]/weight, 2)
v1_idx += 1
v2_idx += 1
return math.sqrt(ss)
def __init__(self, *args, **kwargs):
self.exchangeable_assemblage_comparison_strategy = kwargs.pop("exchangeable_assemblage_comparison_strategy", "joint minimum")
TreeShapeKernel.__init__(self, *args, **kwargs)
AssemblageInducedTreeManager.__init__(self, *args, **kwargs)
def remove_from_cache(self, tree):
for induced_tree in self._tree_assemblage_induced_trees_map[tree]:
TreeShapeKernel.remove_from_cache(self, induced_tree)
TreeShapeKernel.remove_from_cache(self, tree)
AssemblageInducedTreeManager.remove_from_cache(self, tree)
def update_assemblage_induced_tree_cache(self,
tree,
assemblage_leaf_sets):
self.update_cache(tree=tree)
induced_trees = self.generate_induced_trees(tree=tree,
assemblage_leaf_sets=assemblage_leaf_sets)
for induced_tree in induced_trees:
self.update_cache(tree=induced_tree)
def __call__(self,
tree1,
tree2,
tree1_assemblage_leaf_sets,
tree2_assemblage_leaf_sets,
is_tree1_cache_updated=False,
is_tree2_cache_updated=False,
):
main_trees_score = TreeShapeKernel.__call__(self,
tree1=tree1,
tree2=tree2,
is_tree1_cache_updated=is_tree1_cache_updated,
is_tree2_cache_updated=is_tree2_cache_updated,
)
if not is_tree1_cache_updated or tree1 not in self._tree_assemblage_induced_trees_map:
if tree1_assemblage_leaf_sets is None:
raise ValueError("Uncached tree requires specification of 'tree1_assemblage_leaf_sets'")
self.update_assemblage_induced_tree_cache(
tree=tree1,
assemblage_leaf_sets=tree1_assemblage_leaf_sets)
if not is_tree2_cache_updated or tree2 not in self._tree_assemblage_induced_trees_map:
if tree2_assemblage_leaf_sets is None:
raise ValueError("Uncached tree requires specification of 'tree2_assemblage_leaf_sets'")
self.update_assemblage_induced_tree_cache(
tree=tree2,
assemblage_leaf_sets=tree2_assemblage_leaf_sets)
## ++ main tree score
score_table = collections.OrderedDict()
score_table["primary.tree.kernel.trick.distance"] = main_trees_score
induced_trees1 = self._tree_assemblage_induced_trees_map[tree1]
induced_trees2 = self._tree_assemblage_induced_trees_map[tree2]
# assert len(induced_trees1) == len(induced_trees2) == self._num_assemblage_classifications
if not self.is_exchangeable_assemblage_classifications:
if len(induced_trees1) != len(induced_trees2):
raise TypeError("Different numbers of induced trees not supported for non-exchangeable classifications: {} vs. {}".format(len(induced_trees1), len(induced_trees2)))
for idx, (induced_tree1, induced_tree2) in enumerate(zip(induced_trees1, induced_trees2)):
s = TreeShapeKernel.__call__(self,
tree1=induced_tree1,
tree2=induced_tree2,
is_tree1_cache_updated=True,
is_tree2_cache_updated=True,
)
## ++ raw scores direct comparisons of each of the induced trees
score_table["induced.tree.{}.kernel.trick.distance".format(idx+1)] = s
else:
if self.exchangeable_assemblage_comparison_strategy == "joint minimum":
# if lengths are different, we want to fix the smaller set
if len(induced_trees1) > len(induced_trees2):
induced_trees2, induced_trees1 = induced_trees1, induced_trees2
comparison_vector = [0.0] * len(induced_trees1)
current_minimum_distance = None
current_joint_minimum_vector = None
for induced_trees_permutation in itertools.permutations(induced_trees2, len(induced_trees1)):
distances = []
for t2, t1 in zip(induced_trees_permutation, induced_trees1):
distances.append(TreeShapeKernel.__call__(self,
tree1=t1,
tree2=t2,
is_tree1_cache_updated=True,
is_tree2_cache_updated=True,))
euclidean_distance = self._euclidean_distance(distances, comparison_vector)
if current_minimum_distance is None or euclidean_distance < current_minimum_distance:
current_minimum_distance = euclidean_distance
current_joint_minimum_vector = distances
for didx, d in enumerate(distances):
score_table["induced.tree.{}.kernel.trick.distance".format(didx+1)] = d
for didx in range(didx+1, self._num_assemblage_classifications):
score_table["induced.tree.{}.kernel.trick.distance".format(didx+1)] = "NA"
else:
raise NotImplementedError()
return score_table
###############################################################################
## Legacy
[docs]
def robinson_foulds_distance(tree1, tree2, edge_weight_attr="length"):
"""
DEPRECATED: Use :func:``symmetric_difference`` for the common
unweighged Robinson-Fould's distance metric (i.e., the symmetric difference between two trees)
:func:``weighted_robinson_foulds_distance`` or for the RF distance as defined by Felsenstein, 2004.
"""
return weighted_robinson_foulds_distance(tree1, tree2, edge_weight_attr)
[docs]
def mason_gamer_kellogg_score(tree1, tree2, is_bipartitions_updated=False):
"""
Mason-Gamer and Kellogg. Testing for phylogenetic conflict among molecular
data sets in the tribe Triticeae (Gramineae). Systematic Biology (1996)
vol. 45 (4) pp. 524
"""
if tree1.taxon_namespace is not tree2.taxon_namespace:
raise error.TaxonNamespaceIdentityError(tree1, tree2)
if not is_bipartitions_updated:
tree1.encode_bipartitions()
tree2.encode_bipartitions()
else:
if tree1.bipartition_encoding is None:
tree1.encode_bipartitions()
if tree2.bipartition_encoding is None:
tree2.encode_bipartitions()
se1 = tree1.bipartition_encoding
se2 = tree2.bipartition_encoding
bipartitions = sorted(list(set(se1.keys() + se2.keys())))
###############################################################################
## Supporting
def _get_length_diffs(
tree1,
tree2,
edge_weight_attr="length",
value_type=float,
is_bipartitions_updated=False,
bipartition_length_diff_map=False):
"""
Returns a list of tuples, with the first element of each tuple representing
the length of the branch subtending a particular bipartition on ``tree1``, and
the second element the length of the same branch on ``tree2``. If a
particular bipartition is found on one tree but not in the other, a value of zero
is used for the missing bipartition.
"""
length_diffs = []
bipartition_length_diffs = {}
if tree1.taxon_namespace is not tree2.taxon_namespace:
raise error.TaxonNamespaceIdentityError(tree1, tree2)
if not is_bipartitions_updated:
tree1.encode_bipartitions()
tree2.encode_bipartitions()
else:
if tree1.bipartition_encoding is None:
tree1.encode_bipartitions()
if tree2.bipartition_encoding is None:
tree2.encode_bipartitions()
tree1_bipartition_edge_map = dict(tree2.bipartition_edge_map) # O(n*(2*bind + dict_item_cost))
tree2_bipartition_edge_map = tree1.bipartition_edge_map
for bipartition in tree2_bipartition_edge_map: # O n : 2*bind
edge = tree2_bipartition_edge_map[bipartition]
elen1 = getattr(edge, edge_weight_attr) # attr + bind
if elen1 is None:
elen1 = 0 # worst-case: bind
value1 = value_type(elen1) # ctor + bind
try:
e2 = tree1_bipartition_edge_map.pop(bipartition) # attr + dict_lookup + bind
elen2 = getattr(e2, edge_weight_attr) # attr + bind
if elen2 is None:
# allow root edge to have bipartition with no value: raise error if not root edge
if e2.tail_node is None:
elen2 = 0.0
else:
raise ValueError("Edge length attribute is 'None': Tree: %s ('%s'), Split: %s" % (id(tree2), tree2.label, bipartition.leafset_as_newick_string(tree2.taxon_namespace)))
except KeyError: # excep
elen2 = 0.0
value2 = value_type(elen2) # ctor + bind # best case
# if abs(value2-value1) > 1e-5:
# print("{}: {}, {}".format(bipartition.leafset_as_newick_string(tree1.taxon_namespace), value2, value1))
length_diffs.append((value1,value2)) # ctor + listappend
bipartition_length_diffs[bipartition] = length_diffs[-1]
for bipartition in tree1_bipartition_edge_map: # best-case not executed, worst case O(n) : 2*bind
edge = tree1_bipartition_edge_map[bipartition]
elen2 = getattr(edge, edge_weight_attr) # attr + bind
if elen2 is None:
elen2 = 0
value2 = value_type(elen2) # ctor + bind
e1 = tree2_bipartition_edge_map.get(bipartition) # attr + dict_lookup + bind
if e1 is None:
elen1 = 0.0
else:
elen1 = getattr(e1, edge_weight_attr) # attr + bind
if elen1 is None:
# allow root edge to have bipartition with no value: raise error if not root edge
if e1.tail_node is None:
elen1 = 0.0
else:
raise ValueError("Edge length attribute is 'None': Tree: %s ('%s'), Split: %s" % (id(tree1), tree1.label, bipartition))
#elen1 = 0
value1 = value_type(elen1)
length_diffs.append((value1,value2)) # ctor + listappend
bipartition_length_diffs[bipartition] = length_diffs[-1]
# the numbers below do not reflect additions to the code to protect against
# edges with length None
# loops
# best-case:
# O(n * (dict_lookup + 3*attr + 3*ctor + 7*bind + listappend))
# worst-case:
# separated: O(n * (2*dict_lookup + 4*attr + 3*ctor + 8*bind + listappend + excep) + n*(2*dict_lookup + 4*attr + 3*ctor + 8*bind + listappend))
# or:
# O(2n*(2*dict_lookup + 4*attr + 3*ctor + 8*bind + listappend + 0.5*excep))
# total
# best-case:
# O(n * (dict_lookup + 3*attr + 3*ctor + 8*bind + listappend + dict_item_cost))
# worst-case:
# O(2n*(2*dict_lookup + 4*attr + 3*ctor + 9*bind + listappend + 0.5*(dict_item_cost + excep))
if bipartition_length_diff_map:
return length_diffs, bipartition_length_diffs
else:
return length_diffs
def _bipartition_difference(
tree1,
tree2,
dist_fn,
edge_weight_attr="length",
value_type=float,
is_bipartitions_updated=False):
"""
Returns distance between two trees, each represented by a dictionary of
bipartitions (as bipartition_mask strings) to edges, using ``dist_fn`` to calculate the
distance based on ``edge_weight_attr`` of the edges. ``dist_fn`` is a function
that takes a list of pairs of values, where the values correspond to the edge
lengths of a given bipartition on tree1 and tree2 respectively.
"""
length_diffs = _get_length_diffs(
tree1,
tree2,
edge_weight_attr=edge_weight_attr,
value_type=value_type,
is_bipartitions_updated=is_bipartitions_updated)
return dist_fn(length_diffs)
