## Written by Olivier Broennimann. Departement of Ecology and Evolution (DEE). 
## University of Lausanne. Switzerland. October 09.
##
## DESCRIPTION
##
## functions to perform measures of niche overlap and niche equivalency/similarity tests as described in Broennimann et al. (submitted)
## 
## list of functions:
##
## grid.clim(glob,glob1,sp,R) 
## use the scores of an ordination (or SDM predictions) and create a grid z of RxR pixels 
## (or a vector of R pixels when using scores of dimension 1 or SDM predictions) with occurrence densities
## Only scores of one, or two dimensions can be used 
## sp= scores for the occurrences of the species in the ordination, glob = scores for the whole studies areas, glob 1 = scores for the range of sp 
##
## niche.overlap(z1,z2,cor)
## calculate the overlap metrics D and I (see Warren et al 2008) based on two species occurrence density grids z1 and z2 created by grid.clim
## cor=T correct occurrence densities of each species by the prevalence of the environments in their range
##
## niche.equivalency.test(z1,z2,rep)
## runs niche equivalency test(see Warren et al 2008) based on two species occurrence density grids
## compares the observed niche overlap between z1 and z2 to overlaps between random niches z1.sim and z2.sim.
## z1.sim and z2.sim are built from random reallocations of occurences of z1 and z2
## rep is the number of iterations
##
## niche.similarity.test(z1,z2,rep)
## runs niche similarity test(see Warren et al 2008) based on two species occurrence density grids
## compares the observed niche overlap between z1 and z2 to overlaps between z1 and random niches (z2.sim) in available in the range of z2 (z2$Z) 
## z2.sim have the same patterns as z2 but their center are randomly translatated in the availabe z2$Z space and weighted by z2$Z densities
## rep is the number of iterations
##
## plot.niche(z,title,name.axis1,name.axis2)
## plot a niche z created by grid.clim. title,name.axis1 and name.axis2 are strings for the legend of the plot
##
## plot.contrib(contrib,eigen)
## plot the contribution of the initial variables to the analysis. Typically the eigen vectors and eigen values in ordinations
##
## plot.overlap.test(x,type,title)
## plot an histogram of observed and randomly simulated overlaps, with p-values of equivalency and similarity tests. 
## x must be an object created by niche.similarity.test or niche.equivalency.test.
## type is either "D" or "I". title is the title of the plot

##################################################################################################

grid.clim<-function(glob,glob1,sp,R){

# glob: global background dataset for the whole study area, 
# glob1: background for sp1
# sp: occurrence dataset
# R: resolution of the grid
l<-list()

if (ncol(glob)>2) stop("cannot calculate overlap with more than two axes")

if(ncol(glob)==1){ 											#if scores in one dimension (e.g. LDA,SDM predictions,...)
	xmax<-max(glob[,1])
	xmin<-min(glob[,1])
	sp.dens<-density(sp[,1],kernel="gaussian",from=xmin,to=xmax,n=R,cut=0) 		# calculate the density of occurrences in a vector of R pixels along the score gradient
														# using a gaussian kernel density function, with R bins.
	glob1.dens<-density(glob1[,1],kernel="gaussian",from=xmin,to=xmax,n=R,cut=0)	# calculate the density of environments in glob1
	x<-sp.dens$x 											# breaks on score gradient
	z<-sp.dens$y*nrow(sp)/sum(sp.dens$y) 							# rescale density to the number of occurrences in sp
							 							# number of occurrence/pixel
	Z<-glob1.dens$y*nrow(glob)/sum(glob1.dens$y) 						# rescale density to the number of sites in glob1
	z[z<max(z)/1000]<-0 										# remove infinitesimally small number generated by kernel density function
	Z[Z<max(Z)/1000]<-0 										# remove infinitesimally small number generated by kernel density function

	z.uncor<-z/max(z)											# rescale between [0:1] for comparison with other species
	z<-z/Z												# correct for environment prevalence 
	z[is.na(z)]<-0 											# remove n/0 situations
	z[z=="Inf"]<-0 											# remove 0/0 situations
	z.cor<-z/max(z)											# rescale between [0:1] for comparison with other species
	l$x<-x;l$z.uncor<-z.uncor;l$z.cor<-z.cor;l$Z<-Z;l$glob<-glob;l$glob1<-glob1;l$sp<-sp 
}

if(ncol(glob)==2){ #if scores in two dimensions (e.g. PCA)					


	library(adehabitat)
      xmin<-min(glob[,1]);xmax<-max(glob[,1]);ymin<-min(glob[,2]);ymax<-max(glob[,2])			# data preparation	
	glob1r<-data.frame(cbind((glob1[,1]-xmin)/abs(xmax-xmin),(glob1[,2]-ymin)/abs(ymax-ymin)))	# data preparation	
	spr<-data.frame(cbind((sp[,1]-xmin)/abs(xmax-xmin),(sp[,2]-ymin)/abs(ymax-ymin))) 			# data preparation
	mask<-ascgen(cbind((1:R)/R,(1:R)/R),nrcol=R,count=F)								# data preparation
	sp.dens<-kernelUD(spr[,1:2],grid=mask,kern="bivnorm")					# calculate the density of occurrences in a grid of RxR pixels along the score gradients
														# using a gaussian kernel density function, with RxR bins.
	sp.dens<-asc2spixdf(sp.dens[[1]]$UD)							# data manipulation																		
	glob1.dens<-kernelUD(glob1r[,1:2],grid=mask,kern="bivnorm")
	glob1.dens<-asc2spixdf(glob1.dens[[1]]$UD)
	x<-seq(from=min(glob[,1]),to=max(glob[,1]),length.out=R)				# breaks on score gradient 1 
	y<-seq(from=min(glob[,2]),to=max(glob[,2]),length.out=R)				# breaks on score gradient 2 
	z<-matrix(sp.dens$var*nrow(sp)/sum(sp.dens$var),nrow=R,ncol=R,byrow=F) 			#rescale density to the number of occurrences in sp
	Z<-matrix(glob1.dens$var*nrow(glob1)/sum(glob1.dens$var),nrow=R,ncol=R,byrow=F) 	#rescale density to the number of sites in glob1
	z[z<max(z)/1000]<-0 											# remove infinitesimally small number generated by kernel density function
	Z[Z<max(Z)/1000]<-0 											# remove infinitesimally small number generated by kernel density function

	z.uncor<-z/max(z)											# rescale between [0:1] for comparison with other species	
	z<-z/Z												# correct for environment prevalence
	z[is.na(z)]<-0											# remove n/0 situations
	z[z=="Inf"]<-0											# remove n/0 situations
	z.cor<-z/max(z)											# rescale between [0:1] for comparison with other species	
	l$x<-x;l$y<-y;l$z.uncor<-z.uncor;l$z.cor<-z.cor;l$Z<-Z;l$glob<-glob;l$glob1<-glob1;l$sp<-sp
}
return(l)
}

