[Picante-commits] r145 - in branches/gsoc: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 23 04:33:01 CEST 2008
Author: mrhelmus
Date: 2008-07-23 04:33:00 +0200 (Wed, 23 Jul 2008)
New Revision: 145
Added:
branches/gsoc/src/
branches/gsoc/src/picante.c
Modified:
branches/gsoc/R/phylostruct.R
branches/gsoc/R/randomizeSample.R
Log:
not working C code for trialswap so Steve can see it.
Modified: branches/gsoc/R/phylostruct.R
===================================================================
--- branches/gsoc/R/phylostruct.R 2008-07-23 00:47:55 UTC (rev 144)
+++ branches/gsoc/R/phylostruct.R 2008-07-23 02:33:00 UTC (rev 145)
@@ -16,10 +16,10 @@
} else {
nulls<-switch(metric,
- psv = replicate(runs,mean(psv(randomizeSample(samp,null.model=null.model),tree,compute.var=FALSE)[,1])),
- psr = replicate(runs,mean(psr(randomizeSample(samp,null.model=null.model),tree,compute.var=FALSE)[,1])),
- pse = replicate(runs,mean(pse(randomizeSample(samp,null.model=null.model),tree)[,1])),
- psc = replicate(runs,mean(psc(randomizeSample(samp,null.model=null.model),tree)[,1])))
+ psv = replicate(runs,mean(psv(as.matrix(randomizeSample(samp,null.model=null.model)),tree,compute.var=FALSE)[,1])),
+ psr = replicate(runs,mean(psr(as.matrix(randomizeSample(samp,null.model=null.model)),tree,compute.var=FALSE)[,1])),
+ pse = replicate(runs,mean(pse(as.matrix(randomizeSample(samp,null.model=null.model)),tree)[,1])),
+ psc = replicate(runs,mean(psc(as.matrix(randomizeSample(samp,null.model=null.model)),tree)[,1])))
quantiles.null<-quantile(nulls,probs=c(alpha/2,1-(alpha/2)))
mean.null<-mean(nulls)
mean.obs<-switch(metric,
Modified: branches/gsoc/R/randomizeSample.R
===================================================================
--- branches/gsoc/R/randomizeSample.R 2008-07-23 00:47:55 UTC (rev 144)
+++ branches/gsoc/R/randomizeSample.R 2008-07-23 02:33:00 UTC (rev 145)
@@ -1,40 +1,50 @@
+.First.lib <- function(lib,pkg) {
+ library.dynam("picante",pkg,lib)
+}
+
`randomizeSample` <-
-function(samp, null.model=c("frequency","richness","both")) {
- null.model <- match.arg(null.model)
+function(samp, null.model=c("frequency","richness","both"),it=100) {
+ samp<-as.matrix(samp)
+
+ null.model <- match.arg(null.model)
if (identical(null.model,"frequency")) {
return(data.frame(apply(samp,2,sample),row.names=row.names(samp)))
}
if (identical(null.model,"richness")) {
return(t(data.frame(apply(samp,1,sample),row.names=colnames(samp))))
}
- if (identical(null.model,"both")) {
- #check for presence-absence and warn until abundance implemented
- x <- decostand(samp, "pa")
- if (!identical(x,samp)) stop("Null model currently requires a presence-absence matrix.")
- sppFreq <- apply(x, 2, sum)/ncol(x)
- siteRichness <- apply(x, 1, sum)
- sampleList <- vector()
- sppList <- vector()
- for (siteNum in 1:length(siteRichness)) {
- sampleList <- c(sampleList, rep(names(siteRichness)[siteNum],
- siteRichness[siteNum]))
- sppList <- c(sppList, sample(names(sppFreq), siteRichness[siteNum],
- replace = FALSE, prob = sppFreq))
- }
- shuffledList <- data.frame(sample = sampleList, species = sppList,
- p = rep(1, length(sppList)))
- shuffledMatrix <- tapply(shuffledList$p, list(shuffledList$sample,
- shuffledList$species), sum)
- shuffledMatrix[is.na(shuffledMatrix)] <- 0
- if (identical(dim(shuffledMatrix),dim(x))) {
- return(shuffledMatrix[rownames(x), ])
- }
- else
- {
- mergedFrame <- merge(shuffledMatrix, x, all.y = TRUE)
- mergedFrame[is.na(mergedFrame)] <- 0
- return(mergedFrame[rownames(x),colnames(x)])
- }
+ if (identical(null.model,"both"))
+ {
+ ret1<-.C("trialswap",m=as.integer(samp),as.integer(it),as.integer(nrow(samp)),as.integer(ncol(samp)),PACKAGE="picante")
+ return(matrix(ret1$m,nrow=nrow(samp)))
+
+# #check for presence-absence and warn until abundance implemented
+# x <- decostand(samp, "pa")
+# if (!identical(x,samp)) stop("Null model currently requires a presence-absence matrix.")
+# sppFreq <- apply(x, 2, sum)/ncol(x)
+# siteRichness <- apply(x, 1, sum)
+# sampleList <- vector()
+# sppList <- vector()
+# for (siteNum in 1:length(siteRichness)) {
+# sampleList <- c(sampleList, rep(names(siteRichness)[siteNum],
+# siteRichness[siteNum]))
+# sppList <- c(sppList, sample(names(sppFreq), siteRichness[siteNum],
+# replace = FALSE, prob = sppFreq))
+# }
+# shuffledList <- data.frame(sample = sampleList, species = sppList,
+# p = rep(1, length(sppList)))
+# shuffledMatrix <- tapply(shuffledList$p, list(shuffledList$sample,
+# shuffledList$species), sum)
+# shuffledMatrix[is.na(shuffledMatrix)] <- 0
+# if (identical(dim(shuffledMatrix),dim(x))) {
+# return(shuffledMatrix[rownames(x), ])
+# }
+# else
+# {
+# mergedFrame <- merge(shuffledMatrix, x, all.y = TRUE)
+# mergedFrame[is.na(mergedFrame)] <- 0
+# return(mergedFrame[rownames(x),colnames(x)])
+# }
}
}
Added: branches/gsoc/src/picante.c
===================================================================
--- branches/gsoc/src/picante.c (rev 0)
+++ branches/gsoc/src/picante.c 2008-07-23 02:33:00 UTC (rev 145)
@@ -0,0 +1,70 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <Rmath.h>
+
+int intrand(int n) {
+ double u;
+ u = unif_rand();
+ return((int)(u*n));
+}
+
+/* inefficient but probably not important;
+ * allocate and copy in vector to matrix */
+double **vectomat(double *v,int row,int column) {
+ int i,j;
+ double **m;
+
+ m = (double **)R_alloc(row,sizeof(int *));
+ for (i=0; i<row; i++) {
+ m[i] = (double *)R_alloc(column,sizeof(int));
+ for (j=0; j<column; j++) {
+ m[i][j] = v[row*j+i]; /* R uses column-first ordering */
+ }
+ }
+ return(m);
+}
+
+/* copy matrix back into vector */
+void mattovec(double *v,double **m,int row,int column) {
+ int i,j,k;
+ k=0;
+ for (j=0; j<column; j++) {
+ for (i=0; i<row; i++) {
+ v[k++] = m[i][j];
+ }
+ }
+}
+
+
+void trialswap(double *v, int *pintervals, int * prow, int * pcolumn) {
+ long int trial;
+ int i,j,k,l;
+ int row, column;
+ int intervals;
+ double **m;
+
+ row = *prow;
+ column = *pcolumn;
+ intervals = *pintervals;
+
+ m = vectomat(v,row,column);
+
+ GetRNGstate();
+ for(trial=0;trial<intervals;trial++)
+ {
+ i=intrand(row);
+ while((j=intrand(row))==i);
+ k=intrand(column);
+ while((l=intrand(column))==k);
+ if((m[i][k]*m[j][l]==1 && m[i][l]+m[j][k]==0)||(m[i][k]+m[j][l]==0 && m[i][l]*m[j][k]==1))
+ {
+ m[i][k]=1-m[i][k];
+ m[i][l]=1-m[i][l];
+ m[j][k]=1-m[j][k];
+ m[j][l]=1-m[j][l];
+ }
+ }
+ mattovec(v,m,row,column);
+ PutRNGstate();
+}
More information about the Picante-commits
mailing list