ra.png

Code For Ra Timing Tests

To Ra homepage


# time-jit.R

library(jit)

GCTORTURE.FLAG <- 0
TRACE.FLAG <- 0
QUICK.FLAG <- FALSE   # FALSE for full test, TRUE for quick test
NREPEATS <- 5         # for calculating time stddev
if (QUICK.FLAG == 1)
    NREPEATS <- 1

JIT.FLAG <- 1
source("time-jit-lib.R")

if (is.ra) {
    JIT.FLAG <- 2
    source("time-jit-lib.R")
}

# time-jit-lib.R

jit.flag <- 0
gctorture(GCTORTURE.FLAG)

test <- function(f, N)
{
    # na.rm=TRUE is useful when building the framework with very short times
    percent.relative.sd <- function(x) 100 * sd(x, na.rm=TRUE) / mean(x)

    cat(sprintf("%-22.22s %9.0f   ", paste(substitute(f)), N))

    time.jit <- time.no.jit <- double(NREPEATS)

    N <- as.integer(N)  # use an integer index in for loops

    for (i in 1:NREPEATS) {
        if (i %/% 2) { # alternate test order, probably unnecessary
            if (is.ra) {
                jit.flag <<- JIT.FLAG
                time.jit[i] <- system.time(no.jit.result <- f(N))[3]
            }
            jit.flag <<- 0
            time.no.jit[i]  <- system.time(jit.result    <- f(N))[3]
        } else {
            jit.flag <<- 0
            time.no.jit[i]  <- system.time(jit.result    <- f(N))[3]
            if (is.ra) {
                jit.flag <<- JIT.FLAG
                time.jit[i] <- system.time(no.jit.result <- f(N))[3]
            }
        }
    }
    cat(sprintf("%6.2f %3.1f%%      ",
        mean(time.no.jit),
        percent.relative.sd(time.no.jit)))

    if (is.ra)
        cat(sprintf("%6.2f %3.1f%%     %5.3f %3.1f%%",
            mean(time.jit),
            percent.relative.sd(time.no.jit),
            mean(time.jit / time.no.jit),
            percent.relative.sd(time.jit / time.no.jit)))

    cat("\n")

    if (is.ra)
        stopifnot(identical(no.jit.result, jit.result))
}
convolve <- function(N) # from the extending R manual
{
    jit(jit.flag, TRACE.FLAG)
    a <- double(N)
    b <- double(N)
    na <- length(a)
    nb <- length(b)
    ab <- double(na + nb)
    for(i in 1:na)
        for(j in 1:nb)
             ab[i + j] <- ab[i + j] + a[i] * b[j]
    ab
}
# from base/TAOCP.R, a good test of integer arithmetic
.TAOCP1997init <- function(seed)
{
    seed <- as.integer(seed)  # added for jit to prevent type change error
    jit(jit.flag, TRACE.FLAG)
    KK <- 100L; LL <- 37L; MM <- as.integer(2^30)
    KKK <- KK + KK - 1L; KKL <- KK - LL
    ss <- seed - (seed %% 2L) + 2L
    X <- integer(KKK)
    for(j in 1L:KK) {
        X[j] <- ss
        ss <- ss+ss
        if(ss >= MM) ss <- ss - MM + 2L
    }
    X[2L] <- X[2L] + 1L
    ss <- seed
    T <- 69L
    while(T > 0) {
        for(j in KK:2L) X[j + j - 1L] <- X[j]
        for(j in seq(KKK, KKL + 1L, -2L))
            X[KKK - j + 2L] <- X[j] - (X[j] %% 2L)
        for(j in KKK:(KK+1L))
            if(X[j] %% 2L == 1L) {
                X[j - KKL] <- (X[j - KKL] - X[j]) %% MM
                X[j - KK] <- (X[j - KK] - X[j]) %% MM
            }
        if(ss %% 2L == 1L) {
            for(j in KK:1L) X[j + 1L] <- X[j]
            X[1L] <- X[KK + 1L]
            if(X[KK + 1L] %% 2L == 1L)
                X[LL + 1L] <- (X[LL + 1L] - X[KK + 1L]) %% MM
        }
        if(ss) ss <- ss %/% 2L else T <- T - 1L
    }
    rs <- c(X[(LL+1L):KK], X[1L:LL])
    invisible(rs)
}
`base/TAOCP.R` <- function(N)
{
    x <- 0
    for (i in 1:N)
        x <- c(x, .TAOCP1997init(i))
    x
}
# from ROCR package, calculate area under ROC curve
.performance.auc <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, fpr.stop) {

      jit(jit.flag, TRACE.FLAG)

      x <- fp / n.neg
      y <- tp / n.pos

      finite.bool <- is.finite(x) & is.finite(y)
      x <- x[ finite.bool ]
      y <- y[ finite.bool ]
      if (length(x) < 2) {
          stop(paste("Not enough distinct predictions to compute area",
                     "under the ROC curve."))
      }

      if (fpr.stop < 1) {
        ind <- max(which( x <= fpr.stop ))
        tpr.stop <- approxfun( x[ind:(ind+1)], y[ind:(ind+1)] )(fpr.stop)
        x <- c(x[1:ind], fpr.stop)
        y <- c(y[1:ind], tpr.stop)
      }

      ans <- list()
      auc <- 0
      for (i in 2:length(x)) {
          auc <- auc + 0.5 * (x[i] - x[i-1]) * (y[i] + y[i-1])
      }
      ans <- list( c(), auc)
      names(ans) <- c("x.values","y.values")
      return(ans)
  }
