Repository: cogent3/cogent3 (python)

Issue:
add Lin-Rajan-Moret tree distance metric
I had intended to work with Yu Lin on adding this, but now as he has tragically passed away I think it's even more imperative.

The measure, which Yu Lin called the "matching distance" is superior to the Robinson-Foulds measure. The implementation should use [scipy for the assignment problem component](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linear_sum_assignment.html) since we already have scipy as a dependency.

## Proposed API

Adhere to PEP8 names ( see [cogent3 coding style](https://github.com/cogent3/c3dev/wiki/Coding-Style)).

Add a new `cogent3/phylo/tree_distance.py` module.

In that new module, add a function `lin_rajan_moret(tree1: PhyloNode, tree2: PhyloNode): -> float`. Conversion of the trees to the matrices required should be done in separate function(s) within the module that are called within `lin_rajan_moret()`. Those other functions should be semi-private (use a leading underscore, e.g. `_tree_to_matrix()`).

## Proposed docstring

In addition to following the numpydoc style, should include the citation of Yu Lin's paper under the `Notes` section.

## Proposed tests

Use pytest.

Add a `tests/test_tree_distance.py` module (the tests need to be re-organised, so this location will be changed later).

The `lin_rajan_moret()` function should raise a `ValueError` if the trees do not have matching tip sets (use `set(<tree>.get_tip_names())` for the comparison) OR if they are rooted trees (i.e. `len(<tree base>.children) == 2`). Tests should check these conditions raise the `ValueError`.

The `lin_rajan_moret()` should return the expected distance of 0 for identical shapes (but different orientations), e.g. the distance between `(a, b, (c, d))` and `(c, d, (a, b))` should be 0.

For other tests, look at the examples from the original paper.

If any sample trees need to be included for the tests, they should not be large and they should be added to `tests/data`. Check with @GavinHuttley about patterns to use in this case.

Pull request:
ENH: added TreeNode.lin_rajan_moret() method, fixes #1209
[NEW] convenience method on tree objects.

[CHANGED] replace windows line endings in phylo/tree_distance.py

Patch:
diff --git a/src/cogent3/core/tree.py b/src/cogent3/core/tree.py
index 4c570fab2..1c022fe5a 100644
--- a/src/cogent3/core/tree.py
+++ b/src/cogent3/core/tree.py
@@ -1596,6 +1596,27 @@ class TreeNode(object):
         tip = self.get_node_matching_name(outgroup_name)
         return tip.parent.unrooted_deepcopy()
 
+    def lin_rajan_moret(self, tree2) -> float:
+        """return the lin-rajan-moret distance between trees
+
+        float
+            the Lin-Rajan-Moret distance
+
+        Notes
+        -----
+        This is a distance measure that exhibits superior statistical
+        properties compared to Robinson-Foulds. It can only be applied to
+        unrooted trees.
+
+        see: Lin et al. 2012
+        A Metric for Phylogenetic Trees Based on Matching
+        IEEE/ACM Transactions on Computational Biology and Bioinformatics
+        vol. 9, no. 4, pp. 1014-1022, July-Aug. 2012
+        """
+        from cogent3.phylo.tree_distance import lin_rajan_moret
+
+        return lin_rajan_moret(self, tree2)
+
 
 class PhyloNode(TreeNode):
     def __init__(self, *args, **kwargs):
diff --git a/src/cogent3/phylo/tree_distance.py b/src/cogent3/phylo/tree_distance.py
index c9d864cb0..657b42012 100644
--- a/src/cogent3/phylo/tree_distance.py
+++ b/src/cogent3/phylo/tree_distance.py
@@ -1,94 +1,88 @@
-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. 2012
-    A Metric for Phylogenetic Trees Based on Matching
-    IEEE/ACM Transactions on Computational Biology and Bioinformatics
-    vol. 9, no. 4, pp. 1014-1022, July-Aug. 2012
-
-    """
-    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 = numpy.zeros((len(splits), len(tip_names)), dtype=bool)
-    for i, split in enumerate(splits):
-        row = rows[i]
-        # 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
-    return rows
-
-
-def _hamming(vector1: numpy.ndarray, vector2: numpy.ndarray) -> int:
-    return (vector1 != vector2).sum()
-
-
-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()
+from itertools import product
+from typing import List
+
+import numpy
+
+from scipy.optimize import linear_sum_assignment
+
+
+def lin_rajan_moret(tree1: "PhyloNode", tree2: "PhyloNode") -> float:
+    """calculate the lin-rajan-moret distance (matching distance) between trees
+
+    trees should have matching tips and must not be rooted.
+
+    Parameters
+    ----------
+    tree1, tree2: PhyloNode
+        trees to calculate distance between
+
+    Returns
+    -------
+    float
+        the Lin-Rajan-Moret distance
+
+    Notes
+    -----
+    see: Lin et al. 2012 A Metric for Phylogenetic Trees Based on Matching
+    IEEE/ACM Transactions on Computational Biology and Bioinformatics
+    vol. 9, no. 4, pp. 1014-1022, July-Aug. 2012
+    """
+    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 = numpy.zeros((len(splits), len(tip_names)), dtype=bool)
+    for i, split in enumerate(splits):
+        row = rows[i]
+        # 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
+    return rows
+
+
+def _hamming(vector1: numpy.ndarray, vector2: numpy.ndarray) -> int:
+    return (vector1 != vector2).sum()
+
+
+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/tests/test_core/test_tree.py b/tests/test_core/test_tree.py
index 76703702e..6d8d8647c 100644
--- a/tests/test_core/test_tree.py
+++ b/tests/test_core/test_tree.py
@@ -2201,3 +2201,12 @@ def test_get_distances_endpoints(num_tips):
     dists = tree.get_distances(endpoints=names[:num_tips] or None)
     num = 4 if num_tips == 0 else num_tips
     assert len(dists) == (num**2 - num)
+
+
+def test_lrm_method():
+    # this test just exercises the method, the tests on the underlying
+    # function are in test_tree_distance.py
+    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 = a.lin_rajan_moret(b)
+    assert distance == 8