R/phonR.R

In phonR: Tools for Phoneticians and Phonologists

## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## phonR version 1.0-7
## Functions for phoneticians and phonologists
## AUTHOR: Daniel McCloy, drmccloy@uw.edu
## LICENSED UNDER THE GNU GENERAL PUBLIC LICENSE v3.0:
## http://www.gnu.org/licenses/gpl.html
## Development of this package was funded in part by NIH-R01DC006014
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##

#' @export
plotVowels <- function(f1, f2, vowel=NULL, group=NULL,
    ## tokens
    plot.tokens=TRUE, pch.tokens=NULL, cex.tokens=NULL, alpha.tokens=NULL,
    ## means
    plot.means=FALSE, pch.means=NULL, cex.means=NULL, alpha.means=NULL,
    ## hull
    hull.line=FALSE, hull.fill=FALSE, hull.args=NULL,
    ## polygon
    poly.line=FALSE, poly.fill=FALSE, poly.args=NULL, poly.order=NA,
    ## ellipse
    ellipse.line=FALSE, ellipse.fill=FALSE, ellipse.conf=0.6827,
    ellipse.args=NULL,
    ## diphthongs
    diph.arrows=FALSE, diph.args.tokens=NULL, diph.args.means=NULL,
    diph.label.first.only=TRUE, diph.mean.timept=1, diph.smooth=FALSE,
    ## heatmaps
    heatmap=FALSE, heatmap.args=NULL, heatmap.legend=FALSE,
    heatmap.legend.args=NULL,
    ## color, style
    var.col.by=NULL, var.sty.by=NULL, fill.opacity=0.3, label.las=NULL,
    ## legend
    legend.kwd=NULL, legend.args=NULL,
    ## misc
    pretty=FALSE, output="screen", ...)
{
    ## ## ## ## ## ## ## ##
    ## HANDLE EXTRA ARGS ##
    ## ## ## ## ## ## ## ##
    exargs <- list(...)
    font.specified <- "family" %in% names(exargs)
    ## to-be-deprecated items
    if ("var.style.by" %in% names(exargs)) {
        message("[phonR]: Argument 'var.style.by' has been deprecated and ",
                "renamed 'var.sty.by'. In future versions 'var.style.by' will ",
                "no longer work; please update your code.")
        if (!is.null(var.sty.by)) {
            message("Additionally, you have passed both the old 'var.style.by'",
                    " and the new 'var.sty.by' arguments; the old one will be ",
                    "ignored.")
        } else {
            var.sty.by <- exargs$var.style.by
            exargs$var.style.by <- NULL
        }
    }
    ## two arguments get overridden no matter what
    exargs$ann <- FALSE
    exargs$type <- "n"
    ## Some graphical devices only support inches, so we convert here.
    if ("units" %in% names(exargs)) {
        if (!exargs$units %in% c("in", "cm", "mm", "px")) {
            message("[phonR]: Unsupported argument value '", units, "': ",
                    "'units' must be one of 'in', 'cm', 'mm', or 'px'. Using ",
                    "default ('in').")
            exargs$units <- "in"
        }
        if (output %in% c("pdf", "svg", "screen")) {
            if ("width" %in% names(exargs)) {
                if      (exargs$units == "cm") exargs$width <- exargs$width/2.54
                else if (exargs$units == "mm") exargs$width <- exargs$width/25.4
                else if (exargs$units == "px") exargs$width <- exargs$width/72
            }
            if ("height" %in% names(exargs)) {
                if      (exargs$units == "cm") exargs$height <- exargs$height/2.54
                else if (exargs$units == "mm") exargs$height <- exargs$height/25.4
                else if (exargs$units == "px") exargs$height <- exargs$height/72
            }
        }
    }
    ## Some args are only settable by direct par() call (not via plot(), etc).
    ## Not strictly true for "family" or "las", but works better this way.
    par.only <- c("ask", "fig", "fin", "las", "lheight", "mai", "mar",
                "mex", "mfcol", "mfrow", "mfg", "new", "oma", "omd", "omi",
                "pin", "plt", "ps", "pty", "usr", "xlog", "ylog", "ylbias")
    if (output == "screen") {
        par.only <- append(par.only, "family")
    } else {
        file.only <- c("filename", "width", "height", "units", "pointsize",
                    "res", "quality", "compression", "family")
        file.args <- exargs[names(exargs) %in% file.only]
        exargs <- exargs[!(names(exargs) %in% file.only)]
    }
    par.args <- exargs[names(exargs) %in% par.only]
    ## allow separate "las" for axis numbers & axis labels
    if (is.null(label.las)) {
        if ("las" %in% names(par.args)) label.las <- par.args$las
        else                            label.las <- par("las")
    }
    ## split out arguments for annotations. If "pretty", axes get drawn
    ## separately from plot(); other annotation gets drawn in separate mtext()
    ## calls whether "pretty" or not.
    if (pretty) {
        axis.only <- c("cex.axis", "col.axis", "font.axis")
        axis.args <- exargs[names(exargs) %in% axis.only]
    } else {
        axis.only <- c()
    }
    main.only <- c("cex.main", "col.main", "font.main")
    sub.only <- c("cex.sub",  "col.sub",  "font.sub")
    lab.only <- c("cex.lab",  "col.lab",  "font.lab")
    main.args <- exargs[names(exargs) %in% main.only]
    sub.args <- exargs[names(exargs) %in% sub.only]
    lab.args <- exargs[names(exargs) %in% lab.only]
    names(main.args) <- gsub(".main", "", names(main.args))
    names(sub.args) <- gsub(".sub", "", names(sub.args))
    names(lab.args) <- gsub(".lab", "", names(lab.args))
    if ("xlab" %in% names(exargs)) xlab <- exargs$xlab
    else                           xlab <- "F2"
    if ("ylab" %in% names(exargs)) ylab <- exargs$ylab
    else                           ylab <- "F1"
    if ("main" %in% names(exargs)) main <- exargs$main
    else                           main <- ""
    if ("sub" %in% names(exargs))   sub <- exargs$sub
    else                            sub <- ""
    exargs$xlab <- NULL
    exargs$ylab <- NULL
    exargs$main <- NULL
    exargs$sub <- NULL
    ## add "las" into label args
    lab.args <- as.list(c(lab.args, las=label.las))
    ## don't pass args twice!
    args.to.remove <- c(par.only, main.only, axis.only, lab.only, sub.only)
    exargs <- exargs[!names(exargs) %in% args.to.remove]

    ## ## ## ## ## ## ##
    ## OUTPUT PARSING ##
    ## ## ## ## ## ## ##
    output <- tolower(output)
    if (output=="jpeg") output <- "jpg"
    if (output=="tiff") output <- "tif"
    output.types <- c("pdf", "svg", "jpg", "tif", "png", "bmp", "screen")
    output.raster <- c("jpg", "tif", "png", "bmp", "screen")
    if (!(output %in% output.types)) {
        message("[phonR]: Unknown argument value '", output, "': 'output' ",
                "must be one of 'pdf', 'svg', 'png', 'tif', 'bmp', ",
                "'jpg', or 'screen'. Using default ('screen').")
        output <- "screen"
    }

    ## ## ## ## ## ## ## ##
    ## LEGEND KWD CHECK  ##
    ## ## ## ## ## ## ## ##
    legend.kwds <- c("left", "right", "top", "bottom", "center", "topleft",
                    "topright", "bottomleft", "bottomright")
    if (!is.null(legend.kwd)) {
        if (!legend.kwd %in% legend.kwds) {
            stop(paste(c("legend.kwd must be one of '",
                         paste(legend.kwds, collapse="', '"), "'."),
                 collapse=""))
    }   }

    ## ## ## ## ## ## ## ## ## ## ## ## ##
    ##  PRELIMINARY DIPHTHONG HANDLING  ##
    ## ## ## ## ## ## ## ## ## ## ## ## ##
    if (is.vector(f1) && is.vector(f2)) {
        diphthong <- FALSE
    } else if (length(dim(f1)) == 1 && length(dim(f2)) == 1) {
        f1 <- as.vector(f1)
        f2 <- as.vector(f2)
        diphthong <- FALSE
    } else {
        if (!all(dim(f1) == dim(f2))) {
            stop("Unequal dimensions for 'f1' and 'f2'.")
        } else if (length(dim(f1)) > 2) {
            stop("Argument 'f1' has more than two dimensions.")
        } else if (length(dim(f2)) > 2) {
            stop("Argument 'f2' has more than two dimensions.")
        } else if (!is.null(vowel) && length(vowel) != dim(f1)[1]) {
            stop("First axis of 'f1' does not equal length of 'vowel'.")
        }
        diphthong <- TRUE
    }
    if (!diphthong) {
        if (length(f2) != length(f1)) {
            stop("Unequal dimensions for 'f1' and 'f2'.")
        } else if (!is.null(vowel) && length(vowel) != length(f1)) {
            stop("Unequal dimensions for 'f1' and 'vowel'.")
        }
        l <- length(f1)
    } else {
        f1d <- f1
        f2d <- f2
        f1 <- f1d[,diph.mean.timept]
        f2 <- f2d[,diph.mean.timept]
        l <- nrow(f1d)
    }

    ## ## ## ## ## ## ## ## ## ##
    ## DIPHTHONG ARG HANDLING  ##
    ## ## ## ## ## ## ## ## ## ##
    if (diphthong) {
        arrow.only <- c("angle", "length")
        line.only <- c("type")
        if (pretty) {
            type <- ifelse(is.null(pch.means), "l", "o")
            pretty.diph.tokens <- list(length=0.1, angle=20, type=type)
            pretty.diph.means  <- list(length=0.1, angle=20, type=type,
                                       lwd=2*par("lwd"))
            ## user override
            pretty.diph.tokens[names(diph.args.tokens)] <- diph.args.tokens
            pretty.diph.means[names(diph.args.means)] <- diph.args.means
            ## re-unify
            diph.args.tokens <- pretty.diph.tokens
            diph.args.means  <- pretty.diph.means
        }
        if (diph.arrows) {
            diph.arrow.tokens <- diph.args.tokens[!names(diph.args.tokens)
                                                  %in% line.only]
            diph.arrow.means  <- diph.args.means[!names(diph.args.means)
                                                 %in% line.only]
        }
        diph.args.tokens[names(diph.args.tokens) %in% arrow.only] <- NULL
        diph.args.means[names(diph.args.means) %in% arrow.only] <- NULL
    }

    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
    ##  PRELIMINARY HANDLING OF GROUPING FACTOR, COLOR, AND STYLE  ##
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
    if (is.null(vowel)) v <- rep(NA, l)
    else                v <- factor(vowel, levels=unique(vowel))
    if (is.null(group)) gf <- rep("gf", l)
    else                gf <- factor(group, levels=unique(group))
    ## used later to set default polygon color when color varies by vowel
    col.by.vowel <- identical(as.numeric(factor(var.col.by)),
                              as.numeric(factor(vowel)))
    ## var.col.by & var.sty.by
    if (!is.null(var.col.by[1])) {
        if (is.na(var.col.by[1])) {
            legend.col.lab <- NULL
        } else {
            legend.col.lab <- unique(as.character(var.col.by))
        }
        var.col.by <- as.numeric(factor(var.col.by, levels=unique(var.col.by)))
    } else {
        legend.col.lab <- c()
        var.col.by <- rep(1, l)  # default to black
    }
    if (!is.null(var.sty.by[1])) {
        if (is.na(var.sty.by[1])) legend.style.lab <- NULL
        else legend.style.lab <- unique(as.character(var.sty.by))
        var.sty.by <- as.numeric(factor(var.sty.by,
                                        levels=unique(var.sty.by)))
    } else {
        legend.style.lab <- c()
        var.sty.by <- rep(1, l)  # default to solid
    }
    num.col <- length(unique(var.col.by))
    num.sty <- length(unique(var.sty.by))
    ## misc. plotting defaults
    if (is.null(cex.tokens))    cex.tokens <- par("cex")
    if (is.null(cex.means))      cex.means <- par("cex")

    ## ## ## ## ## ## ## ## ## ##
    ##  DEFAULTS FOR "PRETTY"  ##
    ## ## ## ## ## ## ## ## ## ##
    if (pretty) {
        draw.axes <- TRUE
        ## still suppress axes if axes=FALSE
        if ("axes" %in% names(exargs)) {
            draw.axes <- as.logical(exargs$axes)
        }
        ## if no colors specified, use equally spaced HCL values
        ## [-1] avoids duplicate hues 0 and 360
        if (num.col == 1) pretty.col <- "black"
        else {
            hue <- seq(0,  360, length.out=1+num.col)[-1]
            chr <- seq(60, 100, length.out=num.col)
            lum <- seq(60,  40, length.out=num.col)
            pretty.col <- hcl(hue, chr, lum, alpha=1)
        }
        ## PCH: filled / open {circ,tri,squ,diam}, plus, x, inverted open tri
        pretty.pch <- rep(c(16,1,17,2,15,0,18,5,3,4,6), length.out=num.sty)
        ## LTY: custom linetypes more readily distinguishable
        pretty.lty <- c("solid", "44", "F4", "4313", "F3131313", "23F3",
                        "232923", "23258385", "282823B3", "13", "82")
        pretty.args <- list(mgp=c(2,0.5,0), xaxs="i", yaxs="i", axes=FALSE,
                            fg=hcl(0,0,40), tcl=-0.25, xpd=NA,
                            pch=pretty.pch, lty=pretty.lty, col=pretty.col)
        pretty.par.args <- list(mar=c(1,1,5,5), las=1)
        ## legend args
        pretty.legend.args <- list(bty="n", seg.len=1)
        ## let user-specified args override "pretty" defaults
        pretty.args[names(exargs)] <- exargs
        pretty.par.args[names(par.args)] <- par.args
        pretty.legend.args[names(legend.args)] <- legend.args
        ## re-unify to avoid later logic branching
        exargs <- pretty.args
        par.args <- pretty.par.args
        legend.args <- pretty.legend.args
    }

    ## ## ## ## ## ## ##
    ## OTHER DEFAULTS ##
    ## ## ## ## ## ## ##
    vary.col <- !is.na(var.col.by[1])
    vary.sty <- !is.na(var.sty.by[1])
    ## color: use default pallete if none specified and pretty=FALSE
    if (!"col" %in% names(exargs)) exargs$col <- palette()
    ## linetypes & plotting characters
    if (!"lty" %in% names(exargs)) exargs$lty <- seq_len(num.sty)
    if (!"pch" %in% names(exargs)) exargs$pch <- seq_len(num.sty)
    if (!"lwd" %in% names(exargs)) exargs$lwd <- par("lwd")
    ## recycle user-specified colors to the length we need
    if (vary.col) exargs$col <- rep(exargs$col, length.out=num.col)[var.col.by]
    if (vary.sty) exargs$lty <- rep(exargs$lty, length.out=num.sty)[var.sty.by]
    if (vary.sty) exargs$pch <- rep(exargs$pch, length.out=num.sty)[var.sty.by]
    ## set defaults for token and mean plotting characters
    if (is.null(pch.tokens)) pch.t <- exargs$pch
    else                     pch.t <- pch.tokens
    if (is.null(pch.means))  pch.m <- exargs$pch
    else                     pch.m <- pch.means
    ## transparency
    trans.col <- makeTransparent(exargs$col, fill.opacity)
    trans.fg <- makeTransparent(par("fg"), fill.opacity)
    ## ellipse colors
    if (ellipse.line) {
        if ("lty" %in% names(ellipse.args)) {
            ellipse.args$lty <- rep(ellipse.args$lty, length.out=num.col)
            if (vary.sty) ellipse.line.sty <- ellipse.args$lty[var.sty.by]
            else          ellipse.line.sty <- ellipse.args$lty
            ellipse.args$lty <- NULL
        } else {
            ellipse.line.sty <- exargs$lty
        }
        if ("border" %in% names(ellipse.args)) {
            ellipse.args$border <- rep(ellipse.args$border, length.out=num.col)
            if (vary.col) ellipse.line.col <- ellipse.args$border[var.col.by]
            else          ellipse.line.col <- ellipse.args$border
        } else            ellipse.line.col <- exargs$col
    } else {
        ellipse.line.col <- NA[var.col.by]
        ellipse.line.sty <- NA[var.sty.by]
    }
    if (ellipse.fill) {
        if ("col" %in% names(ellipse.args)) {
            if (vary.col) ellipse.fill.col <- ellipse.args$col[var.col.by]
            else          ellipse.fill.col <- ellipse.args$col
        } else            ellipse.fill.col <- trans.col
    } else                ellipse.fill.col <- NA[var.col.by]
    ## hull colors
    if (hull.line) {
        if ("lty" %in% names(hull.args)) {
            if (vary.sty) hull.line.sty <- hull.args$lty[var.sty.by]
            else          hull.line.sty <- hull.args$lty
            hull.args$lty <- NULL
        } else {
            hull.line.sty <- exargs$lty
        }
        if ("border" %in% names(hull.args)) {
            if (vary.col) hull.line.col <- hull.args$border[var.col.by]
            else          hull.line.col <- hull.args$border
            hull.args$border <- NULL
        } else if (col.by.vowel) {
            hull.line.col <- par("fg")
        } else {
            hull.line.col <- exargs$col
        }
    } else {
        hull.line.col <- NA[var.col.by]
        hull.line.sty <- NA[var.sty.by]
    }
    if (hull.fill) {
        if ("col" %in% names(hull.args)) {
            if (vary.col)        hull.fill.col <- hull.args$col[var.col.by]
            else                 hull.fill.col <- hull.args$col
        } else if (col.by.vowel) hull.fill.col <- trans.fg
        else                     hull.fill.col <- trans.col
    } else                       hull.fill.col <- NA[var.col.by]
    ## polygon colors
    if (poly.line) {
        if ("lty" %in% names(poly.args)) {
            poly.args$lty <- rep(poly.args$lty, length.out=num.sty)
            if (vary.sty) poly.line.sty <- poly.args$lty[var.sty.by]
            else          poly.line.sty <- poly.args$lty
            poly.args$lty <- NULL
        } else {
            poly.line.sty <- exargs$lty
        }
        if ("border" %in% names(poly.args)) {
            if (vary.col) poly.line.col <- poly.args$border[var.col.by]
            else          poly.line.col <- poly.args$border
            poly.args$border <- NULL
        } else if (col.by.vowel) {
            poly.line.col <- par("fg")
        } else {
            poly.line.col <- exargs$col
        }
    } else {
        poly.line.col <- NA[var.col.by]
        poly.line.sty <- NA[var.sty.by]
    }
    if (poly.fill) {
        if ("col" %in% names(poly.args)) {
            if (vary.col)        poly.fill.col <- poly.args$col[var.col.by]
            else                 poly.fill.col <- poly.args$col
        } else if (col.by.vowel) poly.fill.col <- trans.fg
        else                     poly.fill.col <- trans.col
    } else                       poly.fill.col <- NA[var.col.by]
    ## handle NAs in linetypes (lty=0 means do not draw)
    ellipse.line.sty[is.na(ellipse.line.sty)] <- 0
    poly.line.sty[is.na(poly.line.sty)] <- 0
    hull.line.sty[is.na(hull.line.sty)] <- 0

    ## ## ## ## ## ## ## ## ## ## ##
    ## TOKEN / MEAN TRANSPARENCY  ##
    ## ## ## ## ## ## ## ## ## ## ##
    if (!is.null(alpha.tokens)) {
        col.tokens <- makeTransparent(exargs$col, alpha.tokens)
    } else {
        col.tokens <- exargs$col
    }
    if (!is.null(alpha.means)) {
        col.means <- makeTransparent(exargs$col, alpha.means)
    } else {
        col.means <- exargs$col
    }

    ## ## ## ## ## ## ## ## ## ## ##
    ## INITIALIZE OUTPUT DEVICES  ##
    ## ## ## ## ## ## ## ## ## ## ##
    if      (output == "pdf") do.call(cairo_pdf, file.args)
    else if (output == "svg") do.call(svg, file.args)
    else if (output == "jpg") do.call(jpeg, file.args)
    else if (output == "tif") do.call(tiff, file.args)
    else if (output == "png") do.call(png, file.args)
    else if (output == "bmp") do.call(bmp, file.args)
    ## font handling for Windows
    is.win <- .Platform$OS.type == "windows"
    if (is.win && font.specified && output %in% output.raster) {
        windowsFonts(phonr=windowsFont(par.args$family))
        par.args$family <- "phonr"
    }
    ## initial call to par()
    op <- par(par.args)

    ## ## ## ## ## ## ## ## ## ##
    ## COLLECT INTO DATAFRAMES ##
    ## ## ## ## ## ## ## ## ## ##
    d <- data.frame(f1=f1, f2=f2, v=v, gf=factor(gf, levels=unique(gf)),
                    col.tokens=col.tokens, col.means=col.means,
                    style=var.sty.by, lty=exargs$lty, lwd=exargs$lwd,
                    ellipse.fill.col=ellipse.fill.col,
                    ellipse.line.col=ellipse.line.col,
                    ellipse.line.sty=ellipse.line.sty,
                    poly.fill.col=poly.fill.col,
                    poly.line.col=poly.line.col,
                    poly.line.sty=poly.line.sty,
                    hull.fill.col=hull.fill.col,
                    hull.line.col=hull.line.col,
                    hull.line.sty=hull.line.sty,
                    pch.m=pch.m, pch.t=pch.t,
                    stringsAsFactors=FALSE)
    if (diphthong) {
        d$f2d <- lapply(apply(f2d, 1, list), unlist, use.names=FALSE)
        d$f1d <- lapply(apply(f1d, 1, list), unlist, use.names=FALSE)
        #d.diph <- d
        #d.diph$f2d <- lapply(apply(f2d, 1, list), unlist, use.names=FALSE)
        #d.diph$f1d <- lapply(apply(f1d, 1, list), unlist, use.names=FALSE)
        #d.diph[names(diph.args.tokens)] <- diph.args.tokens
    }
    if (is.null(vowel) && is.null(group)) {
        byd <- list(d=d)
    } else {
        byd <- by(d, d[c("v", "gf")], identity)
    }
    ## dataframe of means. at this point each element of "byd" should have
    ## exactly 1 vowel and 1 grouping factor (gf) value
    m <- lapply(byd, function(i) {
        if (!is.null(i)) {
            idx <- !duplicated(i$gf)
            fg <- par("fg")
            tfg <- makeTransparent(fg, fill.opacity)
            with(i, data.frame(f2=mean(f2, na.rm=TRUE),
                               f1=mean(f1, na.rm=TRUE),
                               v=v[idx], gf=gf[idx], n=nrow(i),
                               ## mean of a color / style?
                               col.means=uniquify(col.means, fg),
                               style=uniquify(style, 1),
                               poly.fill.col=uniquify(poly.fill.col, tfg),
                               hull.fill.col=uniquify(hull.fill.col, tfg),
                               ellipse.fill.col=uniquify(ellipse.fill.col, tfg),
                               poly.line.col=uniquify(poly.line.col, fg),
                               hull.line.col=uniquify(hull.line.col, fg),
                               ellipse.line.col=uniquify(ellipse.line.col, fg),
                               poly.line.sty=uniquify(poly.line.sty, 1),
                               hull.line.sty=uniquify(hull.line.sty, 1),
                               ellipse.line.sty=uniquify(ellipse.line.sty, 1),
                               pch.m=uniquify(pch.m, 1),
                               lty.means=uniquify(lty, 1),
                               lwd.means=uniquify(lwd, par("lwd")),
                               stringsAsFactors=FALSE))
    }})
    m <- do.call(rbind, m)
    m$gfn <- as.numeric(factor(m$gf, levels=unique(m$gf)))
    if (diphthong) {
        m$f2d <- do.call(rbind, lapply(byd, function(i) if (!is.null(i))
                         with(i, list(colMeans(do.call(rbind, f2d))))))
        m$f1d <- do.call(rbind, lapply(byd, function(i) if (!is.null(i))
                         with(i, list(colMeans(do.call(rbind, f1d))))))
    }
    ## means & covariances for ellipse drawing
    if (ellipse.fill || ellipse.line) {
        mu <- lapply(byd, function(i) { if (!is.null(i)) {
            with(i, list(colMeans(cbind(f2, f1), na.rm=TRUE)))
        }})
        mu <- do.call(rbind, mu)
        sigma <- lapply(byd, function(i) { if (!is.null(i)) {
            with(i, list(var(cbind(f2, f1), na.rm=TRUE)))
        }})
        ## the covariance calculation above still may yield an NA covariance
        ## matrix if a vowel has only 1 token. This is handled later.
        sigma <- do.call(rbind, sigma)
        m$mu <- mu
        m$sigma <- sigma
    }
    ## bym
    bym <- by(m, m$gfn, identity)

    ## ## ## ## ## ## ## ## ## ##
    ##  DETERMINE PLOT BOUNDS  ##
    ## ## ## ## ## ## ## ## ## ##
    ## token extrema
    if (diphthong) {
        plot.bounds <- cbind(range(f2d, finite=TRUE), range(f1d, finite=TRUE))
    } else {
        plot.bounds <- apply(d[,c("f2", "f1")], 2, range, finite=TRUE)
    }
    ## ellipse extrema
    if (ellipse.fill || ellipse.line) {
        ellipse.param <- apply(m, 1, function(i) {
            if (any(is.na(i$sigma))) {
                i$sigma <- matrix(rep(0, 4), nrow=2)
                msg <- ifelse(i$gf == "gf", as.character(i$v),
                              paste0("(", i$gf, ", ", i$v, ")"))
                message("[phonR]: No ellipse drawn for ", msg,
                        " because there is only one token.")
            } else if (i$n == 2) {
                msg <- ifelse(i$gf == "gf", as.character(i$v),
                              paste0("(", i$gf, ", ", i$v, ")"))
                message("[phonR]: No ellipse drawn for ", msg,
                        " because there are only two tokens.")
            }
            list("mu"=i$mu, "sigma"=i$sigma, "n"=i$n, "alpha"=1 - ellipse.conf,
                 "draw"=FALSE)
        })
        ellipse.points <- lapply(ellipse.param,
                                 function(i) do.call(ellipse, i))
        ellipse.bounds <- lapply(ellipse.points,
                                 function(i) apply(i, 2, range, finite=TRUE))
        ellipse.bounds <- apply(do.call(rbind, ellipse.bounds), 2,
                                range, finite=TRUE)
        plot.bounds <- apply(rbind(plot.bounds, ellipse.bounds), 2,
                             range, finite=TRUE)
    }

    ## ## ## ## ## ## ## ##
    ## PREPARE GARNISHES ##
    ## ## ## ## ## ## ## ##
    ## determine plot bounds
    user.set.xlim <- "xlim" %in% names(exargs)
    user.set.ylim <- "ylim" %in% names(exargs)
    if (!user.set.xlim) exargs$xlim <- rev(plot.bounds[,1])
    if (!user.set.ylim) exargs$ylim <- rev(plot.bounds[,2])
    ## calculate nice tickmark intervals
    xticks <- prettyTicks(exargs$xlim)
    yticks <- prettyTicks(exargs$ylim)
    ## ensure axes begin and end at a tickmark (unless xlim/ylim overtly set)
    if (pretty && !user.set.xlim) exargs$xlim <- rev(range(xticks))
    if (pretty && !user.set.ylim) exargs$ylim <- rev(range(yticks))
    ## annotation
    if (pretty) {
        x.args <- list(side=3, line=2)
        y.args <- list(side=4, line=3)
        t.args <- list(side=3, line=4)
        s.args <- list(side=3, line=3)
    } else {
        x.args <- list(side=1, line=par("mgp")[1])
        y.args <- list(side=2, line=par("mgp")[1])
        t.args <- list(side=3, line=1, outer=FALSE)
        s.args <- list(side=4, line=par("mgp")[1] + 1)
    }

    ## ## ## ## ## ## ## ## ## ## ## ##
    ##  PLOT THE AXES AND GARNISHES  ##
    ## ## ## ## ## ## ## ## ## ## ## ##
    do.call(plot, c(list(NA, NA), exargs))
    do.call(mtext, c(xlab, x.args, lab.args))
    do.call(mtext, c(ylab, y.args, lab.args))
    do.call(mtext, c(main, t.args, main.args))
    do.call(mtext, c(sub, s.args, sub.args))
    if (pretty && draw.axes) {
        if (!is.null(xticks[1])) {
            x.axis.args <- c(list(side=3, at=xticks), axis.args)
        } else {
            x.axis.args <- c(list(side=3), axis.args)
        }
        if (!is.null(yticks[1])) {
            y.axis.args <- c(list(side=4, at=yticks), axis.args)
        } else {
            y.axis.args <- c(list(side=4), axis.args)
        }
        do.call(axis, x.axis.args)
        do.call(axis, y.axis.args)
    }

    ## ## ## ## ## ## ##
    ##  PLOT HEATMAP  ##
    ## ## ## ## ## ## ##
    if (heatmap) {
        if (pretty & !"colormap" %in% names(heatmap.args)) {
            heatmap.args$colormap <- plotrix::color.scale(x=0:100,
                                                          cs1=c(0, 180),
                                                          cs2=100,
                                                          cs3=c(25, 100),
                                                          color.spec="hcl")
        }
        if (!"add" %in% names(heatmap.args)) heatmap.args$add <- TRUE
        with(d, do.call(repulsiveForceHeatmap, c(list(f2, f1, type=v),
                                                 heatmap.args)))
    }
    if (heatmap.legend) {
        if (!"x" %in% names(heatmap.legend.args)) {
            heatmap.legend.args$x <- rep(exargs$xlim[1], 2)
            heatmap.legend.args$y <- exargs$ylim - c(0, diff(exargs$ylim) / 2)
        }
        if (!"colormap" %in% heatmap.legend.args) {
            heatmap.legend.args$colormap <- heatmap.args$colormap
        }
        do.call(repulsiveForceHeatmapLegend, heatmap.legend.args)
    }

    ## ## ## ## ## ##
    ##  PLOT HULL  ##
    ## ## ## ## ## ##
    if (hull.fill || hull.line) {
        hull.columns <- c("f2", "f1", "hull.fill.col", "hull.line.col",
                          "hull.line.sty")
        hh <- by(d, d$gf, function(i) i[!is.na(i$f2) & !is.na(i$f1),])
        if (diphthong) {
            ## if diphthong, hull should ignore diph.mean.timept
            hulls <- lapply(hh, function(i) {
                ipts <- with(i, data.frame(f2=do.call(c, f2d),
                                           f1=do.call(c, f1d)))
                hull <- ipts[chull(ipts),]
                hull$hull.fill.col <- unique(i$hull.fill.col)
                hull$hull.line.col <- unique(i$hull.line.col)
                hull$hull.line.sty <- unique(i$hull.line.sty)
                return (hull)
            })
        } else {
            hulls <- lapply(hh, function(i) with(i, i[chull(f2, f1),
                            hull.columns]))
        }
        lapply(hulls, function(i) {
            hull.args$col    <- i$hull.fill.col
            hull.args$border <- i$hull.line.col
            hull.args$lty    <- i$hull.line.sty
            hull.args$x      <- cbind(i$f2, i$f1)
            with(i, do.call(polygon, hull.args))
            })
    }

    ## ## ## ## ## ## ##
    ## PLOT ELLIPSES  ##
    ## ## ## ## ## ## ##
    if (ellipse.fill || ellipse.line) {
        invisible(lapply(seq_along(ellipse.points),
                         function(i) {
                             ellipse.args$border <- m$ellipse.line.col[i]
                             ellipse.args$lty    <- m$ellipse.line.sty[i]
                             ellipse.args$col    <- m$ellipse.fill.col[i]
                             do.call(polygon, c(list(x=ellipse.points[[i]]),
                                                ellipse.args))
                             }))
    }

    ## ## ## ## ## ## ##
    ## PLOT POLYGONS  ##
    ## ## ## ## ## ## ##
    if (!is.na(poly.order[1]) && (poly.fill || poly.line)) {
        if (length(poly.order) != length(unique(poly.order))) {
            message("[phonR]: Duplicate entries in 'poly.order' detected; they",
                    " will be ignored.")
        }
        poly.order <- as.character(poly.order) # as.character in case factor
        v <- unique(as.character(m$v))
        if (length(setdiff(poly.order, v)) > 0) {
            message("[phonR]: There are vowels in 'poly.order' that are not in",
                    " 'vowel'; they will be ignored.")
        }
        poly.order <- intersect(poly.order, v)
        pp <- m
        pp$v <- factor(pp$v, levels=poly.order, ordered=TRUE)
        pp <- pp[order(pp$v),]
        pp <- pp[!is.na(pp$v),]
        pp <- split(pp, pp$gf)
        if (poly.fill) {
            bigenough <- sapply(pp, function(i) nrow(i) > 2)
            lapply(pp[bigenough], function(i) {
                pargs <- poly.args
                pargs$x <- cbind(i$f2, i$f1)
                pargs$col <- i$poly.fill.col
                pargs$border <- NA
                with(i, do.call(polygon, pargs))
            })
        }
        if (poly.line) {
            if (plot.means) type <- "c"
            else type <- "l"
            bigenough <- sapply(pp, function(i) nrow(i) > 1)
            invisible(lapply(pp[bigenough], function(i) {
                pargs <- poly.args
                pargs$x <- i$f2
                pargs$y <- i$f1
                pargs$type <- type
                pargs$cex <- 1.2 * cex.means
                pargs$col <- i$poly.line.col
                pargs$lty <- i$poly.line.sty
                with(i[i$v %in% poly.order,], do.call(points, pargs))
                }))
    }   }

    ## ## ## ## ## ##
    ## PLOT TOKENS ##
    ## ## ## ## ## ##
    if (plot.tokens) {
        if (diphthong) {
            ## setup
            timepts <- length(d$f2d[[1]])
            ## no smoothing splines
            if (!diph.smooth || timepts < 4) {
                if (diph.smooth) message("[phonR]: Cannot smooth diphthong ",
                                         "traces with fewer than 4 timepoints.",
                                         " Plotting connecting segments ",
                                         "instead.")
                ## plot first point
                if (diph.label.first.only) {
                    if (!is.null(pch.tokens)) {
                        with(d, text(f2, f1, labels=pch.t, col=col.tokens,
                                     cex=cex.tokens))
                    } else {
                        with(d, points(f2, f1, pch=pch.t, col=col.tokens,
                                       cex=cex.tokens))
                    }
                    ## if diph.label.first.only, ignore cex and pch from now on
                    diph.args.tokens$type <- "l"
                    if (diph.arrows) diph.arrow.tokens$type <- "l"
                }
                ## prepare tokens args
                d.split <- split(d, seq(nrow(d)))
                d.args <- lapply(d.split, function(i) {
                    with(i, list(f2d[[1]], f1d[[1]], pch=pch.t,
                                 cex=cex.tokens, col=col.tokens, lty=lty))
                })
                ## combine diph.args.means with m.args and plot
                invisible(lapply(d.args, function(i) {
                    i[names(diph.args.tokens)] <- diph.args.tokens
                    do.call(points, i)
                }))
                if (diph.arrows) {
                    ## prepare arrow args
                    d.arr.args <- lapply(d.split, function(i) {
                        xd <- 0.01*diff(i$f2d[[1]][(timepts-1):timepts])
                        yd <- 0.01*diff(i$f1d[[1]][(timepts-1):timepts])
                        with(i, list(x0=f2d[[1]][timepts] - xd,
                                     y0=f1d[[1]][timepts] - yd,
                                     x1=f2d[[1]][timepts], y1=f1d[[1]][timepts],
                                     col=col.tokens, lwd=lwd))
                    })
                    ## combine with diph.arrow.means and plot
                    invisible(lapply(d.arr.args, function(i){
                        i[names(diph.arrow.tokens)] <- diph.arrow.tokens
                        i$lty <- "solid"
                        if ("type" %in% names(i)) i$type <- NULL
                        do.call(arrows, i)
                    }))
                }
            ## diphthong smoothing spline
            } else if (diph.smooth) {  # timepts > 3
                apply(d, 1, function(i) {
                    tryCatch({
                        steep <- with(i,
                                      abs(lm(f1d~f2d)$coefficients["f2d"]) > 1)
                        if (steep) dat <- with(i, cbind(f1d, f2d))
                        else       dat <- with(i, cbind(f2d, f1d))
                        pc <- prcomp(dat, center=FALSE, scale.=FALSE)
                        ss <- smooth.spline(pc$x)
                        ssi <- as.matrix(as.data.frame(predict(ss))) %*%
                            solve(pc$rotation)  # * pc$scale + pc$center
                        end <- nrow(ssi)
                        if (diph.arrows) {
                            curve.range <- 1:(end-1)
                            with(as.data.frame(ssi),
                                 do.call(arrows, c(list(x0=f2[end-1],
                                                        y0=f1[end-1],
                                                        x1=f2[end],
                                                        y1=f1[end],
                                                        col=i$col.tokens),
                                                   diph.arrow.tokens)))
                        } else {
                            curve.range <- 1:end
                        }
                        do.call(lines, c(list(ssi[curve.range]),
                                         diph.args.tokens))
                    },
                    error=function(e){
                        message("[phonR]: Could not plot diphthong smoother. ",
                                "Plotting connecting segments instead.")
                        message(paste(e, ""))
                        if (diph.arrows) {
                            end <- nrow(i)
                            with(i, points(f2[1:end-1], f1[1:end-1],
                                           col=col.tokens, pch=pch.t,
                                           cex=cex.tokens, type="o"))
                            with(i, do.call(arrows, c(list(x0=f2[end-1],
                                                           y0=f1[end-1],
                                                           x1=f2[end],
                                                           y1=f1[end],
                                                           col=col.tokens),
                                                      diph.arrow.args)))
                        } else {
                            with(i, points(f2, f1, col=col.tokens,
                                           pch=pch.t, cex=cex.tokens,
                                           type="o"))
                        }
                    },
                    warning=function(w) message(w),
                    finally={}
                    )
                })
            }
        } else {  # !diphthong
            if (is.null(pch.tokens)) {
                with(d, points(f2, f1, pch=pch.t, cex=cex.tokens,
                               col=col.tokens))
            } else {
                with(d, text(f2, f1, labels=pch.t, cex=cex.tokens,
                             col=col.tokens))
    }   }   }

    ## ## ## ## ## ##
    ## PLOT MEANS  ##
    ## ## ## ## ## ##
    if (plot.means) {
        if (diphthong) {
            ## TODO: implement smoothing splines for means
            ## setup
            timepts <- length(m$f2d[[1]])
            ## plot first point
            if (diph.label.first.only) {
                if (!is.null(pch.means)) {
                    with(m, text(f2, f1, labels=pch.m, col=col.means,
                                 cex=cex.means))
                } else {
                    with(m, points(f2, f1, pch=pch.m, col=col.means,
                                   cex=cex.means))
                }
                ## if diph.label.first.only, ignore cex and pch from now on
                diph.args.means$type <- "l"
                if (diph.arrows) diph.arrow.means$type <- "l"
            }
            ## prepare means args
            m.split <- split(m, seq(nrow(m)))
            m.args <- lapply(m.split, function(i) {
                with(i, list(f2d[[1]], f1d[[1]], pch=pch.m, cex=cex.means,
                             col=col.means, lty=lty.means))
            })
            ## combine diph.args.means with m.args and plot
            invisible(lapply(m.args, function(i) {
                i[names(diph.args.means)] <- diph.args.means
                do.call(points, i)
                }))
            ## plot arrowheads
            if (diph.arrows) {
                ## prepare arrow args
                m.arr.args <- lapply(m.split, function(i) {
                    xd <- 0.01*diff(i$f2d[[1]][(timepts-1):timepts])
                    yd <- 0.01*diff(i$f1d[[1]][(timepts-1):timepts])
                    with(i, list(x0=f2d[[1]][timepts] - xd,
                                 y0=f1d[[1]][timepts] - yd,
                                 x1=f2d[[1]][timepts],
                                 y1=f1d[[1]][timepts],
                                 col=col.means, lwd=lwd.means))
                })
                ## combine with diph.arrow.means and plot
                invisible(lapply(m.arr.args, function(i){
                    i[names(diph.arrow.means)] <- diph.arrow.means
                    i$lty <- "solid"
                    if ("type" %in% names(i)) i$type <- NULL
                    do.call(arrows, i)
                }))
            }
        } else {
            if (is.null(pch.means)) {
                with(m, points(f2, f1, col=col.means, pch=pch.m,
                               cex=cex.means))
            } else {
                with(m, text(f2, f1, labels=pch.m, col=col.means,
                             cex=cex.means))
            }
        }
    }

    ## ## ## ## ##
    ##  LEGEND  ##
    ## ## ## ## ##
    if (!is.null(legend.kwd)) {
        if (is.null(legend.col.lab) && is.null(legend.style.lab)) {
            message("[phonR]: Legend will not be drawn because var.col.by and ",
                    "var.sty.by are both NULL or NA. You will have to use ",
                    "the legend() function.")
        } else {
            legend.merge <- TRUE
            ## legend pch
            legend.pch <- NULL
            if (length(legend.style.lab)) {
                if (plot.means && all(grepl("[[:digit:]]", pch.means))) {
                    legend.pch <- unique(m$pch.m)
                } else if (plot.tokens &&
                               all(grepl("[[:digit:]]", pch.tokens))) {
                    legend.pch <- unique(d$pch.t)
                }
            }
            ## legend col
            legend.col <- NULL
            if (length(legend.col.lab)) {
                if (plot.means) {
                    legend.col <- sapply(bym, function(i) unique(i$col.means))
                } else if (plot.tokens) {
                    legend.col <- sapply(byd, function(i) unique(i$col.tokens))
                }
            }
            ## legend background fill
            legend.bgf <- NULL
            if (hull.fill || poly.fill || ellipse.fill) {
                if (!is.na(m$ellipse.fill.col[1])) {
                    legend.bgf <- sapply(bym, function(i) unique(i$ellipse.fill.col))
                } else if (!is.na(m$poly.fill.col[1])) {
                    legend.bgf <- sapply(bym, function(i) unique(i$poly.fill.col))
                } else {
                    legend.bgf <- unique(hull.fill.col)
                }
            }
            ## legend linteype & border color
            legend.lty <- NULL
            legend.brd <- NULL
            if (hull.line || poly.line || ellipse.line) {
                if (!(length(unique(ellipse.line.sty)) == 1 &&
                          ellipse.line.sty[1] == 0)) {
                    legend.lty <- sapply(bym, function(i) unique(i$ellipse.line.sty))
                } else if (!(length(unique(poly.line.sty)) == 1 &&
                                 poly.line.sty[1] == 0)) {
                    legend.lty <- sapply(bym, function(i) unique(i$poly.line.sty))
                } else {
                    legend.lty <- unique(hull.line.sty)
                }
                if (!is.na(m$ellipse.line.col[1])) {
                    legend.brd <- sapply(bym, function(i) unique(i$ellipse.line.col))
                } else if (!is.na(m$poly.line.col[1])) {
                    legend.brd <- sapply(bym, function(i) unique(i$poly.line.col))
                } else {
                    legend.brd <- unique(hull.line.col)
                }
                if (length(legend.brd) != length(legend.bgf)) {
                    legend.brd <- NULL
                }
            }
            ## handle lty specially; needed for both style & color
            if (!is.null(legend.lty)) {
                if (!length(legend.style.lab)) {
                    legend.lty <- rep(legend.lty,
                                      length.out=length(legend.col.lab))
                } else if (length(legend.col.lab)) {
                    legend.lty <- c(legend.lty, rep(NA, length(legend.col.lab)))
                }
            }
            ## reconcile
            if (identical(legend.style.lab, legend.col.lab)) {
                legend.lab <- legend.col.lab
                legend.merge <- FALSE
            } else {
                legend.lab <- c(legend.style.lab, legend.col.lab)
                legend.pch <- c(legend.pch, rep(NA, length(legend.col.lab)))
                ## handle case: no lines / fills / pchs for color
                if (length(legend.col.lab) && is.null(legend.bgf) &&
                        is.null(legend.brd) && is.null(legend.lty)) {
                    legend.bgf <- c(rep(NA, length(legend.style.lab)),
                                    legend.col)
                    legend.brd <- legend.bgf
                    legend.col <- c(rep(par("fg"), length(legend.style.lab)),
                                    legend.col)
                ## handle case: only hulls (only 1 fill col) but col.by.vowel
                } else if (length(legend.col) != length(legend.bgf) &&
                               !is.null(legend.bgf)) {
                    legend.bgf <- c(rep(NA, length(legend.style.lab)),
                                    legend.col)
                    legend.brd <- legend.bgf
                    legend.col <- c(rep(par("fg"), length(legend.style.lab)),
                                    legend.col)
                ## handle other cases
                } else {
                    nas <- rep(NA, length(legend.style.lab))
                    fgs <- rep(par("fg"), length(legend.style.lab))
                    legend.bgf <- c(nas, legend.bgf)
                    legend.brd <- c(nas, legend.brd)
                    legend.col <- c(fgs, legend.col)
                }
            }
            ## eliminate vacuous args
            if (identical(legend.bgf, logical(0))) legend.bgf <- NULL
            if (identical(legend.brd, logical(0))) legend.brd <- legend.bgf
            ## assemble legend args
            new.legend.args <- list(legend.kwd, legend=legend.lab,
                                    pch=legend.pch, col=legend.col,
                                    lty=legend.lty)
            ## user override & recombine
            new.legend.args[names(legend.args)] <- legend.args
            legend.args <- new.legend.args
            ## can't always pass fill because fill=NULL triggers drawing empty
            ## boxes and border=NULL draws black! :(
            if (!is.null(legend.bgf)) {
                if (!"fill" %in% names(legend.args)) {
                    legend.args$fill <- legend.bgf
                }
                if (!"border" %in% names(legend.args)) {
                    legend.args$border <- legend.brd
                }
            }
            ## avoid warning that merge only works when segments are drawn
            if (!is.null(legend.lty)) {
                if (!"merge" %in% names(legend.args)) {
                    legend.args$merge <- legend.merge
                }
            }
            ## draw legend
            do.call(legend, legend.args)
        }
    }

    ## ## ## ## ##
    ## CLEANUP  ##
    ## ## ## ## ##
    ## close file devices
    if (output != "screen") dev.off()
    ## reset graphical parameters to defaults
    par(op)
}


## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## OMNIBUS NORMALIZATION FUNCTION (convenience function) ##
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
#' @export
normVowels <- function(method, f0=NULL, f1=NULL, f2=NULL, f3=NULL,
                    vowel=NULL, group=NULL, ...) {
    m <- tolower(method)
    methods <- c("bark", "mel", "log", "erb", "zscore", "lobanov",
                "logmean", "shared", "nearey1", "nearey2", "scentroid",
                "wattfabricius")
    if (!(m %in% methods)) {
        stop("Method must be one of: bark, mel, log, erb, ",
             "zscore | lobanov, logmean | nearey1, shared | nearey2, ",
             "scentroid | wattfabricius.")
    }
    f <- cbind(f0=f0, f1=f1, f2=f2, f3=f3)
    if (m %in% "bark") return(normBark(f))
    else if (m %in% "mel") return(normMel(f))
    else if (m %in% "log") return(normLog(f))
    else if (m %in% "erb") return(normErb(f))
    else if (m %in% c("z", "zscore", "lobanov")) return(normLobanov(f, group))
    else if (m %in% c("logmean", "nearey1")) return(normLogmean(f, group, ...))
    else if (m %in% c("shared", "nearey2")) return(normSharedLogmean(f, group, ...))
    else {
        f <- as.matrix(cbind(f1=f1, f2=f2))
        return(normWattFabricius(f, vowel, group))
}    }


