Strain calculation & display (#109)

Fixes incorrect calculation of strain colours

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 = "[email protected]", 
                     role = c("aut", "cre", "cph", "prg"), 

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()`

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)
+}

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)
-}

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()

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)
+})

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)
-})