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

#This function searches the screen logs output by GARLI for the log likelihood of the best tree from multiple replicates
#It will automate the task of extracting log likelihoods from the results of multiple different analyses (e.g., simulated datasets run in GARLI using batch mode)
#If both constrained and unconstrained runs were done, it will compile the scores for both and then calculate the difference between them
#The difference values will be plotted in a histogram, and if an empirical score difference is provided, it will calculate the probability of observing a value in the histogram >= the empirical value
#Takes as arguments: 
#nreps (number of separate GARLI runs to be summarized) 
#matrixname (base name for the data matrices run through GARLI - this is assuming the format 'matrixname, matrix number beginning at 0, .garli.screen.log'
#noconstraint (tells the function that you are including GARLI runs that did not have topology constraints)
#constraint (tells the function that you are including GARLI runs that did have topology constraints) 
#obs.lnldiff (the empirical value for difference in log likelihoods between constrained and unconstrained runs) 

garliBestScoreExtractor<-function(nreps, matrixname, noconstraint=T, constraint=F, obs.lnldiff=(-9999)){
  scores<- list()

  for (i in 1:nreps){
      scores$matrix.num[i] <- paste(matrixname, i-1, sep='')

    if (noconstraint==T){
      outfile.noconst <- scan(file=paste(matrixname, i-1, '.garli.screen.log', sep=''), what='character', sep='\n')

      scores$unconstrainedLL[i] <- as.numeric(strsplit(outfile.noconst[(grep('best', outfile.noconst)[1])],' ')[[1]][4] )
#      cat( scores$unconstrainedLL[i] , '\n', file=paste(matrixname, 'scores.txt'), append=T)
    }
    if (constraint==T){
      outfile.const <- scan(file=paste(matrixname, i-1, '.const.garli.screen.log', sep=''), what='character', sep='\n')
      scores$constrainedLL[i] <- as.numeric(strsplit(outfile.const[(grep('best', outfile.const)[1])],' ')[[1]][4] )
#      cat(scores$constrainedLL[i] , '\n', file=paste(matrixname, 'scores.txt'), append=T)    
      scores$difference[i]<-  (scores$unconstrainedLL[i] - scores$constrainedLL[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')
  }
  if(noconstraint==T & constraint==T){
    hist(scores$difference)
  }
  write.table(scores, file=paste(matrixname,'_besttree_scores.txt', sep=''))


}