## ## ## ## ## ## ## ## ## ## ## ## ## ##
## INDIVIDUAL NORMALIZATION FUNCTIONS  ##
## ## ## ## ## ## ## ## ## ## ## ## ## ##
#' @export
normBark <- function(f) {
    bark <- 26.81 * f / (1960 + f) - 0.53
    bark[bark < 2] <- bark[bark < 2] + 0.15 * (2 - bark[bark < 2])
    bark[bark > 20.1] <- bark[bark > 20.1] + 0.22 * (bark[bark > 20.1] - 20.1)
    return(bark)
}

#' @export
normLog <- function(f) {
    log10(f)
}

#' @export
normMel <- function(f) {
    2595*log10(1+f/700)
}

#' @export
normErb <- function(f) {
    21.4*log10(1+0.00437*f)
}

#' @export
normLobanov <- function(f, group=NULL) {
    if (is.null(group)) {
        return(scale(f))
    } else {
        f <- as.data.frame(f)
        groups <- split(f, group)
        scaled <- lapply(groups, function(x) as.data.frame(scale(x)))
        return(unsplit(scaled, group))
}    }

#' @export
normLogmean <- function(f, group=NULL, exp=FALSE, ...) {
    ## AKA "Nearey1", what Adank confusingly calls "SingleLogmean".
    ## Note that Adank et al 2004 (eq. 8) looks more like a shared logmean,
    ## but in text she says it is applied to each formant separately.
    if (ncol(f) < 2) {
        stop("Missing values: normalization method 'normLogmean' ",
            "requires non-null values for f1 and f2).")
    }
    if (is.null(group)) {
        logmeans <- log(f) - rep(colMeans(log(f), ...), each=nrow(f))
    } else {
        f <- as.data.frame(f)
        groups <- split(f, group)
        logmeans <- lapply(groups,
                           function(x) log(x) - rep(apply(log(x), 2, mean),
                                                    each=nrow(x)))
    logmeans <- unsplit(logmeans, group)
    }
    if (exp) logmeans <- exp(logmeans)
    logmeans
}

