################# Begin Shipley et al.'s function ###############
############### <http://pages.usherbrooke.ca/jshipley/recherche/> ##
Maxent<- function(constraint.means, constraint.matrix, initial.probs = 
 NA)
 {
 # MAxent function in R. Written by Bill Shipley
 # Bill.Shipley@USherbrooke.ca
 # Improved Iterative Scaling Algorithm from Della Pietra, S.
 # Della Pietra, V., Lafferty, J.  (1997). Inducing features of random fields.
 # IEEE Transactions Pattern Analysis And Machine Intelligence
 # 19:1-13
 #.............................
 # constraint.means is a vector of the community-aggregated trait values
 # constraint.matrix is a matrix with each trait defining a row and each
 # species defining a column.  Thus, two traits and 6 species would produce
 # a constraint.matrix with two rows (the two traits) and 6 columns (the
 # six species)
 if(!is.vector(constraint.means)) return(
 "constraint.means must be a vector")
 if(is.vector(constraint.matrix)) {
 constraint.matrix <- matrix(constraint.matrix,
 nrow = 1, ncol = length(
 constraint.matrix))
 }
 n.constraints <- length(constraint.means)
 dim.matrix <- dim(constraint.matrix)
 if(dim.matrix[2] == 1 & dim.matrix[1] > 1) {
 constraint.matrix <- t(constraint.matrix)
 dim.matrix <- dim(constraint.matrix)
 }
 if(n.constraints != dim.matrix[1])
 return("Number of constraint means not equal to number in constraint.matrix"
 )
 n.species <- dim.matrix[2]
 if(is.na(initial.probs)) {
 prob <- rep(1/n.species, n.species)
 }
 if(!is.na(initial.probs)) {
 prob <- initial.probs
 }
 # C.values is a vector containing the sum of each trait value over the species 
 C.values <- apply(constraint.matrix, 1, sum)
 test1 <- test2 <- 1
 while(test1 > 1e-008 & test2 > 1e-008) {
 denom <- constraint.matrix %*% prob
 gamma.values <- log(constraint.means/denom)/
 C.values
 if(sum(is.na(gamma.values)) > 0)
 return("missing values in gamma.values"
 )
 unstandardized <- exp(t(gamma.values) %*% 
 constraint.matrix) * prob
 new.prob <- as.vector(unstandardized/sum(
 unstandardized))
 if(sum(is.na(new.prob)) > 0)
 return("missing values in new.prob")
 test1 <- 100 * sum((prob - new.prob)^2)
 test2 <- sum(abs(gamma.values))
 prob <- new.prob
 }
 list(probabilities = prob, final.moments = denom,
 constraint.means = constraint.means, 
 constraint.matrix = constraint.matrix, entropy
  = -1 * sum(prob * log(prob)))
 }
################# End Shipley et al.'s function ###############

###############################################################
#Analysis of random communities using the method of Shipley et al. Science 314, 812 (2006).
#Code by S. Roxburgh, December 2006
#stephen.roxburgh@ensisjv.com

#Set the number of species, and the number of traits, and program control parameters
NSpecies<-50              # The number of species
NTraits<-8                # The number of traits
Niters<-500               # The number of random communities to generate
BasedOnObservedData<-TRUE  # Whether the 8 random traits are based on observations, or are selected at random from N(25,5)
PlotScatter<-TRUE         # Whether to plot the observed vs predicted scatterplot
ExponentialAbundance<-TRUE  #True if pseudo-observed abundances are exponentially distributed, false if normal N(20,2)

#Set up the graphics
if (PlotScatter) {
  par(mfcol=c(1,2))
  }

#Initialise the trait matrices & vector for saving correlation coefficients
Constraint.matrix<-matrix(0,NTraits,NSpecies)
Community.traits<-matrix(0,NTraits,NSpecies)
Rvector<-matrix(Niters,1)