`ROCR/auc` <- function(N)
{
    .performance.auc(predictions=NULL, labels=NULL,
        cutoffs=NULL, fp=fp, tp=tp,
        fn=NULL, tn=NULL, n.pos=length(tp), n.neg=length(fp),
        n.pos.pred=NULL, n.neg.pred=NULL, fpr.stop=1)
}
# Distribution of determinant of 2x2 matrix
# From V&R S Programming p154

dd.for.jit <- function()
{
    jit(JIT.FLAG, TRACE.FLAG)
    val <- NULL
    nojit(val)                  # allow "c" below to change len of val
    for (a in 0:9)
      for (b in 0:9)
        for (d in 0:9)
          for (e in 0:9)
            val <- c(val, a*b - d*e)
    table(val)
}
dd.for.nojit <- function()
{
    val <- NULL
    for (a in 0:9)
      for (b in 0:9)
        for (d in 0:9)
          for (e in 0:9)
            val <- c(val, a*b - d*e)
    table(val)
}
dd.for <- function(N)
{
    if (jit.flag)
        for (i in 1:N)
            dd.for.jit()
    else
        for (i in 1:N)
            dd.for.nojit()
}
dd.for.prealloc.jit <- function()
{
    jit(JIT.FLAG, TRACE.FLAG)
    val <- double(10000)        # was val <- NULL
    nval <- 0
    for (a in 0:9)
      for (b in 0:9)
        for (d in 0:9)
          for (e in 0:9)
            val[nval <- nval + 1] <- a*b - d*e
                                # was val <- c(val, a*b - d*e)
    table(val)
}
dd.for.prealloc.nojit <- function()
{
    val <- double(10000)        # was val <- NULL
    nval <- 0
    for (a in 0:9)
      for (b in 0:9)
        for (d in 0:9)
          for (e in 0:9)
            val[nval <- nval + 1] <- a*b - d*e
                                # was val <- c(val, a*b - d*e)
    table(val)
}
dd.for.prealloc <- function(N)
{
    if (jit.flag)
        for (i in 1:N)
            dd.for.prealloc.jit()
    else
        for (i in 1:N)
            dd.for.prealloc.nojit()
}
dd.for.tabulate.jit <- function()
{
    jit(JIT.FLAG, TRACE.FLAG)
    val <- double(10000)        # was val <- NULL
    nval <- 0
    for (a in 0:9)
      for (b in 0:9)
        for (d in 0:9)
          for (e in 0:9)
            val[nval <- nval + 1] <- a*b - d*e
                                # was val <- c(val, a*b - d*e)
    tabulate(val)
}
dd.for.tabulate.nojit <- function()
{
    val <- double(10000)        # was val <- NULL
    nval <- 0
    for (a in 0:9)
      for (b in 0:9)
        for (d in 0:9)
          for (e in 0:9)
            val[nval <- nval + 1] <- a*b - d*e
                                # was val <- c(val, a*b - d*e)
    tabulate(val)
}
dd.for.tabulate <- function(N)
{
    if (jit.flag)
        for (i in 1:N)
            dd.for.tabulate.jit()
    else
        for (i in 1:N)
            dd.for.tabulate.nojit()
}
dd.fast.jit <- function()
{
    jit(JIT.FLAG, TRACE.FLAG)
    val <- outer(0:9, 0:9, "*")
    val <- outer(val, val, "-")
    table(val + 82)
}
dd.fast.nojit <- function()
{
    val <- outer(0:9, 0:9, "*")
    val <- outer(val, val, "-")
    table(val + 82)
}
dd.fast <- function(N)
{
    if (jit.flag)
        for (i in 1:N)
            dd.fast.jit()
    else
        for (i in 1:N)
            dd.fast.nojit()
}
dd.fast.tabulate.jit <- function()
{
    jit(JIT.FLAG, TRACE.FLAG)
    val <- outer(0:9, 0:9, "*")
    val <- outer(val, val, "-")
    tabulate(val)
}
dd.fast.tabulate.nojit <- function()
{
    val <- outer(0:9, 0:9, "*")
    val <- outer(val, val, "-")
    tabulate(val)
}
dd.fast.tabulate <- function(N)
{
    if (jit.flag)
        for (i in 1:N)
            dd.fast.tabulate.jit()
    else
        for (i in 1:N)
            dd.fast.tabulate.nojit()
}

