### Source code for required functions for FS CoDisp analyses

##################################
### INSTALL AND LOAD REQUIRED PACKAGES
##################################
# Function to Install and Load R Packages
# Install_And_Load <- function(Required_Packages) {
#     Remaining_Packages <- Required_Packages[!(Required_Packages %in% installed.packages()[,"Package"])];
#     if(length(Remaining_Packages)) {install.packages(Remaining_Packages)}
#     for(package_name in Required_Packages)
#     { library(package_name,character.only=TRUE,quietly=TRUE)  }
#     } # end function  
# 
#   Required_Packages=c("spatstat","geoR","fields","SpatialPack","ggplot2","raster","gstat")
#   Install_And_Load(Required_Packages)

library(spatstat)
library(geoR)
library(fields)
library(SpatialPack)
library(ggplot2)
library(grid)
library(raster)
library(gstat)

# library(ggmap)   # for graphing on a map
# library(ggsubplot) # for embedding graphs
# library(map) # for graphing on a map
# library(mapdata) # for graphing on a map

##################################
### SIMPLE FUNCTIONS
##################################

# basal area function: calculates basal area from DBH values (must be in cm)
basal.area.fn <- function(x){ (pi*(x)^2)/40000 } # calculate basal area in m^2

### Function to draw random values from a truncated log normal distribution
rtlnorm <- function (n, meanlog = 0, sdlog = 1, lower = -Inf, upper = Inf) 
{ 
  ret <- numeric() 
  if (n > 1) 
    n <- n 
  while (length(ret) < n) { 
    y <- rlnorm(n - length(ret), meanlog, sdlog) 
    y <- y[y >= lower & y <= upper] 
    ret <- c(ret, y) 
  } 
  stopifnot(length(ret) == n) 
  ret 
} 

### Function for simulating a bivariate normal distribution
bivariate <- function(x,y){
  mu1 <- 0     # expected value of x
  mu2 <- 0     # expected value of y
  sig1 <- 1    # variance of x
  sig2 <- 1    # variance of y
  rho <- 0.5   # corr(x, y)
  term1 <- 1 / (2 * pi * sig1 * sig2 * sqrt(1 - rho^2))
  term2 <- (x - mu1)^2 / sig1^2
  term3 <- -(2 * rho * (x - mu1)*(y - mu2))/(sig1 * sig2)
  term4 <- (y - mu2)^2 / sig2^2
  z <- term2 + term3 + term4
  term5 <- term1 * exp((-z / (2 *(1 - rho^2))))
  return (term5)
}


##################################
### DATA MANIPULATION
##################################

# List to array function for Co_disp null model output objects
list2ary = function(input.list){  #input a list of lists
  temp.ls <- vector("list")
  for(i in 1:length(input.list)) { temp.ls[i] <- input.list[[i]][1]  } # take the dataframes out of the list and put them in a new list
  rows.cols <- dim(temp.ls[[1]])
  sheets <- length(temp.ls)
  output.ary <- array(unlist(temp.ls), dim = c(rows.cols, sheets))
  colnames(output.ary) <- colnames(temp.ls[[1]])
  row.names(output.ary) <- row.names(temp.ls[[1]])
  return(output.ary)    # output as a 3-D array
}

### Function to extract values from an environmental grid (raster) at point locations in a ppp object. 
# Inputs are the ppp and geoR.env, which is a geoR object holding the environmental data layer.
# Extract works with a buffer around each point (10cm in this case)
# geoR.env <- geo.elev

extract.env.fn <- function(ppp.dat,geoR.env,xmin,xmax,ymin,ymax,quad.size=quad.size){
  
  ppp.df <- data.frame(x=ppp.dat$x,y=ppp.dat$y)   # create a dataframe from the ppp object
  
  X = geoR.env$coords[,1]+quad.size/2
  Y = geoR.env$coords[,2]+quad.size/2
  
  rast <- raster()  # create empty raster to add data to
  extent(rast) <- extent(c(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)) 
  ncol(rast) <- xmax/quad.size 
  nrow(rast) <- ymax/quad.size 
  raster.env <- rasterize(cbind(X,Y), rast, geoR.env$data)
  
  env.dat <- extract(x=raster.env,y=ppp.df,df=TRUE) # extract environmental data at tree locations
  
  env.out <- data.frame(x=ppp.df$x,y=ppp.df$y,z=env.dat$layer)
  env.geo <- as.geodata(env.out,coords.col=1:2,data.col=3)
  return(env.geo)
  
}  
  
#### Function to generate a geodata object (used by packages geoR and the codispersion function) from a ppp object.

# ppp.dat = input ppp object
# xmin, xmax, ymin, ymax = plot dimensions
# method = the measure that is used to generate the 'data' value for the geodata object

ppp.to.geoR.fn <- function(ppp.dat,xmin,xmax,ymin,ymax,quad.size,method=c("abundance","mean.mark","mean.ba","total.ba")){ # function to generate geoR objects with abundance and basal area in 20x20m quadrats. Note that DBH must be measured in cm. Input data= ppp object.
  x <- ppp.dat$x # extract x coordinate
  y <- ppp.dat$y # extract y coordinate
  z <- ppp.dat$marks # extract DBH values
  ba <- (pi*(z)^2)/40000 # calculate basal area in m^2 
  xt <- cut(x,seq(xmin,xmax,quad.size)) # cut x coordinates using 20m spacing
  yt <- cut(y,seq(ymin,ymax,quad.size)) # cut y coordinates using 20m spacing
  coords <- dimnames(table(yt,xt)) # extract quadrat coordinate lists
  qx <- rep(seq(xmin,xmax-quad.size,length=length(coords$xt)),each=length(coords$yt)) # vector of x coordinates for the bottom left corner of the quadrat
  qy <- rep(seq(ymin,ymax-quad.size,length=length(coords$yt)),length(coords$xt)) # vector of y coordinates for the bottom left corner of the quadrat
  if(method=="abundance"){
    out.grid <- table(yt,xt) # count the trees in each quadrat
    out.grid[is.na(out.grid)==T] <- 0 # replace NAs in table with zeros for empty quadrats
  }
  if(method=="mean.mark"){
    out.grid <- tapply(z,list(yt,xt),mean) # calculate mean DBH in each quadrat
    out.grid[is.na(out.grid)==T] <- 0
  }    
  if(method=="mean.ba"){
    out.grid <- tapply(ba,list(yt,xt),mean) # calculate mean ba in each quadrat
    out.grid[is.na(out.grid)==T] <- 0 
  }    
  if(method=="total.ba"){
    out.grid <- tapply(ba,list(yt,xt),sum) # calculate total ba in each quadrat    
    out.grid[is.na(out.grid)==T] <- 0 
  }  
  out.df <- data.frame(qx,qy,as.vector(out.grid))
  out.geo <- as.geodata(out.df,coords.col=1:2,data.col=3)
  return(out.geo)
} # end function

##################################
### CODISP ANALYSIS
##################################

