From 4d0fc9176f30524d15da8f08de4f1dd08573e11b Mon Sep 17 00:00:00 2001 From: ahl27 <30053966+ahl27@users.noreply.github.com> Date: Tue, 1 Sep 2026 09:28:02 +0100 Subject: [PATCH 1/3] Fix issue where trees with no labels in common would silently return without warning --- R/PhyloDistance.R | 25 +++++++++++++++++++++---- man/PhyloDistance-CI.Rd | 4 ++-- man/PhyloDistance.Rd | 5 +++-- src/SEutils.h | 10 +++++----- 4 files changed, 31 insertions(+), 13 deletions(-) diff --git a/R/PhyloDistance.R b/R/PhyloDistance.R index 844ab61..1e7233b 100644 --- a/R/PhyloDistance.R +++ b/R/PhyloDistance.R @@ -99,11 +99,28 @@ 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") diff --git a/man/PhyloDistance-CI.Rd b/man/PhyloDistance-CI.Rd index c4b5c22..699445b 100644 --- a/man/PhyloDistance-CI.Rd +++ b/man/PhyloDistance-CI.Rd @@ -34,8 +34,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. diff --git a/man/PhyloDistance.Rd b/man/PhyloDistance.Rd index 17f61b5..fe46374 100644 --- a/man/PhyloDistance.Rd +++ b/man/PhyloDistance.Rd @@ -46,9 +46,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{ diff --git a/src/SEutils.h b/src/SEutils.h index f5fbf2e..6872109 100644 --- a/src/SEutils.h +++ b/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; From 9df8e7c4d7f790733c72f5658f078847c0bb0ab2 Mon Sep 17 00:00:00 2001 From: ahl27 <30053966+ahl27@users.noreply.github.com> Date: Tue, 1 Sep 2026 09:49:53 +0100 Subject: [PATCH 2/3] Fix issue where comparing tree to itself returns non-zero CI distance This was due to the greedy matching scheme originally implemented, which can cause suboptimal partitions to be paired and return a nonzero distance. Instead, we'll just revert to using the original Hungarian algorithm here for globally optimal pairing (at the expense of computational efficiency). Also updates some math to correct for floating point precision issues. Original greedy matching is maintained via new `Exact` argument. --- R/EvoWeaver-PSPreds.R | 6 +-- R/PhyloDistance.R | 7 ++-- man/PhyloDistance-CI.Rd | 8 ++-- man/PhyloDistance.Rd | 13 ++++-- src/CDend.c | 87 +++++++++++++++++++++++++++++++++++++---- src/CDend.h | 3 +- src/HungarianAlgo.c | 6 +++ src/HungarianAlgo.h | 7 +--- src/R_init_synextend.c | 2 +- src/SynExtend.h | 2 +- 10 files changed, 111 insertions(+), 30 deletions(-) diff --git a/R/EvoWeaver-PSPreds.R b/R/EvoWeaver-PSPreds.R index e4f57b2..6860bde 100644 --- a/R/EvoWeaver-PSPreds.R +++ b/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 diff --git a/R/PhyloDistance.R b/R/PhyloDistance.R index 1e7233b..16edf0b 100644 --- a/R/PhyloDistance.R +++ b/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)) @@ -128,11 +129,11 @@ PhyloDistance <- function(dend1, dend2, Method=c("CI", "RF", "KF", "JRF"), RawSc 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, diff --git a/man/PhyloDistance-CI.Rd b/man/PhyloDistance-CI.Rd index 699445b..a01fdab 100644 --- a/man/PhyloDistance-CI.Rd +++ b/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. @@ -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) } diff --git a/man/PhyloDistance.Rd b/man/PhyloDistance.Rd index fe46374..98f32b2 100644 --- a/man/PhyloDistance.Rd +++ b/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: @@ -83,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") diff --git a/src/CDend.c b/src/CDend.c index 6e9b261..d3441ba 100644 --- a/src/CDend.c +++ b/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 maxSim) maxSim = cursum; + } + } + + int maxdim = nrow > ncol ? nrow : ncol; + double *costMat = malloc(sizeof(double) * maxdim * maxdim); + for (int i=0; i= 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); diff --git a/src/CDend.h b/src/CDend.h index cb24a73..0fab1eb 100644 --- a/src/CDend.h +++ b/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); diff --git a/src/HungarianAlgo.c b/src/HungarianAlgo.c index 15eb41c..f371d44 100644 --- a/src/HungarianAlgo.c +++ b/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: diff --git a/src/HungarianAlgo.h b/src/HungarianAlgo.h index 22eec28..d387e68 100644 --- a/src/HungarianAlgo.h +++ b/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); diff --git a/src/R_init_synextend.c b/src/R_init_synextend.c index a2f1d07..15c1f4e 100644 --- a/src/R_init_synextend.c +++ b/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), diff --git a/src/SynExtend.h b/src/SynExtend.h index 7b531a5..b99a8cd 100644 --- a/src/SynExtend.h +++ b/src/SynExtend.h @@ -25,7 +25,7 @@ SEXP calcScoreHamming(SEXP ov1, SEXP ov2, SEXP NN, SEXP norm); SEXP cladeCollapsePA(SEXP tnPtr, SEXP ANCESTRAL_STATES); // Tree distance -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); SEXP RFDist(SEXP tnPtr1, SEXP tnPtr2, SEXP allLabels); SEXP KFDist(SEXP tnPtr1, SEXP tnPtr2, SEXP allLabels); From 582d33262e0ca622c830a3a64750e4c5d7d30e0f Mon Sep 17 00:00:00 2001 From: ahl27 <30053966+ahl27@users.noreply.github.com> Date: Tue, 1 Sep 2026 10:02:37 +0100 Subject: [PATCH 3/3] Update NEWS.md --- NEWS.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/NEWS.md b/NEWS.md index 9e087d7..a66e74c 100644 --- a/NEWS.md +++ b/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.