##################################################################################################
niche.overlap<-function(z1,z2,cor){ 

# z1 = species 1 occurrence density grid created by grid.clim
# z2 = species 2 occurrence density grid created by grid.clim
# cor=T correct occurrence densities of each species by the prevalence of the environments in their range

l<-list()

if(cor==F){		
p1<-z1$z.uncor/sum(z1$z.uncor) # rescale occurence densities so that the sum of densities is the same for both species
p2<-z2$z.uncor/sum(z2$z.uncor) # rescale occurence densities so that the sum of densities is the same for both species
}

if(cor==T){
p1<-z1$z.cor/sum(z1$z.cor)	# rescale occurence densities so that the sum of densities is the same for both species
p2<-z2$z.cor/sum(z2$z.cor)	# rescale occurence densities so that the sum of densities is the same for both species
}

D <-1-(0.5*(sum(abs(p1-p2))))				# overlap metric D
I <-1-(0.5*(sqrt(sum((sqrt(p1)-sqrt(p2))^2))))	# overlap metric I
l$D<-D
l$I<-I
return(l)
}

##################################################################################################

niche.equivalency.test<-function(z1,z2,rep){

R<-length(z1$x)
l<-list()
x11(2,2,pointsize = 12); par(mar=c(0,0,0,0));

obs.o<-niche.overlap(z1,z2,cor=T)									#observed niche overlap
sim.o<-data.frame(matrix(nrow=rep,ncol=2))							#empty list of random niche overlap
names(sim.o)<-c("D","I")
for (i in 1:rep){
	plot.new(); text(0.5,0.5,paste("runs to go:",rep-i+1))
	
	if(is.null(z1$y)){ #overlap on one axis
	
		occ1.sim<-sample(z1$x,size=nrow(z1$sp),replace=T,prob=z1$z.cor) 		#random sampling of occurrences following the corrected densities distribution
		occ2.sim<-sample(z2$x,size=nrow(z2$sp),replace=T,prob=z2$z.cor)
	
		occ.pool<-c(occ1.sim,occ2.sim)   # pool of random occurrences
		rand.row<-sample(1:length(occ.pool),length(occ1.sim),replace=T) 		# random reallocation of occurrences to datasets
		sp1.sim<-occ.pool[rand.row]
		sp2.sim<-occ.pool[-rand.row]
	
		z1.sim<-grid.clim(z1$glob,z1$glob1,data.frame(sp1.sim),R) # gridding
		z2.sim<-grid.clim(z2$glob,z2$glob1,data.frame(sp2.sim),R)	
	}

	if(!is.null(z1$y)){ 										#overlap on two axes
		coordinates<-which(z1$z.cor>0,arr.ind=T)						# array of cell coordinates ((1,1),(1,2)...)
		weight<-z1$z.cor[z1$z.cor>0] 								#densities in the same format as cells
		coordinates.sim1<-coordinates[sample(1:nrow(coordinates),size=nrow(z1$sp),replace=T,prob=weight),] #random sampling of coordinates following z1$z.cor distribution
		occ1.sim<-cbind(z1$x[coordinates.sim1[,1]],z1$y[coordinates.sim1[,2]]) 	# random occurrences following the corrected densities distribution
	
		coordinates<-which(z2$z.cor>0,arr.ind=T)						# array of cell coordinates ((1,1),(1,2)...)
		weight<-z2$z.cor[z2$z.cor>0] 								#densities in the same format as cells
		coordinates.sim2<-coordinates[sample(1:nrow(coordinates),size=nrow(z2$sp),replace=T,prob=weight),] #random sampling of coordinates following z1$z.cor distribution
		occ2.sim<-cbind(z2$x[coordinates.sim2[,1]],z2$y[coordinates.sim2[,2]]) 	# random occurrences following the corrected densities distribution
	
		occ.pool<-rbind(occ1.sim,occ2.sim) # pool of random occurrences
		rand.row<-sample(1:nrow(occ.pool),nrow(occ1.sim),replace=T) 			# random reallocation of occurrences to datasets
		sp1.sim<-occ.pool[rand.row,]
		sp2.sim<-occ.pool[-rand.row,]
	
		z1.sim<-grid.clim(z1$glob,z1$glob1,data.frame(sp1.sim),R)
		z2.sim<-grid.clim(z2$glob,z2$glob1,data.frame(sp2.sim),R)	
	}

	o.i<-niche.overlap(z1.sim,z2.sim,cor=F)							# overlap between random and observed niches
	sim.o$D[i]<-o.i$D											# storage of overlaps
	sim.o$I[i]<-o.i$I
}

dev.off()
l$sim<-sim.o	# storage
l$obs<-obs.o	# storage
l$p.D<- min((sum(sim.o$D <= obs.o$D ) + 1),(sum(sim.o$D >= obs.o$D ) + 1))*2/(length(sim.o$D) + 1) 	# storage of p-values
l$p.I<- min((sum(sim.o$I <= obs.o$I ) + 1),(sum(sim.o$I >= obs.o$I ) + 1))*2/(length(sim.o$I) + 1)	# storage of p-values

return(l)
}

