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