# Splus Benchmark 2 (8 March 2003)
# version 2, scaled to get 1 +/- 0.1 sec with R 1.6.2
# using the standard ATLAS library (Rblas.dll)
# on a Pentium IV 1.6 Ghz with 1 Gb Ram on Win XP pro

# modified for Splus 6, 8 June 2002 (dropped library(Matrix) that is not provided any more)
# Author : Philippe Grosjean
# eMail  : phgrosjean@sciviews.org
# Web    : http://www.sciviews.org
# License: GPL 2 or above at your convenience (see: http://www.gnu.org)
#
# Several tests are adapted from the Splus Benchmark Test V. 2
# by Stephan Steinhaus (stst@informatik.uni-frankfurt.de)                              # 
# Reference for Escoufier's equivalents vectors (test III.5):
# Escoufier Y., 1970. Echantillonnage dans une population de variables
# aleatoires réelles. Publ. Inst. Statis. Univ. Paris 19 Fasc 4, 1-47.
#
# type source("c:/<dir>/Splus2.ssc") to start the test

runs <- 3			# Number of times the tests are executed
times <- rep(0, 15); dim(times) <- c(5,3)
options(object.size=100000000)

cat("   Splus Benchmark 2\n")
cat("   =================\n")
cat(c("Number of times each test is run__________________________: ", runs))
cat("\n\n")


cat("   I. Matrix calculation\n")
cat("   ---------------------\n")

