boxplots.R

#/usr/bin/env Rscript

args <- commandArgs(trailingOnly = TRUE)

pdf("boxplots.pdf", width=8,height=20)

par(mfrow=c(5,2))
for (a in args) {
    x <- read.table(a)
    boxplot(x, log = 'y', xlab = "Number of Taxa", main = a)
}

effiency.rb

#!/usr/bin/env jruby

require 'java'
require 'beast.jar'

java_import 'dr.inference.trace.TraceAnalysis'

re = /^\d+.\d+ (seconds|minutes|hours)\s*$/.freeze

def parse_time(time)
  time.to_f * (time[/minutes/] ? 60 : 1) * (time[/hours/] ? 3600 : 1)
end

def mean(array)
  array.reduce(:+) / array.size
end

def median(array)
  sorted = array.sort
  len = sorted.length
  return (sorted[(len - 1) / 2] + sorted[len / 2]) / 2.0
end

[2, 4, 8].each do |l|
  ['0.33', '0.66', '1.0'].each do |alpha|
    puts "Working on #{l}-#{alpha}"
    f_post_ess_rel = File.open("#{l}-#{alpha}-post_ess_rel.txt", 'w')
    f_node_ess_rel = File.open("#{l}-#{alpha}-node_ess_rel.txt", 'w')
    f_mcmc_post_ess_abs = File.open("#{l}-#{alpha}-mcmc_post_ess_abs.txt", 'w')
    f_mcmc_node_ess_abs = File.open("#{l}-#{alpha}-mcmc_node_ess_abs.txt", 'w')
    f_hmc_post_ess_abs = File.open("#{l}-#{alpha}-hmc_post_ess_abs.txt", 'w')
    f_hmc_node_ess_abs = File.open("#{l}-#{alpha}-hmc_node_ess_abs.txt", 'w')
    f_epsilon = File.open("#{l}-#{alpha}-epsilon.txt", 'w')
    f_hmc_accept = File.open("#{l}-#{alpha}-hmc_accept.txt", 'w')
    f_mcmc_scale_accept = File.open("#{l}-#{alpha}-mcmc_scale_accept.txt", 'w')
    f_mcmc_uniform_accept = File.open("#{l}-#{alpha}-mcmc_uniform_accept.txt", 'w')
    files = [f_post_ess_rel, f_node_ess_rel, f_mcmc_post_ess_abs, f_mcmc_node_ess_abs, f_hmc_post_ess_abs, f_hmc_node_ess_abs, f_epsilon, f_hmc_accept, f_mcmc_scale_accept, f_mcmc_uniform_accept]
    files.each do |f|
      f.write "4 8 16 32 64\n"
    end
    (1..10).each do |rep|
      files.each do |f|
        f.write "#{rep}"
      end
      [4, 8, 16, 32, 64].each do |taxa|
        dir = "#{taxa}-#{alpha}-#{l}"
        mcmc_trace = TraceAnalysis.analyzeLogFile("#{dir}/#{rep}mcmc.log", 500000)
        hmc_trace = TraceAnalysis.analyzeLogFile("#{dir}/#{rep}hmc.log", 5000)
        mcmc_time, hmc_time = File.foreach("#{dir}/hamilton.#{rep}.out").grep(re).map { |t| parse_time(t) }
        mcmc_post_ess = mcmc_trace.getTrace(0).getTraceStatistics.getESS()
        hmc_post_ess = hmc_trace.getTrace(0).getTraceStatistics.getESS()
	mcmc_node_ess = mcmc_trace.getTrace(4+taxa).getTraceStatistics.getESS
	hmc_node_ess = hmc_trace.getTrace(4+taxa).getTraceStatistics.getESS
        epsilon = File.readlines("#{dir}/#{rep}hmc.ops")[3].split[1]
        hmc_accept = File.readlines("#{dir}/#{rep}hmc.ops")[3].split[5]
        mcmc_scale_accept = File.readlines("#{dir}/#{rep}mcmc.ops")[3].split[5]
        mcmc_uniform_accept = File.readlines("#{dir}/#{rep}mcmc.ops")[4].split[4]
        post_ess_rel = (hmc_post_ess / hmc_time) / (mcmc_post_ess / mcmc_time)
        node_ess_rel = (hmc_node_ess / hmc_time) / (mcmc_node_ess / mcmc_time)
        f_post_ess_rel.write " #{post_ess_rel}"
        f_node_ess_rel.write " #{node_ess_rel}"
        f_mcmc_post_ess_abs.write " #{mcmc_post_ess}"
        f_mcmc_node_ess_abs.write " #{mcmc_node_ess}"
        f_hmc_post_ess_abs.write " #{hmc_post_ess}"
        f_hmc_node_ess_abs.write " #{hmc_node_ess}"
        f_epsilon.write " #{epsilon}"
        f_hmc_accept.write " #{hmc_accept}"
        f_mcmc_scale_accept.write " #{mcmc_scale_accept}"
        f_mcmc_uniform_accept.write " #{mcmc_uniform_accept}"
      end
      files.each do |f|
        f.write "\n"
      end
    end
  end
end