[Vegan-commits] r1819 - pkg/vegan/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 9 23:02:02 CEST 2011
Author: psolymos
Date: 2011-09-09 23:02:02 +0200 (Fri, 09 Sep 2011)
New Revision: 1819
Added:
pkg/vegan/R/permatfull1.R
pkg/vegan/R/permatswap1.R
Modified:
pkg/vegan/R/permatfull.R
pkg/vegan/R/permatswap.R
Log:
new functions return single matrix -> reused in old ones
Modified: pkg/vegan/R/permatfull.R
===================================================================
--- pkg/vegan/R/permatfull.R 2011-09-09 21:00:51 UTC (rev 1818)
+++ pkg/vegan/R/permatfull.R 2011-09-09 21:02:02 UTC (rev 1819)
@@ -2,37 +2,9 @@
`permatfull` <-
function(m, fixedmar="both", shuffle="both", strata=NULL, mtype="count", times=99)
{
-## internal function
-indshuffle <- function(x)
-{
- N <- length(x)
- n <- sum(x)
- out <- numeric(N)
- names(out) <- 1:N
- y <- table(sample(1:N, n, replace = TRUE))
- out[names(out) %in% names(y)] <- y
- names(out) <- NULL
- out
-}
-bothshuffle <- function(x, y=1)
-{
- x[x!=0] <- indshuffle(x[x!=0] - y) + y
- sample(x)
-}
- if (!identical(all.equal(m, round(m)), TRUE))
- stop("function accepts only integers (counts)")
mtype <- match.arg(mtype, c("prab", "count"))
shuffle <- match.arg(shuffle, c("ind", "samp", "both"))
- count <- mtype == "count"
fixedmar <- match.arg(fixedmar, c("none", "rows", "columns", "both"))
- sample.fun <- switch(shuffle,
- "ind"=indshuffle,
- "samp"=sample,
- "both"=bothshuffle)
- m <- as.matrix(m)
- n.row <- nrow(m)
- n.col <- ncol(m)
- if (mtype == "prab") m <- ifelse(m > 0, 1, 0)
if (is.null(strata))
str <- as.factor(rep(1, n.row))
@@ -43,31 +15,12 @@
nstr <- length(unique(str))
if (any(tapply(str,list(str),length) == 1))
stop("strata should contain at least 2 observations")
- perm <- list()
- perm[[1]] <- matrix(0, n.row, n.col)
- for (k in 1:times)
- perm[[k]] <- perm[[1]]
- for (j in 1:nstr) {
- id <- which(str == j)
- if (fixedmar == "none")
- for (i in 1:times)
- if (count) perm[[i]][id,] <- matrix(sample.fun(array(m[id,])), length(id), n.col)
- else perm[[i]][id,] <- commsimulator(m[id,], method="r00")
- if (fixedmar == "rows")
- for (i in 1:times)
- if (count) perm[[i]][id,] <- t(apply(m[id,], 1, sample.fun))
- else perm[[i]][id,] <- commsimulator(m[id,], method="r0")
- if (fixedmar == "columns")
- for (i in 1:times)
- if (count) perm[[i]][id,] <- apply(m[id,], 2, sample.fun)
- else perm[[i]][id,] <- commsimulator(m[id,], method="c0")
- if (fixedmar == "both")
- for (i in 1:times)
- if (count) perm[[i]][id,] <- r2dtable(1, apply(m[id,], 1, sum), apply(m[id,], 2, sum))[[1]]
- else perm[[i]][id,] <- commsimulator(m[id,], method="quasiswap")
- }
+
+ perm <- replicate(times, permatfull1(m, fixedmar, shuffle, strata, mtype), simplify=FALSE)
if (fixedmar == "both")
shuffle <- NA
+ if (mtype == "prab")
+ m <- ifelse(m > 0, 1, 0)
out <- list(call=match.call(), orig=m, perm=perm)
attr(out, "mtype") <- mtype
attr(out, "ptype") <- "full"
Added: pkg/vegan/R/permatfull1.R
===================================================================
--- pkg/vegan/R/permatfull1.R (rev 0)
+++ pkg/vegan/R/permatfull1.R 2011-09-09 21:02:02 UTC (rev 1819)
@@ -0,0 +1,66 @@
+## permatfull1 returns 1 permuted matrix
+## this function is called by permatfull
+`permatfull1` <-
+function(m, fixedmar="both", shuffle="both", strata=NULL, mtype="count")
+{
+ ## internal functions
+ indshuffle <- function(x)
+ {
+ N <- length(x)
+ n <- sum(x)
+ out <- numeric(N)
+ names(out) <- 1:N
+ y <- table(sample(1:N, n, replace = TRUE))
+ out[names(out) %in% names(y)] <- y
+ names(out) <- NULL
+ out
+ }
+ bothshuffle <- function(x, y=1)
+ {
+ x[x!=0] <- indshuffle(x[x!=0] - y) + y
+ sample(x)
+ }
+
+ if (!identical(all.equal(m, round(m)), TRUE))
+ stop("function accepts only integers (counts)")
+ mtype <- match.arg(mtype, c("prab", "count"))
+ shuffle <- match.arg(shuffle, c("ind", "samp", "both"))
+ count <- mtype == "count"
+ fixedmar <- match.arg(fixedmar, c("none", "rows", "columns", "both"))
+ sample.fun <- switch(shuffle,
+ "ind"=indshuffle,
+ "samp"=sample,
+ "both"=bothshuffle)
+ m <- as.matrix(m)
+ n.row <- nrow(m)
+ n.col <- ncol(m)
+ if (mtype == "prab")
+ m <- ifelse(m > 0, 1, 0)
+
+ if (is.null(strata))
+ str <- as.factor(rep(1, n.row))
+ else str <- as.factor(strata)[drop = TRUE]
+
+ levels(str) <- 1:length(unique(str))
+ str <- as.numeric(str)
+ nstr <- length(unique(str))
+ if (any(tapply(str,list(str),length) == 1))
+ stop("strata should contain at least 2 observations")
+ perm <- matrix(0, n.row, n.col)
+ for (j in 1:nstr) {
+ id <- which(str == j)
+ if (fixedmar == "none")
+ if (count) perm[id,] <- matrix(sample.fun(array(m[id,])), length(id), n.col)
+ else perm[id,] <- commsimulator(m[id,], method="r00")
+ if (fixedmar == "rows")
+ if (count) perm[id,] <- t(apply(m[id,], 1, sample.fun))
+ else perm[id,] <- commsimulator(m[id,], method="r0")
+ if (fixedmar == "columns")
+ if (count) perm[id,] <- apply(m[id,], 2, sample.fun)
+ else perm[id,] <- commsimulator(m[id,], method="c0")
+ if (fixedmar == "both")
+ if (count) perm[id,] <- r2dtable(1, apply(m[id,], 1, sum), apply(m[id,], 2, sum))[[1]]
+ else perm[id,] <- commsimulator(m[id,], method="quasiswap")
+ }
+ perm
+}
Modified: pkg/vegan/R/permatswap.R
===================================================================
--- pkg/vegan/R/permatswap.R 2011-09-09 21:00:51 UTC (rev 1818)
+++ pkg/vegan/R/permatswap.R 2011-09-09 21:02:02 UTC (rev 1819)
@@ -3,158 +3,50 @@
function(m, method="quasiswap", fixedmar="both", shuffle="both", strata=NULL,
mtype="count", times=99, burnin = 0, thin = 1)
{
-## internal function
-indshuffle <- function(x)
-{
- N <- length(x)
- n <- sum(x)
- out <- numeric(N)
- names(out) <- 1:N
- y <- table(sample(1:N, n, replace = TRUE))
- out[names(out) %in% names(y)] <- y
- names(out) <- NULL
- out
-}
-bothshuffle <- function(x, y=1)
-{
- x[x!=0] <- indshuffle(x[x!=0] - y) + y
- sample(x)
-}
- if (!identical(all.equal(m, round(m)), TRUE))
- stop("function accepts only integers (counts)")
mtype <- match.arg(mtype, c("prab", "count"))
fixedmar <- match.arg(fixedmar, c("rows", "columns", "both"))
shuffle <- match.arg(shuffle, c("samp", "both"))
count <- mtype == "count"
+ ## evaluating algo type
if (count) {
method <- match.arg(method, c("swap", "quasiswap", "swsh", "abuswap"))
## warning if swapcount is to be used
if (method == "swap") {
warning("quantitative swap method may not yield random null models, use only to study its properties")
isSeq <- TRUE
- if (fixedmar != "both")
- stop("if 'method=\"swap\"', 'fixedmar' must be \"both\"")
} else {
if (method == "abuswap") {
- if (fixedmar == "both")
- stop("if 'method=\"abuswap\"', 'fixedmar' must not be \"both\"")
- direct <- if (fixedmar == "columns")
- 0 else 1
isSeq <- TRUE
} else {
isSeq <- FALSE
- if (fixedmar != "both")
- stop("'fixedmar' must be \"both\"")
}
}
} else {
- if (fixedmar != "both")
- stop("if 'mtype=\"prab\"', 'fixedmar' must be \"both\"")
method <- match.arg(method, c("swap", "quasiswap", "tswap", "backtracking"))
isSeq <- method != "quasiswap"
- }
+ }
- m <- as.matrix(m)
- att <- attributes(m)
- n.row <- nrow(m)
- n.col <- ncol(m)
- if (mtype == "prab") m <- ifelse(m > 0, 1, 0)
-
- if (is.null(strata))
- str <- as.factor(rep(1, n.row))
- else str <- as.factor(strata)[drop = TRUE]
-
- levels(str) <- 1:length(unique(str))
- str <- as.numeric(str)
- nstr <- length(unique(str))
- if (any(tapply(str,list(str),length) == 1))
- stop("strata should contain at least 2 observations")
-
- perm <- list()
- perm[[1]] <- matrix(0, n.row, n.col)
- if (times > 1)
+ ## sequential algos: might have burnin
+ if (isSeq) {
+ perm <- vector("list", times)
+ ## with burnin
+ perm[[1]] <- permatswap1(m, method, fixedmar, shuffle, strata, mtype, burnin+thin)
for (i in 2:times)
- perm[[i]] <- perm[[1]]
-
- for (j in 1:nstr) {
- id <- which(str == j)
- temp <- m[id,]
- nn.row <- nrow(m[id,])
- nn.col <- ncol(m[id,])
- if (isSeq) {
- if (count) {
- if (burnin > 0) {
- if (method == "swap")
- temp <- .C("swapcount", m = as.double(temp),
- as.integer(nn.row), as.integer(nn.col),
- as.integer(burnin), PACKAGE = "vegan")$m
- if (method == "abuswap")
- temp <- .C("abuswap", m = as.double(temp),
- as.integer(nn.row), as.integer(nn.col),
- as.integer(burnin), as.integer(direct), PACKAGE = "vegan")$m
- }
- } else {
- if (burnin > 0)
- temp <- commsimulator(temp, method=method, thin = burnin)
- }
- for (i in 1:times) {
- if (count) {
- if (method == "swap")
- perm[[i]][id,] <- .C("swapcount",
- m = as.double(temp),
- as.integer(nn.row),
- as.integer(nn.col),
- as.integer(thin),
- PACKAGE = "vegan")$m
- if (method == "abuswap")
- perm[[i]][id,] <- .C("abuswap",
- m = as.double(temp),
- as.integer(nn.row),
- as.integer(nn.col),
- as.integer(thin),
- as.integer(direct),
- PACKAGE = "vegan")$m
- } else {
- perm[[i]][id,] <- commsimulator(temp, method=method, thin=thin)
- }
- temp <- perm[[i]][id,]
- } # for i end
- } else {
- if (method != "swsh") {
- r2tabs <- r2dtable(times, rowSums(m[id,]), colSums(m[id,]))
- } else {
- tempPos <- temp[temp > 0]
- }
- for (i in 1:times) {
- if (count) {
- if (method != "swsh") {
- ms <- sum(m[id,] > 0)
- tmp <- r2tabs[[i]]
- ## if fills are equal, no need to restore fill
- if (sum(tmp > 0) != ms) {
- tmp <- .C("rswapcount",
- m = as.double(tmp),
- as.integer(nn.row),
- as.integer(nn.col),
- as.integer(ms),
- PACKAGE="vegan")$m
- }
- perm[[i]][id,] <- matrix(tmp, nrow(perm[[i]][id,]), ncol(perm[[i]][id,]))
- } else { # method == "swsh"
- tmp <- commsimulator(temp, method="quasiswap")
- if (shuffle == "samp") {
- tmp[tmp > 0] <- sample(tempPos)
- } else {
- tmp[tmp > 0] <- bothshuffle(tempPos)
- }
- perm[[i]][id,] <- tmp
- }
- } else perm[[i]][id,] <- commsimulator(temp, method=method)
- }
- burnin <- 0
- thin <- 0
- }
- } # for j end
+ perm[[i]] <- permatswap1(perm[[i-1]], method, fixedmar, shuffle, strata, mtype, thin)
+ ## non-sequential algos: no burnin required
+ } else {
+ if (burnin > 0)
+ warnings("Non sequential algorithm used: 'burnin' argument ignored")
+ burnin <- 0
+ if (thin != 1)
+ warnings("Non sequential algorithm used: 'thin' value set to 1")
+ thin <- 1
+ perm <- replicate(times,
+ permatswap1(m, method, fixedmar, shuffle, strata, mtype, thin),
+ simplify=FALSE)
+ }
+ if (mtype == "prab")
+ m <- ifelse(m > 0, 1, 0)
out <- list(call=match.call(), orig=m, perm=perm)
attr(out, "mtype") <- mtype
attr(out, "ptype") <- "swap"
Added: pkg/vegan/R/permatswap1.R
===================================================================
--- pkg/vegan/R/permatswap1.R (rev 0)
+++ pkg/vegan/R/permatswap1.R 2011-09-09 21:02:02 UTC (rev 1819)
@@ -0,0 +1,136 @@
+## permatswap1 returns 1 permuted matrix
+## this function is called by permatswap
+`permatswap1` <-
+function(m, method="quasiswap", fixedmar="both", shuffle="both", strata=NULL,
+ mtype="count", thin = 1)
+{
+ ## internal functions
+ indshuffle <- function(x)
+ {
+ N <- length(x)
+ n <- sum(x)
+ out <- numeric(N)
+ names(out) <- 1:N
+ y <- table(sample(1:N, n, replace = TRUE))
+ out[names(out) %in% names(y)] <- y
+ names(out) <- NULL
+ out
+ }
+ bothshuffle <- function(x, y=1)
+ {
+ x[x!=0] <- indshuffle(x[x!=0] - y) + y
+ sample(x)
+ }
+ if (!identical(all.equal(m, round(m)), TRUE))
+ stop("function accepts only integers (counts)")
+ mtype <- match.arg(mtype, c("prab", "count"))
+ fixedmar <- match.arg(fixedmar, c("rows", "columns", "both"))
+ shuffle <- match.arg(shuffle, c("samp", "both"))
+ count <- mtype == "count"
+ if (thin < 1)
+ stop("'thin' must be >= 1")
+ if (count) {
+ method <- match.arg(method, c("swap", "quasiswap", "swsh", "abuswap"))
+ if (method == "swap") {
+ isSeq <- TRUE
+ if (fixedmar != "both")
+ stop("if 'method=\"swap\"', 'fixedmar' must be \"both\"")
+ } else {
+ if (method == "abuswap") {
+ if (fixedmar == "both")
+ stop("if 'method=\"abuswap\"', 'fixedmar' must not be \"both\"")
+ direct <- if (fixedmar == "columns")
+ 0 else 1
+ isSeq <- TRUE
+ } else {
+ isSeq <- FALSE
+ if (fixedmar != "both")
+ stop("'fixedmar' must be \"both\"")
+ }
+ }
+ } else {
+ if (fixedmar != "both")
+ stop("if 'mtype=\"prab\"', 'fixedmar' must be \"both\"")
+ method <- match.arg(method, c("swap", "quasiswap", "tswap", "backtracking"))
+ isSeq <- method != "quasiswap"
+ }
+
+ m <- as.matrix(m)
+ att <- attributes(m)
+ n.row <- nrow(m)
+ n.col <- ncol(m)
+ if (mtype == "prab")
+ m <- ifelse(m > 0, 1, 0)
+
+ if (is.null(strata))
+ str <- as.factor(rep(1, n.row))
+ else str <- as.factor(strata)[drop = TRUE]
+
+ levels(str) <- 1:length(unique(str))
+ str <- as.numeric(str)
+ nstr <- length(unique(str))
+ if (any(tapply(str,list(str),length) == 1))
+ stop("strata should contain at least 2 observations")
+
+ perm <- matrix(0, n.row, n.col)
+
+ for (j in 1:nstr) {
+ id <- which(str == j)
+ temp <- m[id,]
+ nn.row <- nrow(m[id,])
+ nn.col <- ncol(m[id,])
+ if (isSeq) {
+ if (count) {
+ if (method == "swap")
+ perm[id,] <- .C("swapcount",
+ m = as.double(temp),
+ as.integer(nn.row),
+ as.integer(nn.col),
+ as.integer(thin),
+ PACKAGE = "vegan")$m
+ if (method == "abuswap")
+ perm[id,] <- .C("abuswap",
+ m = as.double(temp),
+ as.integer(nn.row),
+ as.integer(nn.col),
+ as.integer(thin),
+ as.integer(direct),
+ PACKAGE = "vegan")$m
+ } else {
+ perm[id,] <- commsimulator(temp, method=method, thin=thin)
+ }
+ } else {
+ if (method != "swsh") {
+ r2tabs <- r2dtable(1, rowSums(m[id,]), colSums(m[id,]))[[1]]
+ } else {
+ tempPos <- temp[temp > 0]
+ }
+ if (count) {
+ if (method != "swsh") {
+ ms <- sum(m[id,] > 0)
+ tmp <- r2tabs
+ ## if fills are equal, no need to restore fill
+ if (sum(tmp > 0) != ms) {
+ tmp <- .C("rswapcount",
+ m = as.double(tmp),
+ as.integer(nn.row),
+ as.integer(nn.col),
+ as.integer(ms),
+ PACKAGE="vegan")$m
+ }
+ perm[id,] <- matrix(tmp, nrow(perm[id,]), ncol(perm[id,]))
+ } else { # method == "swsh"
+ tmp <- commsimulator(temp, method="quasiswap")
+ if (shuffle == "samp") {
+ tmp[tmp > 0] <- sample(tempPos)
+ } else {
+ tmp[tmp > 0] <- bothshuffle(tempPos)
+ }
+ perm[id,] <- tmp
+ }
+ } else perm[id,] <- commsimulator(temp, method=method)
+ }
+ } # for j end
+ perm
+}
+
More information about the Vegan-commits
mailing list