Add the SpeedTest suite

Benchmarks/prg/Qlearn.R

@@ -0,0 +1,78 @@
+execute <- function(size = 30000L) {
+  simulate <- function(init, world, gamma, alpha, epsilon, steps) {
+    history = matrix(NA,steps,6)
+    colnames(history) = c("t","s","a","r","rs","sn")
+
+    Q = matrix(0,n.states,n.actions)
+
+    s = init()
+
+    for (t in 1:steps) { 
+      if (runif(1)<epsilon) {a = sample(n.actions,1)}
+      else {a = order(Q[s,])[n.actions]}
+
+      w = world(s,a)
+      r = w$r
+      sn = w$s
+
+      Q[s,a] = (1-alpha) * Q[s,a] + 
+        alpha * (r + gamma * (epsilon*mean(Q[sn,]) + (1-epsilon)*max(Q[sn,])))
+
+      history[t,"t"] = t
+      history[t,"s"] = s
+      history[t,"a"] = a
+      history[t,"r"] = r
+      history[t,"sn"] = sn
+
+      s = sn
+    }
+
+    history[,"rs"] = history[,"r"]
+    for (t in 2:steps){history[t,"rs"] = 0.1*history[t,"rs"] + 0.9*history[t-1,"rs"]}
+
+    list (history=history, Q=Q)
+  }
+
+  init2m = function () { 
+    marks <<- rep(0,n.states)
+    sample(n.states,1)
+  }
+
+  world2m = function (s, a) {
+    bit = as.numeric(s>n.states/2)
+    s = s - bit*n.states/2
+
+    if (runif(1)<0.05)
+    { s = sample(n.states/2,1)
+    }
+    else
+    { if (a>3)
+      { bit = 1-bit
+        a = a-3
+      }
+      s = s + (a-2)
+      if (s<1) s = n.states/2
+      if (s>n.states/2) s = 1
+    }
+
+    r = marks[s] - 10*as.numeric(s==1)
+
+    marks <<- as.numeric (marks>0 | (runif(n.states)<0.3))
+    marks[s] <<- 0
+
+    s = s + bit*n.states/2
+
+    list (s=s, r=r)
+  }
+
+  gamma = 0.95
+  alpha = 0.015
+  epsilon = 0.1
+
+  n.states = 10*2
+  n.actions = 3*2
+  set.seed(1)
+
+  result2m = simulate (init2m, world2m, gamma, alpha, epsilon, size)
+  return (mean(result2m$history[,"r"])) # Average reward
+}

Benchmarks/prg/cholesky.R

@@ -0,0 +1,36 @@
+execute <- function(size = 200L) {
+  cholesky <- function(A) {
+    if (!is.matrix(A) || nrow(A) != ncol(A)) {
+      stop("The argument for cholesky must be a square matrix")
+    }
+
+    p <- nrow(A)
+    U <- matrix(0, p, p)
+
+    for (i in 1:p) {
+      if (i == 1) {U[i,i] <- sqrt(A[i,i])} 
+      else {U[i,i] <- sqrt(A[i,i] - sum(U[1:(i-1),i]^2))}
+
+      if (i<p) {
+        for (j in (i+1):p) {
+          if (i == 1) {U[i, j] <- A[i, j] / U[i, i]} 
+          else {U[i, j] <- (A[i, j] - sum(U[1:(i-1), i] * U[1:(i-1), j])) / U[i, i]}
+        }
+      }
+    }
+    U
+  }
+
+  set.seed(1)
+  A <- matrix(rnorm(size^2), size, size)
+  M <- t(A) %*% A
+
+  U <- vector("list", 10)
+  for (i in 1:10) U[[i]] <- cholesky(M)
+
+  V <- vector("list", 1000)
+  for (i in 1:1000) V[[i]] <- chol(M)
+
+  res <- sum((U[[1]]-V[[1]])^2)
+  return (res)
+}

Benchmarks/prg/cv-basisfun.R

