From 6484729b622148d7861ab7879d6ca7d451cade35 Mon Sep 17 00:00:00 2001
From: Yifan Wang <yifan.wang@aei.mpg.de>
Date: Thu, 26 May 2022 13:28:24 +0000
Subject: [PATCH] add utils code to plot

---
 utils/__pycache__/plot.cpython-37.pyc | Bin 0 -> 3925 bytes
 utils/plot.py                         |  75 ++++++++++++++++++++++++++
 2 files changed, 75 insertions(+)
 create mode 100644 utils/__pycache__/plot.cpython-37.pyc
 create mode 100644 utils/plot.py

diff --git a/utils/__pycache__/plot.cpython-37.pyc b/utils/__pycache__/plot.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..c3c80f2961467f6844fb0ad76b0358ea4bab6608
GIT binary patch
literal 3925
zcmZ?b<>g{vU|?8U-=7r3$H4Fy#DQU61_lNP1_p-WWef}qDGVu$ISf&ZV45k48BDW8
zv4Ux~T=pmqMurroD9#k-D6Uk-1>7kt3mH>bB^XlJq#2@kQaMw(L>QuYQ+ZQ)L>Qv@
zQh8JP(wU<8Q`l2DQaDrCQ@B#tQ`l3uQ@B%jQg~DOQ~2iaL<yv_W(h73N)=eh$jFew
zpCZu05G9-<m?G4|5G9f#ks^{J+QJwmnj)4W-og+imckg!pegwh<a)n!5Xm6W?vNgo
z(Ys?(SdaZnRt5$Jh%mzkdnkj!;iUuv1H&!W;?$zdv?@{O{M>@llGGvv-~5!+qCAD5
z#N5QZmm&-d3?VO57#J8{@-Z+lX!2C?C03*s>*XXSrREf8R;5;Pr{(9B=%pp*X695@
zae+Ay6FDJlkQ)A+)b!N66g`j#SQCFmNoI03k}3R^XabyIYtoA{Q>wU2QY%XIN{drV
zQY)&sixbOI(=yZbQVKGw1adO-Qj7I+6N|D_i@*l+rDdj<7NzQ?Wu}Ap$&BE@U|;~z
zObiSR&Y(CI1H~z03R4S14MP@#3qvzw4MRLb4MQ-4CbJ)03rIN|0|P@aNDC;a(iv(P
zVg+g$YZw=Rlf*)%iA;qo!3>&AzZey6F=#U0VlPT9$Vp62)nvNGq-Su8F?l6J5hx}7
z3f3>rFUr=h%uGwn(<@KROV=;T%u7$nFVE9=4>vS5urxH$FQ^2G>8EFw<Rm8Pg4iH=
z-QvpPlGNP9lFZ~{{nC=moMQcgoct2Kg34c9HaVHaCCT}@1$MF^H$og-rBIZbn_re1
zUtExwoT_J&lb@WJQ*5V)Pz-Wta(+%}ZeFn_<1MC?G_V^X1Spbk@gSQZ12UQg6pV~Q
zj735W3=HrPVrF1q0Lg#@HGqMEp@boev4$aqQJf)#DTTS0sg@Bc!U7jzDq*T&Y-X%w
zEMcl)YG$lus$ooFm1M|bPGOT|s9{2|B|y=^?pI~4Z)9W$52g5kAb%Ibc<<DrqIe@C
z!*~N<eZSP=lGK#=;>6s7oYW$H*Rsr%)V$<Wy^@NOA`VcRWxT~)oLT@j$N+~y1_TU>
zV$R8PDgq@UO_p0s1u3^Uic%9(;**QZz#(#rtt>I8G__ch`4&q-Vo636ONx<^A(%9P
zgexc&6oCVm9i%fSFG&FuuHXWQQHW89k%N(gu?UpPZZR7f8770`5)|ej8WhB!pe~L;
z3E~=XFpD#!FsHEeGIcVvGo&%5u%@uJaFj5WFn2Iyu{1L_GlEig2SYqd2SXMsq{QiE
zVq_>`>tJYR>}2R*Na2!X$YKv>(B$^}#i(Cpfg^yhg>*rFaY<@XW`2=gN+K+%i-eIA
z2*~xIECdoPRsk16H4L#bwTvAMSqwFdDGcHaC17?lqZorE12j7^z61v}Q;`6ud}7W?
z&C_Hm;s&u9Gc_4+F{VI*6cn08{GbR1g(SqPDmhT-#uwxlI~rn3`o9>JLo}fVvuEak
zGiwy9qmhx3p(aNWsIV&%1zE&el$e*ES_De_5Ld9I6y+Cy<GKh`dT25gi8C-TL^0>&
z=7C+K$yg-Ez`y{BVlIRu)In~8#Ih8l93vknmX$$41adDZf<bnI5(hXBePLu^NMQu!
zofO6tmK4SmP=vCku%)n<Fr;wI;i+XRVa#G$z+A%wsuoIE7O*a4sAaBUUci>Zxsb7z
zC6l3+wZfo;y@ml4@v$l~Otoya>@|!jT;dER95rmsjN%Nn93`AJ9L<b1tTk*l4253d
z3?*DOtP8joGW0RTFx7I_a@DY8@zijoaHp_>F;5C_FIz2lg<cJJ3ZEo{3q!0@El&;4
z0^SsUkSl6<7VyDY{3+lHt8h^`Lq$yqdkq)Z9R*-_@RkTH5Uk-{$XLr+BCtTHhBrkp
zopB)(lr02eGcnZi)$*4JED)~Ytl^VjSjbo_0Os-6aDsV)H7r>oHEb!u;tVyMH3Bt)
z=}e%?5mX(5s=^eBUQR}a!b{-{H7pB67cyjtrAVeo^)fOt)UafUFOW!)2DvyzCPfxh
z^n<cStxyfi0?87A1yU)(CDJuQ&5R(k7lOmLhGl_FjSxsj4a)-A6gjB8yg0)`CPsz|
znHrV_atlG>Ay6w^BfLO9MFFIC0%PGOu>A@tiXeUs%L2s|B`8ZNML9(Uq`QV?fpUr}
zNHj$xm_bvms+3DXK|vuPC%;6&(a6Y9p)4~$Co@k0Qk@lpC2gSq%qhvtO|?=;Q~=Qm
ziA6<;m0;nb{QMG731p>El98$a;wgZ*3i)XY0hK|SdFcxIr6mQWCB<A2<)Ct<2voq}
zXaE#}nhqgFpiErE4$5motO4*cpMim)ib=o3Ej<kkp4n@%MsXHr<d^5BCgv4KaYL9X
znQ%^OMM+U=ZekQmS!z+qEw=K^l#-0%TWm@BDVeFoRU9RWNmfZYrK#z&W`e*?`yx<Q
zei;wS{n88!41O=cd{9N1UM;fD{b5y$y(arDmdezkqA0f1<oulcqFbzqIRzPsx7Z*p
zmne?p#DbEH%;fAS_T<C@P{ZXGC#Y!yw)+-yerbs&<1NmV%%arflFa-(kYh?p@^j%8
zF~~+xM1h+spz<jl)KuZCWlUjMz)-`ukP%k0R<TuE=^3TfYBJtpOUp0HO)P;HtwpjR
z_k#)|aB&7M$H0YaDyRShna{>hr3`9$#DkJ2Bz0mdVWF+d3$SM93;UNgppudU*3PM7
zcZ{z#GBT_MSEnyQrDm0mN{Mbsd`xa)Nk&m_bx~?st%630hoP;Zp^=%9sildTk)DZ>
zv5~Q+ib8I&rphgjw9NE)a4<3F7UbOGC`t}UEK1BRzQtBrkdj!EdW*Rrr{oqJB=baZ
zC#Iyt7ndf1Tz`u*HLo-`wFs=JxTL5ER2bi4uF5R9#hO`?np+&j1@Rk5aTITHYDs)r
zVlpJwAl%fHbSRe-Bnyd1ro4h%Y`LkRRCJ50AhjqtHLoNyC-oL{Vp4GwM`}?~eo<0l
zQ50u-YDs)WVnt?g6hBBasLdN+0Csv5D@0rfCSDAx#WF$B9mNF_1GkuqZ?PpJswRjb
zV2K=XdpwE*tTHDvH;Nx@703dp6c0$MD8D2TWaKTd+e<P_a#A6tSAtbR^jCuPL+k{L
zYKq)q0TnY*9H8<9RGi#mgE;yYb6RHlEvCeZC>Btv04Jhbj3wZ@@fKsnEyl_yZfI~O
z7L^p=Vk$@mSCdgJDHZWKV6rrd9n_Y~$u9?oQ+!5ddd4j_u)W1mEFkwo5+pcjf~x#m
z{Gg~oPKBTpngps#IhYu+K@KJ!Mh-?UCKg5^Mgb6HVuP#WVB}-uV&r1vV-{lMVggmr
ze2ij@Ld*h;JdAvdQjA=TBFub@0!$z^Jd8!Ypdyg5ND@Se-(t-x%`K=bvIDgk*a{N!
zQWA^7ZHQZ3xrrqpSLS3U6`6u$1rY*z1(hJFTWk=z$O4o*ct8mnoZI6IN}_l(^O7@C
zi{lGYi{cAPZm|>>7L`PCr|0LSq~^tgimfP4xcVqA5Hlk+Gd-iE$OM#$T|tByh%g5c
zpn?e8ummS{un&+z6cMZ(Ho5sJr8%i~pvG4*$Rtpcl!K9niG>jqa$p`4w(w<Q<Y3}q
L;$Y!m=VAu{J+u;H

literal 0
HcmV?d00001

diff --git a/utils/plot.py b/utils/plot.py
new file mode 100644
index 0000000..ecd1cd6
--- /dev/null
+++ b/utils/plot.py
@@ -0,0 +1,75 @@
+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
-- 
GitLab