#### Modified codispersion function (modified from Cuevas et al. 2013)
## See 'Box 1' in paper for a detailed explanation.
Codisp.Kern<-function(X,Y,k,h,gamma=1)
{
  Kernel<-function(u,gamma)
  {
    v=0
    v=ifelse(abs(u)<=1,(1/beta(0.5,gamma+1))*(1-u^2)^gamma,0)
  }
  ifelse(X$coords==Y$coords,1,
{
  break
  print("The coordinates of X and Y are different")
})

n=length(X$data)
mX <- matrix(X$data,nrow=n,ncol=n,byrow=FALSE)
mY <- matrix(Y$data,nrow=n,ncol=n,byrow=FALSE)
MatriXX <- (mX - t(mX))^2
MatriYY <- (mY - t(mY))^2
MatriXY <- (mX - t(mX))*(mY - t(mY))
mX <- matrix(X$coords[,1],nrow=n,ncol=n,byrow=FALSE)
DesignX <- mX - t(mX)
mY <- matrix(X$coords[,2],nrow=n,ncol=n,byrow=FALSE)
DesignY <- mY - t(mY)

KERNMATRIXX=Kernel((k[1]-DesignX)/h[1],gamma)*Kernel((k[2]-DesignY)/h[1],gamma)

if(h[1]==h[2]&h[1]==h[3]){
  KERNMATRIYY=KERNMATRIXX
  KERNMATRIXY=KERNMATRIXX } else{
  KERNMATRIYY=Kernel((k[1]-DesignX)/h[2],gamma)*Kernel((k[2]-DesignY)/h[2],gamma)
  KERNMATRIXY=Kernel((k[1]-DesignX)/h[3],gamma)*Kernel((k[2]-DesignY)/h[3],gamma) 
}

Numerador=sum(KERNMATRIXY*MatriXY)/(2*sum(KERNMATRIXY))
Denominador1=sum(KERNMATRIYY*MatriYY)/(2*sum(KERNMATRIYY))
Denominador2=sum(KERNMATRIXX*MatriXX)/(2*sum(KERNMATRIXX))
v1=Denominador1
v2=Denominador2
v3=Numerador
v4=Numerador/sqrt(Denominador1*Denominador2)
print(c(v1,v2,v3,v4))
}  

#### Function to run co-dispersion plot-level analyses. Input is two geodata objects (e.g. abundance and basal area)
codisp.coef.fn <- function(geodata1,geodata2){
  out <- vector("list",length=3)
  out[[1]] <- codisp(geodata1$data,geodata2$data,geodata1$coords,nclass=20) # calculate the codispersion coefficent for the whole plot
  out[[2]] <- cor.spatial(geodata1$data,geodata2$data,geodata1$coords) # Tjostheim's Coef
  out[[3]] <- modified.ttest(geodata1$data,geodata2$data,geodata1$coords,nclass=20) # modified t-test
  return(out)
}

#### Function to run co-dispersion window analysis (modified from Cuevas et al. 2013)

# geodata1 = first input data object (a geoR geodata object)
# geodata2 = second input object
# h = c(h1, h2, h3) = a vector of three bandwidth values for X, Y and XY
# grid.cell.size = the size of the cells in the geodata objects
# max.window.size = the maximum lag distance

codisp.fn <- function(geodata1,geodata2,h=h,grid.cell.size=grid.cell.size,max.window.size=max.window.size){
  out <- vector("list",length=2)
  X=geodata1  # input data process 1
  Y=geodata2  # input data process 2
  h=c(h[1],h[2],h[3]) # Set the bandwith for the kernel
  
  k_range <- max.window.size # set the spatial lags over which to calculate codisp
  k1=seq(-k_range,k_range,l=20)       # x-axis values for codispersion graph (lags)
  k2=seq(grid.cell.size,k_range,l=10) # y-axis values for codispersion graph (lags)
  
  MCodisp=matrix(0,ncol=10,nrow=20) # loop through the lags
  for(i in 1:20)     # 'left-right' lags
  {
    for(j in 1:10)   # 'up' lags
    {
      MCodisp[i,j]=Codisp.Kern(X,Y,c(k1[i],k2[j]),h)[4]; # calculate codisp
    }
  }
  Codispersion <- as.numeric(MCodisp) # save codisp object as output
  xx <- rep(k1,length(k2))            # write out values for x-axis
  yy <- rep(k2,each=length(k1))       # write out values for y-axis
  graphing.data <- data.frame(xx,yy,Codispersion) # graphing object
  
  # put both the graphing object and the original object in an output list
  out[[1]] <- graphing.data 
  out[[2]] <- MCodisp
  return(out) 
}


### NULL MODELS
#### Function to generate a list of 'nsim' ppp objects (marked point patterns) under three different null models
ppp.null.fn <- function(ppp.dat,nsim,model=c("RLM","HomP","HetP")){
  #ppp.dat <- ppp.dat[[1]]
  ppp.out <- vector("list",nsim) # create output list object  

  if(model=="RLM"){ # Random labelling model
    for(i in 1:nsim){ # start loop to generate simulations
      ppp.out[[i]] <- rlabel(ppp.dat, labels=marks(ppp.dat), permute=TRUE) # randomise marks
    } # end simulations loop
  } # end RLM loop  

  if(model=="HomP"){ # Homogeneous Poisson model
    for(i in 1:nsim){ # start loop to generate simulations
      ppp.HomP <- rpoint(ppp.dat$n,win=ppp.dat$win) # randomise the observed ppp
      ppp.HomP$marks <- sample(ppp.dat$marks, replace=F) # assign shuffled marks to new ppp
      ppp.out[[i]] <- ppp.HomP # add new marked ppp to output list
    } # end simulations loop    
  } # end HomP loop  
  
  if(model=="HetP"){ # this null model generates random marks based on a lognormal fit to the DBH distribution
    intensity_function <- density.ppp(ppp.dat, bw.diggle) # generate the intensity function  
    LN_params <- fitdistr(ppp.dat$marks,"log-normal") # fit lognormal to DBH distribution
    for(i in 1:nsim){ # start loop to generate simulations
      ppp.HetP <- rpoispp(intensity_function) # generate randomised ppp using intensity function      
      ppp.HetP$marks <- rtlnorm(ppp.HetP$n,meanlog=LN_params$estimate[[1]],sdlog=LN_params$estimate[[2]],1,max(ppp.dat$marks)) # generate marks using parameters of DBH distribution
      ppp.out[[i]] <- ppp.HetP # add new marked ppp to output list
    } # end simulations loop   
  } # end HetP loop  
  return(ppp.out)
} # end function


##################################
### DEALING WITH CODISP OUTPUTS
##################################