#' @export
normNearey1 <- function(f, group=NULL, exp=FALSE, ...) {
    normLogmean(f, group=group, exp=exp, ...)
}

#' @export
normSharedLogmean <- function(f, group=NULL, exp=FALSE, ...) {
    ## AKA "Nearey2"
    if (is.null(group)) {
        ## this is my implementation of Nearey 1978's CLIH (eq. 3.1.10, page 95)
        ## (cf. eqs. 1 & 2 of Morrison & Nearey 2006)
        logmeans <- log(f) - mean(log(unlist(f)), ...)
        ## NOTE: Adank et al 2004 (eq. 9) would suggest this implementation:
        ## logmeans <- log(f) - sum(colMeans(log(f), ...))
    } else {
        f <- as.data.frame(f)
        groups <- split(f, group)
        logmeans <- lapply(groups, function(x) log(x) - mean(log(unlist(x)), ...))
        ## Adank:   lapply(groups, function(x) log(x) - sum(colMeans(log(x), ...)))
        logmeans <- unsplit(logmeans, group)
    }
    if (exp) logmeans <- exp(logmeans)
    logmeans
}

#' @export
normNearey2 <- function(f, group=NULL, exp=FALSE, ...) {
    normSharedLogmean(f, group=group, exp=exp, ...)
}

