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

# run this code below all together - this is the sum of variances function from Hughes et al 2013 PNAS
# disparity functions for Sum of Variances and Sum of Ranges, use all PCO axes

SoV <- function(X) sum(diag(cov(X)))
SoR <- function(X) sum(apply(X,2,range)[2,]-apply(X,2,range)[1,])

# Disparity analysis using SoV and SoR
DtT <- function(MS,PA,boo) {
	skewness.b.sple <- c()			# the skewness of the distribution of the bootstrap resample at everybin
	rar <- min(apply(PA,2,sum))			# the diversity value that the Sum of Ranges is rarefied to
	if (rar==1) rar <- 2
	Variance <- c() ; Range <- c()
	for (tps in 1:ncol(PA)) {			# a loop for computing disparity for each time interval
		occ <- which(PA[,tps]==1)
		if (length(occ)==1) {occ <- c(occ,occ) ; rar <- 2 }
		vari <- c() ; ran <- c()
		for (i in 1:boo) {
			b.sple <- sample(occ,length(occ),replace=T)
			r.sple <- sample(occ,rar,replace=F)
			vari <- c(vari,SoV(MS[b.sple,]))		# records disparity values of boostrap replicates
			ran <- c(ran,SoR(MS[r.sple,]))		# same but for range-based estimates
			}
		outvar <- which(vari==0) ; outran <- which(ran==0)
		ifelse(length(outvar!=0),Variance <- cbind(Variance,rbind(mean(vari[-c(outvar)]),sd(vari[-c(outvar)]))),Variance <- cbind(Variance,rbind(mean(vari),sd(vari)))) # Bootstrap estimates and errors
		ifelse(length(outran!=0),Range <- cbind(Range,rbind(mean(ran[-c(outran)]),sd(ran[-c(outran)]))),Range <- cbind(Range,rbind(mean(ran),sd(ran))))       # id.
		skewness.b.sple <- c(skewness.b.sple,skewness(vari))
		}
	list(Variance=Variance, Range=Range,skewness.b.sple=skewness.b.sple)
	}

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

## run this 
upperTriangle<-function (x, diag = FALSE, byrow = FALSE) 
{
    if (byrow) 
        t(x)[rev(upper.tri(x, diag = diag))]
    else x[upper.tri(x, diag = diag)]
}




## run this 
  cores <- detectCores()
  cl <- makeCluster(cores - 1, outfile = "")
  registerDoParallel(cl)




## run this 

 bootstrapWMPD <- function(dissim, compars) {
  
  dissim <- upperTriangle(dissim)
  compars <- upperTriangle(compars)
  dat <- as.data.frame(cbind(dissim,compars))
  
  weighted.mean <- sum(dissim * compars, na.rm=T) / sum(compars, na.rm=T)
  
  Z <- length(dissim[complete.cases(dissim)])
  
  boot.mean <- vector()
  
  for(i in 1:10000) {
    
    temp.dat <- dat[complete.cases(dissim),][sample.int(Z,Z,replace=T),]
        
    boot.mean[i] <- sum(temp.dat$dissim * temp.dat$compars, na.rm=T) / sum(temp.dat$compars, na.rm=T)
    
  }
  
#   weighted.mean <- mean(boot.mean, na.rm=T)
    
  #Lower 0.05 for the mean
  lower <- sort(boot.mean)[length(boot.mean)*0.05]
  
  #Upper 0.95 for the mean
  upper <- sort(boot.mean)[length(boot.mean)*0.95]
  
  return(cbind(weighted.mean,lower,upper))
  
}


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



MatrixPruner <- function(clad.matrix, taxa2prune = c(), characters2prune = c()) {

# ADD OPTION TO REMOVE CONSTANT CHARACTERS

    # Check that something to prune has been specified:
    if(is.null(taxa2prune) && is.null(characters2prune)) stop("No taxa or characters to prune specified.")

    # Check taxa to prune have names matching those in matrix:
    if(length(setdiff(taxa2prune, rownames(clad.matrix$matrix))) > 0) stop("Taxa specified that are not found in the matrix. Check spelling.")

    # Check character numbers are present in matrix:
    if(length(setdiff(characters2prune, c(1:ncol(clad.matrix$matrix))))) stop("Character numbers specified that are not found in the matrix.")

    # If there are taxa to prune, then prune them:
    if(!is.null(taxa2prune)) clad.matrix$matrix <- clad.matrix$matrix[-match(taxa2prune, rownames(clad.matrix$matrix)), , drop = FALSE]

    # If there are characters to prune:
    if(!is.null(characters2prune)) {
        
        # Remove from matrix
        clad.matrix$matrix <- clad.matrix$matrix[, -characters2prune, drop = FALSE]
        
        # Remove from ordering:
        clad.matrix$ordering <- clad.matrix$ordering[-characters2prune]
        
        # Remove from weights:
        clad.matrix$weights <- clad.matrix$weights[-characters2prune]
        
        # Remove from maximum values:
        clad.matrix$max.vals <- clad.matrix$max.vals[-characters2prune]
        
        # Remove from minimum values:
        clad.matrix$min.vals <- clad.matrix$min.vals[-characters2prune]
        
    }
    
    # Get unique values for each character:
    unique.values <- lapply(lapply(lapply(lapply(lapply(lapply(lapply(apply(clad.matrix$matrix, 2, list), unlist), as.character), strsplit, split = "&"), unlist), sort), unique), as.numeric)
    
    # If any character is now all missing data:
    if(length(which(unlist(lapply(unique.values, length)) == 0)) > 0) {
        
        # For each such character insert a single zero:
        for(i in which(unlist(lapply(unique.values, length)) == 0)) unique.values[[i]] <- c(0)
        
    }
    
    # Update maximum values (post pruning):
    clad.matrix$max.vals <- unlist(lapply(unique.values, max))
    
    # Update minimum values (post pruning):
    clad.matrix$min.vals <- unlist(lapply(unique.values, min))

    # Return pruned matrix:
    return(clad.matrix)

}

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

pairwise.adonis <- function(x, factors, sim.method = 'euclidean', p.adjust.m ='fdr')
{
library(vegan)
co = combn(unique(factors),2)
pairs = c()
F.Model =c()
R2 = c()
p.value = c()

for(elem in 1:ncol(co)){
ad = adonis(x[factors %in% c(co[1,elem],co[2,elem]),] ~ factors[factors %in% c(co[1,elem],co[2,elem])] , method =sim.method);
pairs = c(pairs,paste(co[1,elem],'vs',co[2,elem]));
F.Model =c(F.Model,ad$aov.tab[1,4]);
R2 = c(R2,ad$aov.tab[1,5]);
p.value = c(p.value,ad$aov.tab[1,6])
}
p.adjusted = p.adjust(p.value,method=p.adjust.m)
pairw.res = data.frame(pairs,F.Model,R2,p.value,p.adjusted)
return(pairw.res)
}

#######################################################################################
#################################################################
# RUN THIS FUNCTION ALL AT ONCE - CAN DRAW HULLS AROUND GROUPS

Plot_ConvexHull<-function(xcoord, ycoord, lcolor, bgcolor){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])
  lines(xcoord[hpts], ycoord[hpts], col = lcolor)
  polygon(x = xcoord[hpts], y =  ycoord[hpts], col = adjustcolor(bgcolor, alpha.f = 0.15) , border = 0)

}  

# END OF FUNCTION
###################################################################
