added support for cophenetic distance

src/PhyloNetwork/classes.py

@@ -270,7 +270,7 @@ def CSA(self,tax1,tax2):
     @memoize_method
     def LCSA(self,tax1,tax2):
         csa=self.CSA(tax1,tax2)
-        print self,tax1,tax2,csa
+        #print self,tax1,tax2,csa
         csa.sort(lambda x,y:cmp(self.height(x),self.height(y)))
         return csa[0]        
 
@@ -293,6 +293,18 @@ def nodal_area(self):
         #mat=mat+mat.transpose()
         return sum(abs(mat.flatten()))
 
+    @memoize_method
+    def cophenetic_matrix(self):
+        n=len(self.taxa())
+        matrix=numpy.zeros((n,n),int)
+        for i in range(n):
+            ti=self.taxa()[i]
+            for j in range(i,n):
+                tj=self.taxa()[j]
+                lcsa=self.LCSA(ti,tj)
+                matrix[i,j]=self.depth(lcsa)
+        return matrix
+
     def common_taxa(self,net2):
         common=[]
         taxa1=self.taxa()

src/PhyloNetwork/distances.py

@@ -63,6 +63,21 @@ def nodal_distance_unsplitted(net1,net2,p=1,take_root=False):
         else:
             return (sum(abs(mat.flatten())**p)/2)
 
+def cophenetic_distance(net1,net2,p=1,take_root=False):
+    mat1=net1.cophenetic_matrix()
+    mat2=net2.cophenetic_matrix()
+    try:
+        mat=mat1-mat2
+    except:
+        raise Exception("Networks over different set of taxa")
+    if p==1:
+        return sum(abs(mat.flatten()))/2
+    else:
+        if take_root:
+            return (sum(abs(mat.flatten())**p))**(1.0/p)
+        else:
+            return (sum(abs(mat.flatten())**p))
+
 def transposition_distance(net1,net2):
 #        pi1=permutations.Permutation(net1.matching_permutation())
 #        pi2=permutations.Permutation(net2.matching_permutation())