turning R expressions into functions?

http://seehuhn.de/
-------------- next part --------------
# timeit.R - pairwise comparison for the execution time of R expressions
#
# Copyright (c) 2012  Jochen Voss  <voss at seehuhn.de>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# ----------------------------------------------------------------------
#
# This file provides the R command 'TimeIt' to compare the execution
# time of two R expressions.

FuncIt <- function(k, expr) {
  # Return a function which executes an expression k times.
  #
  # Args:
  #   k: The number of times 'expr' is executed.
  #   expr: An R expression.
  #
  # Returns:
  #   An R function, executing 'expr' in a loop.
  k <- as.numeric(k)
  expr <- eval.parent(substitute(expr))
  fn <- eval(substitute(function() { for (funcit.i in 1:k) { expr } }))
  return(fn)
}

TuneIt <- function(expr, max.seconds=1) {
  # Determine the approximate cost of calling an R expression in a
  # loop.  This function tries loops of different length and the uses
  # linear interpolation to get the result.
  #
  # Args:
  #   expr: The R expression to test.
  #   max.seconds: How much time (approximately) to use for
  #                measurements, in seconds.  This should be much
  #                larger than the resolution of 'system.time'.
  #                Default is 1.
  #
  # Returns:
  #   A vector 'x' of length 2, such that the execution time for 'k'
  #   iterations is approximately 'x[0] + k * x[1]'.
  kk <- c()
  tt <- c()
  k <- 1
  repeat {
    f <- FuncIt(k, expr)
    t <- system.time(f())[1]
    kk <- c(kk, k)
    tt <- c(tt, t)
    if (t > max.seconds / 3) break
    k <- 2 * k
  }
  if (k > 1) {
    fit <- lm(tt ~ kk)
    return(coefficients(fit))
  } else {
    return(c(0, tt))
  }
  return(k)
}

TimeIt <- function(ex1, ex2, total.time=30, verbose=T) {
  # Compare the execution time of two R expressions.
  #
  # Args:
  #   ex1: The first R expression to evaluate
  #   ex2: The second R expression to evaluate
  #   total.time: How much time (approximately) to spend on
  #     measuring, in seconds.  Longer times lead to more
  #     accurate measurements and allow to detect smaller
  #     differences in run time.  Default is 30 seconds.
  #
  # Returns:
  #   An object of class 'TimeIt', summarising the difference in
  #   execution time of the two expressions.
  start <- proc.time()[1]

  ex1 <- substitute(ex1)
  ex2 <- substitute(ex2)

  if (verbose) {
    cat("tuning ...\n")
  }
  # Use at most 20% or 10 seconds (whatever is smaller) of our time
  # budget for tuning.
  tune.time <- min(.2*total.time, 10)
  c1 <- TuneIt(ex1, tune.time / 2)
  c2 <- TuneIt(ex2, tune.time / 2)
  mid <- proc.time()[1] - start
  total.time <- total.time - mid

  block.min <- 1
  block.target <- total.time^(1/4) * block.min^(3/4)
  c <- c1 + c2
  block.k <- max(round((block.target - c[1]) / c[2]), 1)
  f1 <- FuncIt(block.k, ex1)
  f2 <- FuncIt(block.k, ex2)
  ex1.time <- c1[1] + block.k * c1[2]
  ex2.time <- c2[1] + block.k * c2[2]
  pair.time <- ex1.time + ex2.time
  n <- max(round(total.time / pair.time), 2)

  if (verbose) {
    cat("measuring ", n, "*", block.k,
        " samples for each expression ...\n", sep="")
    flush.console()
    progress <- txtProgressBar(min=0, max=mid+n*pair.time, style=3)
  }
  X1 <- c()
  X2 <- c()
  for (i in 1:n) {
    if (verbose) {
      setTxtProgressBar(progress, mid + (i-1)*pair.time)
    }
    X1 <- c(X1, system.time(f1())[1])
    if (verbose) {
      setTxtProgressBar(progress, mid + (i-1)*pair.time + ex1.time)
    }
    X2 <- c(X2, system.time(f2())[1])
  }
  X1 <- X1 / block.k
  X2 <- X2 / block.k
  if (verbose) {
    setTxtProgressBar(progress, mid + n*pair.time)
    close(progress)
    cat("\n")
  }

  name1 <- deparse(ex1)
  name2 <- deparse(ex2)
  l <- max(nchar(name1), nchar(name2))
  pad1 <- rep(" ", l - nchar(name1))
  pad2 <- rep(" ", l - nchar(name2))

  test <- t.test(X1, X2, paired=T)

  res <- list(name=c(name1,name2),
              name.padded=c(paste(name1, pad1, sep=" "),
                            paste(name2, pad2, sep=" ")),
              mean=c(mean(X1),mean(X2)),
              sd=c(sd(X1),sd(X2)),
              ci=test$conf.int,
              p.value=test$p.value)
  class(res) <- "TimeIt"
  return(res)
}

print.TimeIt <- function(x, ...) {
  if (min(x$mean[1], x$mean[2]) >= 0.1) {
    q = 1
    unit = "seconds/call"
  } else {
    q = 1000
    unit = "ms/call"
  }

  old <- options(digits=3)
  cat("execution time comparison:\n")
  cat(x$name.padded[1], "  (", x$mean[1]*q, " ? ", x$sd[1]*q, ") ",
      unit, "\n", sep="")
  cat(x$name.padded[2], "  (", x$mean[2]*q, " ? ", x$sd[2]*q, ") ",
      unit, "\n", sep="")

  cat("CI for difference: [", x$ci[1]*q, ", ", x$ci[2]*q, "] ",
      unit, "\n", sep="")
  cat("\n")
  if (x$ci[1] > 0) {
    cat("'", x$name[2], "' is about ", 100*(x$mean[1]-x$mean[2])/x$mean[1],
        "% faster (p=", x$p.value, ")\n", sep="")
  } else if (x$ci[2] < 0) {
    cat("'", x$name[1], "' is about ", 100*(x$mean[2]-x$mean[1])/x$mean[2],
        "% faster (p=", x$p.value, ")\n", sep="")
  } else {
    cat("difference is less than ",
        50 * (x$ci[2]-x$ci[1]) / mean(x$mean[1], x$mean[2]),
        "% (p=", x$p.value, ")\n", sep="")
  }
  cat("\n")
  options(old)
}