Skip to content
Snippets Groups Projects
Commit 6484729b authored by Yifan Wang's avatar Yifan Wang
Browse files

add utils code to plot

parent 21a7395b
No related branches found
No related tags found
No related merge requests found
File added
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
# PLOTTING OPTIONS
def remove_space(df):
df.columns = [c.replace(' ','') for c in df.columns]
return df
def readlnb(path):
d221 = pd.read_csv(path + "/221/GW150914_PROD1_Kerr_221_0M/Nested_sampler/Evidence.txt",sep=' ')
d220 = pd.read_csv(path + "/220/GW150914_PROD1_Kerr_220_0M/Nested_sampler/Evidence.txt",sep=' ')
d221 = remove_space(d221)
d220 = remove_space(d220)
return d221['lnB'].values[0] - d220['lnB'].values[0]
def read_posA1(path,lmn='221'):
d = pd.read_csv(path+'/'+lmn+'/GW150914_PROD1_Kerr_221_0M/Nested_sampler/posterior.dat',sep=' ')
d.columns = [d.columns[(i+1)%len(d.columns)] for i in range(len(d.columns))]
d.drop(columns='#',inplace=True)
d = remove_space(d)
return d['A2221'].values
def plot_A221violin(time,rootpath):
'''
Plot A221 violin figures
========
time: a time array
rootpath: the root path of PyRing outputs
'''
A221 = {}
for i,t in enumerate(time):
path = rootpath + '/t'+str(i)
A221[i] = read_posA1(path)
x,y = zip(*A221.items())
violinparts = ax.violinplot(y,
showmeans=False,
showmedians=True,
showextrema=False,
vert=True,
widths=0.7)
for pc in violinparts['bodies']:
pc.set_facecolor('tab:blue')
pc.set_edgecolor('tab:blue')
pc.set_alpha(0.4)
#plot error bar:
for i,d in enumerate(y):
x = np.median(d)
#print i,x
dx_l,dx_u = np.percentile(d,[5]),np.percentile(d,[95])
err_low = np.abs(x-dx_l)
err_high = np.abs(dx_u-x)
ax.errorbar(i+1,x,yerr=[err_low,err_high],ecolor='tab:blue',alpha=0.4,
marker=None, capthick=2, capsize=18, linewidth=0.01)
# set x axis
ax.get_xaxis().set_tick_params(direction='out')
ax.xaxis.set_ticks_position('bottom')
labels = ['{:.2f}'.format(i) for i in time]
ax.set_xticks(np.arange(1, len(labels) + 1))
ax.set_xticklabels(labels)
ax.set_xlim(0.25, len(labels) + 0.75)
for tick in ax.get_xticklabels():
tick.set_rotation(60)
ax.set_title('GW150914')
ax.set_ylabel('A_{221}')
ax.set_ylim(0,30)
ax.set_xlabel('$t-t_\mathrm{ref} (tH1=1126259462.42323)$ ms')
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment