[Vegan-commits] r436 - in pkg: R inst man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 23 07:22:37 CEST 2008
Author: jarioksa
Date: 2008-06-23 07:22:37 +0200 (Mon, 23 Jun 2008)
New Revision: 436
Added:
pkg/src/nestedness.c
Modified:
pkg/R/commsimulator.R
pkg/inst/ChangeLog
pkg/man/oecosimu.Rd
Log:
commsimulator: quasiswap in C and *much* faster
Modified: pkg/R/commsimulator.R
===================================================================
--- pkg/R/commsimulator.R 2008-06-18 06:40:29 UTC (rev 435)
+++ pkg/R/commsimulator.R 2008-06-23 05:22:37 UTC (rev 436)
@@ -53,22 +53,9 @@
out <- x
} else if (method == "quasiswap") {
out <- r2dtable(1, rowSums(x), colSums(x))[[1]]
- swp <- matrix(c(-1,1,1,-1), nrow=2)
- bad <- sum(out*out) - sum(out)
- iter <- 0
- while (bad > 0) {
- i <- .Internal(sample(nr, 2, FALSE, NULL))
- j <- .Internal(sample(nc, 2, FALSE, NULL))
- z <- out[i,j]
- if(z[1,1] > 0 && z[2,2] > 0 && sum(-swp*z) >= 2 ) {
- out[i,j] <- z + swp
- bad <- bad - 2*(sum(-swp*z)-2)
- }
- if(z[1,2] > 0 && z[2,1] > 0 && sum(swp*z) >= 2 ) {
- out[i,j] <- z - swp
- bad <- bad - 2*(sum(swp*z)-2)
- }
- }
+ out <- .C("quasiswap", m = as.integer(out), as.integer(nrow(x)),
+ as.integer(ncol(x)), PACKAGE = "vegan")$m
+ dim(out) <- dim(x)
out
}
else if (method == "backtrack") {
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-06-18 06:40:29 UTC (rev 435)
+++ pkg/inst/ChangeLog 2008-06-23 05:22:37 UTC (rev 436)
@@ -2,8 +2,18 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
-Version 1.14-5 (opened June 16, 2008)
+Version 1.14-6 (opened June 23, 2008)
+ * 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
+ sec. Actually, quasiswap is now much faster than ordinary swap
+ which also should be written in C. "Backtracking" is so much more
+ complicated code that it could probably never reach "quasiswap"
+ even if written in C, and it may be dropped in the future.
+
+Version 1.14-5 (closed June 19, 2008)
+
* scores.rda: scaling = 0 returns now unmodified scores from the
object (like documented) without multiplying by the scaling
constant. Gains argument 'const' for user-settable general scaling
Modified: pkg/man/oecosimu.Rd
===================================================================
--- pkg/man/oecosimu.Rd 2008-06-18 06:40:29 UTC (rev 435)
+++ pkg/man/oecosimu.Rd 2008-06-23 05:22:37 UTC (rev 436)
@@ -116,6 +116,9 @@
2001).
The other methods maintain both row and column frequencies.
+ Methods \code{swap} and \code{tswap} implement sequential methods,
+ where the matrix is changed only little in one step, but the changed
+ matrix is used as an input if the next step.
Methods \code{swap} and \code{tswap} inspect random 2x2 submatrices
and if they are checkerboard units, the order of columns is
swapped. This changes the matrix structure, but does not influence
@@ -132,12 +135,25 @@
and uses that ratio to thin the results, so that on average one swap
will be found per step of \code{tswap}. However, the checkerboard
frequency probably changes during swaps, but this is not taken into
- account in estimating the \code{thin}. Both swap methods are sequential,
- and the result of the swap is used as the input of the next swap. One
- swap still changes the matrix only little, and it may be useful to
+ account in estimating the \code{thin}. One swap still changes the
+ matrix only little, and it may be useful to
thin the results so that the statistic is only evaluated after
\code{thin} steps.
+ Methods \code{quasiswap} and \code{backtracking} are not sequential,
+ but each call produces a matrix that is independent of previous
+ matrices, and has the same marginal totals as the original data. The
+ recommended method is \code{quasiswap} which is much faster because
+ it is implemented in C. Method \code{bactkracking} is provided for
+ comparison, but it is so slow that it may be dropped from future
+ releases of \pkg{vegan} (or also implemented in C).
+ Method \code{quasiswap} (\enc{Miklós}{Miklos} & Podani 2004)
+ implements a method where matrix is first filled
+ honouring row and column totals, but with integers that may be larger than
+ one. Then the method inspects random 2x2 matrices and performs a
+ quasiswap on them. Quasiswap is similar to ordinary swap, but it also
+ can reduce numbers above one to ones maintaining marginal
+ totals.
Method \code{backtracking}
implements a filling method with constraints both for row and column
frequencies (Gotelli & Entsminger 2001). The matrix is first filled
@@ -146,14 +162,7 @@
After that begins \dQuote{backtracking}, where some of the
points are removed, and then filling is started again, and this
backtracking is done so may times that all incidences will be filled
- into matrix. Method
- \code{quasiswap} (\enc{Miklós}{Miklos} & Podani 2004) implements a method where
- matrix is first filled
- honouring row and column totals, but with integers that may be larger than
- one. Then the method inspects random 2x2 matrices and performs a
- quasiswap on them. Quasiswap is similar to ordinary swap, but it also
- can reduce numbers above one to ones maintaining marginal
- totals. The \code{quasiswap} method is not sequential, but it produces
+ into matrix. The \code{quasiswap} method is not sequential, but it produces
a random incidence matrix with given marginal totals.
}
@@ -240,6 +249,11 @@
\examples{
data(sipoo)
nestedchecker(sipoo)
+## Matrix temperature
+out <- oecosimu(sipoo, nestedtemp, "r00")
+out
+plot(out)
+plot(out, kind="incid")
## Use the first eigenvalue of correspondence analysis as an index
## of structure: a model for making your own functions.
## This is a minimal structure; fancier functions give fancier results
Added: pkg/src/nestedness.c
===================================================================
--- pkg/src/nestedness.c (rev 0)
+++ pkg/src/nestedness.c 2008-06-23 05:22:37 UTC (rev 436)
@@ -0,0 +1,86 @@
+#include <R.h>
+
+/* Utility functions */
+
+/* Random integer */
+
+#define IRAND(imax) (int) (((double) (imax + 1)) * unif_rand())
+
+/* 2 different random integers */
+
+void i2rand(int *vec, int imax)
+{
+ vec[0] = IRAND(imax);
+ do {
+ vec[1] = IRAND(imax);
+ } while (vec[1] == vec[0]);
+}
+
+
+/*
+ * Quasiswap or sum-of-squares reducing swap of Miklos & Podani. A quasiswap
+ * step takes a random 2x2 submatrix and adds (-1,+1,+1,-1). If the submatrix
+ * was (1,0,0,1) it is swapped to (0,1,1,0), but if it was, say, (2,0,0,1) it
+ * is swapped to (1,1,1,0) which reduces sums-of-squares. We start with a
+ * random matrix with given marginal totals (from R r2dtable) but possibly
+ * some values >1. Then we perform quasiswaps on random 2x2 submatrices so
+ * long that only 1 and 0 are left. The function only does the quasiswaps,
+ * and it assumes that input matrix 'm' (dimensions 'nr', 'nc') was produced
+ * by r2dtable or some other function with given marginal totals, but some
+ * values possibly > 1.
+ */
+
+/* row & col indices to a vector index */
+
+#define INDX(i, j, nr) (i) + (nr)*(j)
+
+void quasiswap(int *m, int *nr, int *nc)
+{
+ int i, n, mtot, ss, row[2], col[2], nr1, nc1, a, b, c, d;
+
+ nr1 = (*nr) - 1;
+ nc1 = (*nc) - 1;
+
+ /* Get matrix total 'mtot' and sum-of-squares 'ss' */
+
+ n = (*nr) * (*nc);
+ for (i = 0, mtot = 0, ss = 0; i < n; i++) {
+ mtot += m[i];
+ ss += m[i] * m[i];
+ }
+
+ /* Get R RNG */
+ GetRNGstate();
+
+ /* Quasiswap while there are entries > 1 */
+
+ while (ss > mtot) {
+ i2rand(row, nr1);
+ i2rand(col, nc1);
+ /* a,b,c,d notation for a 2x2 table */
+ a = INDX(row[0], col[0], *nr);
+ b = INDX(row[0], col[1], *nr);
+ c = INDX(row[1], col[0], *nr);
+ d = INDX(row[1], col[1], *nr);
+ if (m[a] > 0 && m[d] > 0 && m[a] + m[d] - m[b] - m[c] >= 2) {
+ ss -= 2 * (m[a] + m[d] - m[b] - m[c] - 2);
+ m[a]--;
+ m[d]--;
+ m[b]++;
+ m[c]++;
+ } else if (m[b] > 0 && m[c] > 0 &&
+ m[b] + m[c] - m[a] - m[d] >= 2) {
+ ss -= 2 * (m[b] + m[c] - m[a] - m[d] - 2);
+ m[a]++;
+ m[d]++;
+ m[b]--;
+ m[c]--;
+ }
+ }
+
+ /* Set R RNG */
+ PutRNGstate();
+}
+
+#undef IRAND
+#undef INDX
More information about the Vegan-commits
mailing list