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. 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 844ab61..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)) @@ -99,11 +100,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") @@ -111,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 c4b5c22..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. @@ -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) } diff --git a/man/PhyloDistance.Rd b/man/PhyloDistance.Rd index 17f61b5..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: @@ -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") 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/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; 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);