#' @export
normWattFabricius <- function(f, vowel, group=NULL) {
    if (ncol(f) != 2) {
        stop("Wrong dimensions: s-centroid normalization requires an Nx2 ",
             "matrix or data frame of F1 and F2 values.")
    }
    if (is.null(group)) group <- rep("g", nrow(f))
    ## make 2D list (group x vowel) of lists (f1,f2)
    subsets <- by(f, list(vowel, group), identity)
    means <- matrix(sapply(subsets,
                           function(i) ifelse(is.null(i),
                                              data.frame(I(c(f1=NA, f2=NA))),
                                              data.frame(I(colMeans(i))))),
                    nrow=nrow(subsets), dimnames=dimnames(subsets))
    minima <- do.call(rbind, apply(means, 2, function(i) data.frame(t(apply(do.call(rbind, i), 2, min, na.rm=TRUE)))))
    maxima <- do.call(rbind, apply(means, 2, function(i) data.frame(t(apply(do.call(rbind, i), 2, max, na.rm=TRUE)))))
    min.id <- apply(means, 2, function(i) apply(do.call(rbind, i), 2, which.min))
    max.id <- apply(means, 2, function(i) apply(do.call(rbind, i), 2, which.max))
    if (length(unique(min.id["f1",]))>1) {
        message("[phonR]: The vowel with the lowest mean F1 value (usually /i/)",
                "does not match across all speakers/groups. You'll ",
                "have to calculate s-centroid manually.")
        print(data.frame(minF1=minima["f1",],
                vowel=dimnames(means)[[1]][min.id["f1",]],
                group=dimnames(means)[[2]]))
        stop()
    } else if (length(unique(max.id["f1",]))>1) {
        message("[phonR]: The vowel with the highest mean F1 value (usually /a/) ",
                "does not match across all speakers/groups. You'll ",
                "have to calculate s-centroid manually.")
        print(data.frame(maxF1=round(maxima["f1",]),
                vowel=dimnames(means)[[1]][max.id["f1",]],
                group=dimnames(means)[[2]]))
        stop()
    }
    lowvowf2 <- do.call(rbind, means[unique(max.id["f1",]),])[,"f2"]
    centroids <- t(rbind(f1=(2*minima$f1 + maxima$f1) / 3,
                         f2=(minima$f2 + maxima$f2 + lowvowf2) / 3))
    f / centroids[group,]
}


