Repository: cogent3/cogent3
Pull request #1313 (2023-03-13 03:35:38)

NEW: Add Lin Rajan Moret Distance
#1209 

The original source code can be found on [this webpage](https://users.cecs.anu.edu.au/~u1024708/). See: Matching Distance


Patch:
diff --git a/src/cogent3/phylo/tree_distance.py b/src/cogent3/phylo/tree_distance.py
new file mode 100644
index 000000000..a9787673d
--- /dev/null
+++ b/src/cogent3/phylo/tree_distance.py
@@ -0,0 +1,95 @@
+from itertools import product
+from typing import List
+
+import numpy
+
+from scipy.optimize import linear_sum_assignment
+
+from cogent3.core.tree import PhyloNode
+
+
+def lin_rajan_moret(tree1: PhyloNode, tree2: PhyloNode) -> float:
+    """calculate the lin-rajan-moret distance (matchign distance) between two trees
+
+    trees should have matching tips and must not be rooted.
+
+    Parameters
+    ----------
+    tree1: PhyloNode
+    tree2: PhyloNode
+        trees to calculate distance between
+
+    Returns
+    -------
+    float 
+        the lin-rajan-moret distance
+
+    Notes
+    -----
+    see:
+    Lin et al. 2011
+    A Metric for Phylogenetic Trees Based on Matching
+
+    """
+    names = tree1.get_tip_names()
+
+    if set(names) != set(tree2.get_tip_names()):
+        raise ValueError("tree tip names must match")
+    if len(tree1.children) == 2 or len(tree2.children) == 2:
+        raise ValueError("trees must be unrooted")
+
+    names.sort()
+
+    vector1 = _convert_tree_to_vectors(tree1, names)
+    vector2 = _convert_tree_to_vectors(tree2, names)
+
+    matching_distance = _matched_distance(vector1, vector2)
+
+    return float(matching_distance)
+
+
+def _convert_tree_to_vectors(tree: PhyloNode, tip_names: List) -> numpy.ndarray:
+    ref_tip = tip_names[0]
+    name_set = set(tip_names)
+    name_index = {n: i for i, n in enumerate(tip_names)}
+    # we get the tree as a set of splits.
+    splits = tree.subsets()
+    rows = []
+    for split in splits:
+        row = numpy.zeros(len(tip_names), int)
+        # Cogent only returns one side of a
+        # split, so we build the other side
+        if ref_tip in split:
+            names = list(split)
+        else:
+            names = list(name_set - split)
+        indices = [name_index[n] for n in names]
+        row[indices] = True
+        rows.append(row)
+    rows = numpy.array(rows)
+    return rows
+
+
+def _hamming(vector1: numpy.ndarray, vector2: numpy.ndarray) -> int:
+    return numpy.count_nonzero(vector1 != vector2)
+
+
+def _weight(vector1: numpy.ndarray, vector2: numpy.ndarray) -> int:
+    return min(_hamming(vector1, vector2), _hamming(vector1, 1 - vector2))
+
+
+def _bipartite_graph(vector1: numpy.ndarray, vector2: numpy.ndarray) -> numpy.ndarray:
+    if not len(vector1) == len(vector2):
+        raise ValueError("number of edges must be equal")
+
+    B = numpy.empty([len(vector1)] * 2, int)
+    for i, j in product(*[range(len(vector1))] * 2):
+        B[i, j] = _weight(vector1[i], vector2[j])
+    return B
+
+
+def _matched_distance(vector1: numpy.ndarray, vector2: numpy.ndarray) -> int:
+    B = _bipartite_graph(vector1, vector2)
+    matching = linear_sum_assignment(B)
+    return B[matching].sum()
+


Test patch:
diff --git a/.github/workflows/testing_develop.yml b/.github/workflows/testing_develop.yml
index 1b5b9cac6..d3485d35b 100644
--- a/.github/workflows/testing_develop.yml
+++ b/.github/workflows/testing_develop.yml
@@ -2,9 +2,9 @@ name: CI
 
 on:
   push:
-    branches: [ "develop", "data_store_refactor" ]
+    branches: [ "develop", "data_store_refactor", "lin_rajan_moret_feature" ]
   pull_request:
-    branches: [ "develop", "data_store_refactor" ]
+    branches: [ "develop", "data_store_refactor", "lin_rajan_moret_feature" ]
 
 jobs:
   tests:
diff --git a/tests/test_tree_distance.py b/tests/test_tree_distance.py
new file mode 100644
index 000000000..8f8881bb8
--- /dev/null
+++ b/tests/test_tree_distance.py
@@ -0,0 +1,62 @@
+import pytest
+
+from cogent3 import make_tree
+from cogent3.phylo.tree_distance import lin_rajan_moret
+
+
+def test_different_trees():
+    a = make_tree(treestring="(1,(((2,3),4),(5,((6,(7,(8,9))),(10,11)))),12);")
+    b = make_tree(treestring="(1,((((2,3),4),5),((6,7),((8,9),(10,11)))),12);")
+
+    distance = lin_rajan_moret(a, b)
+
+    assert distance == 8
+
+    a = make_tree(treestring="((3:0.00,75:0.00):0.00,(((((((((((((((((((((((1:0.00,57:0.00):0.00,29:0.00):0.00,62:0.00):0.00,43:0.00):0.00,31:0.00):0.00,(59:0.00,89:0.00):0.00):0.00,(67:0.00,100:0.00):0.00):0.00,(((((17:0.00,61:0.00):0.00,44:0.00):0.00,46:0.00):0.00,72:0.00):0.00,82:0.00):0.00):0.00,(((((((((((((4:0.00,15:0.00):0.00,95:0.00):0.00,77:0.00):0.00,54:0.00):0.00,53:0.00):0.00,((12:0.00,64:0.00):0.00,((35:0.00,41:0.00):0.00,(86:0.00,92:0.00):0.00):0.00):0.00):0.00,((((((6:0.00,56:0.00):0.00,32:0.00):0.00,((23:0.00,52:0.00):0.00,81:0.00):0.00):0.00,42:0.00):0.00,10:0.00):0.00,(((((((19:0.00,25:0.00):0.00,39:0.00):0.00,71:0.00):0.00,83:0.00):0.00,(36:0.00,66:0.00):0.00):0.00,94:0.00):0.00,(20:0.00,(68:0.00,79:0.00):0.00):0.00):0.00):0.00):0.00,(7:0.00,38:0.00):0.00):0.00,(13:0.00,(63:0.00,99:0.00):0.00):0.00):0.00,(85:0.00,87:0.00):0.00):0.00,58:0.00):0.00,28:0.00):0.00,8:0.00):0.00):0.00,98:0.00):0.00,47:0.00):0.00,((((30:0.00,34:0.00):0.00,55:0.00):0.00,37:0.00):0.00,(70:0.00,76:0.00):0.00):0.00):0.00,(33:0.00,96:0.00):0.00):0.00,80:0.00):0.00,((((9:0.00,16:0.00):0.00,91:0.00):0.00,(74:0.00,84:0.00):0.00):0.00,27:0.00):0.00):0.00,78:0.00):0.00,(11:0.00,(49:0.00,51:0.00):0.00):0.00):0.00,40:0.00):0.00,(((45:0.00,90:0.00):0.00,97:0.00):0.00,93:0.00):0.00):0.00,48:0.00):0.00,((((5:0.00,21:0.00):0.00,(18:0.00,(26:0.00,(50:0.00,69:0.00):0.00):0.00):0.00):0.00,(14:0.00,22:0.00):0.00):0.00,88:0.00):0.00):0.00,73:0.00):0.00,(24:0.00,65:0.00):0.00):0.00,(2:0.00,60:0.00):0.00);")
+    b = make_tree(treestring="((((((3:0.00,65:0.00):0.00,(28:0.00,57:0.00):0.00):0.00,9:0.00):0.00,38:0.00):0.00,(61:0.00,98:0.00):0.00):0.00,(((((((((((((1:0.00,(((((5:0.00,11:0.00):0.00,84:0.00):0.00,88:0.00):0.00,95:0.00):0.00,56:0.00):0.00):0.00,7:0.00):0.00,68:0.00):0.00,(((20:0.00,55:0.00):0.00,47:0.00):0.00,58:0.00):0.00):0.00,62:0.00):0.00,(78:0.00,93:0.00):0.00):0.00,((4:0.00,64:0.00):0.00,10:0.00):0.00):0.00,(37:0.00,86:0.00):0.00):0.00,97:0.00):0.00,(32:0.00,60:0.00):0.00):0.00,(8:0.00,19:0.00):0.00):0.00,(44:0.00,52:0.00):0.00):0.00,49:0.00):0.00,((((((((((((((((((((((((((((2:0.00,39:0.00):0.00,(31:0.00,69:0.00):0.00):0.00,41:0.00):0.00,29:0.00):0.00,(6:0.00,((83:0.00,96:0.00):0.00,99:0.00):0.00):0.00):0.00,16:0.00):0.00,13:0.00):0.00,27:0.00):0.00,75:0.00):0.00,74:0.00):0.00,(((22:0.00,89:0.00):0.00,36:0.00):0.00,26:0.00):0.00):0.00,((35:0.00,43:0.00):0.00,(46:0.00,67:0.00):0.00):0.00):0.00,((23:0.00,(53:0.00,81:0.00):0.00):0.00,(((34:0.00,71:0.00):0.00,82:0.00):0.00,100:0.00):0.00):0.00):0.00,77:0.00):0.00,(((((14:0.00,72:0.00):0.00,33:0.00):0.00,54:0.00):0.00,51:0.00):0.00,15:0.00):0.00):0.00,(((17:0.00,48:0.00):0.00,25:0.00):0.00,73:0.00):0.00):0.00,80:0.00):0.00,92:0.00):0.00,59:0.00):0.00,40:0.00):0.00,85:0.00):0.00,63:0.00):0.00,((18:0.00,79:0.00):0.00,21:0.00):0.00):0.00,66:0.00):0.00,76:0.00):0.00,(12:0.00,94:0.00):0.00):0.00,(((24:0.00,87:0.00):0.00,(45:0.00,91:0.00):0.00):0.00,((42:0.00,(50:0.00,70:0.00):0.00):0.00,90:0.00):0.00):0.00):0.00,30:0.00):0.00);")
+
+    distance = lin_rajan_moret(a, b)
+
+    assert distance == 1402
+
+    
+
+def test_same_tree():
+    a = make_tree(treestring="(a, b, (c, d));")
+    b = make_tree(treestring="(c, d, (a, b));")
+
+    distance = lin_rajan_moret(a, b)
+
+    assert distance == 0
+
+def test_fails_different_tips():
+    a = make_tree(treestring="(a, b, (c, d));")
+    b = make_tree(treestring="(e, d, (a, b));")
+
+    with pytest.raises(ValueError):
+        lin_rajan_moret(a,b)
+    
+    a = make_tree(treestring="(a, b, (c, d));")
+    b = make_tree(treestring="(a, d, (a, b));")
+
+    with pytest.raises(ValueError):
+        lin_rajan_moret(a,b)
+
+def test_fails_unrooted():
+    a = make_tree(treestring="((a, b), (c, d));")
+    b = make_tree(treestring="(e, d, (a, b));")
+
+    with pytest.raises(ValueError):
+        lin_rajan_moret(a,b)
+    
+    a = make_tree(treestring="(a, b, (c, d));")
+    b = make_tree(treestring="((e, d), (a, b));")
+
+    with pytest.raises(ValueError):
+        lin_rajan_moret(a,b)
+
+    a = make_tree(treestring="((a, b), (c, d));")
+    b = make_tree(treestring="((e, d), (a, b));")
+
+    with pytest.raises(ValueError):
+        lin_rajan_moret(a,b)
\ No newline at end of file