1#!/usr/bin/env python 2import pylab as plt 3from scipy.integrate import simps, trapz 4import numpy as np 5 6class job: 7 def __init__(self, jobname): 8 f = open ( jobname, "r" ) 9 self.ac = [float(x) for x in f.readline().strip('[').strip(']\n').split(',') ] 10 self.ac_mean = [float(x) for x in f.readline().strip('[').strip(']\n').split(',') ] 11 self.ac_std = [float(x) for x in f.readline().strip('[').strip(']\n').split(',') ] 12 line = f.readline() 13 self.time_mean = float(line.strip('[').strip(']\n').split(',')[0]) 14 self.time_std = float(line.strip('[').strip(']\n').split(',')[1]) 15 16 def printing(self): 17 print self.ac 18 print self.ac_mean 19 print self.ac_std 20 print self.time_mean 21 print self.time_std 22 23 def diff_from_ms(self, ms_ac, delta): 24 y = np.array([ np.abs(ms_ac[i] - self.ac[i]) for i in range(len(ms_ac))] ) 25 self.dev = np.abs(simps(y, x = delta)) 26 #print self.dev 27 28delta = range( 0, int(2e4+1), 1000 ) 29 30prefix = "Constantpopsize" 31suffix = "_Processed" 32ms_case = [""] 33fastsimcoal_case = [""] 34scrm_case = [ "window" + `x` for x in [0, 500, 1000, 3000, 5000, 7000, 10000, 30000, 50000, 70000, 100000] ] 35macs_case = [ "retain" + `x` for x in [0, 500, 1000, 3000, 5000, 7000, 10000, 30000, 50000, 70000, 100000] ] 36macs_case.insert(0,"") 37program = ["ms", "fastsimcoal", "scrm", "macs"] 38case = [ms_case, fastsimcoal_case, scrm_case, macs_case] 39joblist = [] 40 41ms = job("Constantpopsizems_Processed") 42#ms.printing() 43 44fig1 = plt.figure(figsize=(9, 6), dpi=80) 45ax1 = fig1.add_subplot(111) 46 47fig2 = plt.figure(figsize=(12, 8), dpi=80) 48ax2 = fig2.add_subplot(111) 49 50fig3 = plt.figure(figsize=(9, 6), dpi=80) 51ax3 = fig3.add_subplot(111) 52 53#ax1.plot ( delta, ms.ac , linewidth=3.0, color = "black") 54#ax1.errorbar ( delta, ms.ac, yerr = [ x/(1000**0.5) *1.96 for x in ms.ac_std ] ) 55 56linestyles = ['-', '-.', '--', ':'] 57colors = [ "blue", "red", "green", "purple", "black", "yellow", "cyan", "magenta", "orange"] 58markers = ["v", "o", "*", ">", "<", "s", "^", "+" , "D", "H", "d","x"] 59legendlist1 = [] 60l1 = [] 61 62legendlist2 = [] 63l2 = [] 64 65legendlist3 = [] 66l3 = [] 67 68for i, program_i in enumerate ( program ): 69 color_j = 0 70 program_dev = [] 71 program_time = [] 72 program_time_err = [] 73 for j, case_j in enumerate ( case[i] ): 74 current_job2 = job( prefix + program_i + case_j + suffix ) 75 current_job2.diff_from_ms (ms.ac, delta) 76 program_dev.append ( current_job2.dev) 77 program_time.append (current_job2.time_mean) 78 program_time_err.append( current_job2.time_std ) 79 #current_dot = ax2.plot ( current_job2.dev, np.log(current_job2.time_mean), markers[j], color = colors[i]) 80 current_dot = ax2.plot ( current_job2.dev, current_job2.time_mean, markers[j], color = colors[i]) 81 l2.append(current_dot) 82 legendlist2.append(program_i + case_j) 83 if j % 3 == 0: 84 current_job = job( prefix + program_i + case_j + suffix ) 85 legendlist1.append( program_i + case_j) 86 current_line = ax1.plot ( delta, current_job.ac, linestyles[i], color = colors[color_j] ) 87# ax1.errorbar ( delta, current_job.ac, yerr = [ x/(1000**0.5) *1.96 for x in current_job.ac_std ], 88# fmt='.', color = colors[color_j] ) 89 l1.append(current_line) 90 relative_ac = [ np.abs(ms.ac[ac_i] - current_job.ac[ac_i]) for ac_i in range(len(ms.ac))] 91 current_line3 = ax3.plot ( delta, relative_ac, linestyles[i], color = colors[color_j] ) 92 # ax3.errorbar ( delta, relative_ac, yerr = [ x/(1000**0.5) *1.96 for x in current_job.ac_std ], 93 # fmt='.', color = colors[color_j] ) 94 l3.append(current_line3) 95 color_j += 1 96 97 #ax2.errorbar ( program_dev, np.log(program_time), yerr = program_time_err, color = colors[i]) 98 ax2.plot ( program_dev, program_time, color = colors[i]) 99 #ax2.errorbar ( program_dev, program_time, yerr = program_time_err, color = colors[i]) 100 101ms_line = ax1.plot ( delta, ms.ac, linewidth=2.0, color = "black") 102#ax1.errorbar ( delta, ms.ac, yerr = [ x/(1000**0.5) *1.96 for x in current_job.ac_std ], 103 #fmt='.', color = colors[color_j] ) 104l1[0] = ms_line 105 106relative_ac = [ float(0) for ac_i in range(len(ms.ac))] 107ms_line3 = ax3.plot ( delta, relative_ac, linewidth=2.0, color = "black") 108#ax3.errorbar ( delta, relative_ac, yerr = [ x/(1000**0.5) *1.96 for x in ms.ac_std ], 109 #fmt='.', color = "black" ) 110l3[0] = ms_line3 111 112ax1.legend ([ x[0] for x in l1], legendlist1, loc = 1) 113ax1.axis([0,20000, 0, 1]) 114ax1.set_xlabel(r'Distance between two sites $\delta$') 115ax1.set_ylabel(r'Autocorrelation $\rho$') 116fig1.savefig("TMRCArhoLD.pdf") 117 118ax3.legend ([ x[0] for x in l1], legendlist1, loc = 1) 119ax3.axis([0,30000, -.01, 0.06]) 120ax3.set_xlabel(r'Distance between two sites $\delta$') 121ax3.set_ylabel(r'Error in Autocorrelation $\rho$') 122fig3.savefig("RelativeTMRCArhoLD.pdf") 123 124 125ax2.set_xlim ([-10, 1300]) 126ax2.set_ylim ([0.1, 14]) 127#ax2.axis([-10, 1300, -2.5, np.log(ms.time_mean)*1.1] ) 128#yticks = ax2.get_yticks() 129#print yticks 130#ylabels = ["%.5g" % (np.exp(float(y))) for y in yticks] 131#ax2.set_yticklabels(ylabels) 132ax2.set_yscale('log') 133ax2.legend ([ x[0] for x in l2], legendlist2, loc = 1, numpoints=1) 134ax2.set_ylabel("Time (sec)") 135ax2.set_xlabel("Deviation") 136fig2.savefig("time_vs_dev.pdf") 137