# This function takes an interspecific dataset of phenotype and performances, calculates an F-matrix 
# and associated specifics or each species, and then assembles this into an F-array, calcualting an 
# eVCV matrix for each performance measure and comparing pairwaise the eVCV matrices using modified 
# mantel tests that include the diagonal.

# Input is a vector identifying each observation to species, a data frame with phenotypic traits 
# where each row corresponds to an individual belonging to one of the species identified in species, 
# perf is a similar data frame with performance traits, tree is a phylogeny relating all the species, 
# and s is the number of randomizations for the mantel tests, by default 10000.

farray <- function(species,pheno,perf,tree,s=10000) {

# calcualte the number of phenotypic traits, number of performance traits and number of species
	nspp <- length(levels(species))
	npheno <- length(pheno)
	nperf <- length(perf)

# Combine species identifiers with the phenotypic and the performance data
	pheno_sp <- cbind(species,pheno)
	perf_sp <- cbind(species, perf)

# Create a 3D array of dimensions to fit the f-matrices of all species
	far3d <- rep(NA,npheno*nperf*nspp)
	dim(far3d) <- c(npheno,nperf,nspp)

# Calculate an F-matrix for each species as well as all associated stats, also output the f-matrix 
# for each species into the 3D array.
	spp_fmats <- NULL
	for (i in 1:length(levels(species))) {
		spp = levels(species)[i]
		temp <- fmat(pheno_sp[species==spp,2:(npheno+1)],perf_sp[species==spp,2:(nperf+1)])
		spp_fmats[[i]] <- temp
		far3d[,,i] <- temp[[3]]}
	names(spp_fmats) <- levels(species)

# Add labels to the far3d object, the interspecific F-array
	dimnames(far3d) <- list(NULL,NULL,levels(species))
	colnames(far3d) <- colnames(perf_sp[,2:(nperf+1)])
	rownames(far3d) <- colnames(pheno_sp[,2:(npheno+1)])

# Produce an object with the F-array decomposed into a matrix for each performance trait
	perf_coefs <- NULL
	for (i in 1:nperf) {
		perf_coefs[[i]] <- t(far3d[,i,])}
	names(perf_coefs) <- colnames(perf_sp[,2:(nperf+1)])

# Calculate the eVCV matrix for each performance trait
	eVCV <- NULL
	for (i in 1:nperf) {
		eVCV[[i]] <- ratematrix(tree,perf_coefs[[i]])}
	names(eVCV) <- colnames(perf_sp[,2:(nperf+1)])

# Calculate a modified Mantel test with s permutations for each pairwise combo of eVCV matrices
	mantel_out <- NULL
	for (i in 1:(length(eVCV)-1)) {
		for (j in (i+1):length(eVCV)) {
			mantel_temp <- mantel(eVCV[[i]],eVCV[[j]],s,type="d_s")
			mantel_temp <- c(names(perf_coefs)[i],names(perf_coefs)[j],mantel_temp)
			mantel_out <- rbind(mantel_out,mantel_temp)}}
	rownames(mantel_out)<-NULL
	if (nrow(mantel_out)==1)
	{mantel_out <- mantel_out} else
	{mantel_out <- cbind(as.factor(mantel_out[,1]),as.factor(mantel_out[,2]),as.data.frame(mantel_out[,3:7]))
	colnames(mantel_out)<-c("v1","v2","R","p","ncol","nrow","nreps")}

# Compile results
	results <- list(spp_fmats,perf_coefs,eVCV,mantel_out)
	names(results) <- c("spp_fmats","perf_coefs","eVCV","mantel")

	results
	}