#! /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 import dendropy from dendropy.utility import error ############################################################################### ## Public Functions 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 -------- :: 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.symmetric_difference(tree1, tree2)) """ t = false_positives_and_negatives( tree1, tree2, is_bipartitions_updated=is_bipartitions_updated) return t[0] + t[1] def unweighted_robinson_foulds_distance(tree1, tree2, is_bipartitions_updated=False): """ Alias for ``symmetric_difference()``. """ return symmetric_difference(tree1, tree2, is_bipartitions_updated) 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) 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) 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) 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 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) 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)