diff --git a/DESCRIPTION b/DESCRIPTION index b6925f4c..5cc0be32 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: TreeDist Type: Package Title: Calculate and Map Distances Between Phylogenetic Trees -Version: 2.6.3.9001 +Version: 2.6.3.9002 Authors@R: c(person("Martin R.", "Smith", email = "martin.smith@durham.ac.uk", role = c("aut", "cre", "cph", "prg"), diff --git a/NEWS.md b/NEWS.md index a512637d..64c1279c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,14 @@ +# TreeDist 2.6.3.9002 (development) + +- Fix calculation error in `StrainCol()` +- App: Display strain in 3D tree space viewer + + # TreeDist 2.6.3.9001 (development) - Support for distances between larger trees + # TreeDist 2.6.3.9000 (development) - Support unrooted trees in `VisualizeMatching()` diff --git a/R/MSTSegments.R b/R/MSTSegments.R new file mode 100644 index 00000000..22f584ee --- /dev/null +++ b/R/MSTSegments.R @@ -0,0 +1,79 @@ +#' Add minimum spanning tree to plot, colouring by stress +#' +#' To identify strain in a multidimensional scaling of distances, it can be +#' useful to plot a minimum spanning tree +#' \insertCite{Gower1966,SmithSpace}{TreeDist}. Colouring each edge of the +#' tree according to its strain can identify areas where the mapping is +#' stretched or compressed. +#' +#' @param mapping Two-column matrix giving _x_ and _y_ coordinates of plotted +#' points. +#' @param mstEnds Two-column matrix identifying rows of `mapping` at end of +#' each edge of the MST, as output by [`TreeTools::MSTEdges()`]. +#' @param distances Matrix or `dist` object giving original distances between +#' each pair of points. +#' @param palette Vector of colours with which to colour edges. +#' @param \dots Additional arguments to [`segments()`]. +#' +#' @examples +#' set.seed(0) +#' library("TreeTools", quietly = TRUE) +#' distances <- ClusteringInfoDist(as.phylo(5:16, 8)) +#' mapping <- cmdscale(distances, k = 2) +#' mstEnds <- MSTEdges(distances) +#' +#' # Set up blank plot +#' plot(mapping, asp = 1, frame.plot = FALSE, ann = FALSE, axes = FALSE, +#' type = "n") +#' # Add MST +#' MSTSegments(mapping, mstEnds, +#' col = StrainCol(distances, mapping, mstEnds)) +#' # Add points at end so they overprint the MST +#' points(mapping) +#' PlotTools::SpectrumLegend( +#' "bottomleft", +#' legend = c("Extended", "Median", "Contracted"), +#' bty = "n", # No box +#' y.intersp = 2, # Expand in Y direction +#' palette = hcl.colors(256L, "RdYlBu", rev = TRUE) +#' ) +#' @template MRS +#' @references \insertAllCited{} +#' @family tree space functions +#' @importFrom graphics segments +#' @export +MSTSegments <- function(mapping, mstEnds, ...) { + segments(mapping[mstEnds[, 1], 1], mapping[mstEnds[, 1], 2], + mapping[mstEnds[, 2], 1], mapping[mstEnds[, 2], 2], ...) +} + +#' @rdname MSTSegments +#' @return `StrainCol()` returns a vector in which each entry is selected from +#' `palette`, with an attribute `logStrain` denoting the logarithm of the +#' mapped over original distance, shifted such that the median value is zero. +#' Palette colours are assigned centred on the median value, with entries +#' early in `palette` assigned to edges in which the ratio of mapped +#' distance to original distance is small. +#' @importFrom grDevices hcl.colors +#' @importFrom TreeTools MSTEdges +#' @export +StrainCol <- function(distances, mapping, mstEnds = MSTEdges(distances), + palette = rev(hcl.colors(256L, "RdYlBu"))) { + distMat <- as.matrix(distances) + logStrain <- apply(mstEnds, 1, function(ends) { + orig <- distMat[ends[1], ends[2]] + mapped <- sum((mapping[ends[1], ] - mapping[ends[2], ]) ^ 2) + ( + log(mapped) / 2 # sqrt + ) - log(orig) # High when mapping extends original distances + }) + strain <- logStrain - median(logStrain[is.finite(logStrain)]) + # Infinite values arise when orig == 0 + maxVal <- max(abs(strain[is.finite(strain)])) + sqrt(.Machine$double.eps) + nCols <- length(palette) + bins <- cut(strain, seq(-maxVal, maxVal, length.out = nCols)) + + # Return: + structure(palette[bins], + logStrain = strain) +} diff --git a/R/plot.R b/R/plot.R index dc8cfe8e..5c99441d 100644 --- a/R/plot.R +++ b/R/plot.R @@ -265,85 +265,3 @@ VisualizeMatching <- function(Func, tree1, tree2, setPar = TRUE, # Return: invisible() } - -#' Add minimum spanning tree to plot, colouring by stress -#' -#' To identify strain in a multidimensional scaling of distances, it can be -#' useful to plot a minimum spanning tree -#' \insertCite{Gower1966,SmithSpace}{TreeDist}. Colouring each edge of the -#' tree according to its strain can identify areas where the mapping is -#' stretched or compressed. -#' -#' @param mapping Two-column matrix giving _x_ and _y_ coordinates of plotted -#' points. -#' @param mstEnds Two-column matrix identifying rows of `mapping` at end of -#' each edge of the MST, as output by [`TreeTools::MSTEdges()`]. -#' @param distances Matrix or `dist` object giving original distances between -#' each pair of points. -#' @param palette Vector of colours with which to colour edges. -#' @param \dots Additional arguments to [`segments()`]. -#' -#' @examples -#' set.seed(0) -#' library("TreeTools", quietly = TRUE) -#' distances <- ClusteringInfoDist(as.phylo(5:16, 8)) -#' mapping <- cmdscale(distances, k = 2) -#' mstEnds <- MSTEdges(distances) -#' -#' # Set up blank plot -#' plot(mapping, asp = 1, frame.plot = FALSE, ann = FALSE, axes = FALSE, -#' type = "n") -#' # Add MST -#' MSTSegments(mapping, mstEnds, -#' col = StrainCol(distances, mapping, mstEnds)) -#' # Add points at end so they overprint the MST -#' points(mapping) -#' PlotTools::SpectrumLegend( -#' "bottomleft", -#' legend = c("Extended", "Median", "Contracted"), -#' bty = "n", # No box -#' y.intersp = 2, # Expand in Y direction -#' palette = hcl.colors(256L, "RdYlBu", rev = TRUE) -#' ) -#' @template MRS -#' @references \insertAllCited{} -#' @family tree space functions -#' @importFrom graphics segments -#' @export -MSTSegments <- function(mapping, mstEnds, ...) { - segments(mapping[mstEnds[, 1], 1], mapping[mstEnds[, 1], 2], - mapping[mstEnds[, 2], 1], mapping[mstEnds[, 2], 2], ...) -} - -#' @rdname MSTSegments -#' @return `StrainCol()` returns a vector in which each entry is selected from -#' `palette`, with an attribute `logStrain` denoting the logarithm of the -#' mapped over original distance, shifted such that the median value is zero. -#' Palette colours are assigned centred on the median value, with entries -#' early in `palette` assigned to edges in which the ratio of mapped -#' distance to original distance is small. -#' @importFrom grDevices hcl.colors -#' @importFrom TreeTools MSTEdges -#' @export -StrainCol <- function(distances, mapping, mstEnds = MSTEdges(distances), - palette = rev(hcl.colors(256L, "RdYlBu"))) { - distMat <- as.matrix(distances) - x <- mapping[, 1] - y <- mapping[, 2] - logStrain <- apply(mstEnds, 1, function(ends) { - orig <- distMat[ends[1], ends[2]] - mapped <- sum((x[ends[1]] - y[ends[2]]) ^ 2) - ( - log(mapped) / 2 # sqrt - ) - log(orig) # High when mapping extends original distances - }) - strain <- logStrain - median(logStrain[is.finite(logStrain)]) - # Infinite values arise when orig == 0 - maxVal <- max(abs(strain[is.finite(strain)])) + sqrt(.Machine$double.eps) - nCols <- length(palette) - bins <- cut(strain, seq(-maxVal, maxVal, length.out = nCols)) - - # Return: - structure(palette[bins], - logStrain = strain) -} diff --git a/inst/treespace/app.R b/inst/treespace/app.R index f4d727f7..b0e9c5d6 100644 --- a/inst/treespace/app.R +++ b/inst/treespace/app.R @@ -1378,9 +1378,16 @@ server <- function(input, output, session) { rgl::text3d(proj[, 1], proj[, 2], proj[, 3], thinnedTrees()) } if (mstSize() > 0) { - apply(mstEnds(), 1, function(segment) - rgl::lines3d(proj[segment, 1], proj[segment, 2], proj[segment, 3], - col = "#bbbbbb", lty = 1)) + rgl::segments3d( + proj[t(mstEnds()), ], + col = if ("mstStrain" %in% input$display) { + rep(StrainCol(distances(), proj[, 1:3]), + each = 2) # each end of each segment + } else { + "#bbbbbb" + }, + lty = 1 + ) } }) rgl::rglwidget() diff --git a/man/MSTSegments.Rd b/man/MSTSegments.Rd index 3df43f76..7a163192 100644 --- a/man/MSTSegments.Rd +++ b/man/MSTSegments.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot.R +% Please edit documentation in R/MSTSegments.R \name{MSTSegments} \alias{MSTSegments} \alias{StrainCol} diff --git a/tests/testthat/_snaps/plot/mst-example-plot.svg b/tests/testthat/_snaps/MSTSegments/mst-example-plot.svg similarity index 99% rename from tests/testthat/_snaps/plot/mst-example-plot.svg rename to tests/testthat/_snaps/MSTSegments/mst-example-plot.svg index e9a970a2..0b407aeb 100644 --- a/tests/testthat/_snaps/plot/mst-example-plot.svg +++ b/tests/testthat/_snaps/MSTSegments/mst-example-plot.svg @@ -25,17 +25,17 @@ - - - - - + + + + + - - - - - + + + + + diff --git a/tests/testthat/test-MSTSegments.R b/tests/testthat/test-MSTSegments.R new file mode 100644 index 00000000..59713c35 --- /dev/null +++ b/tests/testthat/test-MSTSegments.R @@ -0,0 +1,47 @@ +test_that("MST example plots as expected", { + skip_if_not_installed("graphics", "4.3") + skip_if_not_installed("vdiffr", "1.0") + vdiffr::expect_doppelganger("MST example plot", function() { + distances <- structure( + c(3, 2.3, 2.3, 2.3, 3, 1.7, 2.3, 2.3, 2.3, 4.5, 4.5, 1.7, 1.7, 2.3, 3, + 2.3, 1.7, 2.3, 2.3, 3.8, 4.4, 1.6, 3.1, 2.3, 2.6, 3, 1.5, 2.6, 4.7, + 4.6, 2.6, 2.3, 1.5, 3, 2.6, 1.5, 4.3, 4.7, 1.7, 2.6, 1.5, 3, 1.6, 4.6, + 4.7, 2.3, 2.3, 1.7, 1.7, 4.4, 3.8, 3.1, 3.1, 1.5, 4.4, 4.4, 3.1, 2.6, + 4.2, 4.8, 3, 4.8, 4.2, 4.7, 4.3, 1.6), + class = "dist", Size = 12L, Diag = FALSE, Upper = FALSE + ) # dput(round(ClusteringInfoDist(as.phylo(5:16, 8)), 1)) + + mapping <- structure( + c(0.75, 0.36, 0.89, 0.85, 0.89, 0.36, 0.64, 0.6, 0.6, 0.85, -3.4, -3.4, + 0, -0.25, 1.33, 0.39, -1.33, 0.25, 0, -1.52, 1.52, -0.39, -0.58, 0.58), + dim = c(12L, 2L) + ) # dput(round(cmdscale(distances, k = 2), 2)) + + mstEnds <- MSTEdges(distances) + + # Set up blank plot + plot(mapping, asp = 1, frame.plot = FALSE, ann = FALSE, axes = FALSE, + type = "n") + + # Add MST + MSTSegments(mapping, mstEnds, + col = StrainCol(distances, mapping, mstEnds)) + + # Add points at end so they overprint the MST + points(mapping) + PlotTools::SpectrumLegend( + "bottomleft", + legend = c("Extended", "Median", "Contracted"), + bty = "n", + y.intersp = 2, + lend = "square", + palette = hcl.colors(256L, "RdYlBu", rev = TRUE) + ) + }) +}) + +test_that("StrainCol() handles zeroes", { + distances <- dist(c(1, 1, 10, 100)) + mapping <- cmdscale(distances) + expect_equal(attr(StrainCol(distances, mapping), "logStrain")[1], Inf) +}) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R index 05a2d367..0c87fe24 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -168,50 +168,3 @@ test_that("VisualizeMatching() handles unrooted trees", { Plot = TreeDistPlot) }) }) - -test_that("MST example plots as expected", { - skip_if_not_installed("graphics", "4.3") - skip_if_not_installed("vdiffr", "1.0") - vdiffr::expect_doppelganger("MST example plot", function() { - distances <- structure( - c(3, 2.3, 2.3, 2.3, 3, 1.7, 2.3, 2.3, 2.3, 4.5, 4.5, 1.7, 1.7, 2.3, 3, - 2.3, 1.7, 2.3, 2.3, 3.8, 4.4, 1.6, 3.1, 2.3, 2.6, 3, 1.5, 2.6, 4.7, - 4.6, 2.6, 2.3, 1.5, 3, 2.6, 1.5, 4.3, 4.7, 1.7, 2.6, 1.5, 3, 1.6, 4.6, - 4.7, 2.3, 2.3, 1.7, 1.7, 4.4, 3.8, 3.1, 3.1, 1.5, 4.4, 4.4, 3.1, 2.6, - 4.2, 4.8, 3, 4.8, 4.2, 4.7, 4.3, 1.6), - class = "dist", Size = 12L, Diag = FALSE, Upper = FALSE - ) # dput(round(ClusteringInfoDist(as.phylo(5:16, 8)), 1)) - - mapping <- structure( - c(0.75, 0.36, 0.89, 0.85, 0.89, 0.36, 0.64, 0.6, 0.6, 0.85, -3.4, -3.4, - 0, -0.25, 1.33, 0.39, -1.33, 0.25, 0, -1.52, 1.52, -0.39, -0.58, 0.58), - dim = c(12L, 2L) - ) # dput(round(cmdscale(distances, k = 2), 2)) - - mstEnds <- MSTEdges(distances) - - # Set up blank plot - plot(mapping, asp = 1, frame.plot = FALSE, ann = FALSE, axes = FALSE, - type = "n") - - # Add MST - MSTSegments(mapping, mstEnds, - col = StrainCol(distances, mapping, mstEnds)) - - # Add points at end so they overprint the MST - points(mapping) - PlotTools::SpectrumLegend( - "bottomleft", - legend = c("Extended", "Median", "Contracted"), - bty = "n", - y.intersp = 2, - palette = hcl.colors(256L, "RdYlBu", rev = TRUE) - ) - }) -}) - -test_that("StrainCol() handles zeroes", { - distances <- dist(c(1, 1, 10, 100)) - mapping <- cmdscale(distances) - expect_equal(attr(StrainCol(distances, mapping), "logStrain")[1], Inf) -})