R/svg.R

In TGuillerme/Inapp: Reconstruction of Inapplicable Discrete Characters on Phylogenetic Trees

################################################################################
##                                                                            ##
## This file contains functions that allow the visualization of the           ##
## inapplicable algorithm to be rendered as an SVG element in HTML documents. ##
##                                                                            ##
## This makes it possible to build interactive JavaScript functions to allow  ##
## users to engage with the visualization in Javascript, and to style the     ##
## graphical output using CSS.                                                ##
##                                                                            ##
################################################################################

#' SVG canvas
#'
#' Creates an object of class `SVGCanvas`.
#'
#' `SVGCanvas` objects are used to visualize character reconstructions on
#' phylogenetic trees.
#' The vizualizations are rendered as SVG elements that can be saved
#' as stand-alone SVG files, or included in HTML documents.
#'
#' If incorporated in an HTML vignette, Javascript and CSS can be used to
#' style the output and enable users to interact with the visualizations.
#'
#' @param trees A list of trees, perhaps of class multiPhylo
#' @param uniqueTrees (optional) A list generated by unique(trees).  Can be
#'                              prepopulated to avoid time-intensive recalculation.
#' @template outgroupTipsParam
#' @param analysisNames Character vector of length `length(trees)`, to be used to
#' specify the analytical method under which each tree in `trees` was generated
#' (e.g. equal weights parsimony, maximum likelihood).  Will be displayed as a
#' legend beside each tree.  A list of methods will be displayed beside an output
#' tree as a legend; analyses associated with a tree will be given the CSS class
#' `thisAn`; other analyses will bear the class `notThisAn`.
#' Defaults to `character(0)`, which suppresses this legend.
#' @param width,height Numeric specifying the size of the plot in pixels
#' @param xMargins,yMargins Numeric vectors of length 2 specifying, in pixels,
#'  the margin to be left at each side of the plot
#' @param textY Numeric specifying how far to adjust the vertical position of text
#'
#' @return An object of class SVGCanvas
#' @exportClass SVGCanvas
#' @importFrom phangorn Ancestors Children
#' @importFrom TreeTools SplitFrequency RootTree
#'
#' @author Martin R. Smith
#'
#' @examples {
#' \dontrun{
#'  # Generate multiPhylo object containing four trees, with first and last identical:
#'    trees <- structure(list(
#'      ape::rtree(10, br=NULL),
#'      ape::rtree(10, br=NULL),
#'      ape::rtree(10, br=NULL),
#'      NULL), class='multiPhylo')
#'    trees[[4]] <- trees[[1]]
#'
#'    treesGeneratedBy <- c('rtree', 'rtree', 'rtree', 'duplication')
#'
#'    canvas <- SVGCanvas(trees=trees,
#'                        outgroupTips='t1',
#'                        analysisNames=treesGeneratedBy,
#'                        width=300, height=300)
#'
#'    PlotCharacterMapping(char = '0011----11',
#'                         stateLabels = c('Absent', 'Present'),
#'                         singleTree = trees[[1]],
#'                         legendText = 'This is printed on PNGs',
#'                         canvas = canvas,
#'                         svgFilename = 'tree_number_%s.svg')
#'
#'   }
#' }
SVGCanvas <- function (trees, outgroupTips=NULL, analysisNames=character(0),
                       width = 682, height = 682,
                       uniqueTrees=unique(trees),
                       xMargins=c(10L, 10L), yMargins=c(10L, 10L), textY = 3L) {

    rootedTrees <- lapply(uniqueTrees, RootTree, outgroupTips=outgroupTips)
    nTree <- length(uniqueTrees)
    eachTree <- seq_len(nTree)
    treeIndex <- vapply(trees,
                     function (tr1)
                       which(vapply(uniqueTrees, all.equal, logical(1), tr1,
                                    use.tip.label=TRUE, use.edge.length=FALSE)),
                     integer(1))
    nTip <- length(uniqueTrees[[1]]$tip.label)
    nNode <- nTip - 1L
    internal <- nTip + seq_len(nNode)
    terminal <- seq_len(nTip)
    tipLabel <- matrix('', nTree, nTip)
    nEdge <- dim(uniqueTrees[[1]]$edge)[1]
    eachEdge <- seq_len(nEdge)
    nodeSupport <- matrix(0, nTree, nNode - 2L) # Two root nodes irrelevant
    parents <- children <- matrix(0, nTree, nEdge)
    xStep <- nodeX <- nodeY <- vector('list', nTree)
    yStep <- round((height - yMargins[2] - yMargins[1]) / nTip, 1)
    nAnalyses <- length(analysisNames)

    for (i in eachTree) {
        tree <- rootedTrees[[i]]
        tipLabel[i, ] <- tree$tip.label
        treeEdge <- tree$edge
        parents[i, ] <- treeEdge[, 1]
        children[i, ] <- treeEdge[, 2]
        ancestors <- Ancestors(tree, type='all')
        nAncestors <- vapply(ancestors, length, 1)
        xStep[[i]] <- round(seq(xMargins[1], to=width-xMargins[2],
                                len=max(nAncestors) + 1L), 1)
        nodeX[[i]] <- xStep[[i]][nAncestors + 1L]
        tmpNodeY <- double(length(nAncestors))
        tmpNodeY[children[i, children[i, ] <= nTip]] <- seq(yMargins[1], height-yMargins[2], len=nTip)

        for (nodeI in rev(internal)) {
            tmpNodeY[nodeI] <- mean(tmpNodeY[Children(tree, nodeI)])
        }
        nodeY[[i]] <- round(tmpNodeY, 1)
        nodeSupport[i, ] <- round(SplitFrequency(tree, trees) / length(trees), 2)
    }

    structure( # For a S3 object
        list(height=height, width=width,
             rootedTrees = rootedTrees,
             analysisNames = analysisNames,
             nTree = nTree, eachTree = eachTree, treeIndex=treeIndex,
             parents = parents, children = children, tipLabel = tipLabel,
             internal = internal, terminal = terminal,
             nodeX = nodeX, nodeY = nodeY, textY = textY, yStep = yStep,
             tipLabel = tipLabel,
             nodeSupport = nodeSupport),
        class = 'SVGCanvas')
}

