[Picante-commits] r154 - in branches/gsoc: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 25 01:06:10 CEST 2008
Author: skembel
Date: 2008-07-25 01:06:10 +0200 (Fri, 25 Jul 2008)
New Revision: 154
Modified:
branches/gsoc/R/randomizeSample.R
branches/gsoc/src/picante.c
Log:
Implement simple richness-frequency null code in C, cleaning up randomizeSample R code to allow comparison of different nulls.
Modified: branches/gsoc/R/randomizeSample.R
===================================================================
--- branches/gsoc/R/randomizeSample.R 2008-07-24 08:42:30 UTC (rev 153)
+++ branches/gsoc/R/randomizeSample.R 2008-07-24 23:06:10 UTC (rev 154)
@@ -3,51 +3,77 @@
}
`randomizeSample` <-
-function(samp, null.model=c("frequency","richness","both"),it=10000) {
- samp<-as.matrix(samp)
-
- null.model <- match.arg(null.model)
+function(samp, null.model=c("frequency","richness","richnessindex","both","trialswap"), lang=c("C","R"), it=1000) {
+
+ samp <- as.matrix(samp)
+ lang <- match.arg(lang)
+ 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(lang,"R")) {
+ return(data.frame(apply(samp,2,sample),row.names=row.names(samp)))
+ } else {
+ ret1 <- .C("frequency", m=as.numeric(samp), as.integer(nrow(samp)), as.integer(ncol(samp)), PACKAGE="picante")
+ return(matrix(ret1$m,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(samp))))
+ }
}
+
if (identical(null.model,"richness")) {
+ if (identical(lang,"R")) {
+ return(t(data.frame(apply(samp,1,sample),row.names=colnames(samp))))
+ } else {
+ ret1 <- .C("richness", m=as.numeric(samp), as.integer(nrow(samp)), as.integer(ncol(samp)), PACKAGE="picante")
+ return(matrix(ret1$m,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(samp))))
+ }
+ }
+
+ if (identical(null.model,"richnessindex")) {
ret1 <- .C("richnessindex", m=as.numeric(samp), as.integer(nrow(samp)), as.integer(ncol(samp)), PACKAGE="picante")
- return(matrix(ret1$randm,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(samp))))
- #return(t(data.frame(apply(samp,1,sample),row.names=colnames(samp))))
+ return(matrix(ret1$randm,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(samp))))
}
+
if (identical(null.model,"both"))
{
- ret1 <- .C("trialswap", m=as.numeric(samp), as.integer(it), as.integer(nrow(samp)), as.integer(ncol(samp)), PACKAGE="picante")
- return(matrix(ret1$m,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(samp))))
+ if (identical(lang,"R")) {
+ #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)])
+ }
+ } else {
+ ret1 <- .C("independentswap", m=as.numeric(samp), as.integer(it), as.integer(nrow(samp)), as.integer(ncol(samp)), PACKAGE="picante")
+ return(matrix(ret1$m,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(samp))))
+ }
}
+
+ if (identical(null.model,"trialswap"))
+ {
+ if (identical(lang,"R")) {
+ stop("Trial swap implemented in C only")
+ } else {
+ ret1 <- .C("trialswap", m=as.numeric(samp), as.integer(it), as.integer(nrow(samp)), as.integer(ncol(samp)), PACKAGE="picante")
+ return(matrix(ret1$m,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(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)])
-# }
-
}
-
Modified: branches/gsoc/src/picante.c
===================================================================
--- branches/gsoc/src/picante.c 2008-07-24 08:42:30 UTC (rev 153)
+++ branches/gsoc/src/picante.c 2008-07-24 23:06:10 UTC (rev 154)
@@ -2,11 +2,11 @@
#include <stdlib.h>
#include <math.h>
#include <R.h>
-#include <Rmath.h>
+#include <Rmath.h>
/*libraries to include*/
-
+
/*This takes an integer n and returns one random integer*/
-int intrand(int n)
+int intrand(int n)
{
double u;
u = unif_rand();
@@ -15,16 +15,16 @@
/* inefficient but probably not important;
* allocate and copy in vector to matrix */
-double **vectomat(double *v,int row,int column)
+double **vectomat(double *v,int row,int column)
{
int i,j;
double **m;
m = (double **)R_alloc(row,sizeof(double *));
- for (i=0; i<row; i++)
+ for (i=0; i<row; i++)
{
m[i] = (double *)R_alloc(column,sizeof(double));
- for (j=0; j<column; j++)
+ for (j=0; j<column; j++)
{
m[i][j] = v[row*j+i]; /* R uses column-first ordering */
}
@@ -33,13 +33,13 @@
}
/* copy matrix back into vector */
-void mattovec(double *v,double **m,int row,int column)
+void mattovec(double *v,double **m,int row,int column)
{
int i,j,k;
k=0;
- for (j=0; j<column; j++)
+ for (j=0; j<column; j++)
{
- for (i=0; i<row; i++)
+ for (i=0; i<row; i++)
{
v[k++] = m[i][j];
}
@@ -47,7 +47,7 @@
}
-void trialswap(double *v, int *pintervals, int *prow, int *pcolumn)
+void trialswap(double *v, int *pintervals, int *prow, int *pcolumn)
{
long int trial;
int i,j,k,l;
@@ -84,72 +84,171 @@
mattovec(v,m,row,column);
PutRNGstate();
}
-
-void richnessindex(double *v, int *prow, int *pcolumn)
-{
- int i, j, k, q, p, b, sum=0, we;
- int row, column;
- double **m, **randm;
- row = *prow;
- column = *pcolumn;
- q = column;
- double *randrow[column];
-
- m = vectomat(v,row,column);//Change the vector v passed from R back to a Matrix
-
- /*make list of ONES THE LENGTH OF THE COLUMNS*/
- int s[column];
- int t;
- for (t=0;t<column;t++)
- {
- s[t]=1;
- }
- //////////////////////////////////
-
- randm = (double **)R_alloc(row,sizeof(double *)); //Initiate the Random Matrix
- //randrow = (double *)R_alloc(column,sizeof(double)); //Initiate a vector to hold the community that is being randomized
-
- /*Begin looping to create random matrix*/
- for (i=0; i<row; i++)
- {
- randm[i] = (double *)R_alloc(column,sizeof(double)); //Initiate the Random community at row i
-
- for (p=0;p<column;p++) //Get obs community i
- {
- randrow[i][p]=m[i][p];
- }
-
- for (j=0; j<column-1; j++) // randomly sample one species at a time across all rows of community i
- {
- for(b=0; b<q; b++)// Sum vector s to give the number to be passed to intrand
- {
- sum = sum + s[i];
- }
- k=intrand(sum);
- randm[i][j]=randrow[i][k];//Assign element randm i,j a random element from the community
- double *randrow[q-1];//Reinitialize randrow THIS WILL PROB BE A BUG IN THE CODE
-
- /*Fill randrow with all the elements of community i that have not been sampled yet
- A bug may also come up here if q=k...*/
- we=0;
- while(we<q-1)
- {
- for (p=0;p<q-k;p++)
- {
- randrow[i][we]=m[i][p];
- we++;
- }
- for (p=k+1;p<q;p++)
- {
- randrow[i][we]=m[i][p];
- we++;
- }
- }
- q=q-1;
- }
- }
-
- mattovec(v,randm,row,column);
-
- PutRNGstate();
-}
+
+void richnessindex(double *v, int *prow, int *pcolumn)
+{
+ int i, j, k, q, p, b, sum=0, we;
+ int row, column;
+ double **m, **randm;
+ row = *prow;
+ column = *pcolumn;
+ q = column;
+ double *randrow[column];
+
+ m = vectomat(v,row,column);//Change the vector v passed from R back to a Matrix
+
+ /*make list of ONES THE LENGTH OF THE COLUMNS*/
+ int s[column];
+ int t;
+ for (t=0;t<column;t++)
+ {
+ s[t]=1;
+ }
+ //////////////////////////////////
+
+ randm = (double **)R_alloc(row,sizeof(double *)); //Initiate the Random Matrix
+ //randrow = (double *)R_alloc(column,sizeof(double)); //Initiate a vector to hold the community that is being randomized
+
+ /*Begin looping to create random matrix*/
+ for (i=0; i<row; i++)
+ {
+ randm[i] = (double *)R_alloc(column,sizeof(double)); //Initiate the Random community at row i
+
+ for (p=0;p<column;p++) //Get obs community i
+ {
+ randrow[i][p]=m[i][p];
+ }
+
+ for (j=0; j<column-1; j++) // randomly sample one species at a time across all rows of community i
+ {
+ for(b=0; b<q; b++)// Sum vector s to give the number to be passed to intrand
+ {
+ sum = sum + s[i];
+ }
+ k=intrand(sum);
+ randm[i][j]=randrow[i][k];//Assign element randm i,j a random element from the community
+ double *randrow[q-1];//Reinitialize randrow THIS WILL PROB BE A BUG IN THE CODE
+
+ /*Fill randrow with all the elements of community i that have not been sampled yet
+ A bug may also come up here if q=k...*/
+ we=0;
+ while(we<q-1)
+ {
+ for (p=0;p<q-k;p++)
+ {
+ randrow[i][we]=m[i][p];
+ we++;
+ }
+ for (p=k+1;p<q;p++)
+ {
+ randrow[i][we]=m[i][p];
+ we++;
+ }
+ }
+ q=q-1;
+ }
+ }
+
+ mattovec(v,randm,row,column);
+
+ PutRNGstate();
+}
+
+void richness(double *v, int *prow, int *pcolumn)
+{
+
+ int i,j,k,l;
+ int row, column;
+ double tmp;
+ double **m;
+
+ row = *prow;
+ column = *pcolumn;
+
+ m = vectomat(v,row,column);
+
+ GetRNGstate();
+ for(i=0;i<row;i++)
+ {
+ for (j=0;j<column;j++)
+ {
+ while((k=intrand(column))==j);//choose another column (species) at random
+ tmp = m[i][j];
+ m[i][j] = m[i][k];
+ m[i][k] = tmp;
+ }
+ }
+ mattovec(v,m,row,column);
+ PutRNGstate();
+}
+
+void frequency(double *v, int *prow, int *pcolumn)
+{
+
+ int i,j,k,l;
+ int row, column;
+ double tmp;
+ double **m;
+
+ row = *prow;
+ column = *pcolumn;
+
+ m = vectomat(v,row,column);
+
+ GetRNGstate();
+ for(i=0;i<column;i++)
+ {
+ for (j=0;j<row;j++)
+ {
+ while((k=intrand(row))==j);//choose another row (sample) at random
+ tmp = m[j][i];
+ m[j][i] = m[k][i];
+ m[k][i] = tmp;
+ }
+ }
+ mattovec(v,m,row,column);
+ PutRNGstate();
+}
+
+void independentswap(double *v, int *pintervals, int *prow, int *pcolumn)
+{
+ long int swap;
+ int swapped;
+ int i,j,k,l;
+ int row, column;
+ int intervals;
+ double tmp;
+ double **m;
+
+ row = *prow;
+ column = *pcolumn;
+ intervals = *pintervals;
+
+ m = vectomat(v,row,column);
+
+ GetRNGstate();
+ for(swap=0;swap<intervals;swap++)
+ {
+ swapped = 0;
+ while(swapped==0) {
+ i=intrand(row);//Choose a random row
+ while((j=intrand(row))==i);//make sure that you do not randomly choose the same row as before
+ k=intrand(column);//Choose a random column
+ while((l=intrand(column))==k);//make sure that you do not randomly choose the same column as before
+ if((m[i][k]>0.0 && m[j][l]>0.0 && m[i][l]+m[j][k]==0.0)||(m[i][k]+m[j][l]==0.0 && m[i][l]>0.0 && m[j][k]>0.0))
+ {
+ //currently swaps abundances within columns (=species)
+ //should have a switch to swap abundances within rows, columns, or random
+ tmp = m[i][k];
+ m[i][k] = m[j][k];
+ m[j][k] = tmp;
+ tmp = m[i][l];
+ m[i][l] = m[j][l];
+ m[j][l] = tmp;
+ swapped = 1;
+ }
+ }
+ }
+ mattovec(v,m,row,column);
+ PutRNGstate();
+}
More information about the Picante-commits
mailing list