##################################################################################################

niche.similarity.test<-function(z1,z2,rep){

R<-length(z1$x)
x11(2,2,pointsize = 12); par(mar=c(0,0,0,0));						# countdown window
l<-list()
obs.o<-niche.overlap(z1,z2,cor=T) 								#observed niche overlap
sim.o<-data.frame(matrix(nrow=rep,ncol=2))						#empty list of random niche overlap
names(sim.o)<-c("D","I")

for (k in 1:rep){
	plot.new(); text(0.5,0.5,paste("similarity tests:","\n","runs to go:",rep-k+1))	# countdown
	
	if(is.null(z2$y)){
		center<-which(z2$z.cor==1,arr.ind=T) 					# define the centroid of the observed niche
		Z<-z2$Z/max(z2$Z)
		rand.center<-sample(1:R,size=1,replace=F,prob=Z)				# randomly (weighted by environment prevalence) define the new centroid for the niche
		
		xshift<-rand.center-center							# shift on x axis
		z2.sim<-z2
		z2.sim$z.cor<-rep(0,R)								# set intial densities to 0
		for(i in 1:R){
			i.trans<-i+xshift
			if(i.trans>R|i.trans<0)next()						# densities falling out of the env space are not considered
			z2.sim$z.cor[i.trans]<-z2$z.cor[i]					# shift of pixels
		}
		z2.sim$z.cor<-(z2$Z!=0)*1*z2.sim$z.cor 					# remove densities out of existing environments
	}

	if(!is.null(z2$y)){
		centroid<-which(z2$z.cor==1,arr.ind=T)[1,] 				# define the centroid of the observed niche
		Z<-z2$Z/max(z2$Z)
		rand.centroids<-which(Z>0,arr.ind=T)					# all pixels with existing environments in the study area
		weight<-Z[Z>0]
		rand.centroid<-rand.centroids[sample(1:nrow(rand.centroids)
			,size=1,replace=F,prob=weight),]					# randomly (weighted by environment prevalence) define the new centroid for the niche
		xshift<-rand.centroid[1]-centroid[1]					# shift on x axis
		yshift<-rand.centroid[2]-centroid[2]					# shift on y axis
		z2.sim<-z2
		z2.sim$z.cor<-matrix(rep(0,R*R),ncol=R,nrow=R)				# set intial densities to 0
		for(i in 1:R){
			for(j in 1:R){
				i.trans<-i+xshift
				j.trans<-j+yshift
				if(i.trans>R|i.trans<0)next()					# densities falling out of the env space are not considered
				if(j.trans>R|j.trans<0)next()
				z2.sim$z.cor[i.trans,j.trans]<-z2$z.cor[i,j]		# shift of pixels
			}
		}
		z2.sim$z.cor<-(z2$Z!=0)*1*z2.sim$z.cor 					# remove densities out of existing environments
	}
	o.i<-niche.overlap(z1,z2.sim,cor=T)							# overlap between random and observed niches
	sim.o$D[k]<-o.i$D										# storage of overlaps
	sim.o$I[k]<-o.i$I
}
dev.off()
l$sim<-sim.o											# storage
l$obs<-obs.o											# storage
l$p.D<- min((sum(sim.o$D <= obs.o$D ) + 1),(sum(sim.o$D >= obs.o$D ) + 1))*2/(length(sim.o$D) + 1)	# storage of p-values
l$p.I<- min((sum(sim.o$I <= obs.o$I ) + 1),(sum(sim.o$I >= obs.o$I ) + 1))*2/(length(sim.o$I) + 1)	# storage of p-values

return(l)
}

