tree.pyx: Add t.distance_matrix() and reimplement/deprecate cophenetic_matrix().

We use now the much faster implementation of distance_matrix() in
operations.pyx. And since cophenetic_matrix() is less clear, has less
options, and the array of leaf names that it returns is unnecessary,
mark it as deprecated. We may want to remove it completely in the future.

Author: jordibc

ete4/core/tree.pyx

@@ -1747,108 +1747,32 @@ cdef class Tree:
         """
         ops.resolve_polytomy(self, descendants)
 
-    def cophenetic_matrix(self):
-        """Return a cophenetic distance matrix of the tree.
-
-        The `cophenetic matrix
-        <https://en.wikipedia.org/wiki/Cophenetic>`_ is a matrix
-        representation of the distance between each node.
-
-        If we have a tree like::
-
-                 ╭╴A
-             ╭╴y╶┤
-             │   ╰╴B
-          ╴z╶┤
-             │   ╭╴C
-             ╰╴x╶┤
-                 │   ╭╴D
-                 ╰╴w╶┤
-                     ╰╴E
-
-        where w, x, y, z are internal nodes, then::
-
-          d(A,B) = d(y,A) + d(y,B)
-
-        and::
-
-          d(A,E) = d(z,A) + d(z,E)
-                 = (d(z,y) + d(y,A)) + (d(z,x) + d(x,w) + d(w,E))
-
-        To compute it, we use an idea from
-        https://gist.github.com/jhcepas/279f9009f46bf675e3a890c19191158b
-
-        First, for each node we find its path to the root. For example::
-
-          A -> A, y, z
-          E -> E, w, x, z
-
-        and make these orderless sets. Then we XOR the two sets to
-        only find the elements that are in one or other sets but not
-        both. In this case A, E, y, x, w.
-
-        The distance between the two nodes is the sum of the distances
-        from each of those nodes to the parent
-
-        One more optimization: since the distances are symmetric, and
-        distance to itself is zero we user itertools.combinations
-        rather than itertools.permutations. This cuts our computes
-        from theta(n^2) 1/2n^2 - n (= O(n^2), which is still not
-        great, but in reality speeds things up for large trees).
+    def distance_matrix(self, selector=None, topological=False, squared=False):
+        """Return a matrix of paired distances between all the selected leaves.
 
-        For this tree, we will return the two dimensional array::
-
-                    A                B                C                D                E
-          A         0         d(A,y) + d(B,y)  d(A,z) + d(C,z)  d(A,z) + d(D,z)  d(A,z) + d(E,z)
-          B  d(B,y) + d(A,y)         0         d(B,z) + d(C,z)  d(B,z) + d(D,z)  d(B,z) + d(E,z)
-          C  d(C,z) + d(A,z)  d(C,z) + d(B,z)         0         d(C,x) + d(D,x)  d(C,x) + d(E,x)
-          D  d(D,z) + d(A,z)  d(D,z) + d(B,z)  d(D,x) + d(C,x)         0         d(D,w) + d(E,w)
-          E  d(E,z) + d(A,z)  d(E,z) + d(B,z)  d(E,x) + d(C,x)  d(E,w) + d(D,w)         0
-
-        We will also return the one dimensional array with the leaves
-        in the order in which they appear in the matrix (i.e. the
-        column and/or row headers).
+        :param selector: Function that returns True for the selected leaves.
+            If None, all leaves will be selected.
+        :param topological: If True, the distance between nodes will be the
+            number of nodes between them (instead of the sum of branch lenghts).
+        :param squared: If True, the output matrix will be squared and symmetrical.
+            Otherwise, only the upper triangle is returned (to save memory).
         """
-        # TODO: Consider naming this distance_matrix(), and also writing it
-        #       more clearly (and much faster, and with more options) as:
-        #   return ops.distance_matrix(self, ...)
+        return ops.distance_matrix(self, selector, topological, squared)
 
-        leaves = list(self.leaves())
-        paths = {x: set() for x in leaves}
-
-        # get the paths going up the tree
-        # we get all the nodes up to the last one and store them in a set
-
-        for n in leaves:
-            if n.is_root:
-                continue
-            movingnode = n
-            while not movingnode.is_root:
-                paths[n].add(movingnode)
-                movingnode = movingnode.up
-
-        # now we want to get all pairs of nodes using itertools combinations. We need AB AC etc but don't need BA CA
-
-        leaf_distances = {x.name: {} for x in leaves}
-
-        for (leaf1, leaf2) in itertools.combinations(leaves, 2):
-            # figure out the unique nodes in the path
-            uniquenodes = paths[leaf1] ^ paths[leaf2]
-            distance = sum(x.dist for x in uniquenodes)
-            leaf_distances[leaf1.name][leaf2.name] = leaf_distances[leaf2.name][leaf1.name] = distance
+    def cophenetic_matrix(self):
+        """Return a cophenetic matrix, and an array with the leaf names.
 
-        allleaves = sorted(leaf_distances.keys()) # the leaves in order that we will return
+        The `cophenetic matrix <https://en.wikipedia.org/wiki/Cophenetic>`_ is
+        a matrix where each element is the distance between two leaves.
 
-        output = [] # the two dimensional array that we will return
+        The matrix elements correspond to the order of the leaf names array.
+        """
+        from warnings import warn
+        warn('Use distance_matrix() instead', DeprecationWarning)
 
-        for i, n in enumerate(allleaves):
-            output.append([])
-            for m in allleaves:
-                if m == n:
-                    output[i].append(0) # distance to ourself = 0
-                else:
-                    output[i].append(leaf_distances[n][m])
-        return output, allleaves
+        matrix = ops.distance_matrix(self, squared=True)
+        names = [leaf.name for leaf in self.leaves()]
+        return matrix, names
 
     # TODO: The next two functions should really be parsers (which may
     # be also used when calling __init__().