# Comparing observed values to null model results from output arrays and graphing average null model results
# inputs are the null model input array object, the observed CoDisp result list, and a choice of null model
# This function deals with abundance and basal area values
comparison.abba.fn = function(null.input.ary,CoDisp_obs,spe=spe,model=c("RLM","HomP","HetP")){
  out.df <- CoDisp_obs[[1]] # extract the observed Codispersion result as a dataframe
  for(i in 1:length(null.input.ary[,1,1])){ # loop through each cell
    nsims <- length(null.input.ary[1,1,])
    obser <- out.df$Codispersion[i] # observed codispersion value
    expec <- null.input.ary[i,3,]
    prop.greater.than <- length(which(expec>obser))/nsims
    prop.less.than <- length(which(expec<obser))/nsims
    out.df$P.value[i]<-min(prop.greater.than,prop.less.than)
  }
  
  out.df$null_mean <- apply(null.input.ary[,3,],MARGIN=1,mean) # calculate mean codispersion value for each cell from the array of null model results
  out.df$diff <- out.df$Codispersion-out.df$null_mean # observed minus expected
  out.df$P.value.cat <- factor(ifelse(out.df$P.value<0.025,"Sig.","Non-sig.")) # significance at alpha=0.05
  
  ## Graph results
  if(model=="RLM"){gtitle="Random labelling model"}
  if(model=="HomP"){gtitle="Homogeneous Poisson model"}
  if(model=="HetP"){gtitle="Heterogeneous Poisson model"}
  
  # observed graphs
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Codispersion of abundance and \n total basal area for",spe,"observed")))
  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7),limits=c(0, 1))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Codispersion of abundance and \n total basal area for",spe,"observed")))
  
  # null model graphs  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=null_mean))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Mean codispersion of abundance and \n total basal area for",spe,gtitle)))
  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=null_mean))+scale_fill_gradientn(colours=rainbow(7),limits=c(0, 1))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Mean codispersion of abundance and \n total basal area for",spe,gtitle)))
  
  # Significance and direction graph
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=diff))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Observed minus mean expected \n codispersion values for",spe,gtitle)))
  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=diff))+scale_fill_gradientn(colours=rainbow(7),limits=c(-1, 1))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Observed minus mean expected \n codispersion values for",spe,gtitle)))  
  
  # P-value graph
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Significance of observed Codispersion values \n for",spe,gtitle)))
  
  print( ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value))+scale_fill_gradientn(colours=c("#FF6666","#0000FF"),limits=c(0,0.5))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Significance of observed Codispersion values \n for",spe,gtitle)))
  
  # P-value category graph
  my.cols <- c("steelblue3","firebrick3")
  if(levels(out.df$P.value.cat)[1]=="Sig."){my.cols <- c("firebrick3") }
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value.cat))+scale_fill_manual(values=my.cols)+scale_colour_discrete(name="P.value.cat",limits=c(0,1))+coord_fixed(ratio=1)+t1+xlab(expression(k[1]))+ylab(expression(k[2]))+ggtitle(paste("Significance of observed Codispersion values \n for",spe,gtitle,"\n (alpha=5%, two sided, one sample t-test)")))  
}

# Comparing observed values to null model results from output arrays and graphing average null model results
# inputs are the null model input array object, the observed CoDisp result list, and a choice of null model
# This function deals with species' co-occurrences  
# 'spe' is text to identify the two species, e.g. spe="tsugca and querru"
comparison.abab.fn = function(null.input.ary,CoDisp_obs,model=c("HomP","HetP"),spe=spe,binwidth=binwidth){
  out.df <- CoDisp_obs[[1]] # extract the observed Codispersion result as a dataframe
  for(i in 1:length(null.input.ary[,1,1])){ # loop through each cell
    nsims <- length(null.input.ary[1,1,])
    obser <- out.df$Codispersion[i] # observed codispersion value
    expec <- null.input.ary[i,3,]
    prop.greater.than <- length(which(expec>obser))/nsims
    prop.less.than <- length(which(expec<obser))/nsims
    out.df$P.value[i]<-min(prop.greater.than,prop.less.than)
  }
  
  out.df$null_mean <- apply(null.input.ary[,3,],MARGIN=1,mean) # calculate mean codispersion value for each cell from the array of null model results
  out.df$diff <- out.df$Codispersion-out.df$null_mean # observed minus expected
  out.df$P.value.cat <- factor(ifelse(out.df$P.value<0.025,"Sig.","Non-sig.")) # significance at alpha=0.05
  
  ## Graph results
  if(model=="HomP"){gtitle="Homogeneous Poisson model"}
  if(model=="HetP"){gtitle="Heterogeneous Poisson model"}
  
  # observed graphs
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Codispersion of abundance of",spe,"observed")))
    
  print( ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Codispersion of abundance of",spe,"observed"))  )
  
  # null model graphs  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=null_mean))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Mean codispersion of abundance of",spe,gtitle)))
  
  print( ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=null_mean))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Mean codispersion of abundance of",spe,gtitle))  )
  
  # Significance and direction graph
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=diff))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Observed minus mean expected \n codispersion values for",spe,gtitle)))
  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=diff))+scale_fill_gradientn(colours=rainbow(7),limits=c(-1, 1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Observed minus mean expected \n codispersion values for",spe,gtitle)))  
  
  # P-value graph
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Significance of observed Codispersion values \n for",spe,gtitle)))
  
  print( ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value))+scale_fill_gradientn(colours=c("#FF6666","#0000FF"),limits=c(0,0.5))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Significance of observed Codispersion values \n for",spe,gtitle)))
    
  # P-value category graph
  my.cols <- c("steelblue3","firebrick3")
  if(levels(out.df$P.value.cat)[1]=="Sig."){my.cols <- c("firebrick3") }
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value.cat))+scale_fill_manual(values=my.cols)+scale_colour_discrete(name="P.value.cat",limits=c(0,1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Significance of observed Codispersion values \n for",spe,gtitle,"\n (alpha=5%, two sided, one sample t-test)")))  
}


# Comparing observed values to null model results from output arrays and graphing average null model results
# inputs are the null model input array object, the observed CoDisp result list, and a choice of null model
# This function deals with situations where the datasets being compared are specified by the user
# dat1 = text describing the first data object, e.g. "abundance of adults"
# dat2 = text describing the second data object
# dataset = text describing the dataset, e.g. "Plot 1"
comparison.fn = function(null.input.ary,CoDisp_obs,model=c("HomP","RLM","HetP"),dat1=dat1,dat2=dat2,dataset=dataset,binwidth=binwidth){
  out.df <- CoDisp_obs[[1]] # extract the observed Codispersion result as a dataframe
  for(i in 1:length(null.input.ary[,1,1])){ # loop through each cell
    nsims <- length(null.input.ary[1,1,])
    obser <- out.df$Codispersion[i] # observed codispersion value
    expec <- null.input.ary[i,3,]
    prop.greater.than <- length(which(expec>obser))/nsims
    prop.less.than <- length(which(expec<obser))/nsims
    out.df$P.value[i]<-min(prop.greater.than,prop.less.than)
    #t.test(x=expec,mu=obser,alternative="two.sided")$p.value # calculate probability for cell    
  }
  
  out.df$null_mean <- apply(null.input.ary[,3,],MARGIN=1,mean) # calculate mean codispersion value for each cell from the array of null model results
  out.df$diff <- out.df$Codispersion-out.df$null_mean # observed minus expected
  out.df$P.value.cat <- factor(ifelse(out.df$P.value<0.025,"Sig.","Non-sig.")) # significance at alpha=0.05
  
  ## Graph results
  if(model=="HomP"){gtitle="Homogeneous Poisson model"}
  if(model=="HetP"){gtitle="Heterogeneous Poisson model"}
  if(model=="RLM"){gtitle="Random labelling model"}
  
  # observed graphs
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Codispersion of",dat1,",",dat2,"observed")))
  
  print( ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Codispersion of",dat1,",",dat2,"observed"))  )
  
  # null model graphs  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=null_mean))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Mean codispersion of",dat1,",",dat2,gtitle)))
  
  print( ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=null_mean))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Mean codispersion of",dat1,",",dat2,gtitle))  )
  
  # Significance and direction graph
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=diff))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Observed minus mean expected \n codispersion values for",dat1,",",dat2,gtitle)))
  
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=diff))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1, 1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Observed minus mean expected \n codispersion values for",dat1,",",dat2,gtitle)))  
  
  # P-value graph
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Significance of observed Codispersion values \n for",dat1,",",dat2,gtitle)))
  
  print( ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value))+scale_fill_gradientn(colours=c("#FF6666","#0000FF"),limits=c(0,0.5))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Significance of observed Codispersion values \n for",dat1,",",dat2,gtitle))  )
  
  # P-value category graph
  my.cols <- c("steelblue3","firebrick3")
  if(levels(out.df$P.value.cat)[1]=="Sig."){my.cols <- c("firebrick3") }
  print(ggplot(out.df,aes(xx,yy))+geom_tile(aes(fill=P.value.cat))+scale_fill_manual(values=my.cols)+scale_colour_discrete(name="P.value.cat",limits=c(0,1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Significance of observed Codispersion values \n for",dat1,",",dat2,gtitle,"\n (alpha=5%, two sided, one sample t-test)")))  
}

# Function to return a data frame with the null model comparison results
codisp.compare <- function(null.input.ary,CoDisp_obs){
  out.df <- CoDisp_obs # observed Codispersion result df
  for(i in 1:length(null.input.ary[,1,1])) { # loop through each cell
    nsims <- length(null.input.ary[1,1,])
    obser <- out.df$Codispersion[i] # observed codispersion value
    expec <- null.input.ary[i,3,]
    prop.greater.than <- length(which(expec>obser))/nsims
    prop.less.than <- length(which(expec<obser))/nsims
    out.df$P.value[i]<-min(prop.greater.than,prop.less.than)
  } # end cell loop
  
  out.df$null_mean <- apply(null.input.ary[,3,],MARGIN=1,mean) # calculate mean codispersion value for each cell from the array of null model results
  out.df$diff <- out.df$Codispersion-out.df$null_mean # observed minus expected
  out.df$P.value.cat <- factor(ifelse(out.df$P.value<0.05,"Sig.","Non-sig.")) # significance at alpha=0.05
  return(out.df)
}

##################################
### GRAPHING
##################################

#### Graphing function for ViewPort Grid graphics
vplayout <- function(x,y) { viewport(layout.pos.row=x,layout.pos.col=y) }

## gglot options
t1<-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text = element_text(colour="black",size=20,angle=0),
  axis.title = element_text(colour="black",size=20),
  legend.key = element_blank(),
  legend.title = element_text(colour="black",size=14),
  legend.text = element_text(colour="black",size=14),
  plot.margin = unit(c(.2,.2,.2,.2),"lines"),
  panel.margin = unit(.2,"lines"),
  plot.background = element_rect(fill=NA)
  )

t1.no.leg_lab <-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text = element_text(colour="black",size=20,angle=0),
  axis.title = element_blank(),
  legend.position="none",
  plot.margin = unit(c(.5,.2,.2,.2),"lines"),
  panel.margin = unit(.2,"lines"),
  plot.background = element_rect(fill=NA)
  )