for (loop in 1:Niters) {  # main program loop
 if (BasedOnObservedData) {
   #Observed trait means and SD'ds, in the order they appear in Fig 1
   #For the first trait (lifespan) the mean value is the probability of perreniality.
   #All traits are selected from the Normal distribution, except traits 5-7 (AG, stem and leaf biomass) which are Exponential. 
   #Summary statistics for species traits in intermediate successional field sites, taken from D. Vile et al. Ecology 87, 504 (2006)
   ObsMeans<-c(0.416666667,2.417792258,212.1666667,21.49166667,1.89,2.77,1.5675,62.51666667)
   ObsStdev<-c(0,0.811281016,40.21721326,4.710039246,2.241687676,3.733275432,3.269668109,30.34042708)

   #Fill trait matrix with random values
   #First do lifespan
   for (i in 1:NSpecies) {
     val<-runif(1,0,1)
     if (val <= ObsMeans[1]) {
       Constraint.matrix[1,i]<-1   
     }
     else {
       Constraint.matrix[1,i]<-0
     }
   }
   #If all elements of lifespan are zero, make the first one 1, to allow estimation to proceed
   if (sum(Constraint.matrix[1,]) == 0) {
     Constraint.matrix[1,1]<-1
     }

   #Then do the remainder of the traits
   Constraint.matrix[2,]<-rnorm(NSpecies,ObsMeans[2],ObsStdev[2]) # log Seed number
   Constraint.matrix[3,]<-rnorm(NSpecies,ObsMeans[3],ObsStdev[3]) # Seed maturation date
   Constraint.matrix[4,]<-rnorm(NSpecies,ObsMeans[4],ObsStdev[4]) # SLA
   Constraint.matrix[5,]<-rexp(NSpecies)*ObsMeans[5] # Vm1
   Constraint.matrix[6,]<-rexp(NSpecies)*ObsMeans[6] # Stm
   Constraint.matrix[7,]<-rexp(NSpecies)*ObsMeans[7] # Lm
   Constraint.matrix[8,]<-rnorm(NSpecies,ObsMeans[8],ObsStdev[8]) # Height
  }
  else #if not basing analysis on observed data
  {
    Constraint.matrix<-rnorm(NSpecies*NTraits,20,5)
    dim(Constraint.matrix)<-c(NTraits,NSpecies)
  }

 #Create a vector of pseudo-observed abundances, from exponential or normal distribution
 if (ExponentialAbundance) {
   Abund<-rexp(NSpecies)
   }
   else {
   Abund<-rnorm(NSpecies,20,1)
   }
 RelAbund<-Abund/sum(Abund)

 #Create the community level traits, based on pseudo-observed relative abundances
 for (i in 1:NSpecies) {
  for (j in 1:NTraits) {
   Community.traits[j,i]<- Constraint.matrix[j,i] * RelAbund[i]
 }}

 #Calculate community-aggregated traits - i.e summed over all species for each trait
 Constraint.means<-apply(Community.traits,1,sum)

 #Run the entropy algorithm
 out<-Maxent(Constraint.means, Constraint.matrix); out
 #save the correlation coefficient
 Rvector[loop]<-cor(RelAbund,out$probabilities)
 #fit a linear regression line for each optimisation
 outlinear<-lm(out$probabilities~RelAbund) 
 #Only update the scatter plot if required
 if (PlotScatter) {
    #If this is the first iteration, then plot, otherwise append
   if (loop == 1) { 
    plot(RelAbund,out$probabilities,xlim=c(0,0.5),ylim=c(0,0.5),xlab = "Pseudo-observed relative abundance", ylab = "Predicted relative abundance",pch=20)
   }
   else 
   {
    points(RelAbund,out$probabilities,main = cat(Rvector[loop],loop),xlim=c(0,0.5),ylim=c(0,0.5),pch=20)
   }
   #Only plot the first 25 regression lines
   if (loop <= 25) { 
     abline(outlinear)
   }
  }
}

#Print the results to textfile - [NOTE, saved as r2]
cat(NSpecies, NTraits,"\n",file = "Results.txt", sep = " ", Fill = NULL, labels = NULL, append = 0)
for (i in 1:Niters) {
  cat(Rvector[i]*Rvector[i],"\n",file = "Results.txt", sep = " ", Fill = NULL, labels = NULL, append = 1)
  }

#Plot frequency histogram [NOTE, plotted as r], and reset graphics
hist(Rvector)
par(mfrow=c(1,1))