#' @export
is.SVGCanvas <- function (x) inherits(x, 'SVGCanvas')
#' @export
length.SVGCanvas <- function (x) x$nTree

#' SVG Tree
#'
#' Plots one of the character reconstructions embodied within an [SVGCanvas]
#' object on a specified phylogenetic tree, returning SVG text that can be
#' saved to file or included directly in an HTML document.
#'
#' @param treeNo Integer specifying which of the trees on the SVGForest to plot
#' @template canvasParam
#' @template treeNamesParam
#' @param char Character string specifying character to optimise on the tree
#' @param charIndex Character string specifying the number of the character,
#'        for inclusion under the SVG element's `data-char` attribute.
#' @param stateLabels Character vector specifying label for each applicable state
#'        of the character
#' @param analysisLabels Character vector specifying label for each analysis
#'
#' @return Character string describing an SVG object that depicts the tree,
#'   which can be written to file or included in markdown destined for HTML
#'   output.
#'
#' @importFrom TreeTools SupportColour
#' @export
#'
#' @author Martin R. Smith
SVGTree <- function (treeNo, canvas, char, stateLabels,
                     treeNames, charIndex = 0, analysisLabels=character(0)) {
    # Cache details from canvas
    tree <- canvas$rootedTree[[treeNo]]
    nodeX <- canvas$nodeX[[treeNo]]
    nodeY <- canvas$nodeY[[treeNo]]
    parents <- canvas$parents
    children <- canvas$children
    internal <- canvas$internal
    terminal <- canvas$terminal
    tipLabel <- canvas$tipLabel
    textY <- canvas$textY
    yStep <- canvas$yStep
    nodeSupport <- canvas$nodeSupport[treeNo, ]

    statesMatrix <- apply.reconstruction(tree, char[tree$tip.label])
    fitchStates <- apply.reconstruction(tree, char[tree$tip.label],
                                        method='Fitch', inapplicable=1)
    matrixData <- MatrixData(statesMatrix, fitchStates, state.labels=stateLabels)
    legendLabels <- matrixData$legend
    legendCol <- matrixData$legend_col
    edgeCol <- matrixData$edge_col
    flagCol <- matrixData$tips_colours
    tipFlag <- matrixData$tips_labels

    LineStyle <- function (x) ifelse(x == 'lightgrey', '" class="inapplicable', '')
    edges <- paste0(vapply(unique(edgeCol), function(colour) {
        theseEdges <- edgeCol == colour
        paste0('<path d="', paste0(
            'M', nodeX[parents[treeNo, theseEdges]],
            ',', nodeY[parents[treeNo, theseEdges]],
            'V', nodeY[children[treeNo, theseEdges]],
            'H', nodeX[children[treeNo, theseEdges]],
            collapse=''),
            '" stroke="', colour, LineStyle(colour), '"></path>')
    }, character(1)), collapse='')
    tips <- paste0('<text x="', (nodeX[terminal] + 0L),
                   '" y="', nodeY[terminal] + textY,
                   '" fill="', flagCol[terminal],
                   '" class="flag">[', tipFlag[terminal], ']</text>',
                   '<text x="', (nodeX[terminal] + 22L),
                   '" y="', nodeY[terminal] + textY,
                   '" fill="', ifelse(flagCol == 'lightgrey', '#ccc',
                                      ifelse(flagCol == 'grey', '#999', '#000')),
                   '" class="taxonLabel">',
                   gsub('_', ' ', tipLabel[treeNo, terminal], fixed=TRUE), '</text>',
                   collapse='')
    nodes <- paste0('<text x="', (nodeX[internal][-1] + 2L),
                    '" y="', nodeY[internal][-1] + textY,
                    '" class="node" fill="', SupportColour(nodeSupport),
                    '">', nodeSupport, '</text>', collapse='')
    nRegions <- length(statesMatrix$regions)
    if (statesMatrix$score == fitchStates$score) {
        fitchNote1 <- fitchNote2 <- ''
    } else {
        fitchNote1 <- paste0(' (<tspan class="score fitch">', fitchStates$score,
                             '</tspan> with Fitch)')
        fitchNote2 <- ' (not counted by Fitch)'
    }
    ciCaption <- paste0('<text x="', canvas$width - 4L,
                        '" y="', yStep + textY, '" text-anchor="end" class="stepsLabel">',
                        'Character adds <tspan class="score">', statesMatrix$score,
                        '</tspan> to tree score', fitchNote1, '</text><text x="',
                        canvas$width - 4L, '" y="', (yStep * 2) + textY,
                        '" text-anchor="end" class="stepsLabel">',
                        '<tspan class="score">', nRegions,
                        '</tspan> additional region', ifelse(nRegions == 1, '', 's'),
                        fitchNote2, '</text><text x="', canvas$width - 4L,
                        '" y="', (yStep * 3) + textY, '" text-anchor="end" class="stepsLabel">',
                        # ci_text[i, treeNo],
                        '</text><text x="', canvas$width - 4L,
                        '" y="', (yStep * 5) + textY, '" text-anchor="end" class="stepsLabel">',
                        '<tspan dy="16" fill="#34caaf">An optimal tree under:</tspan>',
                        paste0('<tspan x="', canvas$width - 2L,
                               '" style="font-style:italic" class="',
                               ifelse(analysisLabels %in% treeNames[canvas$treeIndex==treeNo],
                                      'this', 'notThis'),
                               'An" dy="1.2em">', analysisLabels, '</tspan>', collapse=''),
                        '</text>')
    dudSteps <- matrixData$dud_steps

    svgSource <- paste0('<svg xmlns="http://www.w3.org/2000/svg" version="1.1',
                        '" viewBox="0 0 ', canvas$width, ' ', canvas$height,
                        '" preserveAspectRatio="xMidYMid meet',
                        '" class="tree" data-char="', charIndex, '">',
                        ciCaption, tips, edges, nodes,
                        '</svg>')
    # Return:
    svgSource
}

#' Matrix data
#'
#' Generates and summarises data from a `states_matrix` that has been output
#' by `apply.reconstruction`.
#'
#' @param states_matrix The output of `apply.reconstruction`
#' @param fitch_states The output of `apply.reconstruction(method='Fitch')`
#' @param state.labels Character vector translating tokens into characer states,
#'   e.g. `c('0: Absent', '1: Present')`
#'
#' @return A list summarising details of the matrix, with named entries;
#'  `edge_col` (edge colours); `edge_col_array` (formatted as a Javascript array),
#'   `legend` (detailing what colours mean); `legend_col` (colour for legend entries);
#'  `tips_labels`; `tips_colours`; `dud_steps` (Steps reconstructed by Fitch,
#'  but not the inapplicable algorithm)
#' @export
#' @keywords internal
#' @author Martin R. Smith
MatrixData <- function (states_matrix, fitch_states = states_matrix,
                        state.labels) {
    tree <- states_matrix$tree
    regions <- states_matrix$regions
    changes <- states_matrix$changes
    steps <- fitch_states$changes
    dud_steps <- steps[!steps %in% changes]
    n_tip <- states_matrix$n_tip
    plot.convert.state <- function(char, missing = FALSE) {
        plot.convert.inappli <- function(X) {
            return(ifelse(X == -1, "-", X))
        }
        plot.convert.missing <- function(X, all_states) {
            if (length(all_states) > 1 && length(X) == length(all_states) &&
                all(sort(X) == sort(all_states))) {
                return("?")
            }
            else {
                return(X)
            }
        }
        if (missing) {
            all_states <- unique(unlist(char))
            char <- lapply(char, plot.convert.missing,
                                all_states)
        }
        char <- lapply(char, plot.convert.inappli)
        return(unlist(lapply(char, function(X) paste(as.character(X),
                                                          collapse = ""))))
    }
    get.NA.edges <- function(states_matrix, tree, pass = 4) {
        check.applicable <- function(nodes, states_matrix, pass) {
            node1 <- states_matrix[[pass + 1]][nodes[1]][[1]]
            node2 <- states_matrix[[pass + 1]][nodes[2]][[1]]
            all_char <- sort(unique(unlist(states_matrix$Char)))
            options(warn = -1)
            node2 <- ifelse(all(node2 == all_char), node1, node2)
            options(warn = 0)
            return(ifelse(all(c(node1, node2) == -1), 0, 1))
        }
        return(apply(tree$edge, 1, check.applicable, states_matrix,
                     pass))
    }

    edge_col <- "black"
    tips_labels <- plot.convert.state(states_matrix[[1]][1:n_tip], missing = TRUE)

    palette <- generate.palette(tips_labels)
    tips_colours <- palette[[1]]
    state_colours <- palette[[2]]
    edge_palette <- palette[[3]]

    if (!is.null(unlist(states_matrix$Up2))) {
        na_edges <- get.NA.edges(states_matrix, tree, pass = 4) ==
            1
        edge_final <- ifelse(na_edges, "0", "-")
        edge_col <- ifelse(na_edges, "black", "grey")
    } else {
        edge_final = 0
    }

    if (!is.null(unlist(states_matrix$Up2))) {
        final_state <- states_matrix$Up2
    } else {
        final_state <- states_matrix$Up1
    }
    max_final <- max(unlist(final_state))
    all_states <- -1:max_final
    col_states <- c("-", seq_len(max_final + 1L) - 1L)
    colour.edge <- function(edge) {
        parent <- all_states %in% final_state[[edge[1]]]
        child <- all_states %in% final_state[[edge[2]]]
        common <- parent & child
        if (sum(common) == 1) {
            col_states[common]
        }
        else if (sum(child) == 1) {
            col_states[child]
        }
        else if (sum(parent) == 1 && !identical(parent,
                                                (col_states == "-"))) {
            col_states[parent]
        }
        else "?"
    }
    edge_final <- apply(tree$edge, 1, colour.edge)
    edge_col <- as.character(edge_palette[edge_final])

    if (length(state.labels) == length(edge_palette) - 2) {
        state.labels <- c(state.labels, "Ambiguous", "Inapplicable")
    } else if (length(state.labels) == length(edge_palette) - 1) {
        state.labels <- c(state.labels, "Ambiguous")
    } else if (length(state.labels) != length(edge_palette)) {
        warning("State labels do not seem to match states.  You need to label all states from 0 to the maximum observed.")
    }

    edge_col_array <- paste0('["', paste0(edge_col, collapse='", "'), '"]')

    state_labels <- paste(names(edge_palette), gsub("^['\"]|['\"]$",
                                                    "", state.labels), sep = ": ")
    observed <- names(edge_palette) %in% edge_final
    list (edge_col = edge_col,
          edge_col_array = edge_col_array,
          legend = state_labels[observed],
          legend_col = edge_palette[observed],
          tips_labels = tips_labels,
          tips_colours = as.character(state_colours[tips_colours]),
          dud_steps = dud_steps)
}

#' Plot character mapping
#'
#' Plot a tree depicting a character reconstruction
#'
#' This function will plot a reconstruction of the distribution of character `char`
#' on a set of phylogenetic trees.
#' The function is intended for use in R Markdown (`.Rmd`) documents, which might
#' be rendered using the `knitr`, `rmarkdown` or `bookdown` pacakges.
#' These packages can produce output in PDF or HTML format; this package will
#' generate SVG files for HTML output, and a static PNG otherwise.
#'
#' The SVG output will generate a separate plot for each tree in `canvas`;
#' JavaScript can be specified using [`PrintSwitcher()`] or
#' [`PrintJavascript()`] in order
#' to allow the user to switch between trees.  Due to the space considerations
#' implicit in PDF documents, the PNG images include only a single tree,
#' specified by `singleTree`.
#'
#' HTML output can be forced by setting `options('localInstance' = TRUE)`.
#'
#' @param char Character vector specifying distribution of tokens among (named)
#'   taxa
#' @param charIndex Character vector providing a reference for the character,
#'   typically its number, to be sent to [`SVGTree()`]
#' @param stateLabels Character vector specifying the labels to apply to each
#'   state
#' @param singleTree A single tree to be plotted in Latex output
#' @param legendText Character giving legend text to print in Latex output
#' @template canvasParam
#' @template treeNamesParam
#' @param analysisLabels Character vector specifying names of each analysis,
#'   to be printed on plot
#' @param svgFilename Character string specifying location to save each file,
#'   containing the expression \code{\%s}, which will be replaced with the number
#'   of the tree.
#' @param SetPar Graphical parameters to set before plotting PNG tree
#'   in Latex output.  Specify `NULL` to retain existing parameters.
#'
#' @return Prints the tree in an appropriate markdown format
#' @importFrom knitr is_html_output
#' @examples
#'   # An example is given in the help page for SVGCanvas; type ?SVGCanvas
#' @export
#' @author Martin R. Smith
PlotCharacterMapping <- function (char, stateLabels, singleTree,
                                  legendText = '',
                                  SetPar = par(mar = rep(0.2, 4), cex = 0.7),
                                  canvas = NULL, treeNames = NULL,
                                  analysisLabels = canvas$analysisNames,
                                  charIndex = character(0),
                                  svgFilename = 'tree_%s.svg') {
    if (char[1] == '?' && length(unique(char)) == 1) {
        cat("<p>All taxa are coded as ambiguous for this character.</p>")
        legendLabels <- "?: Not scored"
        legendCol <- "darkgrey"
    } else if (all (char %in% c('?', '-'))) {
        cat("<p>All taxa are coded as ambiguous or inapplicable for this character.</p>")
        legendLabels <- c("?: Not scored", "-: Inapplicable")
        legendCol <- "darkgrey"
    } else if (!is.null(canvas) && (
        !is.null(getOption('localInstance')) || is_html_output())) {
        for (treeNo in canvas$eachTree) {
            svgSource <- SVGTree(treeNo=treeNo, canvas=canvas,
                                 char=char, charIndex=charIndex,
                                 stateLabels=stateLabels, treeNames=treeNames,
                                 analysisLabels=analysisLabels)
            write(svgSource, file=sprintf(svgFilename, treeNo))
        }
        # Just write a single tree to HTML output
        cat(svgSource)
    } else {
        origPar <- SetPar
        plot.states.matrix(apply.reconstruction(singleTree, char,
                                                match.tip.char = TRUE),
                           passes = 0, counts = 1:2, show.labels = 1,
                           col.states = TRUE, state.labels = stateLabels,
                           use.edge.length = TRUE, legend.pos = 'topright')
        legend('bottomleft',  bty = 'n', legendText)
        if(!is.null(origPar)) par(origPar)
    }
}


#' Print Switcher
#'
#' Prints an input box that allows the SVG tree displayed to be toggled,
#' in HTML output.  To be used alongside [`PlotCharacterMapping()`].
#'
#' @param nTrees Integer specifying the number of trees that can be switched,
#'               used to enforce the highest allowed value in the switcher box
#'
#' @importFrom knitr is_html_output
#' @export
#'
#' @author Martin R. Smith
PrintSwitcher <- function (nTrees) {
    if (!is.null(getOption('localInstance')) || is_html_output()) {
        cat(paste0('<div class="switcher">',
                   '<span class="selectTree">Tree number:</span>',
                   '<input class="switcherNumber" type="number" min="1" max="',
                   nTrees, '"  oninput="switchTree(this)" value="', nTrees, '"></input>',
                   '</div><div class="toggleDetails">[Show details]</div>'))
    }
}

#' Print Javascript
#'
#' Allows arbitrary JavaScript content to be imported into an HTML document
#' produced with `bookdown` or `rmarkdown`.  Suggested for use in Rmd documents
#' that will be rendered in HTML; useful in supoprting user interaction with
#' SVG elements produced by [`PlotCharacterMapping()`].
#'
#' Note that the content will be ignored when files are compiled for PDF output.
#'
#' @param filepath Path to Javascript template
#'
#' @return Prints the contents of the Javascript file to to stdout,
#'  executing any R code included in r blocks in the file.
#' @export
#'
#' @author Martin R. Smith
PrintJavascript <- function (filepath) {
    javaLines <- c('<script>', readLines(filepath), '</script>')
    REGEXP_inlineR <- '(.*)(`r (.*)`)(.*)'
    inlineR <- grep(REGEXP_inlineR, javaLines)
    javaLines[inlineR] <- vapply(javaLines[inlineR], function (line) {
        rCode <- gsub(REGEXP_inlineR, "\\3", line)
        gsub(REGEXP_inlineR,
             sprintf("\\1%s\\4", eval(parse(text=rCode))),
             line)
    }, character(1))
    writeLines(javaLines)
}