[Vegan-commits] r547 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Nov 2 21:23:56 CET 2008
Author: psolymos
Date: 2008-11-02 21:23:56 +0100 (Sun, 02 Nov 2008)
New Revision: 547
Modified:
pkg/R/permatfull.R
pkg/R/permatswap.R
Log:
permatswap has "quasiswap(count)" option
Modified: pkg/R/permatfull.R
===================================================================
--- pkg/R/permatfull.R 2008-11-02 16:05:36 UTC (rev 546)
+++ pkg/R/permatfull.R 2008-11-02 20:23:56 UTC (rev 547)
@@ -32,7 +32,7 @@
m <- as.matrix(m)
n.row <- nrow(m)
n.col <- ncol(m)
- if (mtype == "prab") m <- matrix(as.numeric(m > 0), n.row, n.col)
+ if (mtype == "prab") m <- ifelse(m > 0, 1, 0)
if (is.null(reg) && is.null(hab)) str <- as.factor(rep(1, n.row))
if (!is.null(reg) && is.null(hab)) str <- as.factor(reg)
if (is.null(reg) && !is.null(hab)) str <- as.factor(hab)
Modified: pkg/R/permatswap.R
===================================================================
--- pkg/R/permatswap.R 2008-11-02 16:05:36 UTC (rev 546)
+++ pkg/R/permatswap.R 2008-11-02 20:23:56 UTC (rev 547)
@@ -1,19 +1,88 @@
## permatswap function
`permatswap` <-
-function(m, reg=NULL, hab=NULL, mtype="count", method="swap", times=100, burnin = 10000, thin = 1000)
+function(m, method="quasiswap", reg=NULL, hab=NULL, mtype="count", times=100, burnin = 10000, thin = 1000)
{
+## temporary internal function
+## quasiswapcount is based on the idea of Carsten Dorman from function swap.web::bipartite
+quasiswapcount <-
+function(m)
+{
+
+## internal fun
+incRease <- function(x, increase)
+{
+ x<- as.vector(x)
+ X <- as.numeric(x>0)
+ ## sX: number of non-zero cells
+ sX <- sum(X)
+ ## Smallest antidiagonal(1) and diagonal(2) element
+ choose <- c(min(x[c(2,3)]), min(x[c(1,4)]))
+ ## if fill should be increased
+ if (increase) {
+ if (identical(X, c(0,1,1,0)) || identical(X, c(0,1,1,1)) || identical(X, c(1,1,1,0))) {
+ if (choose[1] > 1)
+ ## this step to modify choose[] is needed to increase fill
+ return(choose[1] - sample(choose[1]-1, 1))
+ else return(0)
+ } else {
+ if (identical(X, c(1,0,0,1)) || identical(X, c(1,0,1,1)) || identical(X, c(1,1,0,1))) {
+ if (choose[2] > 1)
+ return(-(choose[2] - sample(choose[2]-1, 1)))
+ else return(0)
+ } else return(0)
+ }
+ ## if fill should be decreased
+ } else {
+
+ if (identical(X, c(1,1,1,1)) || identical(X, c(0,1,1,1)) || identical(X, c(1,1,1,0))) {
+ return(choose[1])
+ } else {
+ if (identical(X, c(1,0,1,1)) || identical(X, c(1,1,0,1)))
+ return(-choose[2])
+ else return(0)
+ }
+ }
+} ## end of internal fun
+
+ x <- as.matrix(m)
+ n.col <- ncol(x)
+ n.row <- nrow(x)
+ fill.orig <- sum(x > 0)
+ ## marginal totals are retained, but not matrix fill
+ x <- r2dtable(1, apply(x, 1, sum), apply(x, 2, sum))[[1]]
+ ## if fill is o.k., no need for further steps (this is rarely the case)
+ if (sum(x > 0) == fill.orig)
+ return(x)
+ ## if fill is not o.k.
+ else {
+ ## while loop as long as fill is o.k.
+ fill.changed <- sum(x > 0)
+ while(fill.changed != fill.orig) {
+ rr <- sample(n.row, 2)
+ rc <- sample(n.col, 2)
+ ev <- incRease(x[rr, rc], increase = fill.changed < fill.orig)
+ if (ev != 0)
+ x[rr, rc] <- x[rr, rc] + matrix(c(ev,-ev,-ev,ev), 2, 2)
+ fill.changed <- sum(x > 0)
+ }
+ return(x)
+ }
+}
+
+
+
if (!identical(all.equal(m, round(m)), TRUE))
stop("function accepts only integers (counts)")
mtype <- match.arg(mtype, c("prab", "count"))
count <- mtype == "count"
if (count) {
- method <- match.arg(method, "swap")
- } else {method <- match.arg(method, c("swap", "tswap"))}
+ method <- match.arg(method, c("swap", "quasiswap"))
+ } else {method <- match.arg(method, c("swap", "quasiswap", "tswap"))}
m <- as.matrix(m)
n.row <- nrow(m)
n.col <- ncol(m)
- if (mtype == "prab") m <- matrix(as.numeric(m > 0), n.row, n.col)
+ if (mtype == "prab") m <- ifelse(m > 0, 1, 0)
if (is.null(reg) && is.null(hab)) str <- as.factor(rep(1, n.row))
if (!is.null(reg) && is.null(hab)) str <- as.factor(reg)
if (is.null(reg) && !is.null(hab)) str <- as.factor(hab)
@@ -26,30 +95,38 @@
perm <- list()
for (i in 1:times)
perm[[i]] <- matrix(0, n.row, n.col)
-
for (j in 1:nstr) {
id <- which(str == j)
temp <- m[id,]
- if (count)
- for (k in 1:burnin)
- temp <- .C("swapcount", m = as.double(temp),
- as.integer(n.row), as.integer(n.col),
- as.integer(1), PACKAGE = "vegan")$m
- else
- for (k in 1:burnin)
- temp <- commsimulator(temp, method=method)
- for (i in 1:times) {
+ if (method != "quasiswap") {
if (count)
- perm[[i]][id,] <- .C("swapcount",
- m = as.double(temp),
- as.integer(n.row),
- as.integer(n.col),
- as.integer(thin),
- PACKAGE = "vegan")$m
- else perm[[i]][id,] <- commsimulator(temp, method=method, thin=thin)
+ for (k in 1:burnin)
+ temp <- .C("swapcount", m = as.double(temp),
+ as.integer(n.row), as.integer(n.col),
+ as.integer(1), PACKAGE = "vegan")$m
+ else
+ for (k in 1:burnin)
+ temp <- commsimulator(temp, method=method)
+ for (i in 1:times) {
+ if (count)
+ perm[[i]][id,] <- .C("swapcount",
+ m = as.double(temp),
+ as.integer(n.row),
+ as.integer(n.col),
+ as.integer(thin),
+ PACKAGE = "vegan")$m
+ else perm[[i]][id,] <- commsimulator(temp, method=method, thin=thin)
temp <- perm[[i]][id,]
+ } # for i end
+ } else {
+ for (i in 1:times) {
+ if (count)
+ perm[[i]][id,] <- quasiswapcount(temp) ## this should be replaced by .C()
+ else perm[[i]][id,] <- commsimulator(temp, method=method)
+ }
+ thin <- burnin <- 0
}
- }
+ } # for j end
specs <- list(reg=reg, hab=hab, burnin=burnin, thin=thin)
out <- list(call=match.call(), orig=m, perm=perm, specs=specs)
attr(out, "mtype") <- mtype
@@ -60,3 +137,5 @@
class(out) <- c("permat", "list")
return(out)
}
+
+
More information about the Vegan-commits
mailing list