fix: PhyloDistance bugs

NEWS.md

@@ -1,3 +1,12 @@
+# SynExtend 1.25.2
+* `PhyloDistance` now throws a warning if two trees with no labels in common are
+compared.
+* `PhyloDistance` now uses the Hungarian algorithm for `Method="CI"` so that
+comparing a tree to itself correctly returns 0 distance.
+* `PhyloDistance` now has a `Exact` argument, defaulting to `TRUE` to use the
+Hungarian algorithm for `Method="CI"`. `Exact=FALSE` will use the previous
+greedy matching algorithm for backwards compatibility.
+
 # SynExtend 1.25.1
 * Large rewrite of gff import, `SummarizePairs` and other functions, added `EvaluatePairs` function and some other ancillary functions.
 

R/EvoWeaver-PSPreds.R

@@ -242,7 +242,7 @@ TreeDistance.EvoWeaver <- function(ew, Subset=NULL, Verbose=TRUE,
         # GRF/CI
         if (bitmask[1]){
           # GRF/CI
-          s <- .Call("GRFInfo", p1, p2, interlabs, FALSE, 0, PACKAGE="SynExtend")
+          s <- .Call("GRFInfo", p1, p2, interlabs, FALSE, 0, FALSE, PACKAGE="SynExtend")
           normval <- 0.5*(s[2] + s[3])
 
           if (is.na(normval) || normval == 0){
@@ -290,7 +290,7 @@ TreeDistance.EvoWeaver <- function(ew, Subset=NULL, Verbose=TRUE,
         }
         # JRF
         if (bitmask[3]){
-          s <- .Call("GRFInfo", p1, p2, interlabs, TRUE, JRFk, PACKAGE="SynExtend")
+          s <- .Call("GRFInfo", p1, p2, interlabs, TRUE, JRFk, FALSE, PACKAGE="SynExtend")
           normval <- (s[2] + s[3])
           if (is.na(normval) || normval == 0)
             pairscoresList$JRF[ctr+1] <- 0
@@ -300,7 +300,7 @@ TreeDistance.EvoWeaver <- function(ew, Subset=NULL, Verbose=TRUE,
         }
         # Nye
         if (bitmask[4]){
-          s <- .Call("GRFInfo", p1, p2, interlabs, TRUE, 1, PACKAGE="SynExtend")
+          s <- .Call("GRFInfo", p1, p2, interlabs, TRUE, 1, FALSE, PACKAGE="SynExtend")
           normval <- (s[2] + s[3])
           if (is.na(normval) || normval == 0)
             pairscoresList$Nye[ctr+1] <- 0

R/PhyloDistance.R

@@ -77,9 +77,10 @@ JRFDist <- function(val, RawScore=FALSE){
   return(retval)
 }
 
-PhyloDistance <- function(dend1, dend2, Method=c("CI", "RF", "KF", "JRF"), RawScore=FALSE, JRFExp=2){
+PhyloDistance <- function(dend1, dend2, Method=c("CI", "RF", "KF", "JRF"), RawScore=FALSE, JRFExp=2, Exact=TRUE){
   Method <- match.arg(Method)
   stopifnot("inputs must both be dendrograms!"=is(dend1, 'dendrogram') && is(dend2, 'dendrogram'))
+  stopifnot("Exact must be logical"=is.logical(Exact) && length(Exact) == 1)
   if (is.integer(JRFExp)) JRFExp <- as.numeric(JRFExp)
   stopifnot("ExpVal must be numeric or integer"=is.numeric(JRFExp))
 
@@ -99,23 +100,40 @@ PhyloDistance <- function(dend1, dend2, Method=c("CI", "RF", "KF", "JRF"), RawSc
   }
   incommonLabs <- intersect(labels(dend1), labels(dend2))
   if (length(incommonLabs) == 0){
-    if(RawScore)
-      val <- c(0, NA, NA)
-    else
-      val <- 1
+    warning("Input dendrograms have no shared tip labels.")
+    if(RawScore){
+      if (Method == 'CI') {
+        val <- c("Similarity"=0, "dend1.Entropy"=NA_real_, "dend2.Entropy"=NA_real_, "p.value"=NA_real_)
+      } else if (Method == 'RF') {
+        val <- c("UniqueSplits"=0, "dend1.Splits"=NA_real_, "dend2.Splits"=NA_real_)
+      } else if (Method == 'JRF') {
+        val <- c("Distance"=0, "dend1.NumSplits"=NA_real_, "dend2.NumSplits"=NA_real_)
+      } else {
+        val <- c("Similarity"=0, "dend1"=NA_real_, "dend2"=NA_real_)
+      }
+    } else {
+      val <- NA_real_
+    }
   } else {
+    if (length(incommonLabs) < length(labels(dend1)) || length(incommonLabs) < length(labels(dend2))){
+      warning("Input dendrograms have incomplete label overlap. Pruning to ",
+              length(incommonLabs), " shared tip labels.")
+    }
+    if (length(incommonLabs) < 4){
+      warning("Fewer than 4 shared tip labels between dendrograms. Unrooted trees with < 4 leaves have no non-trivial internal splits.")
+    }
     tree1ptr <- .Call("initCDend", dend1, PACKAGE="SynExtend")
     on.exit(rm(tree1ptr))
     tree2ptr <- .Call("initCDend", dend2, PACKAGE="SynExtend")
     on.exit(rm(tree2ptr))
 
     if (Method == 'CI'){
       val <- .Call("GRFInfo", tree1ptr, tree2ptr,
-                   incommonLabs, FALSE, 0, PACKAGE="SynExtend")
+                   incommonLabs, FALSE, 0, Exact, PACKAGE="SynExtend")
       return(CIDist(val, incommonLabs, RawScore))
     } else if (Method == 'JRF'){
       val <- .Call("GRFInfo", tree1ptr, tree2ptr,
-                   incommonLabs, TRUE, JRFExp, PACKAGE="SynExtend")
+                   incommonLabs, TRUE, JRFExp, Exact, PACKAGE="SynExtend")
       return(JRFDist(val, RawScore))
     } else if (Method == 'RF'){
       val <- .Call("RFDist", tree1ptr, tree2ptr,

man/PhyloDistance-CI.Rd

@@ -12,10 +12,7 @@ information of branch partitions.
 
 \details{
 This function is called as part of \code{\link{PhyloDistance}} and calculates tree distance using the clustering information
-approach first described in Smith (2020). This function iteratively pairs
-internal tree branches of a phylogeny based on their similarity, then scores
-overall similarity as the sum of these measures. The similarity score is then
-converted to a distance by normalizing by the average entropy of the two trees.
+approach first described in Smith (2020). By default (\code{Exact=TRUE}), internal tree branches are paired via the Hungarian algorithm. Setting \code{Exact=FALSE} uses greedy matching to increase speed at the cost of precision. The total similarity score is converted to a distance by normalizing by the average entropy of the two trees.
 This metric has been demonstrated to outperform numerous other metrics in capabilities;
 see the original publication cited in References for more information.
 
@@ -34,8 +31,8 @@ with different branch lengths may return a distance of 0.
 If \code{RawScore=TRUE}, returns a named length 4 vector with the first entry the similarity
 score, subsequent entries the entropy values for each tree, and the last entry the approximate p-value for the result based on simulations.
 
-If the trees have no leaves in common, the function will return \code{1} if
-\code{RawScore=FALSE}, and \code{c(0, NA, NA, NA)} if \code{TRUE}.
+If the trees have no leaves in common, the function will issue a warning and return \code{1} if
+\code{RawScore=FALSE}, and a named length 4 vector \code{c(Similarity=0, dend1.Entropy=NA, dend2.Entropy=NA, p.value=NA)} if \code{TRUE}.
 }
 \references{
 Smith, Martin R. \emph{Information theoretic generalized Robinson–Foulds metrics for comparing phylogenetic trees.} Bioinformatics, 2020. \bold{36}(20):5007-5013.
@@ -68,6 +65,9 @@ tree2 <- as.dendrogram(hclust(dm2))
 # get CI distance
 PhyloDistance(tree1, tree2, Method="CI")
 
+# get CI distance with fast greedy pairing
+PhyloDistance(tree1, tree2, Method="CI", Exact=FALSE)
+
 # get similarity score with individual entropies
 PhyloDistance(tree1, tree2, Method="CI", RawScore=TRUE)
 }

man/PhyloDistance.Rd

@@ -10,9 +10,8 @@ Calculates distance between two unrooted phylogenies using a variety of metrics.
 \usage{
 PhyloDistance(dend1, dend2,
               Method=c("CI", "RF", "KF", "JRF"),
-              RawScore=FALSE, JRFExp=2)
+              RawScore=FALSE, JRFExp=2, Exact=TRUE)
 }
-%- maybe also 'usage' for other objects documented here.
 \arguments{
   \item{dend1}{
   An object of class \code{dendrogram}, representing an unrooted bifurcating phylogenetic
@@ -31,6 +30,9 @@ PhyloDistance(dend1, dend2,
 \item{JRFExp}{
 \code{k}-value used in calculation of JRF Distance. Unused if \code{Method} is not \code{"JRF"}.
 }
+\item{Exact}{
+Logical; Determines if exact optimal bipartite matching via the Hungarian algorithm should be used for \code{Method="CI"}. Has no effect for other \code{Method} arguments. Defaults to \code{TRUE}; setting \code{FALSE} will use a greedy algorithm for higher speed at the cost of accuracy. See \link[=CIDist]{Clustering Information Distance}} for more information.
+}
 }
 \details{
 This function implements a variety of tree distances, specified by the value of \code{Method}. The following values are supported, along with links to documentation pages for each function:
@@ -46,9 +48,10 @@ Information on each of these algorithms, how scores are calculated, and referenc
 }
 \value{
 Returns a normalized distance, with 0 indicating identical trees and 1 indicating
-maximal difference. If the trees have no leaves in common, the function will return \code{1} if \code{RawScore=FALSE}, or \code{c(0,NA,NA)} if \code{RawScore=TRUE}.
+maximal difference. If the trees have no leaves in common, the function issues a warning and returns \code{NA} if \code{RawScore=FALSE}, or a named vector of \code{NA} component values if \code{RawScore=TRUE}.
+If the trees have incomplete leaf overlap, both trees are pruned to the shared leaf set with a warning.
 
-If \code{RawScore=TRUE}, returns a vector of the components used to calculate the distance. This is typically a length 3 vector, but specific details can be found on the description for each algorithm linked above.
+If \code{RawScore=TRUE}, returns a named vector of the components used to calculate the distance. Specific details can be found on the description for each algorithm linked above.
 }
 
 \author{
@@ -82,8 +85,11 @@ tree2 <- as.dendrogram(hclust(dm2))
 # Robinson-Foulds Distance
 PhyloDistance(tree1, tree2, Method="RF")
 
-# Clustering Information Distance
-PhyloDistance(tree1, tree2, Method="CI")
+# Clustering Information Distance (fast greedy pairing)
+PhyloDistance(tree1, tree2, Method="CI", Exact=FALSE)
+
+# Clustering Information Distance (exact Hungarian matching)
+PhyloDistance(tree1, tree2, Method="CI", Exact=TRUE)
 
 # Kuhner-Felsenstein Distance
 PhyloDistance(tree1, tree2, Method="KF")

src/CDend.c

@@ -174,10 +174,11 @@ SEXP calcScoreHamming(SEXP ov1, SEXP ov2, SEXP NN, SEXP norm){
 
 /**** Tree Distances ****/
 // RF Distance with information-theoretic scoring (clustering info)
-SEXP GRFInfo(SEXP tnPtr1, SEXP tnPtr2, SEXP allLabels, SEXP shouldUseJRF, SEXP JRFExp){
+SEXP GRFInfo(SEXP tnPtr1, SEXP tnPtr2, SEXP allLabels, SEXP shouldUseJRF, SEXP JRFExp, SEXP exactMatch){
   treeNode *tree1 = checkPtrExists(tnPtr1);
   treeNode *tree2 = checkPtrExists(tnPtr2);
   bool useJRF = LOGICAL(shouldUseJRF)[0];
+  bool exactVal = LOGICAL(exactMatch)[0];
   double jaccardExp = 0;
   if (useJRF)
     jaccardExp = REAL(JRFExp)[0];
@@ -223,7 +224,7 @@ SEXP GRFInfo(SEXP tnPtr1, SEXP tnPtr2, SEXP allLabels, SEXP shouldUseJRF, SEXP J
     entropy1 = (double) t1pln;
     entropy2 = (double) t2pln;
   } else {
-    RFscore = scorePMs(part1, part2, t1pln, t2pln, numLabels);
+    RFscore = scorePMs(part1, part2, t1pln, t2pln, numLabels, exactVal);
     entropy1 = calcEntropy(part1, numLabels, t1pln);
     entropy2 = calcEntropy(part2, numLabels, t2pln);
   }
@@ -815,7 +816,72 @@ int reallocPartitionMap(bool **pSets, int lh, int plen){
   return ctr;
 }
 
-double scorePMs(bool **pm1, bool **pm2, int pm1l, int pm2l, int lh){
+static double scorePMsHungarian(bool **pm1, bool **pm2, int pm1l, int pm2l, int lh){
+  int nrow = pm1l;
+  int ncol = pm2l;
+  double *simMat = malloc(sizeof(double) * nrow * ncol);
+  double maxSim = 0.0;
+
+  int counts[8];
+  bool *curS, *curL, v1, v2;
+
+  for (int i=0; i<nrow; i++){
+    curS = pm1[i];
+    for (int j=0; j<ncol; j++){
+      curL = pm2[j];
+      memset(counts, 0, sizeof(counts));
+      for (int k=0; k<lh; k++){
+        v1 = curS[k];
+        v2 = curL[k];
+        counts[0] += v1;
+        counts[1] += v2;
+        counts[2] += !v1;
+        counts[3] += !v2;
+        counts[4] += v1 && v2;
+        counts[5] += v1 && !v2;
+        counts[6] += !v1 && v2;
+        counts[7] += !v1 && !v2;
+      }
+      double cursum = 0.0;
+      cursum += PclDist(counts[0], counts[1], counts[4], lh); // A1 A2
+      cursum += PclDist(counts[0], counts[3], counts[5], lh); // A1 B2
+      cursum += PclDist(counts[2], counts[1], counts[6], lh); // B1 A2
+      cursum += PclDist(counts[2], counts[3], counts[7], lh); // B1 B2
+      simMat[i * ncol + j] = cursum;
+      if (cursum > maxSim) maxSim = cursum;
+    }
+  }
+
+  int maxdim = nrow > ncol ? nrow : ncol;
+  double *costMat = malloc(sizeof(double) * maxdim * maxdim);
+  for (int i=0; i<maxdim; i++){
+    for (int j=0; j<maxdim; j++){
+      if (i < nrow && j < ncol){
+        costMat[i * maxdim + j] = maxSim - simMat[i * ncol + j];
+      } else {
+        costMat[i * maxdim + j] = maxSim;
+      }
+    }
+  }
+
+  int *assignment = hungarian(costMat, maxdim);
+
+  double retval = 0.0;
+  for (int i=0; i<nrow; i++){
+    int assigned_col = assignment[i];
+    if (assigned_col >= 0 && assigned_col < ncol){
+      retval += simMat[i * ncol + assigned_col];
+    }
+  }
+
+  free(simMat);
+  free(costMat);
+  free(assignment);
+
+  return retval;
+}
+
+static double scorePMsGreedy(bool **pm1, bool **pm2, int pm1l, int pm2l, int lh){
   bool firstlonger = pm1l > pm2l;
   bool **longPm = firstlonger ? pm1 : pm2;
   bool **shortPm = firstlonger ? pm2 : pm1;
@@ -872,6 +938,15 @@ double scorePMs(bool **pm1, bool **pm2, int pm1l, int pm2l, int lh){
   return retval;
 }
 
+double scorePMs(bool **pm1, bool **pm2, int pm1l, int pm2l, int lh, bool exactMatch){
+  if (pm1l == 0 || pm2l == 0) return 0.0;
+  if (exactMatch){
+    return scorePMsHungarian(pm1, pm2, pm1l, pm2l, lh);
+  } else {
+    return scorePMsGreedy(pm1, pm2, pm1l, pm2l, lh);
+  }
+}
+
 double calcEntropy(bool **pm, int lh, int pml){
   double res = 0.0;
   double p1, p2;
@@ -882,10 +957,8 @@ double calcEntropy(bool **pm, int lh, int pml){
       p1 += pm[i][j];
       p2 += !pm[i][j];
     }
-    p1 /= lh;
-    p2 /= lh;
-    res += p1 == 0 ? 0 : (-1 * p1 * log2(p1));
-    res += p2 == 0 ? 0 : (-1 * p2 * log2(p2));
+    res += PclDist(p1, p1, p1, lh);
+    res += PclDist(p2, p2, p2, lh);
   }
 
   return(res);

src/CDend.h

@@ -13,6 +13,7 @@
 
 #include "SynExtend.h"
 #include "SEutils.h"
+#include "HungarianAlgo.h"
 
 typedef unsigned long ulong;
 
@@ -45,7 +46,7 @@ void findMapping(treeNode *node, int *mapping, unsigned int *hashvals, int lenHa
 /* Tree Distance */
 void internalPartitionMap(treeNode *node, bool **pSets, unsigned int *hvs, int lh, int rootv);
 int reallocPartitionMap(bool **pSets, int lh, int plen);
-double scorePMs(bool **pm1, bool **pm2, int pm1l, int pm2l, int lh);
+double scorePMs(bool **pm1, bool **pm2, int pm1l, int pm2l, int lh, bool exactMatch);
 double calcEntropy(bool **pm, int lh, int pml);
 ulong RFHashMap(treeNode *node, ulong *htable, ulong *keys, unsigned int *hvs, int lh, int rootv);
 ulong KFHashMap(treeNode *node, ulong *htable, double *dists, ulong *keys, unsigned int *hvs, int lh, int rootv);

src/HungarianAlgo.c

@@ -1,5 +1,11 @@
 #include "HungarianAlgo.h"
 
+static bool has_alloced_mem = false, has_alloced_vec = false, has_alloced_assign = false;
+static uint8_t *sa = NULL;
+static bool *cc = NULL, *cr = NULL;
+static int *av = NULL;
+static double *vec = NULL;
+
 SEXP HungarianAssignment(SEXP MATVEC, SEXP DIM){
   /*
    * Arguments:

src/HungarianAlgo.h

@@ -14,12 +14,7 @@
 
 #include "SynExtend.h"
 
-// Global variables for freeing
-bool has_alloced_mem, has_alloced_vec, has_alloced_assign;
-uint8_t *sa;
-bool *cc, *cr;
-int *av;
-double *vec;
+
 
 // Steps
 int* hungarian(double *costMatrix, int n);

src/R_init_synextend.c

@@ -39,7 +39,7 @@ static const R_CallMethodDef callMethods[] = { // method call, pointer, num args
   CALLDEF(cladeCollapsePA, 2),
   CALLDEF(calcAllTreeLengths, 1),
   CALLDEF(printTree, 1),
-  CALLDEF(GRFInfo, 5),
+  CALLDEF(GRFInfo, 6),
   CALLDEF(RFDist, 3),
   CALLDEF(KFDist, 3),
   CALLDEF(calcDValue, 2),

src/SEutils.h

@@ -21,7 +21,7 @@ void *safe_calloc(size_t nitems, size_t size);
 void *safe_realloc(void *ptr, size_t new_size);
 
 /*** Other Utility Functions ***/
-inline void *void_deref(void *v, int i, size_t size){
+static inline void *void_deref(void *v, int i, size_t size){
   return i ? (void *)((char *)v + (i*size)) : v;
 }
 
@@ -31,13 +31,13 @@ inline void *void_deref(void *v, int i, size_t size){
 /*** Random numbers ***/
 
 // random float number in range [0, 1]
-double inline frand(){ return unif_rand(); }
+static inline double frand(){ return unif_rand(); }
 
 // random integer
-int inline irand(){ return (int) floor(frand() * INT_MAX);}
+static inline int irand(){ return (int) floor(frand() * INT_MAX);}
 
 // random number from normal distribution with mean mu and standard deviation sd
-double inline rnorm(double mu, double sd){ return sd * (norm_rand()) + mu; }
+static inline double rnorm(double mu, double sd){ return sd * (norm_rand()) + mu; }
 
 
 /*** Random permutations using Fisher-Yates shuffling ***/
@@ -84,7 +84,7 @@ void seedRNGState32(struct RNGstate32 *r, uint64_t seed);
 uint32_t xorshift32b(struct RNGstate32 *r);
 
 /*** Random math functions ***/
-long inline doubleFactorial(int n){
+static inline long doubleFactorial(int n){
   long retval = 1;
   while (n > 0) retval *= n--;
   return retval;