t1.no.leg_lab.18 <-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text = element_text(colour="black",size=18,angle=0),
  axis.title = element_blank(),
  legend.position="none",
  plot.margin = unit(c(.5,.2,.2,.2),"lines"),
  panel.margin = unit(.2,"lines"),
  plot.background = element_rect(fill=NA)
  )

t1.no.leg <-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text = element_text(colour="black",size=20,angle=0),
  axis.title = element_text(colour="black",size=20),
  legend.text = element_text(colour="black",size=18),
  legend.position="none",
  plot.margin = unit(c(.5,.2,.2,.2),"lines"),
  panel.margin = unit(.2,"lines"),
  plot.background = element_rect(fill=NA)
  #axis.title.x = element_blank(),
  #axis.title.y = element_blank()
)

t1.no.leg_lab.side <-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text.y = element_text(colour="black",size=14,angle=0),
  axis.text.x = element_text(colour="black",size=14,angle=45,hjust=1),
  axis.title = element_text(colour="black",size=14),
  legend.key = element_blank(),
  legend.title = element_blank(),
  legend.text = element_text(colour="black",size=14),
  plot.margin = unit(c(.5,.2,.2,.2),"lines"),
  panel.margin = unit(.2,"lines"),
  plot.background = element_rect(fill=NA)
  )

t1.no.lab <-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text = element_text(colour="black",size=20,angle=0),
  axis.title = element_text(colour="black",size=20),
  legend.key = element_blank(),
  legend.title = element_text(colour="black",size=20),
  legend.text = element_text(colour="black",size=20),
  plot.margin = unit(c(.2,.2,.2,.2),"lines"),
  panel.margin = unit(.2,"lines"),
  plot.background = element_rect(fill=NA),
  axis.title.x = element_blank(),
  axis.title.y = element_blank()
)

t1.no.lab.20pt <-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text = element_text(colour="black",size=20,angle=0),
  axis.title = element_text(colour="black",size=20),
  legend.key = element_blank(),
  legend.title = element_text(colour="black",size=20),
  legend.text = element_text(colour="black",size=20),
  plot.margin = unit(c(.2,.2,.2,.2),"lines"),
  panel.margin = unit(.2,"lines"),
  plot.background = element_rect(fill=NA),
  axis.title.x = element_blank(),
  axis.title.y = element_blank()
)

t1.fat.margins <-theme(                              
  plot.background = element_blank(), 
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  panel.border = element_blank(), 
  panel.background = element_blank(),
  axis.line = element_line(size=.4),
  axis.text = element_text(colour="black",size=20,angle=0),
  axis.title = element_text(colour="black",size=20),
  legend.key = element_blank(),
  legend.title = element_text(colour="black",size=14),
  legend.text = element_text(colour="black",size=14),
  plot.margin = unit(c(1,1,1,4),"lines"),  # top, right, bottom, and left 
  panel.margin = unit(c(1,1,1,1),"lines"),
  plot.background = element_rect(fill=NA)
  )