## ## ## ## ## ## ## ## ## ## ## ##
## SECONDARY ANALYSIS FUNCTIONS  ##
## ## ## ## ## ## ## ## ## ## ## ##

#' @export
vowelMeansPolygonArea <- function(f1, f2, vowel, poly.order, group=NULL) {
    if (length(poly.order) != length(unique(poly.order))) {
        message("[phonR]: Duplicate entries in 'poly.order' detected; they ",
                "will be ignored.")
    }
    poly.order <- unique(as.character(poly.order)) # as.character in case factor
    v <- unique(as.character(vowel))
    if (length(setdiff(poly.order, v)) > 0) {
        message("[phonR]: There are vowels in 'poly.order' that are not in ",
                "'vowel'; they will be ignored.")
        poly.order <- intersect(poly.order, v)
    }
    if (is.null(group))  group <- "all.points"
    df <- data.frame(f2=f2, f1=f1, v=factor(vowel, levels=poly.order), g=group)
    df <- df[order(df$v),]
    bygrouparea <- c(by(df, df$g,
                        function(i) splancs::areapl(cbind(tapply(i$f2, i$v, mean),
                                                          tapply(i$f1, i$v, mean)))))
}

#' @export
convexHullArea <- function(f1, f2, group=NULL) {
    if (is.null(group))  group <- "all.points"
    df <- data.frame(x=f2, y=f1, g=group, stringsAsFactors=FALSE)
    bygrouppts <- by(df, df$g, function(i) i[chull(i$x, i$y), c("x", "y")])
    bygrouparea <- sapply(bygrouppts, function(i) {
        splancs::areapl(as.matrix(data.frame(x=i$x, y=i$y, stringsAsFactors=FALSE)))
    })
}

