[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