# Multiplot function for making graph panels
# e.g., matrix(c(1,2,3,3), nrow=2, byrow=TRUE)
  
  multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
    require(grid)
    
    # Make a list from the ... arguments and plotlist
    plots <- c(list(...), plotlist)
    
    numPlots = length(plots)
    
    # If layout is NULL, then use 'cols' to determine layout
    if (is.null(layout)) {
      # Make the panel
      # ncol: Number of columns of plots
      # nrow: Number of rows needed, calculated from # of cols
      layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                       ncol = cols, nrow = ceiling(numPlots/cols))
    }
    
    if (numPlots==1) {
      print(plots[[1]])
      
    } else {
      # Set up the page
      grid.newpage()
      pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
      
      # Make each plot, in the correct location
      for (i in 1:numPlots) {
        # Get the i,j matrix positions of the regions that contain this subplot
        matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
        
        print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                        layout.pos.col = matchidx$col))
      }
    }
  } 


### Plot a PCA biplot
# scale.factor - how much to scale the Vectors by
PCbiplot <- function(PC, x="PC1", y="PC2",site.labels=c("TRUE","FALSE"),scale.factor=scale.factor) {
    # PC being a prcomp object
    data <- data.frame(obsnames=1:length(PC$x), PC$x)
    
    if(site.labels=="TRUE"){plot <- ggplot(data, aes(x=PC1, y=PC2)) + geom_text(alpha=.4, size=3, aes(label=obsnames))}
    if(site.labels=="FALSE"){plot <- ggplot(data, aes(x=PC1, y=PC2)) + geom_point(alpha=.4, size=3)}
   
    plot <- plot + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2)
   
    datapc <- data.frame(varnames=rownames(PC$rotation), PC$rotation)
    
    mult <- min( (max(data[,"PC2"]) - min(data[,"PC2"])/(max(datapc[,"PC2"])-min(datapc[,"PC2"]))), (max(data[,"PC1"]) - min(data[,"PC1"])/(max(datapc[,"PC1"])-min(datapc[,"PC1"])))      )
    
    datapc <- transform(datapc, PC1 = scale.factor * mult * (get("PC1")), PC12 = scale.factor * mult * (get("PC2"))  )
    plot <- plot + coord_equal() + geom_text(data=datapc, aes(x=PC1, y=PC2, label=varnames), size = 5, vjust=1, color="red")
  
    print(plot <- plot + geom_segment(data=datapc, aes(x=0, y=0, xend=PC1, yend=PC2), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75, color="red") + t1)
}
#PCbiplot(pca1,site.labels="FALSE", scale.factor=2.5)
#fit <- prcomp(USArrests, scale=T)
#PCbiplot(fit)

## Function to generate variograms and cross variograms for the two geo.data objects used in codispersion analysis (observed patterns)

# labels is a two element vector used for labelling the graphs
# e.g. labels=c("species1","species2")

cross.variog.fn <- function(geodata1,geodata2,lab=missing(lab)){
  Obs_graphs <- vector(mode="list",length=3) # create empty object to store graphs
  D1.dat <- data.frame(X=geodata1$coords[,1],Y=geodata1$coords[,2],D1=geodata1$data) # put geodata object into a dataframe
  D2.dat <- data.frame(X=geodata2$coords[,1],Y=geodata2$coords[,2],D2=geodata2$data)
  # Plot the observed raster patterns
  g1 <- ggplot(D1.dat, aes(x=X, y=Y, size=D1))+geom_point(colour="black", fill="#4dac26", shape=21)+coord_fixed(ratio=1)
  g2 <- ggplot(D2.dat, aes(x=X, y=Y, size=D2))+geom_point(colour="black", fill="steelblue2", shape=21)+coord_fixed(ratio=1)

    ## Plot the variograms and cross variogram
    v.dat <- data.frame(geodata1$coords,dat1=scale(geodata1$data),dat2=scale(geodata2$data))
     g <- gstat(id="D1", formula=dat1~1, locations=~qx+qy, data = v.dat)    
     g <- gstat(g, id="D2", formula=dat2~1, locations=~qx+qy, data = v.dat)
     v <- variogram(g, cutoff=(min(xmax,ymax)*0.67), cross=TRUE) # 2/3 the min. of the two plot dimensions
    g3 <- ggplot(v,aes(x=dist,y=gamma,group=id,colour=id))+geom_line(lwd=2) + labs(x="Distance (m)",y = "Semivariance") 

    if(missing(lab)==FALSE){  # put labels on the graphs
    Obs_graphs[[1]] <- g1+t1.no.leg_lab+ggtitle(lab[1])
    Obs_graphs[[2]] <- g2+t1.no.leg_lab+ggtitle(lab[2])
    Obs_graphs[[3]] <- g3+t1.no.leg     
    }

    if(missing(lab)==TRUE){  # don't put a label on the legend
    Obs_graphs[[1]] <- g1+t1.no.leg_lab
    Obs_graphs[[2]] <- g2+t1.no.leg_lab
    Obs_graphs[[3]] <- g3+t1.no.leg     
    }
  
    return(Obs_graphs) 
  
  }  # end of function


# Function to print a codispersion graph using the CoDisp output object
print.CoDisp <- function(CoDisp.obj=CoDisp.obj,scaled=c("TRUE","FALSE"),contours=c("TRUE","FALSE"),binwidth=binwidth,input=input,gtitle=gtitle){
  
  if(scaled=="FALSE"){
#  print(ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab(expression(h[1]))+ylab(expression(h[2]))+ggtitle(paste("Codispersion of",input,gtitle)))
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")+ggtitle(paste("Codispersion of",input,gtitle))  
  }
  
  if(scaled=="TRUE"){
  if(contours=="TRUE"){
#    print(ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+t1+xlab(expression(h[1]))+ylab(expression(h[2]))
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+t1+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+ggtitle(paste("Codispersion of",input,gtitle))    
  }
  
  if(contours=="FALSE"){
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+t1+xlab(expression(h[1]))+ylab(expression(h[2]))+ggtitle(paste("Codispersion of",input,gtitle))    
  } 
  } # end of scaled

  return(g1)

  } # end of function

# Function to print a codispersion graph using the CoDisp output object
# With plain output (no labels)
print.CoDisp.plain <- function(CoDisp.obj=CoDisp.obj,scaled=TRUE,contours=TRUE,labels=TRUE,legend=TRUE,binwidth=binwidth){
  
  if(labels=="TRUE"){
  if(scaled=="FALSE"){
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1)+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)")  
  }
  
  if(scaled=="TRUE"){
  if(contours=="TRUE"){
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+xlab("Spatial lag in X (m)")+ylab("Spatial lag in Y (m)") + stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)    
  }
  
  if(contours=="FALSE"){
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+xlab(expression(h[1]))+ylab(expression(h[2]))    
  } 
  } # end of scaled
  }  # end of labels == TRUE
  
  
  if(labels=="FALSE"){
  if(scaled=="FALSE"){
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=rainbow(7))+coord_fixed(ratio=1) +xlab(NULL) +ylab(NULL)
  }
  
  if(scaled=="TRUE"){
  if(contours=="TRUE"){
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+ stat_contour(aes(x=xx,y=yy,z=Codispersion),binwidth=binwidth)+xlab(NULL) +ylab(NULL)    
  }
  
  if(contours=="FALSE"){
    g1 <- ggplot(CoDisp.obj,aes(xx,yy))+geom_tile(aes(fill=Codispersion))+scale_fill_gradientn(colours=c("#0000FF","#FFFFFF","#FF6666"),limits=c(-1,1))+coord_fixed(ratio=1)+xlab(NULL) +ylab(NULL)  
  } 
  } # end of scaled
  
  }  # end of labels == FALSE  
  
  if(legend=="TRUE")  {  g1 <- g1 + t1.no.leg_lab }
  if(legend=="FALSE") {  g1 <- g1 + t1.no.leg }

  return(g1)

  } # end of function

