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

notebooks

parent bda75d0e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
notebook/4-inj-after-making-GW190521.ipynb 0 → 100644
+ 353
− 0
View file @ d6442691
%% Cell type:markdown id:5e728de6 tags:
# Plotting Options
%% Cell type:code id:furnished-potato tags:
``` python
import numpy as np
import h5py,glob,os
from pycbc import conversions
from scipy import stats
from tqdm import tqdm
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
# PLOTTING OPTIONS
fig_width_pt = 3*246.0 # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27 # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt # width in inches
fig_height = fig_width*golden_mean+0.5 # height in inches
fig_size = [fig_width,fig_height]
# PLOTTING OPTIONS
fig_width_pt = 3*246.0 # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27 # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt # width in inches
fig_height = fig_width*golden_mean # height in inches
fig_size = [fig_width,fig_height]
params = { 'axes.labelsize': 24,
'font.family': 'serif',
'font.serif': 'Computer Modern Raman',
'font.size': 24,
'legend.fontsize': 20,
'xtick.labelsize': 24,
'ytick.labelsize': 24,
'axes.grid' : True,
'text.usetex': True,
'savefig.dpi' : 100,
'lines.markersize' : 14,
'figure.figsize': fig_size}
mpl.rcParams.update(params)
from pesummary.core.plots.bounded_1d_kde import bounded_1d_kde
```
%% Cell type:code id:e3b939a3 tags:
``` python
binmin = 0
binmax = 1000
binnum = int(int(binmax-binmin) + 1)
print(binnum)
bins = np.linspace(binmin,binmax,binnum)
```
%% Output
1001
%% Cell type:markdown id:01419cad tags:
# GW190521 real data run
%% Cell type:code id:5cbd2e49 tags:
``` python
f = h5py.File('/work/yifan.wang/testingGR-o3b/t4-mpi-massive/GW190521-mpv1000/extract.hdf','r')
```
%% Cell type:code id:a3f7f975 tags:
``` python
d = bounded_1d_kde(f['samples']['parity_mpvinverse'][()])
```
%% Cell type:code id:da49b929 tags:
``` python
r = d(bins[:-1])
r = r/np.sum(r)/(bins[2]-bins[1])
```
%% Cell type:code id:8f9bcae6 tags:
``` python
combine = np.zeros(len(bins[:-1]))
event = []
for i in range(100):
txt = np.loadtxt('./INJDATA/GW190521/'+str(i)+'.txt')
event.append(txt[:,1])
combine += txt[:,1]
```
%% Cell type:code id:cdc06bac tags:
``` python
plt.plot(bins[:-1],combine/100,label='injection average')
plt.plot(bins[:-1],r,label='GW190521')
plt.plot(100,100,color='gray',zorder=-1,alpha=1,label='injections')
for i in range(100):
plt.plot(bins[:-1],event[i],color='gray',zorder=-1,alpha=0.1)
plt.legend()
plt.ylim(0,0.004)
plt.xlim(0,1000)
plt.title('GW190521')
plt.ylabel('Probability Density')
plt.xlabel('$M_\mathrm{PV}^{-1}$ [GeV$^{-1}$]')
plt.savefig('GW190521-inj.png',bbox_inches='tight')
```
%% Output
%% Cell type:markdown id:97916621 tags:
# Bayes factor
%% Cell type:code id:753ebecb tags:
``` python
0.001/r[0]
```
%% Output
10.997900661198125
%% Cell type:code id:584ed5f1 tags:
``` python
np.log(0.001/r[0])
```
%% Output
2.3977044056023975
%% Cell type:code id:13cc1368 tags:
``` python
np.median(f['samples']['parity_mpvinverse'][:])/np.std(f['samples']['parity_mpvinverse'][:])
```
%% Output
2.11694748406522
%% Cell type:code id:a7a900e7 tags:
``` python
r0 = [event[i][0] for i in range(100)]
```
%% Cell type:code id:80fec17a tags:
``` python
np.sort(0.001/np.array(r0))
```
%% Output
array([ 0.10760022, 0.12231412, 0.1491996 , 0.17254167, 0.19083087,
0.19924071, 0.20795018, 0.20839901, 0.20857025, 0.21132103,
0.21489632, 0.21980263, 0.22463276, 0.23087591, 0.2325678 ,
0.23780638, 0.24164496, 0.24436879, 0.25928547, 0.27801257,
0.28535016, 0.29337103, 0.30878492, 0.31968284, 0.31973142,
0.32378121, 0.32481493, 0.34276225, 0.35061967, 0.35403893,
0.36487152, 0.36761182, 0.37567379, 0.41706441, 0.43634365,
0.44148127, 0.44879488, 0.48232095, 0.48260404, 0.48874957,
0.49077087, 0.49716043, 0.4997546 , 0.50248692, 0.50269499,
0.5054598 , 0.50973671, 0.51876029, 0.52082413, 0.5305977 ,
0.53718161, 0.53774539, 0.54849621, 0.54964462, 0.57738073,
0.58006832, 0.59143698, 0.59951194, 0.60704105, 0.61505825,
0.61837905, 0.66972387, 0.69290913, 0.70589659, 0.71326893,
0.73171696, 0.75443428, 0.76732157, 0.77484612, 0.77513832,
0.78708929, 0.79322574, 0.79375416, 0.82795317, 0.83535967,
0.8373832 , 0.86032494, 0.87910256, 0.89052968, 0.93234874,
0.93450264, 0.97997868, 0.98925307, 0.99747858, 1.03073057,
1.06253297, 1.17100283, 1.18475713, 1.19467888, 1.22474887,
1.27758112, 1.36423179, 1.37571616, 1.39121304, 1.42592869,
1.51822558, 1.79046334, 1.82822198, 2.34268803, 25.30624655])
%% Cell type:code id:18494e04 tags:
``` python
np.argsort(0.001/np.array(r0))
```
%% Output
array([46, 84, 56, 75, 13, 64, 90, 31, 48, 78, 58, 86, 73, 76, 82, 72, 74,
61, 33, 43, 59, 99, 37, 69, 18, 5, 38, 23, 10, 65, 70, 44, 77, 51,
36, 28, 15, 60, 57, 85, 2, 39, 16, 62, 97, 79, 94, 49, 91, 12, 81,
8, 21, 54, 11, 71, 29, 24, 45, 52, 30, 93, 42, 67, 34, 68, 53, 1,
22, 98, 96, 20, 32, 47, 3, 55, 89, 92, 17, 83, 4, 35, 25, 95, 27,
7, 26, 40, 80, 6, 87, 0, 19, 88, 9, 66, 14, 41, 63, 50])
%% Cell type:code id:f9d5ef68 tags:
``` python
```
notebook/4-inj-after-making-GW191109.ipynb 0 → 100644
+ 381
− 0
View file @ d6442691
%% Cell type:markdown id:5e728de6 tags:
# Plotting Options
%% Cell type:code id:furnished-potato tags:
``` python
import numpy as np
import h5py,glob,os
from pycbc import conversions
from scipy import stats
from tqdm import tqdm
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
# PLOTTING OPTIONS
fig_width_pt = 3*246.0 # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27 # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt # width in inches
fig_height = fig_width*golden_mean+0.5 # height in inches
fig_size = [fig_width,fig_height]
# PLOTTING OPTIONS
fig_width_pt = 3*246.0 # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27 # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt # width in inches
fig_height = fig_width*golden_mean # height in inches
fig_size = [fig_width,fig_height]
params = { 'axes.labelsize': 24,
'font.family': 'serif',
'font.serif': 'Computer Modern Raman',
'font.size': 24,
'legend.fontsize': 20,
'xtick.labelsize': 24,
'ytick.labelsize': 24,
'axes.grid' : True,
'text.usetex': True,
'savefig.dpi' : 100,
'lines.markersize' : 14,
'figure.figsize': fig_size}
mpl.rcParams.update(params)
from pesummary.core.plots.bounded_1d_kde import bounded_1d_kde
```
%% Cell type:code id:e3b939a3 tags:
``` python
binmin = 0
binmax = 1000
binnum = int(int(binmax-binmin) + 1)
print(binnum)
bins = np.linspace(binmin,binmax,binnum)
```
%% Output
1001
%% Cell type:markdown id:01419cad tags:
# GW191109 real data run
%% Cell type:code id:5cbd2e49 tags:
``` python
#f = h5py.File('/work/yifan.wang/testingGR-o3b/t2-o3bevents/o3btgr_output/posterior_files/H1L1V1-EXTRACT_POSTERIOR_GW191109_010717-1126259200-400.hdf','r')
```
%% Cell type:code id:3a61f1c2 tags:
``` python
f = h5py.File('/work/yifan.wang/testingGR-o3b/t4-mpi-massive/GW191109-xphm-mpv1000/extract.hdf','r')
```
%% Cell type:code id:a3f7f975 tags:
``` python
d = bounded_1d_kde(f['samples']['parity_mpvinverse'][()])
```
%% Cell type:code id:da49b929 tags:
``` python
r = d(bins[:-1])
r = r/np.sum(r)/(bins[2]-bins[1])
```
%% Cell type:code id:a95404f5 tags:
``` python
combine = np.zeros(len(bins[:-1]))
event = []
for i in range(100):
try:
txt = np.loadtxt('./INJDATA/GW191109/'+str(i)+'.txt')
event.append(txt[:,1])
combine += txt[:,1]
except OSError:
print(i, 'is missing')
pass
```
%% Output
93 is missing
95 is missing
%% Cell type:code id:cdc06bac tags:
``` python
plt.plot(bins[:-1],combine/len(event),label='injection average')
plt.plot(bins[:-1],r,label='GW191109')
plt.plot(100,100,color='gray',zorder=-1,alpha=1,label='injections')
for i in range(len(event)):
plt.plot(bins[:-1],event[i],color='gray',zorder=-1,alpha=0.1)
plt.legend()
plt.ylim(0,0.008)
plt.xlim(0,1000)
plt.title('GW191109')
plt.ylabel('Probability Density')
plt.xlabel('$M_\mathrm{PV}^{-1}$ [GeV$^{-1}]$')
plt.savefig('GW191109-inj.png',bbox_inches='tight')
```
%% Output
%% Cell type:markdown id:72dcda91 tags:
# Bayes factor
%% Cell type:code id:e579a76c tags:
``` python
0.001/r[0]
```
%% Output
22.91526523853397
%% Cell type:code id:1d74a4f7 tags:
``` python
np.log(0.001/r[0])
```
%% Output
3.1318032927810275
%% Cell type:code id:4a0a3e72 tags:
``` python
np.median(f['samples']['parity_mpvinverse'][:])/np.std(f['samples']['parity_mpvinverse'][:])
```
%% Output
2.868612528466472
%% Cell type:code id:d438dbaa tags:
``` python
r0 = [event[i][0] for i in range(len(event))]
```
%% Cell type:code id:5faba4ec tags:
``` python
np.sort(0.001/np.array(r0))
```
%% Output
array([8.83795087e-02, 1.02076437e-01, 1.10617837e-01, 1.30700778e-01,
1.31428371e-01, 1.33774799e-01, 1.36146538e-01, 1.37394732e-01,
1.38391179e-01, 1.38965426e-01, 1.43365367e-01, 1.45384636e-01,
1.48743947e-01, 1.50077646e-01, 1.52467160e-01, 1.57079810e-01,
1.59191364e-01, 1.60117123e-01, 1.60406757e-01, 1.63495208e-01,
1.66470644e-01, 1.67045033e-01, 1.68924678e-01, 1.69367639e-01,
1.75360979e-01, 1.80799326e-01, 1.80954840e-01, 1.83642663e-01,
1.84350308e-01, 1.87689601e-01, 1.89200947e-01, 1.89867948e-01,
1.94188976e-01, 1.95842505e-01, 1.99637842e-01, 2.00523467e-01,
2.03587028e-01, 2.06710778e-01, 2.09104127e-01, 2.10247681e-01,
2.10286886e-01, 2.10675623e-01, 2.12720803e-01, 2.20862808e-01,
2.28109817e-01, 2.28441933e-01, 2.31354921e-01, 2.34007312e-01,
2.35457025e-01, 2.38924365e-01, 2.42519147e-01, 2.48013964e-01,
2.60268112e-01, 2.73629212e-01, 2.73833444e-01, 2.75720919e-01,
2.76792374e-01, 2.78969480e-01, 2.83135542e-01, 2.84783716e-01,
2.92218514e-01, 3.06362287e-01, 3.15467735e-01, 3.32061435e-01,
3.32647636e-01, 3.33390023e-01, 3.36321519e-01, 3.50452398e-01,
3.61283899e-01, 3.75044368e-01, 3.77143838e-01, 3.90488979e-01,
4.13463999e-01, 4.25031735e-01, 4.26539346e-01, 4.37707141e-01,
4.38064153e-01, 4.87927565e-01, 4.95877889e-01, 7.60147971e-01,
8.67196798e-01, 9.64104866e-01, 1.21248718e+00, 1.35817430e+00,
1.62178447e+00, 2.00076594e+00, 2.12834552e+00, 2.21873177e+00,
2.99720682e+00, 3.62726461e+00, 4.36987900e+00, 5.33267802e+00,
6.23817245e+00, 8.02823099e+00, 3.89530315e+01, 2.95794359e+02,
1.15620162e+03, 6.69388930e+89])
%% Cell type:code id:3fbdd43f tags:
``` python
np.argsort(0.001/np.array(r0))
```
%% Output
array([89, 68, 42, 69, 79, 78, 61, 57, 21, 94, 53, 95, 32, 46, 52, 13, 7,
93, 54, 88, 29, 11, 16, 2, 23, 38, 85, 28, 87, 62, 76, 24, 92, 34,
41, 22, 40, 45, 27, 77, 14, 97, 19, 59, 70, 5, 86, 51, 37, 49, 90,
0, 75, 55, 17, 48, 71, 1, 44, 9, 67, 25, 47, 12, 64, 80, 82, 58,
91, 63, 31, 84, 20, 35, 8, 10, 72, 66, 73, 96, 26, 50, 65, 30, 74,
15, 18, 6, 33, 60, 3, 56, 83, 81, 39, 36, 43, 4])
%% Cell type:code id:f9d5ef68 tags:
``` python
```
notebook/4-inj.ipynb 0 → 100644
+ 384
− 0
View file @ d6442691
%% Cell type:markdown id:5e728de6 tags:
# Plotting Options
%% Cell type:code id:furnished-potato tags:
``` python
import numpy as np
import h5py,glob,os
from pycbc import conversions
from scipy import stats
from tqdm import tqdm
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
# PLOTTING OPTIONS
fig_width_pt = 3*246.0 # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27 # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt # width in inches
fig_height = fig_width*golden_mean+0.5 # height in inches
fig_size = [fig_width,fig_height]
params = { 'axes.labelsize': 22,
'font.family': 'serif',
'font.serif': 'Computer Modern Raman',
'font.size': 24,
'legend.fontsize': 20,
'xtick.labelsize': 24,
'ytick.labelsize': 24,
'axes.grid' : False,
'text.usetex': True,
'savefig.dpi' : 100,
'lines.markersize' : 14,
'figure.figsize': fig_size}
mpl.rcParams.update(params)
from pesummary.core.plots.bounded_1d_kde import bounded_1d_kde
```
%% Cell type:markdown id:3bc713d8 tags:
# Collecting all the posterior files
%% Cell type:code id:2fe054f7 tags:
``` python
posdir = glob.glob('/work/yifan.wang/testingGR-o3b/injstudy/GW190521/injections-GW190521_output/posterior_files/*.hdf')
len(posdir)
path = '/work/yifan.wang/testingGR-o3b/injstudy/GW190521/injections-GW190521_output/posterior_files/'
```
%% Cell type:code id:smart-program tags:
``` python
mpvinv = []
for i in range(100):
try:
f = h5py.File(path+'H1L1V1-EXTRACT_POSTERIOR_INJ'+str(i)+'-1126259200-400.hdf','r')
mpvinv.append(f['samples']['parity_mpvinverse'][()])
f.close()
except FileNotFoundError:
print(i,' is missing.')
pass
```
%% Cell type:code id:e3b939a3 tags:
``` python
binmin = 0
binmax = 1000
binnum = int(int(binmax-binmin) + 1)
print(binnum)
bins = np.linspace(binmin,binmax,binnum)
```
%% Output
1001
%% Cell type:code id:21b585cd tags:
``` python
combine = np.zeros(len(bins[:-1]))
k = 0
for i in tqdm(range(len(mpvinv))):
#i is the sorted index for Mpvinv
density = bounded_1d_kde(mpvinv[i], xlow=0., xhigh=1000, method="Reflection")
distribution = density(bins[:-1])
norm = np.sum(distribution) * (bins[1]-bins[0])
np.savetxt('./INJDATA/GW190521/'+str(i)+'.txt',np.transpose([bins[:-1],distribution/norm]))
#combine += distribution/norm
```
%% Output
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [1:16:44<00:00, 46.04s/it]
%% Cell type:markdown id:01419cad tags:
# GW190521 real data run
%% Cell type:code id:5cbd2e49 tags:
``` python
f = h5py.File('/work/yifan.wang/testingGR-o3b/t4-mpi-massive/GW190521-mpv1000/extract.hdf','r')
```
%% Cell type:code id:a3f7f975 tags:
``` python
d = bounded_1d_kde(f['samples']['parity_mpvinverse'][()])
```
%% Cell type:code id:da49b929 tags:
``` python
r = d(bins[:-1])
r = r/np.sum(r)/(bins[2]-bins[1])
```
%% Cell type:code id:cdc06bac tags:
``` python
plt.plot(bins[:-1],combine/len(mpvinvplus),label='injection')
plt.plot(bins[:-1],r,label='GW190521')
plt.legend()
plt.savefig('GW190521-inj.png',bbox_inches='tight')
```
%% Output
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
/local/user/yifan.wang/ipykernel_3238571/3463940405.py in <module>
----> 1 plt.plot(bins[:-1],combine/len(mpvinvplus),label='injection')
2 plt.plot(bins[:-1],r,label='GW190521')
3 plt.legend()
4 plt.savefig('GW190521-inj.png',bbox_inches='tight')
NameError: name 'mpvinvplus' is not defined
%% Cell type:markdown id:14be1bbe tags:
# GW191109
%% Cell type:code id:e80548bc tags:
``` python
path2 = '/work/yifan.wang/testingGR-o3b/injstudy/GW191109-rerun/injections-GW191109_output/posterior_files/'
for i in tqdm(range(100)):
try:
f = h5py.File(path2+'H1L1V1-EXTRACT_POSTERIOR_INJ'+str(i)+'-1126259200-400.hdf','r')
density = bounded_1d_kde(f['samples']['parity_mpvinverse'][()], xlow=0., xhigh=1000, method="Reflection")
distribution = density(bins[:-1])
norm = np.sum(distribution) * (bins[1]-bins[0])
np.savetxt('./INJDATA/GW191109/'+str(i)+'.txt',np.transpose([bins[:-1],distribution/norm]))
except FileNotFoundError:
print(i,' is missing.')
pass
```
%% Output
93%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 93/100 [08:39<00:36, 5.23s/it]
93 is missing.
95%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 95/100 [08:43<00:18, 3.80s/it]
95 is missing.
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [09:09<00:00, 5.49s/it]
%% Cell type:code id:793f0491 tags:
``` python
f = h5py.File('/work/yifan.wang/testingGR-o3b/t2-o3bevents/o3btgr_output/posterior_files/H1L1V1-EXTRACT_POSTERIOR_GW191109_010717-1126259200-400.hdf','r')
```
%% Cell type:code id:18d4b73e tags:
``` python
d = bounded_1d_kde(f['samples']['parity_mpvinverse'][()])
```
%% Cell type:code id:fb1c3d34 tags:
``` python
combine = np.zeros(len(bins[:-1]))
k = 0
for i in tqdm(range(100)):
try:
f = np.loadtxt('./INJDATA/GW191109/'+str(i)+'.txt')
combine += f[:,1]
except OSError:
pass
```
%% Output
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [00:00<00:00, 138.42it/s]
%% Cell type:code id:bcd8caf2 tags:
``` python
plt.plot(bins[:-1],combine/98)
plt.plot(bins[:-1],d(bins[:-1]),label='GW191109')
```
%% Output
[<matplotlib.lines.Line2D at 0x7f3494f854e0>]
%% Cell type:code id:f9d5ef68 tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint