###############################################################################
# Functions for plotting birth-death skylinesdev
#' Get range and pad by some percentage of the range on either side
#' Good for plot limits
paddedrange <- function(y, pad=0.1) {
yrange <- c(min(y), max(y))
return( c(yrange[1]-diff(yrange)*pad*0.5, yrange[2]+diff(yrange)*pad*0.5) )
}
#' Plot a Skyline.
#'
#' Assume each column is a different time point (ncol(skylinemat) == length(times) must be fulfilled)
#'
#' @param times The time points to draw the skyline at. length(times) should be equal to ncol(skyline_mat) + 1.
#' (Because there are length(times) intervals and the last interval needs an end-time).
#' If length(times) == ncol(skyline_mat) the last interval is arbitrarily given the same length as
#' the second last interval.
#'
#' @param skyline_mat
#'
#' @param type Type of skyline to plot.
#'
#' "smooth" : If skyline_mat contains the HPDs plot a smooth line and polygon (type='l')
#'
#' "step" : If skyline_mat contains the HPDs plot a stepped line and polygon (type='S')
#'
#' "lines" : Plot every row in the matrix as a smooth line (type='l')
#'
#' "steplines" : Plot every row in the matrix as a stepped line (type='S')
#'
#' @param traces Number of traces to draw if type="traces"
#'
#' @param new Create a new set of axes (the skyline is fitted to the plotting devie)
#'
#' @param add Add to the current plot (do not create a new plotting device)
#'
#' @param ... Parameters passed to plotting function
#'
#' @export
plotSkyline <- function(times, skyline_mat, type="smooth", traces=1000, col=pal.dark(cblack), fill=pal.dark(cgray, 0.25), lwd=1, lty=1,
new=TRUE, add=FALSE, xlims=NULL, ylims=NULL, ...) {
##################
# Initialization #
##################
intervals <- ncol(skyline_mat)
skyline_mat <- cbind(skyline_mat, skyline_mat[,intervals])
# Check dimensions
if (length(times) == intervals)
times <- c(times, times[intervals]+(times[intervals]-times[intervals-1]))
else
if (length(times) < intervals || length(times) > (intervals+1))
stop("Dimension mismatch between times and skyline_mat!")
# Do not open a new device (add to the current plot)
if (add == TRUE) par(new=TRUE)
# Create a new set of axes that fits this skyline
if (new == TRUE) {
# Always slightly pad the y-limits (even if provided) to plot values at the boundary properly
if (is.null(ylims))
ylims <- paddedrange(skyline_mat)
else
ylims <- paddedrange(ylims,0.01)
if (is.null(xlims)) xlims <- range(times)
# Create new plot
plot(1, type='n', ylim=ylims, xlim=xlims, ...)
# Get and set clipping area
usr <- par("usr")
clip(xlims[1], xlims[2], ylims[1], ylims[2])
}
#######################
# Plot actual skyline #
#######################
# Plot traces
if (type == "lines" || type == "steplines") {
if (traces > 0 && traces < nrow(skyline_mat))
ind <- sample(nrow(skyline_mat), size=traces, replace=FALSE)
else
ind <- 1:nrow(skyline_mat)
for (i in 1:length(ind)) {
lines(times, skyline_mat[ind[i],], col=col, lwd=lwd, lty=lty, type=if(type=="lines") 'l' else 's')
}
} else
# Plot HPDs smooth
if (type == "smooth" && nrow(skyline_mat) == 3) {
polygon(c(times, rev(times)), c(skyline_mat[1,], rev(skyline_mat[3,])), col=fill, border=NA)
lines(times, skyline_mat[2,], col=col, lwd=lwd, lty=lty)
} else
# Plot HPDs stepped
if (type == "step" && nrow(skyline_mat) == 3) {
for (i in 2:(intervals+1)) {
rect(times[i-1], skyline_mat[1,i-1], times[i], skyline_mat[3,i-1], col=fill, border=NA)
}
lines(times, skyline_mat[2,], col=col, lwd=lwd, lty=lty, type='s')
} else
stop("Invalid type parameter for Skyline plot!")
##################
# Reset settings #
##################
if (add == TRUE) par(new=FALSE)
if (new == TRUE) do.call("clip", as.list(usr))
}
#' Plot a pretty skyline. Contains default options for drawing axes to the plot
#'
#' When plotting the y-axis on the right need to ensure that the plot margins are big enough!
#' See plotSkyline help for documentation on other parameters!
#'
#' @param side Side to draw the y-axis
#' @param xline Line to draw x-axis label on
#' @param yline Line to draw y-axis label on
#'
#' @todo Plot dates on x-axis
#' Plot tree in background
#'
#'
#' @export
plotSkylinePretty <- function(times, skyline_mat, type="smooth", traces=1000, col=pal.dark(cblack), fill=pal.dark(cgray, 0.25), col.axis=NA, lwd=1, lty=1,
xaxis=TRUE, yaxis=TRUE, xlab="", ylab="", xline=1, yline=1, xticks=NULL, yticks=NULL, xticklabels=NULL,
xlims=NULL, ylims=NULL, xgrid=FALSE, ygrid=FALSE, gridcol=pal.light(cgray),
cex.label=1.4, cex.axis=1, axispadding=0, side=2, ...) {
# Plot the basic skyline first
plotSkyline(times, skyline_mat, type=type, traces=traces, col=col, fill=fill, lwd=lwd, lty=lty, axes=FALSE, xlab=NA, ylab=NA, xlims=xlims, ylims=ylims, ...)
# General settings
if (is.null(xticks)) xticks <- axTicks(1)
if (is.null(yticks)) yticks <- axTicks(2)
if (is.null(ylims)) ylims <- paddedrange(skyline_mat)
if (is.null(xlims)) {
intervals <- ncol(skyline_mat)
if (length(times) == intervals)
times <- c(times, times[intervals]+(times[intervals]-times[intervals-1]))
xlims <- range(times)
}
ypos <- min(ylims, yticks)-axispadding*diff(ylims)
xpos2 <- min(xlims, xticks)-(axispadding*diff(xlims))
xpos4 <- max(xlims, xticks)+(axispadding*diff(xlims))
# x-grid
if (xgrid == TRUE) {
for (tick in yticks) lines(c(xpos2,xpos4), rep(tick,2), col=gridcol, lty=3, lwd=0.5)
}
# y-grid
if (ygrid == TRUE) {
for (tick in xticks) lines(rep(tick,2), c(ypos,ylims[2]), col=gridcol, lty=3, lwd=0.5)
}
# x-axis
if (xaxis == TRUE) {
if (is.null(xticklabels))
xticklabels <- TRUE
axis(1, at=xticks, labels=xticklabels, pos=ypos, lwd=0, lwd.ticks=lwd, las=1, cex.axis=cex.axis)
lines(range(xlims, xticks), rep(ypos,2), lwd=lwd)
mtext(xlab, side=1, line=xline, cex=cex.label)
}
# y-axis
if (yaxis == TRUE) {
if (side == 2)
xpos <- xpos2
else
xpos <- xpos4
if (is.na(col.axis)) col.axis <- pal.dark(cblack)
axis(side, at=yticks, pos=xpos, lwd=0, lwd.ticks=lwd, las=1, col=col.axis, col.axis=col.axis, cex.axis=cex.axis)
lines(rep(xpos,2), range(ylims, yticks), col=col.axis, lwd=lwd)
mtext(ylab, side=side, line=yline, col=col.axis, cex=cex.label)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.