##################################
### SIMULATING PATTERNS
##################################


##############
##### Function to simulate co-occurrence patterns on grids

#   grid.points=20     # scale: grain of grid
#   sp1.pattern="CSR"  # The desired pattern for species 1: CSR,decreasing.x,increasing.x,decreasing.xy,bivariate.normal
#   sp2.pattern="CSR"  # The desired pattern for species 2: CSR,decreasing.x,increasing.x,decreasing.xy,bivariate.normal
#   xmin=0             # Dimensions of the plot area, e.g. 300 x 300m 
#   xmax=300
#   ymin=0
#   ymax=300
#   Print=TRUE         # Whether you want plots of the spatial patterns or not

copp.fn <- function(grid.points = grid.points,sp1.pattern = c("CSR","decreasing.x","increasing.x","decreasing.xy","increasing.xy","bivariate.normal"),sp2.pattern = c("CSR","decreasing.x","increasing.x","decreasing.xy","increasing.xy","bivariate.normal"),xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,Print=c("TRUE","FALSE")){
  
  # 1. Create an empty list to add output geodata objects
  copp.sim <- vector("list",2)
  
  # 2. Set up underlying grid coordinates
  X <- seq(from=xmin,to=xmax-grid.points,by=grid.points)
  Y <- seq(from=ymin,to=ymax-grid.points,by=grid.points)
  gridxy <- expand.grid(x=X,y=Y)
  
  # 3a. Create a set of quadrat abundance values for sp1 based on the 'sp1.pattern' argument
  
  if(sp1.pattern=="CSR"){Z <- rnorm(length(gridxy$x),mean=30,sd=10) }
  
  if(sp1.pattern=="decreasing.x"){Z <- 1+(rev(2*gridxy$x+5))/10 }
  
  if(sp1.pattern=="decreasing.x"){Z <- 1+(rev(2*gridxy$x+5))/10 }
  
  if(sp1.pattern=="increasing.x"){Z <- 1+(2*gridxy$x+5)/10 }
  
  if(sp1.pattern=="decreasing.xy"){  
    Z <- 1+rev(((gridxy$x+1)^2+(gridxy$y+1)^2)/3000) # (x-u)^2+(y-v)^2
  }
  
  if(sp1.pattern=="increasing.xy"){  
    Z <- 1+((gridxy$x+2)^2+(gridxy$y+1)^2)/3000 # (x-u)^2+(y-v)^2
  }
  
  if(sp1.pattern=="bivariate.normal"){ 
    Z <- 300*bivariate(((gridxy$x-min(gridxy$x))/(max(gridxy$x)-min(gridxy$x))*4)-2,((gridxy$y-min(gridxy$y))/(max(gridxy$y)-min(gridxy$y))*4)-2)  
  } #  bivariate normal

  # 3b. Add data to the output list as a geodata object
  copp.sp1.df <- data.frame(x=gridxy$x,y=gridxy$y,Z=jitter(Z,mean(Z)/5))
  copp.sim[[1]] <- as.geodata(copp.sp1.df,coords.col=1:2,data.col=3)
  
  # 4a. Create a set of quadrat abundance values for sp2 based on the 'sp1.pattern' argument
  
  if(sp2.pattern=="CSR"){Z <- rnorm(length(gridxy$x),mean=30,sd=10) }
  
  if(sp2.pattern=="decreasing.x"){Z <- 1+(rev(2*gridxy$x+5))/10 }
  
  if(sp2.pattern=="decreasing.x"){Z <- 1+(rev(2*gridxy$x+5))/10 }
  
  if(sp2.pattern=="increasing.x"){Z <- 1+(2*gridxy$x+5)/10 }
  
  if(sp2.pattern=="decreasing.xy"){  
    Z <- 1+rev(((gridxy$x+1)^2+(gridxy$y+1)^2)/3000) # (x-u)^2+(y-v)^2
  }
  
  if(sp2.pattern=="increasing.xy"){  
    Z <- 1+((gridxy$x+2)^2+(gridxy$y+1)^2)/3000 # (x-u)^2+(y-v)^2
  }
  
  if(sp2.pattern=="bivariate.normal"){ 
    Z <- 300*bivariate(((gridxy$x-min(gridxy$x))/(max(gridxy$x)-min(gridxy$x))*4)-2,((gridxy$y-min(gridxy$y))/(max(gridxy$y)-min(gridxy$y))*4)-2)  
  } #  bivariate normal
  
  # 4b. Add data to the output list as a geodata object
  copp.sp2.df <- data.frame(x=gridxy$x,y=gridxy$y,Z=jitter(Z,mean(Z)/10))
  copp.sim[[2]] <- as.geodata(copp.sp2.df,coords.col=1:2,data.col=3) 
  
  # 5. Print map of points if desired
  if(Print=="TRUE"){
    print(qplot(x, y, data=copp.sp1.df, size=Z,main=paste("sp1.pattern = ",sp1.pattern))+ theme_bw())
    print(qplot(x, y, data=copp.sp2.df, size=Z,main=paste("sp2.pattern = ",sp2.pattern))+ theme_bw())   
    } # end Print loop    
  
  # 6. Output the list of geodata objects
  return(copp.sim)  
  } # end function

##############
##### Function to simulate anisotropic point patterns

# # dimensions of the plot
# xmin=0
# xmax=200
# ymin=0
# ymax=200
# grid.points=5 # the distance between points on the underlying grid (must divide evenly into the plot dimensions)
# env.func = "unimodal"
# ppp.model="Thomas" # or "CSR"
# kappa=20
# sigma=0.5
# mu=10
# lambda=200
# pattern.method="abundance" # or "quant.marks"
# marks.method="unimodal"
# minmark=1
# maxmark=80
# sp.pattern= "random"  # "decreasing.x","increasing.x","decreasing.xy","increasing.xy","bivariate.normal" # the distribution pattern of the species
# sp.maxab=15
# ntrees = the number of trees you want 

