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

#this script reads in a score file from a batch Lscore run in PAUP in which Lscores were calculated in two different ways (e.g., under 2 substitution models) for a set of trees (e.g., trees from simulated datasets generated for parametric bootstrapping). 
#It calculates the difference in Lscore for each set of scores and plots a histogram of all the differences
#it also saves a file with the scores and differences
#it takes matrixname (if desired - this is the base matrix name for each simulated dataset), scorefile (the name of the file that holds the Lscore results - e.g., file produced by a PAUP run of the batch file generated by garliLscorebat2) and nreps (the number of matrices) as arguments
#user can also include a value for obs.lnldiff (difference in log likelihood scores from empirical data) and the script will calculate the proportion of the distribution that is greater than or equal to the empirical result


paupScoreHistMaker<- function(matrixname='matrix', scorefile, nreps, obs.lnldiff=(-9999)){

  scores<- list()
  score.file<- strsplit(scan(file=scorefile, what='integer', sep='\n'),'\t')
  counter<- 2
  for (i in 1:nreps){

    scores$matrix.num[i]<- paste(matrixname, i-1, sep='')

    scores$null[i]<-as.numeric(score.file[[counter]][2])
    counter<- counter+2
    scores$alternative[i]<-as.numeric(score.file[[counter]][2])
    counter<- counter+2
    scores$difference[i]<-  (scores$null[i] - scores$alternative[i])
  }

  if(obs.lnldiff > (-9999)){
    cat('The probability of observing a simulated value >= to the observed value is ', sum(scores$difference >= obs.lnldiff)/nreps, '\n')
  }

  hist(scores$difference)
  write.table(scores, file=paste(scorefile,'_scores.txt', sep=''))

}