mNNrand  <- function(loc, label, nrand, nnid=NULL, splancs=T) {
 
 # Estimate randomization distribution of cell counts 
 #  given specific pattern of locations and labels
 
 # output: (k^2) x nrand matrix
 
 if (is.null(nnid)) {
   nnid <- NNid(loc, splancs=splancs);  # find the nn of each pt
   }
   
 k <- length(unique(label));
 
 rand <- matrix(0,nrow=k*k, ncol=nrand);
 
 for (i in 1:nrand) {
    l <- sample(label);		# permute the labels
    rand[,i] <- as.vector(t(table(l,l[nnid])))
   }

temp <- names(table(l));
dimnames(rand) <- list(t(outer(temp,temp,paste,sep='')), NULL);
 rand;
 }
 
mNNrandtest <- function(rand, info) {
# compute Chi-square tests based on Expected counts and Var, using
#  each simulated randomization

# input: rand: k^2 x Nrand matrix of counts under random labelling
#        info: list with expected counts and VarN from mNNinfo

# output: (k+1) x Nrand matrix of k individual spp Chi-square test
#         and then the overall test for each column of counts in rand

nrand <- ncol(rand);
k <- sqrt(nrow(rand));	# number of types of points

expNij <- as.vector(t(info$EN));

randtest <- matrix(0,nrow=k+1,ncol=nrand);

for (j in 1:k) {
    # calculate Chi-sq test for each type of label
    # easy to deal with gen. inv. here by dropping last column
    
    # get indices for this block of counts
  i1 <- 1+(j-1)*k;	# 1 when j=1, k+1 when j=2, 
  i2 <- i1 + k-2;	# k-1 when j=1, 2k-1 when j=2
    
  e <- as.matrix(expNij[i1:i2]);
  vinv <- solve(info$VarN[i1:i2,i1:i2]);

# now for each observed randomization
  for (i in 1:nrand) {
    o <- as.matrix(rand[i1:i2,i]);
    randtest[j,i] <- t(o-e) %*% vinv %*% (o-e);
    }
  }

# now do omnibus test
e <- as.matrix(expNij);
vinv <- ginv(info$VarN);

for (i in 1:nrand) {
   o <- as.matrix(rand[,i]);
   randtest[k+1,i] <- t(o-e) %*% vinv %*% (o-e);
   }

randtest
}


mNNsim <- function(n, nrand) {

# simulate standardized distribution of cell counts in multivariate NN 
#  given only vector of # of labels of each type

# location process assumed to be CSR
#  labels randomly assigned to points

k <- length(n);		# number of types of points
npt <- sum(n);		# total number of points 

l <- rep(1:k, n);	# vector of labels

rand <- matrix(0,nrow=k*k, ncol=nrand);
randtest <- matrix(0,nrow=k+1,ncol=nrand);
	# row rows 1:k are Chi-sq tests for type i, k+1 = omnibus Chi-sq
	
for (i in 1:nrand) {
  loc <- matrix(runif(2*npt),ncol=2);	# CSR locations
  	# because these aren't sorted, labels are randomized w.r.t locs
  temp <- mNNinfo(loc,l);
  	# get observed, expected counts
  obsNij <- as.vector(t(temp$ON));
  expNij <- as.vector(t(temp$EN))
  rand[,i] <- (obsNij - expNij)/sqrt(diag(temp$VarN));
  
  for (j in 1:k) {
    # calculate Chi-sq test for each type of label
    # easy to deal with gen. inv. here by dropping last column
    
    # get indices for this block of counts
    i1 <- 1+(j-1)*k;	# 1 when j=1, k+1 when j=2, 
    i2 <- i1 + k-2;	# k-1 when j=1, 2k-1 when j=2
    
    
    o <- as.matrix(obsNij[i1:i2]);
    e <- as.matrix(expNij[i1:i2]);
    v <- temp$VarN[i1:i2,i1:i2];
    randtest[j,i] <- t(o-e) %*% solve(v) %*% (o-e);
    }
    # now do omnibus test
    o <- as.matrix(obsNij);
    e <- as.matrix(expNij);
    
    randtest[k+1,i] <- t(o-e) %*% ginv(temp$VarN) %*% (o-e);
    }
    list(rand=rand, randtest=randtest);
    
}