app.sim.fn <- function(grid.points = grid.points,env.func = c("uniform","CSR","decreasing.x","increasing.x","decreasing.xy","increasing.xy","bivariate.normal"),pattern.method=c("quant.marks","abundance"),ppp.model=c("CSR","Thomas"),marks.method=c("random","decreasing.x","increasing.x","decreasing.xy","increasing.xy","bivariate.normal"),sp.pattern=c("random","decreasing.x","increasing.x","decreasing.xy","increasing.xy","bivariate.normal"),ntrees=ntrees,sp.maxab=sp.maxab,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,minmark=minmark,maxmark=maxmark,kappa=kappa,sigma=sigma,mu=mu,lambda=lambda,Print=c("TRUE","FALSE")){ # begin function
  
  # 1. Set up underlying grid coordinates
  X <- seq(from=xmin,to=xmax-grid.points,by=grid.points)
  Y <- seq(from=ymin,to=ymax-grid.points,by=grid.points)
  gridxy <- expand.grid(x=X,y=Y)
  
  # 2. Create a set of marks to use as values for the environmental variable based on the 'env.func' argument
 
  if(env.func=="uniform"){Z <- jitter(rep(50,(length(X)*length(Y)))) }
  
  if(env.func=="CSR"){Z <- rnorm(n=length(gridxy$x),mean=50,sd=15) }
  
  if(env.func=="decreasing.x"){Z <- 1+(rev(2*gridxy$x+5))/10}
  
  if(env.func=="increasing.x"){Z <- 1+(2*gridxy$x+5)/10}
  
  if(env.func=="decreasing.xy"){  
    Z <- 1+rev(((gridxy$x+1)^2+(gridxy$y+1)^2)/3000) # (x-u)^2+(y-v)^2
  }
  
  if(env.func=="increasing.xy"){  
    Z <- 1+((gridxy$x+2)^2+(gridxy$y+1)^2)/3000 # (x-u)^2+(y-v)^2
  }
  
 if(env.func=="bivariate.normal"){ 
   Z <- bivariate(((gridxy$x-min(gridxy$x))/(max(gridxy$x)-min(gridxy$x))*4)-2,((gridxy$y-min(gridxy$y))/(max(gridxy$y)-min(gridxy$y))*4)-2)  
   } #  bivariate normal

 epp.df <- data.frame(x=gridxy$x,y=gridxy$y,Z=Z)
 epp.sim <- as.ppp(epp.df,marks=Z,W=owin(c(xmin,xmax),c(ymin,ymax))) 
 #plot(epp.sim)
 
 # 3. Marked point pattern
 # 3a. Create a ppp of trees using the selected model 
if(pattern.method=="quant.marks"){
 if(ppp.model=="CSR"){ temp <- rpoispp(lambda=lambda,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))  
   mpp.sim <- temp[1:ntrees]  } # generate 2000 trees in the plot
 
 if(ppp.model=="Thomas"){ temp <- rThomas(kappa=kappa,sigma=sigma,mu=mu,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))
   mpp.sim <- temp[1:ntrees]  }# generate 2000 trees in the plt
 
 # 3b. Assign the marks to the point pattern using the selected method of codispersion
 if(marks.method=="random"){ mrks <- rtlnorm(mpp.sim$n,meanlog=log(maxmark/2),sdlog=log(maxmark/15),lower=minmark,upper=maxmark) } # generate a random set of marks drawn from a lognormal distribution 
 
 if(marks.method=="decreasing.x"){ temp <- -18*mpp.sim$x+5
     mrks <- temp+abs(min(temp)) }
 
 if(marks.method=="increasing.x"){ mrks <- 18*mpp.sim$x+5 }
 
 if(marks.method=="decreasing.xy"){  
   temp <- -((mpp.sim$x+1)^2+(mpp.sim$y+1)^2)
   mrks <- temp+abs(min(temp))   } # (x-u)^2+(y-v)^2
 
 if(marks.method=="increasing.xy"){  
   mrks <- ((mpp.sim$x+1)^2+(mpp.sim$y+1)^2) } # (x-u)^2+(y-v)^2
 
 if(marks.method=="bivariate.normal"){ mrks <- bivariate(((mpp.sim$x-min(mpp.sim$x))/(max(mpp.sim$x)-min(mpp.sim$x))*3)-2,((mpp.sim$y-min(mpp.sim$y))/(max(mpp.sim$y)-min(mpp.sim$y))*4)-2)  } #  bivariate.normal
 
  mrks1 <- (mrks-min(mrks))/(max(mrks)-min(mrks)) # scale marks between 0 and 1
  mpp.sim$marks <- minmark+mrks1/max(mrks1)*(maxmark-minmark) # spread marks between max and min mark 
  mpp.sim$window <- owin(c(xmin,xmax),c(ymin,ymax)) # rescale window to metres
  mpp.sim$x <- mpp.sim$x*100 # rescale x and y values to metres
  mpp.sim$y <- mpp.sim$y*100
      } # end quant.marks
 
  # 4. Generate a species point pattern using the selected method 
     # 4a. First generate a grid of values in the selected pattern.

  if(pattern.method=="abundance"){
    if(sp.pattern=="random"){ ab <- runif(n=length(c(gridxy$x)),min=0,max=sp.maxab) }
    
    if(sp.pattern=="decreasing.x"){ ab <- 1+(rev(2*gridxy$x+5))/10 }
    
    if(sp.pattern=="increasing.x"){ ab <- 1+(2*gridxy$x+5)/10 }
    
    if(sp.pattern=="decreasing.xy"){ ab <- 1+rev(((gridxy$x+1)^2+(gridxy$y+1)^2)/3000) } # (x-u)^2+(y-v)^2
    
    if(sp.pattern=="increasing.xy"){ ab <- 1+((gridxy$x+2)^2+(gridxy$y+1)^2)/3000 } # (x-u)^2+(y-v)^2
    
    if(sp.pattern=="bivariate.normal"){ 
      ab <- bivariate(((gridxy$x-min(gridxy$x))/(max(gridxy$x)-min(gridxy$x))*4)-2,((gridxy$y-min(gridxy$y))/(max(gridxy$y)-min(gridxy$y))*4)-2)  } #  bivariate normal

    AB <- round(ab/max(ab)*sp.maxab,0) # scale abundance to maximum number of individuals per grid cell   
    mpp.df <- data.frame(x=gridxy$x,y=gridxy$y,ab=AB)
    mpp.sim <- as.ppp(mpp.df,marks=ab,W=owin(c(xmin,xmax),c(ymin,ymax))) 
    
  } # end abundance loop

  # 5. Put both the environment ppp object and the species2 ppp object into an output list object
    app.sim <- vector("list")
    app.sim[[1]] <- epp.sim
    app.sim[[2]] <- mpp.sim

  # 6. Print map of points if desired
    if(Print=="TRUE"){
      par(mfrow=c(1,2))
      (plot(epp.sim,main=paste("env.func =",env.func),cex.main=0.7))
      if(pattern.method=="quant.marks"){
      (plot(mpp.sim,main=paste("mrks =",marks.method,mpp.sim$n),cex.main=0.7))  }

      if(pattern.method=="abundance"){
        (plot(mpp.sim,main=paste("species =",sp.pattern,mpp.sim$n),cex.main=0.7))  }
    } # end Print loop    

  return(app.sim)

  } # end function


#app.sim <- app.sim.fn(grid.points=5,env.func="increasing.x",pattern.method="abundance",sp.pattern="increasing.x",sp.maxab=20,xmin=0,xmax=200,ymin=0,ymax=200,Print="TRUE")
 

##############
##### Function to simulate marked point patterns

