[Vegan-commits] r1930 - in pkg/vegan: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 4 21:42:48 CEST 2011
Author: jarioksa
Date: 2011-10-04 21:42:47 +0200 (Tue, 04 Oct 2011)
New Revision: 1930
Modified:
pkg/vegan/DESCRIPTION
pkg/vegan/R/permutest.cca.R
pkg/vegan/inst/ChangeLog
pkg/vegan/man/anova.cca.Rd
Log:
premature implementation of a rough version of parallelized permutest.cca
Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION 2011-10-04 09:05:04 UTC (rev 1929)
+++ pkg/vegan/DESCRIPTION 2011-10-04 19:42:47 UTC (rev 1930)
@@ -1,7 +1,7 @@
Package: vegan
Title: Community Ecology Package
-Version: 2.1-2
-Date: September 29, 2011
+Version: 2.1-3
+Date: October 4, 2011
Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre,
Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos,
M. Henry H. Stevens, Helene Wagner
Modified: pkg/vegan/R/permutest.cca.R
===================================================================
--- pkg/vegan/R/permutest.cca.R 2011-10-04 09:05:04 UTC (rev 1929)
+++ pkg/vegan/R/permutest.cca.R 2011-10-04 19:42:47 UTC (rev 1930)
@@ -7,11 +7,51 @@
`permutest.cca` <-
function (x, permutations = 99,
model = c("reduced", "direct", "full"), first = FALSE,
- strata, ...)
+ strata = NULL, parallel = 1, ...)
{
model <- match.arg(model)
isCCA <- !inherits(x, "rda")
isPartial <- !is.null(x$pCCA)
+ ## Function to get the F statistics in one loop
+ getF <- function (R, ...)
+ {
+ mat <- matrix(0, nrow = R, ncol = 3)
+ for (i in seq_len(R)) {
+ take <- permuted.index(N, strata)
+ Y <- E[take, ]
+ if (isCCA)
+ wtake <- w[take]
+ if (isPartial) {
+ if (isCCA) {
+ XZ <- .C("wcentre", x = as.double(Z), as.double(wtake),
+ as.integer(N), as.integer(Zcol),
+ PACKAGE = "vegan")$x
+ dim(XZ) <- c(N, Zcol)
+ QZ <- qr(XZ)
+ }
+ Y <- qr.resid(QZ, Y)
+ }
+ if (isCCA) {
+ XY <- .C("wcentre", x = as.double(X), as.double(wtake),
+ as.integer(N), as.integer(Xcol),
+ PACKAGE = "vegan")$x
+ dim(XY) <- c(N, Xcol)
+ Q <- qr(XY)
+ }
+ tmp <- qr.fitted(Q, Y)
+ if (first)
+ cca.ev <- La.svd(tmp, nv = 0, nu = 0)$d[1]^2
+ else cca.ev <- sum(tmp * tmp)
+ if (isPartial || first) {
+ tmp <- qr.resid(Q, Y)
+ ca.ev <- sum(tmp * tmp)
+ }
+ else ca.ev <- Chi.tot - cca.ev
+ mat[i,] <- cbind(cca.ev, ca.ev, (cca.ev/q)/(ca.ev/r))
+ }
+ mat
+ }
+ ## end getF()
if (first) {
Chi.z <- x$CCA$eig[1]
q <- 1
@@ -20,7 +60,8 @@
Chi.z <- x$CCA$tot.chi
names(Chi.z) <- "Model"
q <- x$CCA$qrank
- }
+ }
+ ## Set up
Chi.xz <- x$CA$tot.chi
names(Chi.xz) <- "Residual"
r <- nrow(x$CA$Xbar) - x$CCA$QR$rank - 1
@@ -30,9 +71,6 @@
if (!isCCA)
Chi.tot <- Chi.tot * (nrow(x$CCA$Xbar) - 1)
F.0 <- (Chi.z/q)/(Chi.xz/r)
- F.perm <- numeric(permutations)
- num <- numeric(permutations)
- den <- numeric(permutations)
Q <- x$CCA$QR
if (isCCA) {
w <- x$rowsum # works with any na.action, weights(x) won't
@@ -62,41 +100,18 @@
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
- for (i in 1:permutations) {
- take <- permuted.index(N, strata)
- Y <- E[take, ]
- if (isCCA)
- wtake <- w[take]
- if (isPartial) {
- if (isCCA) {
- XZ <- .C("wcentre", x = as.double(Z), as.double(wtake),
- as.integer(N), as.integer(Zcol),
- PACKAGE = "vegan")$x
- dim(XZ) <- c(N, Zcol)
- QZ <- qr(XZ)
- }
- Y <- qr.resid(QZ, Y)
- }
- if (isCCA) {
- XY <- .C("wcentre", x = as.double(X), as.double(wtake),
- as.integer(N), as.integer(Xcol),
- PACKAGE = "vegan")$x
- dim(XY) <- c(N, Xcol)
- Q <- qr(XY)
- }
- tmp <- qr.fitted(Q, Y)
- if (first)
- cca.ev <- La.svd(tmp, nv = 0, nu = 0)$d[1]^2
- else cca.ev <- sum(tmp * tmp)
- if (isPartial || first) {
- tmp <- qr.resid(Q, Y)
- ca.ev <- sum(tmp * tmp)
- }
- else ca.ev <- Chi.tot - cca.ev
- num[i] <- cca.ev
- den[i] <- ca.ev
- F.perm[i] <- (cca.ev/q)/(ca.ev/r)
+ ## permutations
+ if (parallel > 1 && getRversion() >= "2.14" && require(parallel)
+ && .Platform$OS.type == "unix") {
+ R <- ceiling(permutations/parallel)
+ tmp <- do.call(rbind, mclapply(seq_len(parallel), getF, R = R,
+ mc.cores = parallel))
+ } else {
+ tmp <- getF(R = permutations)
}
+ num <- tmp[1:permutations,1]
+ den <- tmp[1:permutations,2]
+ F.perm <- tmp[1:permutations,3]
## Round to avoid arbitrary ordering of statistics due to
## numerical inaccuracy
F.0 <- round(F.0, 12)
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2011-10-04 09:05:04 UTC (rev 1929)
+++ pkg/vegan/inst/ChangeLog 2011-10-04 19:42:47 UTC (rev 1930)
@@ -2,8 +2,33 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
-Version 2.1-2 (opened September 29, 20119
+Version 2.1-3 (opened October 4, 2011)
+ * permutest.cca: First attempt of setting 'parallel' processing in
+ permutest.cca. Currently the parallelization only works in R
+ 2.14.0 (alpha) and later with the 'parallel' package, and in
+ unix-like operating systems (Linux and MacOS X were
+ tested). Function permutest.cca gets a new argument 'parallel'
+ (defaults 1) that gives the number of desired parallel
+ processes. The argument is silently ignored if the system is not
+ capable of parallel processing (missing 'parallel' package,
+ Windows). The argument may be bassed to permutest.cca() from
+ anova.cca(), but currently setting the random number generator
+ seed will fail, and the results probably will be wrong. This
+ feature is only for testing. The functionality cannot be included
+ cleanly: it depends on the package 'parallel', but suggesting
+ 'parallel' fails R CMD check in the current R release (2.13.2)
+ which does not yet have 'parallel'. So we get warnings:
+ 'library' or 'require' call not declared from: parallel, and
+ permutest.cca: no visible global function definition for
+ ‘mclapply’.
+ Perhaps we delay adding this feature, and cancel this submission
+ later. However, with these warnings, the function passes tests in
+ R 2.13.2. (It fails in R 2.14.0 alpha since it suggests 'rgl', and
+ that package fails in R 2.14.0.)
+
+Version 2.1-2 (opened October 4, 2011)
+
* permutest.cca could not be update()d, because "permutest.cca"
was not exported from NAMESPACE -- only "permutest" was
exported. Another buglet (and this calls for checking other 'call'
Modified: pkg/vegan/man/anova.cca.Rd
===================================================================
--- pkg/vegan/man/anova.cca.Rd 2011-10-04 09:05:04 UTC (rev 1929)
+++ pkg/vegan/man/anova.cca.Rd 2011-10-04 19:42:47 UTC (rev 1930)
@@ -25,7 +25,7 @@
\method{permutest}{cca}(x, permutations = 99,
model = c("reduced", "direct", "full"),
- first = FALSE, strata, ...)
+ first = FALSE, strata, parallel = 1, ...)
}
\arguments{
@@ -53,6 +53,13 @@
\item{strata}{An integer vector or factor specifying the strata for
permutation. If supplied, observations are permuted only within the
specified strata.}
+
+ \item{parallel}{Number of parallel processes. The parallel
+ processing is only possible in \R version 2.14.x and later, and
+ currently only works in unix-like operating systems, such as Linux
+ and MacOS X. The argument is silently ignored if the system is not
+ capable of parallel processing. }
+
}
\details{
Functions \code{anova.cca} and \code{permutest.cca} implement an ANOVA
More information about the Vegan-commits
mailing list