##################################################################################################

plot.niche<-function(z,title,name.axis1,name.axis2,cor=F){

	if(is.null(z$y)){
		R<-length(z$x)
		x<-z$x
		xx<-sort(rep(1:length(x),2))

		if(cor==F)y1<-z$z.uncor/max(z$z.uncor)
		if(cor==T)y1<-z$z.cor/max(z$z.cor)
		Y1<-z$Z/max(z$Z)
		yy1<-sort(rep(1:length(y1),2))[-c(1:2,length(y1)*2)]
		YY1<-sort(rep(1:length(Y1),2))[-c(1:2,length(Y1)*2)]

		plot(x,y1,type="n",xlab=name.axis1,ylab="density of occurrence")
		polygon(x[xx],c(0,y1[yy1],0,0),col="grey")
		lines(x[xx],c(0,Y1[YY1],0,0))
	}

	if(!is.null(z$y)){
		if(cor==F)image(z$x,z$y,z$z.uncor,col = gray(100:0 / 100),zlim=c(0.000001,max(z$z.uncor)),xlab=name.axis1,ylab=name.axis2)
		if(cor==T)image(z$x,z$y,z$z.cor,col = gray(100:0 / 100),zlim=c(0.000001,max(z$z.cor)),xlab=name.axis1,ylab=name.axis2)
		contour(z$x,z$y,z$Z,add=T,levels=quantile(z$Z[z$Z>0],c(0,0.5)),drawlabels=F,lty=c(1,2))
	}
	title(title)
}

plot.contrib<-function(contrib,eigen){

	if(ncol(contrib)==1){
				h<-c(unlist(contrib))
				n<-row.names(contrib)
				barplot(h,space=0,names.arg=n)
				title(main="variable contribution")
		}
	if(ncol(contrib)==2){
					s.corcircle(contrib[,1:2]/max(abs(contrib[,1:2])),grid = F)
					title(main="correlation circle", sub=paste("axis1 = ",round(eigen[1]/sum(eigen)*100,2),"%","axis2 = ",round(eigen[2]/sum(eigen)*100,2),"%"))
		}
}


plot.overlap.test<-function (x,type,title) 
{
	if(type=="D") {
		obs <- x$obs$D
		sim <- x$sim$D
		p<-x$p.D
	}
	if(type=="I") {
		obs <- x$obs$I
		sim <- x$sim$I
		p<-x$p.I
	}
	r0 <- c(sim, obs)
    	l0 <- max(sim) - min(sim)
    	w0 <- l0/(log(length(sim), base = 2) + 1)
    	xlim0 <- range(r0) + c(-w0, w0)
    	h0 <- hist(sim, plot = FALSE, nclass = 10)
    	y0 <- max(h0$counts)
    	hist(sim, plot = TRUE, nclass = 10, xlim = xlim0, col = grey(0.8), main= title, xlab=type,sub = paste("p.value = ",round(p,5)))
    	lines(c(obs, obs), c(y0/2, 0),col="red")
    	points(obs, y0/2, pch = 18, cex = 2,col="red")
    	invisible()
}

##################################################################################################