@@ -0,0 +1,58 @@
+execute <- function(reps = 10L) {
+  Phi <- function(xv, s) {
+    m <- seq(-2*s, 1+2*s, by=s)
+    Phi <- matrix(1, length(xv), length(m)+1)
+    for (j in 1:length(m)) {
+      Phi[,1+j] <- exp(-0.5*(xv-m[j])^2/s^2)
+    }
+    Phi
+  }
+
+  penalized.least.squares <- function(Phi, tv, lambda=0) {
+    S <- diag(ncol(Phi)) 
+    S[1,1] <- 0
+    as.vector(solve(lambda*S + t(Phi)%*%Phi) %*% t(Phi) %*% tv)
+  }
+
+  predictions <- function(Phi, w) {
+    as.vector(Phi %*% w)
+  }
+
+  validation.error <- function(xv, tv, s, lambda, vix) { 
+    w <- penalized.least.squares(Phi(xv[-vix],s), tv[-vix], lambda)
+    p <- predictions(Phi(xv[vix],s), w)
+    mean((tv[vix]-p)^2)
+  }
+
+  val.array <- function(xv, tv, try.s, try.lambda, vix) {
+    V <- matrix(NA, length(try.s), length(try.lambda))
+    rownames(V) <- paste("s=",try.s,sep="")
+    colnames(V) <- paste("lambda=",try.lambda,sep="")
+
+    for (i in 1:length(try.s)) {
+      for (j in 1:length(try.lambda)) {
+        V[i,j] <- validation.error(xv, tv, try.s[i], try.lambda[j], vix)
+      }
+    }
+    V
+  }
+
+  cross.val.array <- function(xv, tv, try.s, try.lambda, S) {
+    n <- length(tv)
+    for (r in 1:S) {
+      vix <- floor(1 + n*(r-1)/S) : floor(n*r/S)
+      V <- val.array(xv, tv, try.s, try.lambda, vix)
+      V.sum <- if (r==1) V else V.sum + V
+    }
+    V.sum / S
+  }
+
+  set.seed(2)
+  xt <- runif(50)
+  tt <- sin(1+xt^2) + rnorm(50, 0, 0.03)
+  try.lambda <- c(0.001,0.01,0.1,1,10)
+  try.s <- c(0.02,0.1,0.5,2.5)
+
+  for (i in 1:reps) r <- round(sqrt(cross.val.array(xt,tt,try.s,try.lambda,5)),4)
+  return(r)
+}

Benchmarks/prg/data.R

@@ -0,0 +1,37 @@
+execute <- function(size = 10L) {
+  sim <- function() {
+    N <- size*1000; out <- list()
+    df <- data.frame(matrix(NA, N, size))
+    colnames(df) <- paste("VAR", 1:size, sep="")
+    for (i in 1:size) {
+      if (i %% 4 == 0) df[,i] <- rnorm(N, 0.1*(i+(1:N)), 25)
+      else if (i %% 4 == 1) df[,i] <- runif(N) < 0.7
+      else if (i %% 4 == 2) df[,i] <- paste("xyz", rpois(N,3), sep="")
+      else df[,i] <- 0.01*(1:N) + rpois(N,3)
+    }
+    df
+  }
+
+  analyse <- function(df,pr=FALSE) {
+    s <- summary(df)
+    if (pr) c(out, s)
+    df2 <- df[,!sapply(df,is.character)]
+    m <- lm (VAR4 ~ ., data = df2)
+    if (pr) c(out, m)
+    sm <- summary(m)
+    if (pr) c(out, sm)
+    coef(m)
+  }
+
+  set.seed(1)
+  out <- list()
+  test_data_file <- tempfile()
+  for (i in 1:50) data <- sim()
+  for (i in 1:25) write.table(data,test_data_file,quote=FALSE)
+  for (i in 1:50) rdata <- read.table(test_data_file,head=TRUE,stringsAsFactors=FALSE)
+  for (i in 1:150) res <- analyse(rdata)
+  analyse(rdata,TRUE)
+
+  unlink(test_data_file)
+  list(out, res)
+}

Benchmarks/prg/em.R

@@ -0,0 +1,24 @@
+execute <- function(iter = 10L) {
+  EM.censored.poisson <- function(n, m, c, lambda0 = (n * m + c) / (n + c), iterations) {
+    log.likelihood <- function(lambda) 
+    {n * m * log(lambda) - (n + c) * lambda + c * log(1 + lambda)}
+
+    lambda <- lambda0
+    old.ll <- log.likelihood(lambda)
+
+    for (i in 1:iter) {
+      p1 <- lambda / (1+lambda)
+      lambda <- (n*m + c*p1) / (n+c)
+      new.ll <- log.likelihood(lambda)
+
+      if (new.ll - old.ll < -1e-6) {stop("Log likelihood decreased!")}
+      old.ll <- new.ll
+    }
+    lambda
+  }
+
+  for (i in 1:3) r1 <- EM.censored.poisson(5,6.1,20,(iter*1000))
+  for (i in 1:30000) r2 <- EM.censored.poisson(5,6.1,20,iter)
+
+  list(r1,r2)
+}

Benchmarks/prg/gcd.R

