PermutationTestMeans _ function(data, B = 999, 
				treatment,
				alternative = "two.sided",
				combine=NULL,
				combinationFunction = combinePValues.Fisher,
				seed=.Random.seed){

if(any(!is.element(alternative, c("two.sided", "greater", "less"))))
stop("alternative was not two.sided, greater, or less")

n _ numRows(data)
if(!is.matrix(data))
        dim(data) _ c(n, length(data)/n)
u _ unique(treatment)
if(length(u) != 2)
	stop("PermutationTestMeans compare two means, 
		requires that treatment have two unique values")
set1 _ which(treatment == u[1])
set2 _ which(treatment == u[2])
seed.start _ seed

if(!missing(seed)) set.seed(seed)

observed _ colMeans(data[set1,]) - colMeans(data[set2,])
p _ length(observed)
replicates _ unlist(lapply(1:B,function(j, indices1, data) 
                            colSums(data[indices1[,j],,drop=F]),
                   indices1 = samp.permute(n,B)[set1,,drop=F],data = data))
replicates _ (replicates/length(set1) -(colSums(data)-replicates)/length(set2))
dim(replicates) _ c(p, B)
replicates _ t(replicates)

alternative _ rep(alternative, length=p)

if(!length(combine)){
# Find individual p-value(s), careful about ties.
Pvalues _ rep(NA, p)
for(i in 1:p){
               X _ c(observed[i], replicates[,i])
		counts _ switch(alternative[i],
		less = sum(observed[i] >= replicates[,i]),
		greater = sum(observed[i] <=replicates[,i]),
		two.sided = sum( abs(observed[i])<= abs(replicates[,i]))
                               )
Pvalues[i] _ (counts+1)/(B+1)
	     }
                    }
else {
# Find individual p-value(s), careful about ties.
# And get p-values for every replication, then combine them.
  lambda _ matrix(NA, 1 + B, p)
  for(i in 1:p){
  X _ c(observed[i], replicates[,i])
  counts _ switch(alternative[i],
                less = ceiling(rank(X)),
                greater = ceiling(rank(-X)),
                two.sided = (B+2)-ceiling(rank(abs(X)))
               )
  lambda[,i] _ counts/(B+1)
               }

Pvalues _ lambda[1,]
comblambda _ numeric(B+1)

if(!is.list(combine)) combine _ list(combine)
combP _ lapply(combine,function(cols, lambda, f){
                                            comblambda _ f(lambda[ , cols, drop=F])
                                            sum(comblambda[1] <= comblambda)
                                                },
				            lambda=lambda, f=combinationFunction
               )
     }


obj _ list(call = match.call(),
           observed=observed,
	   replicates = replicates,
	   estimate = data.frame(alternative = alternative,
				   "p-value" = Pvalues),
	   B = B,
	   n = n,
 	   dim.obs = NULL,
 	   "p-value" = Pvalues,
	   "combined-p-value" = if(length(combine)) unlist(combP) / (B+1),
	   seed.start = seed.start,
	   seed.end = .Random.seed)

oldClass(obj) _ c("PermutationTestMeans", "resamp")
obj
}