rcpp/html/benchmark_8r_source.html

00001 #!/usr/bin/r -t
00002 #
00003 # Comparison benchmark
00004 #
00005 # This shows how Armadillo improves on the previous version using GNU GSL,
00006 # and how both are doing better than lm.fit()
00007 #
00008 # Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
00009 #
00010 # This file is part of Rcpp.
00011 #
00012 # Rcpp is free software: you can redistribute it and/or modify it
00013 # under the terms of the GNU General Public License as published by
00014 # the Free Software Foundation, either version 2 of the License, or
00015 # (at your option) any later version.
00016 #
00017 # Rcpp is distributed in the hope that it will be useful, but
00018 # WITHOUT ANY WARRANTY; without even the implied warranty of
00019 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00020 # GNU General Public License for more details.
00021 #
00022 # You should have received a copy of the GNU General Public License
00023 # along with Rcpp.  If not, see <http://www.gnu.org/licenses/>.
00024
00025 source("lmArmadillo.R")
00026 source("lmGSL.R")
00027
00028 set.seed(42)
00029 n <- 25000
00030 k <- 9
00031 X <- cbind( rep(1,n), matrix(rnorm(n*k), ncol=k) )
00032 truecoef <- 1:(k+1)
00033 y <- as.numeric(X %*% truecoef + rnorm(n))
00034
00035 N <- 100
00036
00037 lmgsl <- lmGSL()
00038 lmarma <- lmArmadillo()
00039
00040 tlm <- mean(replicate(N, system.time( lmfit <- lm(y ~ X - 1) )["elapsed"]), trim=0.05)
00041 tlmfit <- mean(replicate(N, system.time(lmfitfit <- lm.fit(X, y))["elapsed"]), trim=0.05)
00042 tlmgsl <- mean(replicate(N, system.time(lmgsl(y, X))["elapsed"]), trim=0.05)
00043 tlmarma <- mean(replicate(N, system.time(lmarma(y, X))["elapsed"]), trim=0.05)
00044
00045 res <- c(tlm, tlmfit, tlmgsl, tlmarma)
00046 data <- data.frame(results=res, ratios=tlm/res)
00047 rownames(data) <- c("lm", "lm.fit", "lmGSL", "lmArma")
00048 cat("For n=", n, " and k=", k, "\n", sep="")
00049 print(t(data))
00050 print(t(1/data[,1,drop=FALSE])) # regressions per second
00051