# (1)
cumulate <- 0; a <- 0; b <- 0
for (i in 1:runs) {
  timing <- dos.time({
    a <- abs(matrix(rnorm(1500*1500, sd=0.1), ncol=1500, nrow=1500));
    b <- t(a);
    dim(b) <- c(750, 3000);
    a <- t(b)
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[1, 1] <- timing
cat(c("Creation, transp., deformation of a 1500x1500 matrix (sec): ", timing, "\n"))
remove(c("a", "b"))

# (2)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- abs(matrix(rnorm(800*800, sd=0.5), ncol=800, nrow=800));
  timing <- dos.time({ 
    b <- a^1000 
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[2, 1] <- timing
cat(c("800x800 normal distributed random matrix ^1000______ (sec): ", timing, "\n"))
remove(c("a", "b"))

# (3)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- rnorm(2000000)
  timing <- dos.time({
    b <- sort(a)
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[3, 1] <- timing
cat(c("Sorting of 2,000,000 random values__________________ (sec): ", timing, "\n"))
remove(c("a", "b"))

# (4)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- rnorm(700*700); dim(a) <- c(700, 700)
  timing <- dos.time({
    b <- crossprod(a, a)		# equivalent to: b <- t(a) %*% a
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[4, 1] <- timing
cat(c("700x700 cross-product matrix (b = a' * a)___________ (sec): ", timing, "\n"))
remove(c("a", "b"))

# (5)
cumulate <- 0; c <- 0; qra <-0
for (i in 1:runs) {
  a <- rnorm(600*600); dim(a) <- c(600,600)
  b <- 1:600
  timing <- dos.time({
    qra <- qr(a, tol = 1e-10);
    c <- qr.coef(qra, b)
    #Rem: a little faster than c <- lsfit(a, b)$coef
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[5, 1] <- timing
cat(c("Linear regression over a 600x600 matrix (c = a \\ b') (sec): ", timing, "\n"))
remove(c("a", "b", "c", "qra"))

times[ , 1] <- sort(times[ , 1])
cat("                      --------------------------------------------\n")
cat(c("                 Trimmed geom. mean (2 extremes eliminated): ", exp(mean(log(times[2:4, 1]))), "\n\n"))

cat("   II. Matrix functions\n")
cat("   --------------------\n")

# (1)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- rnorm(800000)
  timing <- dos.time({
    b <- fft(a)
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[1, 2] <- timing
cat(c("FFT over 800,000 random values______________________ (sec): ", timing, "\n"))
remove(c("a", "b"))

# (2)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- matrix(rnorm(320*320), nrow=320, ncol=320)
  timing <- dos.time({
    # b <- eigen.Matrix(a, vectors = F)$values not found any more in Splus 6 for Windows!
    b <- eigen.default(a, symmetric=F, only.values=T)$values
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[2, 2] <- timing
cat(c("Eigenvalues of a 320x320 random matrix______________ (sec): ", timing, "\n"))
remove(c("a", "b"))

# (3)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- matrix(rnorm(650*650), nrow=650, ncol=650)
  timing <- dos.time({
    b <- det(a)
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[3, 2] <- timing
cat(c("Determinant of a 650x650 random matrix______________ (sec): ", timing, "\n"))
remove(c("a", "b"))

# (4)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- rnorm(900*900); dim(a) <- c(900, 900)
  a <- crossprod(a, a)
  timing <- dos.time({
    b <- chol(a)
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[4, 2] <- timing
cat(c("Cholesky decomposition of a 900x900 matrix__________ (sec): ", timing, "\n"))
remove(c("a", "b"))

# (5)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  a <- matrix(rnorm(400*400), nrow=400, ncol=400)
  timing <- dos.time({
    # b <- solve.Matrix(a) is not available any more in Splus 6 for Windows
    b <- solve(a)
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[5, 2] <- timing
cat(c("Inverse of a 400x400 random matrix__________________ (sec): ", timing, "\n"))
remove(c("a", "b"))

times[ , 2] <- sort(times[ , 2])
cat("                      --------------------------------------------\n")
cat(c("                      Trimmed mean (2 extremes eliminated): ", mean(times[2:4, 2]), "\n\n"))

cat("   III. Programmation\n")
cat("   ------------------\n")

# (1)
cumulate <- 0; a <- 0; b <- 0; phi <- 1.6180339887498949
for (i in 1:runs) {
  a <- floor(runif(750000, max=1000))
  timing <- dos.time({
    b <- (phi^a - (-phi)^(-a))/sqrt(5)
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[1, 3] <- timing
cat(c("750,000 Fibonacci numbers calculation (vector calc)_ (sec): ", timing, "\n"))
remove(c("a", "b", "phi"))

# (2)
cumulate <- 0; b <- 0
for (i in 1:runs) {
	a <- 2250
  timing <- dos.time({
    a <- 1:a; b <- 1 / outer(a - 1, a, "+")
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[2, 3] <- timing
cat(c("Creation of a 2250x2250 Hilbert matrix (matrix calc) (sec): ", timing, "\n"))
remove(c("a", "b"))

# (3)
cumulate <- 0; c <- 0
gcd2 <- function(x, y) {if (sum(y > 1.0E-4) == 0) x else {y[y == 0] <- x[y == 0]; gcd2(y, x %% y)}}
for (i in 1:runs) {
  a <- ceiling(runif(70000, max=1000))
  b <- ceiling(runif(70000, max=1000))
  timing <- dos.time({	  
    c <- gcd2(a, b)                            # gcd2 is a recursive function
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[3, 3] <- timing
cat(c("Grand common divisors of 70,000 pairs (recursion)___ (sec): ", timing, "\n"))
remove(c("a", "b", "c", "gcd2"))

# (4)
cumulate <- 0; b <- 0
for (i in 1:runs) {
  b <- rep(0, 220*220); dim(b) <- c(220, 220)
  timing <- dos.time({
  	# Rem: there are faster ways to do this
  	# but here we want to time loops (220*220 'for' loops)! 
    for (j in 1:220) {
      for (k in 1:220) {
        b[k,j] <- abs(j - k) + 1
      }
    }
  })
  cumulate <- cumulate + timing
}
timing <- cumulate/runs
times[4, 3] <- timing
cat(c("Creation of a 220x220 Toeplitz matrix (loops)_______ (sec): ", timing, "\n"))
remove("b")

# (5)
cumulate <- 0; p <- 0; vt <- 0; vr <- 0; vrt <- 0; rvt <- 0; RV <- 0; j <- 0; k <- 0;
x2 <- 0; R <- 0; Rxx <- 0; Ryy <- 0; Rxy <- 0; Ryx <- 0; Rvmax <- 0
# Calculate the trace of a matrix (sum of its diagonal elements)
Trace <- function(y) {sum(c(y)[1 + 0:(min(dim(y)) - 1) * (dim(y)[1] + 1)], na.rm=FALSE)}
for (i in 1:runs) {
  x <- abs(rnorm(37*37)); dim(x) <- c(37, 37)
  timing <- dos.time({
    # Calculation of Escoufier's equivalent vectors
    p <- ncol(x)
    vt <- 1:p                                  # Variables to test
    vr <- NULL                                 # Result: ordered variables
    RV <- 1:p                                  # Result: correlations
    vrt <- NULL
    for (j in 1:p) {                           # loop on the variable number
      Rvmax <- 0
      for (k in 1:(p-j+1)) {                   # loop on the variables
        x2 <- cbind(x, x[,vr], x[,vt[k]])
        R <- cor(x2)                           # Correlations table
        Ryy <- R[1:p, 1:p]
        Rxx <- R[(p+1):(p+j), (p+1):(p+j)]
        Rxy <- R[(p+1):(p+j), 1:p]
        Ryx <- t(Rxy)
        rvt <- Trace(Ryx %*% Rxy) / sqrt(Trace(Ryy %*% Ryy) * Trace(Rxx %*% Rxx)) # RV calculation
        if (rvt > Rvmax) {
          Rvmax <- rvt                         # test of RV
          vrt <- vt[k]                         # temporary held variable
        }
      }
      vr[j] <- vrt                             # Result: variable
      RV[j] <- Rvmax                           # Result: correlation
      vt <- vt[vt!=vr[j]]                      # reidentify variables to test
    }
  })
  cumulate <- cumulate + timing
}
times[5, 3] <- timing
cat(c("Escoufier's method on a 37x37 matrix (mixed)________ (sec): ", timing, "\n"))
remove(c("x", "p", "vt", "vr", "vrt", "rvt", "RV"))
remove(c("x2", "R", "Rxx", "Ryy", "Rxy", "Ryx", "Rvmax", "Trace")) 

times[ , 3] <- sort(times[ , 3])
cat("                      --------------------------------------------\n")
cat(c("                      Trimmed mean (2 extremes eliminated): ", mean(times[2:4, 3]), "\n\n"))

cat(c("Total time for all 15 tests_________________________ (sec): ", sum(times), "\n"))
cat(c("Overall mean (sum of I, II and III trimmed means/3)_ (sec): ", exp(mean(log(times[2:4, ]))), "\n"))
remove(c("cumulate", "timing", "times", "runs"))
cat("                      --- End of test ---\n\n")