[Vegan-commits] r439 - in pkg: R inst src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 28 07:36:08 CEST 2008
Author: jarioksa
Date: 2008-06-28 07:36:08 +0200 (Sat, 28 Jun 2008)
New Revision: 439
Modified:
pkg/R/permat.R
pkg/R/swapcount.R
pkg/inst/ChangeLog
pkg/src/nestedness.c
Log:
Big permat update: C translation of swapcount, fixing matrix update and number of thins
in swapcount, changing vector operators (&,|) to logical (&&,||), and minor cleanup.
Modified: pkg/R/permat.R
===================================================================
--- pkg/R/permat.R 2008-06-28 05:15:59 UTC (rev 438)
+++ pkg/R/permat.R 2008-06-28 05:36:08 UTC (rev 439)
@@ -1,179 +1,195 @@
-## permatfull function
-`permatfull` <-
-function(m, fixedmar="both", reg=NULL, hab=NULL, mtype="count", times=100)
-{
- mtype <- match.arg(mtype, c("prab", "count"))
- count <- if (mtype == "count") TRUE else FALSE
- fixedmar <- match.arg(fixedmar, c("none", "rows", "columns", "both"))
- 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 (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)
- if (!is.null(reg) & !is.null(hab)) str <- interaction(reg, hab, 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()
- for (k in 1:times)
- perm[[k]] <- matrix(0, n.row, n.col)
- for (j in 1:nstr) {
- id <- which(str == j)
- if (fixedmar == "none")
- for (i in 1:times)
- if (count) perm[[i]][id,] <- matrix(sample(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,] <- apply(m[id,], 2, sample)
- else perm[[i]][id,] <- commsimulator(m[id,], method="r0")
- if (fixedmar == "cols")
- for (i in 1:times)
- if (count) perm[[i]][id,] <- t(apply(m[id,], 1, sample))
- 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")
- }
- specs <- list(reg=reg, hab=hab)
- out <- list(call=match.call(), orig=m, perm=perm, specs=specs)
- attr(out, "mtype") <- mtype
- attr(out, "ptype") <- "full"
- attr(out, "fixedmar") <- fixedmar
- attr(out, "times") <- times
- class(out) <- c("permat", "list")
- return(out)
-}
-
-## permatswap function
-`permatswap` <-
-function(m, reg=NULL, hab=NULL, mtype="count", method="swap", times=100, burnin = 10000, thin = 1000)
-{
- mtype <- match.arg(mtype, c("prab", "count"))
- count <- if (mtype == "count") TRUE else FALSE
- if (count) {
- method <- match.arg(method, "swap")
- } else {method <- match.arg(method, c("swap", "tswap", "backtrack"))}
-
- 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 (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)
- if (!is.null(reg) & !is.null(hab)) str <- interaction(reg, hab, 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()
- 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 <- swapcount(temp)
- else for (k in 1:burnin) temp <- commsimulator(temp, method=method)
- for (i in 1:times)
- if (count) perm[[i]][id,] <- swapcount(temp, thin=thin)
- else perm[[i]][id,] <- commsimulator(temp, method=method, thin=thin)
- temp <- perm[[i]][id,]}
- 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
- attr(out, "ptype") <- "swap"
- attr(out, "fixedmar") <- "both"
- attr(out, "times") <- times
- class(out) <- c("permat", "list")
- return(out)
-}
-
-## S3 plot method for permat
-`plot.permat` <-
-function(x, ...)
-{
- n <- attr(x, "times")
- bray <- numeric(n)
- for (i in 1:n) bray[i] <- sum(abs(x$orig-x$perm[[i]]))/sum(x$orig+x$perm[[i]])
- plot(bray,type="n",ylim=c(0,1),ylab="Bray-Curtis dissimilarity",xlab="Runs", ...)
- for (i in c(0.4, 0.6)) abline(i,0, lty=2, col="grey")
- lines(bray,col="red")
- lines(lowess(bray),col="blue",lty=2)
- title(sub=paste("(mean = ", substitute(z, list(z=round(mean(bray),3))),
- ", min = ", substitute(z, list(z=round(min(bray),3))),
- ", max = ", substitute(z, list(z=round(max(bray),3))), ")", sep=""))
- invisible(NULL)
-}
-
-## S3 print method for permat
-`print.permat` <-
-function(x, digits=3, ...)
-{
- if (attr(x, "ptype") != "sar" & !is.null(x$specs$reg) | !is.null(x$specs$hab))
- restr <- TRUE else restr <- FALSE
- cat("Object of class 'permat'\n\nCall: ")
- print(x$call)
- cat("Matrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
- cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
- cat("\n\nMatrix dimensions:", nrow(x$orig), "rows,", ncol(x$orig), "columns")
- cat("\nSum of original matrix:", sum(x$orig))
- cat("\nFill of original matrix:", round(sum(x$orig>0)/(nrow(x$orig)*ncol(x$orig)),digits))
- cat("\nNumber of permuted matrices:", length(x$perm),"\n")
-}
-
-## S3 summary method for permat
-`summary.permat` <-
-function(object, digits=2, ...)
-{
- x <- object
- n <- attr(x, "times")
- if (attr(x, "ptype") != "sar" & !is.null(x$specs$reg) | !is.null(x$specs$hab))
- restr <- TRUE else restr <- FALSE
- if (restr) {
- if (!is.null(x$specs$reg) & is.null(x$specs$hab)) int <- x$specs$reg
- if (is.null(x$specs$reg) & !is.null(x$specs$hab)) int <- x$specs$hab
- if (!is.null(x$specs$reg) & !is.null(x$specs$hab))
- int <- interaction(x$specs$reg, x$specs$hab, drop=TRUE)
+## permatfull function
+`permatfull` <-
+function(m, fixedmar="both", reg=NULL, hab=NULL, mtype="count", times=100)
+{
+ mtype <- match.arg(mtype, c("prab", "count"))
+ count <- mtype == "count"
+ fixedmar <- match.arg(fixedmar, c("none", "rows", "columns", "both"))
+ 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 (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)
+ if (!is.null(reg) && !is.null(hab)) str <- interaction(reg, hab, 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()
+ for (k in 1:times)
+ perm[[k]] <- matrix(0, n.row, n.col)
+ for (j in 1:nstr) {
+ id <- which(str == j)
+ if (fixedmar == "none")
+ for (i in 1:times)
+ if (count) perm[[i]][id,] <- matrix(sample(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,] <- apply(m[id,], 2, sample)
+ else perm[[i]][id,] <- commsimulator(m[id,], method="r0")
+ if (fixedmar == "cols")
+ for (i in 1:times)
+ if (count) perm[[i]][id,] <- t(apply(m[id,], 1, sample))
+ 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")
+ }
+ specs <- list(reg=reg, hab=hab)
+ out <- list(call=match.call(), orig=m, perm=perm, specs=specs)
+ attr(out, "mtype") <- mtype
+ attr(out, "ptype") <- "full"
+ attr(out, "fixedmar") <- fixedmar
+ attr(out, "times") <- times
+ class(out) <- c("permat", "list")
+ return(out)
+}
+
+## permatswap function
+`permatswap` <-
+function(m, reg=NULL, hab=NULL, mtype="count", method="swap", times=100, burnin = 10000, thin = 1000)
+{
+ mtype <- match.arg(mtype, c("prab", "count"))
+ count <- mtype == "count"
+ if (count) {
+ method <- match.arg(method, c("swap", "Cswap"))
+ } else {method <- match.arg(method, c("swap", "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 (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)
+ if (!is.null(reg) && !is.null(hab)) str <- interaction(reg, hab, 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()
+ 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 <- switch(method,
+ swap = swapcount(temp),
+ Cswap = .C("swapcount", m = as.double(temp),
+ as.integer(n.row), as.integer(n.col),
+ as.integer(1))$m)
+ else
+ for (k in 1:burnin)
+ temp <- commsimulator(temp, method=method)
+ for (i in 1:times) {
+ if (count)
+ perm[[i]][id,] <- switch(method,
+ swap = swapcount(temp, thin=thin),
+ Cswap = .C("swapcount",
+ m = as.double(temp),
+ as.integer(n.row),
+ as.integer(n.col),
+ as.integer(thin))$m) else perm[[i]][id,] <- commsimulator(temp, method=method, thin=thin)
+ temp <- perm[[i]][id,]
+ }
+ }
+ 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
+ attr(out, "ptype") <- "swap"
+ attr(out, "fixedmar") <- "both"
+ attr(out, "times") <- times
+ class(out) <- c("permat", "list")
+ return(out)
+}
+
+## S3 plot method for permat
+`plot.permat` <-
+function(x, ...)
+{
+ n <- attr(x, "times")
+ bray <- numeric(n)
+ for (i in 1:n) bray[i] <- sum(abs(x$orig-x$perm[[i]]))/sum(x$orig+x$perm[[i]])
+ plot(bray,type="n",ylim=c(0,1),ylab="Bray-Curtis dissimilarity",xlab="Runs", ...)
+ for (i in c(0.4, 0.6)) abline(i,0, lty=2, col="grey")
+ lines(bray,col="red")
+ lines(lowess(bray),col="blue",lty=2)
+ title(sub=paste("(mean = ", substitute(z, list(z=round(mean(bray),3))),
+ ", min = ", substitute(z, list(z=round(min(bray),3))),
+ ", max = ", substitute(z, list(z=round(max(bray),3))), ")", sep=""))
+ invisible(NULL)
+}
+
+## S3 print method for permat
+`print.permat` <-
+function(x, digits=3, ...)
+{
+ if (attr(x, "ptype") != "sar" & !is.null(x$specs$reg) | !is.null(x$specs$hab))
+ restr <- TRUE else restr <- FALSE
+ cat("Object of class 'permat'\n\nCall: ")
+ print(x$call)
+ cat("Matrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
+ cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
+ cat("\n\nMatrix dimensions:", nrow(x$orig), "rows,", ncol(x$orig), "columns")
+ cat("\nSum of original matrix:", sum(x$orig))
+ cat("\nFill of original matrix:", round(sum(x$orig>0)/(nrow(x$orig)*ncol(x$orig)),digits))
+ cat("\nNumber of permuted matrices:", length(x$perm),"\n")
+}
+
+## S3 summary method for permat
+`summary.permat` <-
+function(object, digits=2, ...)
+{
+ x <- object
+ n <- attr(x, "times")
+ if (attr(x, "ptype") != "sar" && !is.null(x$specs$reg) || !is.null(x$specs$hab))
+ restr <- TRUE else restr <- FALSE
+ if (restr) {
+ if (!is.null(x$specs$reg) && is.null(x$specs$hab)) int <- x$specs$reg
+ if (is.null(x$specs$reg) && !is.null(x$specs$hab)) int <- x$specs$hab
+ if (!is.null(x$specs$reg) && !is.null(x$specs$hab))
+ int <- interaction(x$specs$reg, x$specs$hab, drop=TRUE)
nlev <- length(unique(int))
- ssum <- numeric(n)}
- bray <- psum <- pfill <- vrow <- vcol <- numeric(n)
- for (i in 1:n) {
- bray[i] <- sum(abs(x$orig-x$perm[[i]]))/sum(x$orig+x$perm[[i]])
- psum[i] <- sum(x$orig) == sum(x$perm[[i]])
- pfill[i] <- sum(x$orig > 0) == sum(x$perm[[i]] > 0)
- vrow[i] <- sum(rowSums(x$orig) == rowSums(x$perm[[i]])) == nrow(x$orig)
- vcol[i] <- sum(colSums(x$orig) == colSums(x$perm[[i]])) == ncol(x$orig)
- if (restr) ssum[i] <- {sum(rowSums(aggregate(x$orig,list(int),sum)[,-1]) ==
- rowSums(aggregate(x$perm[[i]],list(int),sum)[,-1])) == nlev}
- }
- strsum <- if (restr) sum(ssum)/n else NA
- outv <- c(sum=sum(psum)/n, fill=sum(pfill)/n, rowsums=sum(vrow)/n, colsums=sum(vcol)/n, strsum=strsum)
-
- cat("Summary of object of class 'permat'\n\nCall: ")
- print(x$call)
- cat("Matrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
- cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
- cat("\n\nMatrix dimensions:", nrow(x$orig), "rows,", ncol(x$orig), "columns")
- cat("\nSum of original matrix:", sum(x$orig))
- cat("\nFill of original matrix:", round(sum(x$orig>0)/(nrow(x$orig)*ncol(x$orig)),digits))
- cat("\nNumber of permuted matrices:", length(x$perm),"\n")
-
- cat("\nMatrix sums retained:", round(100*outv[1], digits), "%")
- cat("\nMatrix fill retained:", round(100*outv[2], digits), "%")
- cat("\nRow sums retained: ", round(100*outv[3], digits), "%")
- cat("\nColumn sums retained:", round(100*outv[4], digits), "%")
- if (restr) cat("\nSums within strata retained:", round(100*outv[5], digits), "%")
- cat("\n\nBray-Curtis dissimilarities among original and permuted matrices:\n")
- print(summary(bray))
- out <- list(bray=bray, restr=outv)
- class(out) <- c("summary.permat", "list")
- invisible(out)
-}
+ ssum <- numeric(n)}
+ bray <- psum <- pfill <- vrow <- vcol <- numeric(n)
+ for (i in 1:n) {
+ bray[i] <- sum(abs(x$orig-x$perm[[i]]))/sum(x$orig+x$perm[[i]])
+ psum[i] <- sum(x$orig) == sum(x$perm[[i]])
+ pfill[i] <- sum(x$orig > 0) == sum(x$perm[[i]] > 0)
+ vrow[i] <- sum(rowSums(x$orig) == rowSums(x$perm[[i]])) == nrow(x$orig)
+ vcol[i] <- sum(colSums(x$orig) == colSums(x$perm[[i]])) == ncol(x$orig)
+ if (restr) ssum[i] <- {sum(rowSums(aggregate(x$orig,list(int),sum)[,-1]) ==
+ rowSums(aggregate(x$perm[[i]],list(int),sum)[,-1])) == nlev}
+ }
+ strsum <- if (restr) sum(ssum)/n else NA
+ outv <- c(sum=sum(psum)/n, fill=sum(pfill)/n, rowsums=sum(vrow)/n, colsums=sum(vcol)/n, strsum=strsum)
+
+ cat("Summary of object of class 'permat'\n\nCall: ")
+ print(x$call)
+ cat("Matrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
+ cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
+ cat("\n\nMatrix dimensions:", nrow(x$orig), "rows,", ncol(x$orig), "columns")
+ cat("\nSum of original matrix:", sum(x$orig))
+ cat("\nFill of original matrix:", round(sum(x$orig>0)/(nrow(x$orig)*ncol(x$orig)),digits))
+ cat("\nNumber of permuted matrices:", length(x$perm),"\n")
+
+ cat("\nMatrix sums retained:", round(100*outv[1], digits), "%")
+ cat("\nMatrix fill retained:", round(100*outv[2], digits), "%")
+ cat("\nRow sums retained: ", round(100*outv[3], digits), "%")
+ cat("\nColumn sums retained:", round(100*outv[4], digits), "%")
+ if (restr) cat("\nSums within strata retained:", round(100*outv[5], digits), "%")
+ cat("\n\nBray-Curtis dissimilarities among original and permuted matrices:\n")
+ print(summary(bray))
+ out <- list(bray=bray, restr=outv)
+ class(out) <- c("summary.permat", "list")
+ invisible(out)
+}
Modified: pkg/R/swapcount.R
===================================================================
--- pkg/R/swapcount.R 2008-06-28 05:15:59 UTC (rev 438)
+++ pkg/R/swapcount.R 2008-06-28 05:36:08 UTC (rev 439)
@@ -4,7 +4,6 @@
## internal, is the 2x2 matrix diagonal or anti-diagonal
isDiag <- function(x) {
x<- as.vector(x)
- x<- as.vector(x)
X <- as.numeric(x>0)
sX <- sum(X)
choose <- c(min(x[c(2,3)]), min(x[c(1,4)]))
@@ -13,14 +12,14 @@
d <- choose[ch]
if (ch == 2) ch <- -1
return(d * ch)}
- if (identical(X, c(0,1,1,0)) | identical(X, c(0,1,1,1)) | identical(X, c(1,1,1,0)))
+ if (identical(X, c(0,1,1,0)) || identical(X, c(0,1,1,1)) || identical(X, c(1,1,1,0)))
return(choose[1])
- if (identical(X, c(1,0,0,1)) | identical(X, c(1,0,1,1)) | identical(X, c(1,1,0,1)))
+ if (identical(X, c(1,0,0,1)) || identical(X, c(1,0,1,1)) || identical(X, c(1,1,0,1)))
return(-choose[2])
- if (sX < 2 | identical(X, c(0,0,1,1)) | identical(X, c(1,1,0,0)) |
- identical(X, c(0,1,0,1)) | identical(X, c(1,0,1,0)))
+ if (sX < 2 | identical(X, c(0,0,1,1)) || identical(X, c(1,1,0,0)) ||
+ identical(X, c(0,1,0,1)) || identical(X, c(1,0,1,0)))
return(0)
- }
+ }
x <- as.matrix(m)
n.col <- ncol(x)
n.row <- nrow(x)
@@ -31,9 +30,11 @@
ev <- isDiag(x[ran.row, ran.col])
if (ev != 0) {
if (identical(sum(x[ran.row, ran.col] > 0),
- sum(x[ran.row, ran.col] + matrix(c(ev,-ev,-ev,ev), 2, 2) > 0)))
- x[ran.row, ran.col] <- x[ran.row, ran.col] + matrix(c(ev,-ev,-ev,ev), 2, 2)
- changed <- changed + 1}
+ sum(x[ran.row, ran.col] + matrix(c(ev,-ev,-ev,ev), 2, 2) > 0))) {
+ x[ran.row, ran.col] <- x[ran.row, ran.col] + matrix(c(ev,-ev,-ev,ev), 2, 2)
+ changed <- changed + 1
+ }
}
+ }
return(x)
}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-06-28 05:15:59 UTC (rev 438)
+++ pkg/inst/ChangeLog 2008-06-28 05:36:08 UTC (rev 439)
@@ -4,6 +4,21 @@
Version 1.14-6 (opened June 23, 2008)
+ * permatswap (nestedness.c): translated Peter Solymos's
+ swapcount.R to C. This is still experimental code, and the user
+ interface is undocumented, except here: use method = "Cswap" in
+ permatswap.
+
+ * permatswap: was not updating swap matrix but using the same
+ starting matrix after burnin for every swap.
+
+ * swapcount.R: was advancing 'thin' counter even when a swap was
+ rejected so that fewer than requesed 'thins' were done.
+
+ * permat.R, swapcount.R: genereal cleanup, most importantly
+ replacing vector operators & and | with logical operators && and
+ ||.
+
* commsimulator: "quasiswap" written in C and *much* faster
now. Times are for MacBook Intel 1.86 GHz and 100 matrices:
"sipoo" from 7 min to 4 sec, "BCI" from 2+ hrs to 45
Modified: pkg/src/nestedness.c
===================================================================
--- pkg/src/nestedness.c 2008-06-28 05:15:59 UTC (rev 438)
+++ pkg/src/nestedness.c 2008-06-28 05:36:08 UTC (rev 439)
@@ -1,8 +1,9 @@
#include <R.h>
+#include <Rmath.h>
/* Utility functions */
-/* Random integer */
+/* Random integer 0..imax */
#define IRAND(imax) (int) (((double) (imax + 1)) * unif_rand())
@@ -143,8 +144,8 @@
m[b] = 1;
m[c] = 1;
break;
- } else if (m[c] == 1 && m[b] == 1 && m[d] == 0 &&
- m[a] == 0) {
+ }
+ if (m[c] == 1 && m[b] == 1 && m[d] == 0 && m[a] == 0) {
m[a] = 1;
m[d] = 1;
m[b] = 0;
@@ -156,5 +157,94 @@
PutRNGstate();
}
+
+/* 'swapcount' is a C translation of Peter Solymos's R code. It is
+ * similar to 'swap', but can swap > 1 values and so works for
+ * quantitative (count) data.
+ */
+
+
+/* 'isDiag' is a utility function for 'swapcount' to find the largest
+ * value that can be swapped and whether in diagonal or antidiagonal
+ * way. The input is a 2x2 submatrix 'sm'.
+*/
+
+double isDiag(double *sm)
+{
+ int i, sX;
+ double choose[2];
+
+ /* sX: number of non-zero cells */
+ for (i = 0, sX = 0; i++; i < 4)
+ if (sm[i] > 0)
+ sX++;
+
+ /* Smallest diagonal and antidiagonal element */
+ choose[0] = (sm[1] < sm[2]) ? sm[1] : sm[2];
+ choose[1] = (sm[0] < sm[3]) ? -sm[0] : -sm[3];
+
+ if (sX == 4) {
+ /* Either choose could be returned, but RNG is not needed,
+ * because sm already is in random order, and we always return
+ * choose[0] */
+ return choose[0];
+ }
+ if ((sm[0] == 0 && sm[1] > 0 && sm[2] > 0 && sm[3] == 0) ||
+ (sm[0] == 0 && sm[1] > 0 && sm[2] > 0 && sm[3] > 0) ||
+ (sm[0] > 0 && sm[1] > 0 && sm[2] > 0 && sm[3] == 0))
+ return choose[0];
+ if ((sm[0] > 0 && sm[1] == 0 && sm[2] == 0 && sm[3] > 0) ||
+ (sm[0] > 0 && sm[1] == 0 && sm[2] > 0 && sm[3] > 0) ||
+ (sm[0] > 0 && sm[1] > 0 && sm[2] == 0 && sm[3] > 0))
+ return choose[1];
+ if (sX < 2 ||
+ (sm[0] == 0 && sm[1] == 0 && sm[2] > 0 && sm[3] > 0) ||
+ (sm[0] > 0 && sm[1] > 0 && sm[2] == 0 && sm[3] == 0) ||
+ (sm[0] == 0 && sm[1] > 0 && sm[2] == 0 && sm[3] > 0) ||
+ (sm[0] > 0 && sm[1] == 0 && sm[2] > 0 && sm[3] == 0))
+ return 0;
+}
+
+void swapcount(double *m, int *nr, int *nc, int *thin)
+{
+ int row[2], col[2], k, ij[4], changed, oldn, newn,
+ pm[4] = {1, -1, -1, 1} ;
+ double sm[4], ev;
+
+ GetRNGstate();
+
+ changed = 0;
+ while (changed < *thin) {
+ /* Select a random 2x2 matrix*/
+ i2rand(row, *nr - 1);
+ i2rand(col, *nc - 1);
+ ij[0] = INDX(row[0], col[0], *nr);
+ ij[1] = INDX(row[1], col[0], *nr);
+ ij[2] = INDX(row[0], col[1], *nr);
+ ij[3] = INDX(row[1], col[1], *nr);
+ for (k = 0; k < 4; k ++)
+ sm[k] = m[ij[k]];
+ /* The largest value that can be swapped */
+ ev = isDiag(sm);
+ if (ev != 0) {
+ /* Check that the fill doesn't change*/
+ for (k = 0, oldn = 0, newn = 0; k < 4; k++) {
+ if(sm[k] > 0)
+ oldn++;
+ if (sm[k] + pm[k]*ev > 0)
+ newn++;
+ }
+ /* Swap */
+ if (oldn == newn) {
+ for (k = 0; k < 4; k++)
+ m[ij[k]] += pm[k]*ev;
+ changed++;
+ }
+ }
+ }
+
+ PutRNGstate();
+}
+
#undef IRAND
#undef INDX
More information about the Vegan-commits
mailing list