#' @export
repulsiveForce <- function(x, y, type, xform=log, exclude.inf=TRUE) {
    dmat <- as.matrix(dist(cbind(x, y)))
    if (exclude.inf) dmat[dmat == 0] <- min(dmat[dmat>0]) / 2
    force <- sapply(seq_along(type), function(i) {
        sum(1 / dmat[i, !(type %in% type[i])] ^ 2)
    })
    if (!is.null(xform)) force <- xform(force)
    force
}

#' @export
repulsiveForceHeatmap <- function(x, y, type=NULL, xform=log, exclude.inf=TRUE,
                                  resolution=10, colormap=NULL, fast=FALSE, ...) {
    ## default to grayscale
    if (is.null(colormap)) colormap <- plotrix::color.scale(x=0:100, cs1=0, cs2=0,
                                                            cs3=c(25,100),
                                                            color.spec="hcl")
    ## create grid encompassing vowel space
    gridlist <- createGrid(x, y, resolution)
    xgrid <- gridlist$x
    ygrid <- gridlist$y
    grid <- gridlist$g
    grid$z <- NA
    grid$v <- NA
    ## if using fast method, pre-calculate force of vowel tokens
    force <- repulsiveForce(x, y, type, xform, exclude.inf)
    if (fast) vertices <- data.frame(x=x, y=y, z=force)[!is.na(x) & !is.na(y),]
    else      vertices <- data.frame(x=x, y=y, z=type)[!is.na(x) & !is.na(y),]
    if (fast) {
        ## segregate grid points based on Delaunay triangulation of vowel space
        triang.obj <- with(vertices, deldir::deldir(x, y, z=z, suppressMsge=TRUE))
        triangs <- deldir::triang.list(triang.obj)
        sg.idxs <- lapply(triangs, function(i) splancs::inpip(grid, i[c("x", "y")],
                                                              bound=TRUE))
        subgrids <- lapply(sg.idxs, function(i) grid[i,])
        ## ignore triangles too small to contain any grid points
        vacant <- sapply(subgrids, function(i) dim(i)[1] == 0)
        triangs <- triangs[!vacant]
        subgrids <- subgrids[!vacant]
        sg.idxs <- sg.idxs[!vacant]
        ## assign a force value to each grid point
        sg.force <- mapply(function(i, j) fillTriangle(i[,"x"], i[,"y"],
                                                       j[, c("x", "y", "z")]),
                           subgrids, triangs, SIMPLIFY=FALSE)
        grid[do.call(c, sg.idxs), "z"] <- do.call(c, sg.force)
    } else {
        if (is.null(type)) stop("More accurate force calculation method (fast=FALSE) ",
                                "requires non-null values for 'type'.")
        hull <- vertices[chull(vertices),]  # polygon of hull
        sg.idxs <- splancs::inpip(grid[,1:2], hull)
        subgrid <- grid[sg.idxs,]
        ## which vowel is closest to each grid point?
        dmat <- apply(subgrid, 1, function(i) apply(as.matrix(vertices[c("x", "y")]),
                                                    1, function(j) dist(rbind(i, j))))
        ## still need to exclude inf, in case duplicate pts have diff vowels
        if (exclude.inf) dmat[dmat == 0] <- min(dmat[dmat>0]) / 2
        ## if there is a tie of which vowel is closest, pick first one (arbitrarily)
        which.min <- apply(dmat, 2, function(i) which(i == min(i))[1])
        #how.many.min <- apply(dmat, 2, function(i) length(which(i == min(i))))
        ## TODO: sensible way to do tiebreaker? wherever how.many.min > 1, look
        ## at next closest vowel iteratively until you find one that matches one
        ## of the vowels tied as closest... could get ugly.
        sg.vowel <- vertices$z[which.min]
        sg.force <- sapply(seq_along(sg.vowel),
                           function(i) sum(1/dmat[!(vertices$z %in% sg.vowel[i]), i] ^ 2))
        if (!is.null(xform)) sg.force <- xform(sg.force)
        grid[sg.idxs, "z"] <- sg.force
    }
    image(xgrid, ygrid, matrix(grid$z, nrow=length(xgrid)),
          col=colormap, ...)
}