@@ -0,0 +1,47 @@
+execute <- function(size = 500L) {
+  gcd1 <- function(a, b) if (a > b) gcd1(a-b,b) else if (b > a) gcd1(a,b-a) else a
+
+  gcd2 <- function(a, b) {
+    if (a == 0) b 
+    else if (b == 0) a 
+    else if (a > b) gcd2(a%%b,b) 
+    else gcd2(a,b%%a) 
+  }
+
+  gcd3 <- function(a, b) {
+    repeat { 
+      if (a > b) { a <- a-b; next }
+      if (b > a) { b <- b-a; next }
+      return (a);
+    }
+  }
+
+  gcd4 <- function(a, b) {
+    repeat {
+      if (a == 0) return (b)
+      if (b == 0) return (a)
+      if (a > b) a <- a%%b
+      else b <- b%%a
+    }
+  }
+
+  gcd5 <- function(a, b) {
+    while (a != b) {
+      if (a > b) a <- a - b
+      else b <- b - a
+    }
+    a
+  }
+
+  gcd_table <- function(n, gcd) {
+    tbl <- matrix(integer(), n, n)
+    for (i in 1:n)
+      for (j in 1:n)
+        tbl[i, j] <- gcd(i, j)
+    tbl
+  }
+  r1 <- gcd_table(size,gcd1); r2 <- gcd_table(size,gcd2); r3 <- gcd_table(size,gcd3)
+  r4 <- gcd_table(size,gcd4); r5 <- gcd_table(size,gcd5)
+  res <- all(sapply(list(r2, r3, r4, r5), identical, r1))
+  return(res)
+}

Benchmarks/prg/gp.R

@@ -0,0 +1,79 @@
+execute <- function(size = 1L) {
+  gp_cov <- function(x1, x2, s, r) {
+    n <- length(x1)
+    m <- length(x2)
+    C <- matrix(nrow=n,ncol=m)
+
+    for (j in seq_len(m)) {C[,j] <- s^2 * exp(-(r*(x1-x2[j]))^2)}
+    C
+  }
+
+  gp_log_likelihood <- function(x, y, s, r, v) {
+    C <- gp_cov(x, x, s, r)
+    diag(C) <- diag(C) + v^2
+    U <- chol(C)
+    u <- backsolve(U, y, transpose = TRUE)
+    as.vector(-sum(u^2)/2 - sum(log(diag(U))) - length(y)*log(2*pi)/2)
+  }
+
+  gp_search <- function(x, y, s.vec, r.vec, v.vec) {
+    best.ll <- -Inf
+    for (s in s.vec) {
+      for (r in r.vec) {
+        for (v in v.vec) {
+          ll <- gp_log_likelihood(x, y, s, r, v)
+          if (ll >= best.ll) {
+            best.s <- s
+            best.r <- r
+            best.v <- v
+            best.ll <- ll
+          }
+        }
+      }
+    }
+    c(s = best.s, r = best.r, v = best.v)
+  }
+
+  gp_predict <- function(x, y, s, r, v, x.test) {
+    C <- gp_cov(x,x,s,r)
+    diag(C) <- diag(C) + v^2
+    U <- chol(C)
+    u <- backsolve(U,y,transpose=TRUE)
+    u <- backsolve(U,u)
+    K <- gp_cov(x,x.test,s,r)
+    as.vector (t(K) %*% u)
+  }
+
+  tf <- function(x) 1.2 * sin(0.3 + 0.2 * x^2 + 2.7 * sin(2 * x + 0.2))
+
+  f <- function (N, n) { 
+    set.seed(1)
+    hlist <<- list()
+
+    for (i in 1:N) {
+      x <- rnorm(n)
+      x.test <- seq(-2,2,length=1000)
+      y <- tf(x) + rnorm(n,0,0.11)
+
+      h <- gp_search (x, y, s.vec = c(0.7,1.0,1.2,1.4,2.0), 
+                            r.vec = c(1.0,1.4,2.0,2.8), 
+                            v.vec = c(0.1,0.14,0.2,0.28))
+      hlist[[i]] <<- h
+
+      p <- gp_predict (x, y, h["s"], h["r"], h["v"], x.test)
+    }
+
+    cbind (x=x.test, y=tf(x.test), p=p)
+  }
+
+  R1 <- f((size*10), 100)
+  res <- hlist
+  R2 <- f(size, 350)
+
+  list(res, R1[c(1,500,1000),"p"], R2[c(1,500,1000),"p"])
+
+  # Original benchmark result checking
+  #
+  # plot(R[,"x"],R[,"y"],type="l")
+  # points(R[,"x"],R[,"p"],pch=20)
+}