rcpp/html/benchmarkLongley_8r_source.html

00001 #!/usr/bin/r -t
00002 #
00003 # Comparison benchmark -- using old and small Longley data set
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 suppressMessages(library(utils))
00026 suppressMessages(library(Rcpp))
00027 suppressMessages(library(inline))
00028 suppressMessages(library(datasets))
00029
00030 source("lmArmadillo.R")
00031 source("lmGSL.R")
00032
00033 data(longley)
00034
00035 longleydm <- data.matrix(data.frame(intcp=1, longley))
00036 X <- longleydm[,-8]
00037 y <- as.numeric(longleydm[,8])
00038
00039 N <- 1000
00040
00041 lmgsl <- lmGSL()
00042 lmarma <- lmArmadillo()
00043
00044 tlm <- mean(replicate(N, system.time( lmfit <- lm(y ~ X - 1) )["elapsed"]), trim=0.05)
00045 tlmfit <- mean(replicate(N, system.time(lmfitfit <- lm.fit(X, y))["elapsed"]), trim=0.05)
00046 tlmgsl <- mean(replicate(N, system.time(lmgsl(y, X))["elapsed"]), trim=0.05)
00047 tlmarma <- mean(replicate(N, system.time(lmarma(y, X))["elapsed"]), trim=0.05)
00048
00049 res <- c(tlm, tlmfit, tlmgsl, tlmarma)
00050 data <- data.frame(results=res, ratios=tlm/res)
00051 rownames(data) <- c("lm", "lm.fit", "lmGSL", "lmArma")
00052 cat("For Longley\n")
00053 print(t(data))
00054 print(t(1/data[,1,drop=FALSE])) # regressions per second
00055