#' @export
repulsiveForceHeatmapLegend <- function (x, y, labels=c("low", "high"), pos=c(1, 3),
                                         colormap=NULL, smoothness=50, lend=2,
                                         lwd=12, ...) {
    if (is.null(colormap)) {  # default to grayscale
        colormap <- plotrix::color.scale(x=0:100, cs1=0, cs2=0, cs3=c(0,100),
                                         alpha=1, color.spec="hcl")
    }
    xvals <- seq(x[1], x[2], length.out=smoothness)
    yvals <- seq(y[1], y[2], length.out=smoothness)
    cols <- colormap[round(seq(1, length(colormap), length.out=smoothness))]
    invisible(plotrix::color.scale.lines(xvals, yvals, col=cols, lend=lend,
                                         lwd=lwd, ...))
    if (!is.null(labels)) {
        text(x, y, labels=labels, pos=pos, xpd=TRUE)
    }
}


## ## ## ## ## ## ## ## ## ## ## ## ##
## UTILITY FUNCTIONS: NOT EXPORTED  ##
## ## ## ## ## ## ## ## ## ## ## ## ##
prettyTicks <- function(lim) {
    axrange <- abs(diff(lim))
    step <- 10^(floor(log(axrange,10)))
    coef <- ifelse(axrange/step < 1, 0.1,
                   ifelse(axrange/step < 2, 0.2,
                          ifelse(axrange/step < 5, 0.5, 1)))
    step <- step*coef
    lims <- c(ceiling(max(lim)/step)*step, floor(min(lim)/step)*step)
    if (diff(lims) < 0) {step <- -step}
    seq(lims[1],lims[2],step)
}

ellipse <- function(mu, sigma, n, alpha=0.05, npoints=250, draw=TRUE, ...) {
    if (all(sigma == matrix(rep(0, 4), nrow=2))) return(rbind(mu, mu))
    else if (n < 3) return(rbind(mu, mu))
    p <- length(mu)
    es <- eigen(sigma)
    e1 <- es$vectors %*% diag(sqrt(es$values))
    theta <- seq(from=0, to=2 * pi, length.out=npoints)
    unit.circle <- cbind(cos(theta), sin(theta))
    ## for small n, confidence ellipses for multivariate sample means are
    ## distributed as Hotelling's T^2, which is (with appropriate scale factor)
    ## equivalent to an F distribution (hence the qf() function). For large n,
    ## this asymptotically approaches qchisq(1-alpha, p)
    scale.factor <- p * (n - 1) / (n - p)
    critical.radius <- sqrt(scale.factor * qf(1 - alpha, df1=p, df2=n - p))
    ## if we needed it, this would be the t-squared statistic
    # tsquared <- n * t(mu) %*% solve(sigma) %*% mu
    pts <- t(mu - (e1 %*% t(critical.radius * unit.circle)))
    if (draw) {
        colnames(pts) <- c("x", "y")
        polygon(pts, ...)
    }
    invisible(pts)
}

makeTransparent <- function (color, opacity) {
    rgba <- t(col2rgb(color, alpha=TRUE))
    rgba[,4] <- round(255 * opacity)
    colnames(rgba) <- c("red", "green", "blue", "alpha")
    trans.color <- do.call(rgb, c(as.data.frame(rgba), maxColorValue=255))
}

uniquify <- function(x, default.val) {
    ## handle factors smartly
    y <- suppressWarnings(as.numeric(as.character(x)))
    if (any(is.na(y))) x <- as.character(x)
    else x <- y
    ux <- unique(x)
    ifelse(length(ux) == 1, ux, default.val)
}

fillTriangle <- function(x, y, vertices) {
    ## pineda's triangle filling algorithm
    x0 <- vertices[1,1]
    x1 <- vertices[2,1]
    x2 <- vertices[3,1]
    y0 <- vertices[1,2]
    y1 <- vertices[2,2]
    y2 <- vertices[3,2]
    z0 <- vertices[1,3]
    z1 <- vertices[2,3]
    z2 <- vertices[3,3]
    e0xy <- (x-x0)*(y1-y0)-(y-y0)*(x1-x0)
    e1xy <- (x-x1)*(y2-y1)-(y-y1)*(x2-x1)
    e2xy <- (x-x2)*(y0-y2)-(y-y2)*(x0-x2)
    e0x2 <- (x2-x0)*(y1-y0)-(y2-y0)*(x1-x0)
    e1x0 <- (x0-x1)*(y2-y1)-(y0-y1)*(x2-x1)
    e2x1 <- (x1-x2)*(y0-y2)-(y1-y2)*(x0-x2)
    f0 <- e0xy / e0x2
    f1 <- e1xy / e1x0
    f2 <- e2xy / e2x1
    z <- f0*z2 + f1*z0 + f2*z1
}

createGrid <- function(x, y, resolution) {
    ## create grid encompassing all pts (x, y) with `resolution` pts along
    ## short dimension
    vertices <- data.frame(x=x, y=y)
    vertices <- vertices[!is.na(vertices$x) & !is.na(vertices$y),]
    bounding.rect <- apply(vertices[c("x", "y")], 2, range, na.rm=TRUE)
    xr <- abs(diff(bounding.rect[,"x"]))
    yr <- abs(diff(bounding.rect[,"y"]))
    if (xr > yr) {
        xres <- round(resolution * xr / yr)
        yres <- resolution
    } else {
        xres <- resolution
        yres <- round(resolution * yr / xr)
    }
    xgrid <- seq(floor(bounding.rect[1, 1]), ceiling(bounding.rect[2, 1]),
                 length.out=xres)
    ygrid <- seq(floor(bounding.rect[1, 2]), ceiling(bounding.rect[2, 2]),
                 length.out=yres)
    grid <- expand.grid(x=xgrid, y=ygrid)
    list(g=grid, x=xgrid, y=ygrid)
}