#Created by Kurt Galbreath (kurt.galbreath@wwu.edu); copyright 2009

#For a given set of widespread (i.e., present in more than one population) alleles, generate multiple random samples of matching characteristics (i.e., same number of widespread alleles, each found in same number of populations) and produce a histogram showing the frequency distribution for the number of alleles that are distributed across a specific subdivision between populations. Compare the distribution to the empirical result to determine the probability that widespread alleles could have been partitioned as they are by chance alone.

#arguments are as follows:
#alleles = a vector in which each element represents one allele and the value is the number of populations in which it is found (e.g., c(2,2,5,3,3,2,2,4,4,4)  )

#popns = the total number of populations from which to sample (e.g., if you have 4 populations on one side of the subdivision and 3 on the other, then you would input 7 here) 

#partition = the popn number that marks the dividing line (i.e., counting from one, the first popn that is in the second partition; in the example above where subdivision 1 has 4 popns and subdivision 2 has 3, the value to input here would be either 5 or 4 depending on which subdivision was listed first)

#nreps = the number of iterations (samples) to run (default=10000)      



alleleRandomizer<- function(alleles, popns, partition, nreps=10000, obs.amphi){

  prop.amphi <- rep(NA, nreps)
#prop.amphi holds results of iterations, which will be the proportion of alleles that are located in popns on both sides of the geographic partition


  for (i in 1:nreps){
    temp.counter <- rep(1, length(alleles))

    for (j in 1:length(alleles)){
      temp.sample <- sample(popns, size=alleles[j], replace=F)
        print(temp.sample)
        if ((alleles[j]) == sum(temp.sample < partition) |  (sum(temp.sample < partition) == 0)){
          temp.counter[j] <- 0
        }
    }
    cat(temp.counter, '\n')
    prop.amphi[i] <- sum(temp.counter)
  }

  hist(prop.amphi)
  cat('The chance probability of seeing the observed number of widespread alleles (or fewer) spanning the subdivision:',  sum(prop.amphi <= obs.amphi) / nreps, '\n')
}