###########################################################################
#  "slicedice.R" A program used to do Sex Randomization and Count Cells of Morphospace 
#  published in:
#
#  Butler, Sawyer, and Losos (2007) Sexual Dimorphism and Adaptive Radiation in Anolis lizards, 
#  Nature. in press.
#                  written by Marguerite Butler, 06/06/06
#                  
#  Input file: a tab-delimited file with multivariate data expressed as canonical variates or 
# PCA scores. These data should be labeled "can1", "can2", "can3". This analysis only considers 
# the first three axes.
# Input file should also contain "sex" as "M" or "F", and " 
###########################################################################

s0 <- read.delim("input.txt", header=T)
attach(s0)


#---------------------------------------------------------
#
# 3D: slices = 1.0, can1, can2 and can3
#
# This code bins the data into units defined by slice width (1.0), 
# assigns a unique index to each 3D cube of morphospace
# then counts up how much morphospace is occupied by all, M, and F only 
##########################################################################
unit=1.0
min1 <- floor(min(can1))
min2 <- floor(min(can2))
min3 <- floor(min(can3))
max1 <- ceiling(max(can1))
max2 <- ceiling(max(can2))
max3 <- ceiling(max(can3))

r1 <- seq(min1, max1, by=unit)
r2 <- seq(min2, max2, by=unit)
r3 <- seq(min3, max3, by=unit)

cut(s0$can1, br=r1, labels=FALSE ) -> bin1
cut(s0$can2, br=r2, labels=FALSE ) -> bin2
cut(s0$can3, br=r3, labels=FALSE ) -> bin3

nrows <- max(bin1)
ncols <- max(bin2) 

index <- bin1 + (bin2-1)*nrows + (bin3-1)*(nrows*ncols)
length(unique(index)) -> ncubesall			   # 289  Observed score, ALL
length(unique(index[sex=="F"])) -> ncubesFobs  # 182                  F only
length(unique(index[sex=="M"])) -> ncubesMobs  # 154                  M only
											   # 47 overlaps, 107 unique M, 135 unique F
obsscore <- ncubesall - ncubesMobs             # 135 Amt. added by adding females
#---------------------------------------------------------

#  Do randomizations over SEX, compute Male morphospace occupancy
#  number of successes = cases where male morphospace is less than or equal to ncubesMobs
#  this is because total morphospace is the same, therefore if this is true, random "females" 
#  occupy more unique morphospace than the original data. 

ns = 0
for (i in 1:20000) 
{
  by(s0$sex, s0$species, function(x) sample(x)) -> rsex
  unlist(rsex, use.names=FALSE) -> rsex

  length(unique(index[rsex=="M"])) -> Mperm  
  ns <- ns + (Mperm <= ncubesMobs) 
}
print(paste("ns = ", ns)
print(paste("P-value = ", pval <- ns/20000))
# The overall P-value is ns/#randomizations. 
#> ns  
#[1] 236
#> pval
#[1] 0.0118
#> (117+119)/10000
#[1] 0.0118
#> 



count_cubes <-function(eco, ss) {
  cubes <- length(unique(index[s0$ecomorph==eco & s0$sex==ss])) 
  others <- length(unique(index[!(s0$ecomorph==eco & s0$sex==ss)]))
  unique_cubes <- ncubesall-others
  nindiv <- sum(s0$ecomorph==eco & s0$sex==ss) 

  results<- data.frame(eco, ss, nindiv, cubes, unique_cubes)
  return(results)
}

eco <- unique(ecomorph)
ss <- unique(sex)
rr <- count_cubes(eco[1], ss[1])

for (i in 1: length(ss))
{
	for (j in i:length(eco))
	{
		
		rr <- rbind(rr, count_cubes(eco[j], ss[i]))
	}	
}

print(rr[-1,])  # This prints out the following output

#   eco ss nindiv cubes unique_cubes
#2   TG  F     87    57           32
#3   CG  F     15    13            9
#4   TC  F     87    47           37
#5   GB  F     45    38           26
#6   TW  F     39    33           27
#7   CG  M     16    16           12
#8   TC  M     63    38           26
#9   GB  M     43    36           23
#10  TW  M     19    16           10