# kappa=20
# sigma=0.05
# mu=10
# xmin=0
# xmax=200
# ymin=0
# ymax=200
# nsim=1
# n=500
# 
# 
# minmark=1
# maxmark=80
# meanmark=30
# sdmark=5
# 
# ppp.model="rThomas"
# 
# 
# lambda=100
# neigh_radius=0 # all the trees are the same size (0.5 times the max diameter)
# neigh_radius=5
# neigh_radius=10
# neigh_radius=20 # large neighbourhoods mean small trees are near small trees and large trees are very isolated
# Print=TRUE

mpp.sim.fn <- function(ppp.model=c("HomP","rThomas"),marks.method=c("random","inv_neigh","prop_neigh"),xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,minmark=minmark,maxmark=maxmark,meanmark=meanmark,sdmark=sdmark,kappa=kappa,sigma=sigma,mu=mu,lambda=lambda,neigh_radius=10,Print=c("TRUE","FALSE")){ # begin function

  # 1. Generate the point pattern either by CSR or a Thomas process
  if(ppp.model=="HomP"){ mpp.sim <- rpoispp(lambda=lambda,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))  } 
  
  if(ppp.model=="rThomas"){ mpp.sim <- rThomas(kappa=kappa,sigma=sigma,mu=mu,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100))) } 
  
  # 2. Generate the marks via a truncated log-normal distribution
  mpp.sim$marks <- rtlnorm(mpp.sim$n,meanlog=log(meanmark),sdlog=log(sdmark),lower=minmark,upper=maxmark) # if marks.method = "random"
  
  # 3. Calculate the number of nearest neighbours
  nneigh <- applynbd(mpp.sim,R=neigh_radius/100,FUN=function(Y,...){(Y$n-1)}) # count the number of neighbours within 'neigh_radius' distance of each point (distance is converted to decimeters)
  
  # 4. Assign marks to the point pattern under the selected 'marks.method' rule  
  if(marks.method=="inv_neigh"){ # dbh is the inverse of the number of neighbours within distance 'neigh_radius_dist' plus one
    mpp.sim$marks <- (1/(nneigh+1))*maxmark
  }

  if(marks.method=="prop_neigh"){ # dbh is proportional to the number of neighbours within distance 'neigh_radius_dist' plus one
    mpp.sim$marks <- (nneigh+1)/(max(nneigh)+1)*maxmark
  }
  
  # 5. Create the point pattern in a window scaled in metres
  mpp.sim$window <- owin(c(xmin,xmax),c(ymin,ymax))
  mpp.sim$x <- mpp.sim$x*100
  mpp.sim$y <- mpp.sim$y*100
  
  # 6. Print map of points if desired
  if(Print=="TRUE"){(plot(mpp.sim,main=paste("ppp.model =",ppp.model,"\n", "marks distribution = trunc.log.normal","\n marks model =",marks.method,"\n N =",mpp.sim$n),cex.main=0.7))}

  return(mpp.sim)
  
  } # end function

#plot(mpp.sim)

##############
##### Function to simulate paired point patterns

# xmin=0
# xmax=200
# ymin=0
# ymax=200
# 
# n=500
# 
# kappa=20
# sigma=0.05
# mu=10
# 
# sp1.pattern="CSR"
# sp1.pattern="Thomas"
# 
# association="independent.sp2.random"
# association="independent.sp2.Thomas"
# association="aggregated"
# association="segregated"
# 
# pattern.scale=5 # this is the distance away that neighbours must be either segregated or aggregated
# 
# par(mfrow=c(1,2))
# plot(ppp.pairs.sim[[1]])
# plot(ppp.pairs.sim[[2]])


ppp.pairs.sim.fn <- function(n,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,sp1.pattern=c("CSR","Thomas"),kappa=kappa,sigma=sigma,mu=mu,association=c("independent.sp2.random","independent.sp2.Thomas","aggregated","segregated"),pattern.scale=pattern.scale,Print=c("TRUE","FALSE")){
  
  # 1. Make an empty list to hold the two resultant point patterns
  ppp.pairs.sim <- vector("list",2) # create an empty list to hold simulated ppp objects
  
  # 2. Generate the first species' point pattern either by CSR or a Thomas process
  if(sp1.pattern=="CSR"){
    ppp.pairs.sim[[1]] <- rpoint(n,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))  
    pp <- rpoint(n*100,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100))) 
  }
  
  if(sp1.pattern=="Thomas"){
    x <- rThomas(kappa=kappa,sigma=sigma,mu=mu,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))
    z <- n * density(x, sigma=sigma)
    ppp.pairs.sim[[1]] <- rpoint(n, z,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100))) 
    x <- rThomas(kappa=kappa,sigma=sigma,mu=mu,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))
    z <- n * density(x, sigma=sigma)
    pp <- rpoint(n*100, z,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100))) 
    } # end sp1.pattern==Thomas
  
  # 3. Generate the second species' point pattern according to the association rule selected by the user
  if(association=="independent.sp2.random"){ppp.pairs.sim[[2]] <- rpoint(n,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))  }
    
  if(association=="independent.sp2.Thomas"){
    x <- rThomas(kappa=kappa,sigma=sigma,mu=mu,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))
    z <- n * density(x, sigma=sigma)
    ppp.pairs.sim[[2]] <- rpoint(n, z,win=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100))) 
     } 
      
  if(association=="aggregated"){
    nndists.rows <- sample(which(nncross(pp,ppp.pairs.sim[[1]],what="dist")<(pattern.scale/100)),n)
    new.x <- pp$x[nndists.rows]
    new.y <- pp$y[nndists.rows]
    ppp.pairs.sim[[2]] <- as.ppp(cbind(new.x,new.y),W=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100))) 
      } # end aggregated associaton 
    
  if(association=="segregated"){   
    nndists.rows <- sample(which(nncross(pp,ppp.pairs.sim[[1]],what="dist")>(pattern.scale/100)),n)
    new.x <- pp$x[nndists.rows]
    new.y <- pp$y[nndists.rows]
    ppp.pairs.sim[[2]] <- as.ppp(cbind(new.x,new.y),W=owin(c(xmin/100,xmax/100),c(ymin/100,ymax/100)))    
  } # end segregated associaton 
    
  # 4. Print map of points if desired
  if(Print=="TRUE"){
    par(mfrow=c(1,2))
    plot(ppp.pairs.sim[[1]],main=paste("sp 1\n pattern =",sp1.pattern),cex.main=0.7)
    plot(ppp.pairs.sim[[2]],main=paste("sp 2\n association =\n",association),cex.main=0.7) 
    
  # 5. Rescale the point patterns in metres
  for(i in 1:2){
    ppp.pairs.sim[[i]]$window <- owin(c(xmin,xmax),c(ymin,ymax))
    ppp.pairs.sim[[i]]$x <- ppp.pairs.sim[[i]]$x*100
    ppp.pairs.sim[[i]]$y <- ppp.pairs.sim[[i]]$y*100
  }
  
        
  }
  return(ppp.pairs.sim)
  
} # end function