looped.dnorm <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    mu <- 0
    sigma <- 1
    x <- 0
    for (i in 1:N) # from one of Luke's compiler documents
        x <- x + (1/sqrt(2 * pi)) * exp(-0.5 * ((x - mu)/sigma)^2) / sigma
    x
}

`while  x <- x + 1` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    x <- 0
    while (x < N)
        x <- x+1
    x
}
`while  x <- x + 1i` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    i <- 0i
    Ni <- N + 0i
    while (i != Ni)
        i <- i + 1
    i
}
`repeat x <- x + 1` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    x <- 0
    repeat {
        x <- x+1
        if (x >= N)
            break;
    }
    x
}
`repeat x <- x + 1i` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    x <- 0i
    repeat {
        x <- x+1
        if (Re(x) >= N)
            break;
    }
    x
}

`for.if` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    iA <- seq(2,N); x <- double(N)
    for (i in iA) {
        if (i %% 2)
            x <- x + 1
        else
            x <- x + 100
    }
    x
}

# Tests from Vadim Ogranovich post.
# See http://tolstoy.newcastle.edu.au/R/devel/05/04/0678.html
# Expressions are the same as Luke's email reply except
# that x and iA are local.

`vadim1 1` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    iA <- seq(2,N); x <- double(N)
    for (i in iA)
        1
    x
}
`vadim2 i` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    iA <- seq(2,N);
    for (i in iA)
        i
    i
}
`vadim3 i-1` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    iA <- seq(2,N);
    for (i in iA)
        i-1
    i
}
`vadim4 x[i-1]` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    iA <- seq(2,N); x <- double(N); x[1] <- 1; x[2] <- 2
    for (i in iA)
        x[i-1]
    x
}
`add1   x <- x + 1` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    x <- 0
    for(i in 1:N)
        x <- x+1
    x
}
`vadim5 x[i] <- 1.0` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    iA <- seq(2,N); x <- double(N); x[1] <- 1; x[2] <- 2
    for (i in iA)
        x[i] <- 1.0
    x
}
`vadim6 x[i] <- x[i-1]` <- function(N)
{
    jit(jit.flag, TRACE.FLAG)
    iA <- seq(2,N); x <- double(N); x[1] <- 1; x[2] <- 2
    for (i in iA)
        x[i] <- x[i-1]
    x
}
# Loop counts are chosen so the jitted times are greater
# than about a second (with Ra version 5 on a 3G Pentium).
# This is necessary for plausible timing results.
# We adjust N below so each jitted test takes roughly the same time.

N <- if (QUICK.FLAG) 2e6 else 2e7

cat("is.ra", is.ra, "NREPEATS", NREPEATS, "QUICK.FLAG", QUICK.FLAG, "JIT.FLAG", JIT.FLAG, "\n\n")
cat("testname                       N     time rsd     jit-time rsd    reltime rsd\n\n")

test(convolve,       if (QUICK.FLAG) 500 else 1600)
test(`base/TAOCP.R`, if (QUICK.FLAG) 20 else 80)
test(looped.dnorm,   if (QUICK.FLAG) 1e5 else 8e5)

set.seed(1)   # for reproducibility
fp = c(0, cumsum(runif(N / 10) > .5)) # cumulative false positives for `ROCR/auc`
tp = c(0, cumsum(runif(N / 10) > .5)) # cumulative true positives
test(`ROCR/auc`, N / 10)

cat("\n")
test(dd.for,           if (QUICK.FLAG) 10 else N / 30000)
test(dd.for.prealloc,  if (QUICK.FLAG) 10 else N / 30000)
test(dd.for.tabulate,  if (QUICK.FLAG) 10 else N / 30000)
test(dd.fast,          if (QUICK.FLAG) 10 else N / 30000)
test(dd.fast.tabulate, if (QUICK.FLAG) 10 else N / 30000)
cat("\n")
test(`while  x <- x + 1`, N / 5)
test(`repeat x <- x + 1`, N / 5)
test(`for.if`, N / 500)
test(`while  x <- x + 1i`, N / 5)
test(`repeat x <- x + 1i`, N / 5)
cat("\n")
test(`vadim1 1`, N)
test(`vadim2 i`, N)
test(`vadim3 i-1`, N)
test(`vadim4 x[i-1]`, N)
test(`add1   x <- x + 1`, N)
test(`vadim5 x[i] <- 1.0`, N)
test(`vadim6 x[i] <- x[i-1]`, N)
cat("\n")

To Ra homepage