from pylab import *
import arviz as az
import h5py
import pandas as pd
import seaborn as sns
import ringdown
sns.set_context('notebook')
sns.set_palette('colorblind')
def read_strain(file, dname):
with h5py.File(file, 'r') as f:
t0 = f['meta/GPSstart'][()]
T = f['meta/Duration'][()]
h = f['strain/Strain'][:]
dt = T/len(h)
raw_strain = ringdown.Data(h, index=t0 + dt*arange(len(h)), ifo=dname)
return raw_strain
h_raw_strain = read_strain('H-H1_GWOSC_16KHZ_R1-1126259447-32.hdf5', 'H1')
l_raw_strain = read_strain('L-L1_GWOSC_16KHZ_R1-1126259447-32.hdf5', 'L1')
h_raw_strain.get_psd(nperseg=int(h_raw_strain.fsamp)).iloc[1:].plot(label=h_raw_strain.ifo)
l_raw_strain.get_psd(nperseg=int(l_raw_strain.fsamp)).iloc[1:].plot(label=l_raw_strain.ifo)
xlabel(r'$f / \mathrm{Hz}$');
ylabel(r'$P\left( f \right) / \mathrm{Hz}^{-1}$');
xscale('log');
yscale('log');
legend(loc='best');
fit = ringdown.Fit(model='mchi_aligned', modes=[(1, -2, 2, 2, 0), (1, -2, 2, 2, 1)])
fit.result = az.from_netcdf('./srate16384-0p2s-32sdata.nc')
fit.add_data(h_raw_strain)
fit.add_data(l_raw_strain)
T=0.2
srate=16384
fit.set_target(1126259462.4083147, ra=1.95, dec=-1.27, psi=0.82, duration=T)
fit.condition_data(ds=int(round(h_raw_strain.fsamp/srate)), flow=20)
fit.compute_acfs()
az.plot_trace(fit.result, var_names=['A', 'M', 'chi', 'f', 'tau'], compact=True);
L = [a.iloc[:fit.n_analyze].cholesky for a in fit.acfs.values()]
hs = fit.result.posterior.h_det.stack(samples=('chain', 'draw'))
data_pd = fit.analysis_data
data_dict = list(data_pd.values())
def rd_snr(L,hs,data=None,optimal=True, network=True) -> np.ndarray:
'''
L: cholesky factor
hs: h strain
data: raw data
'''
whs = linalg.solve(L, hs)
opt_ifo_snrs = linalg.norm(whs, axis=1)
if optimal:
snrs = opt_ifo_snrs
else:
# whiten it with the Cholesky factors, so shape will remain (ifo, time)
wdata = linalg.solve(L, data)
# take inner product between whitened template and data, and normalize
snrs = einsum('ijk,ij->ik', whs, wdata)/opt_ifo_snrs
if network:
# take norm across detectors
return linalg.norm(snrs, axis=0)
else:
return snrs
optsnr = rd_snr(L,hs,data=None,optimal=True,network=True)
mfsnr = rd_snr(L,hs,data_dict,optimal=False,network=True)
plot(optsnr,label='opt snr')
plot(mfsnr,label='mf snr')
legend()
<matplotlib.legend.Legend at 0x7fa7a280c160>
from pyRing import pyRing
from pyRing import likelihood
from pyRing import noise
from tqdm import tqdm
WARNING:cpnest.utils:Setting verbosity to 0
/work/yifan.wang/virtualenv/ringdown/lib/python3.7/site-packages/pyRing/utils.py:16: UserWarning: surfinBH is not automatically installed due to possible conflicts. If you wish to use its functionalities, it needs to be installed separately.
warnings.warn("surfinBH is not automatically installed due to possible conflicts. If you wish to use its functionalities, it needs to be installed separately.")
# stacking the posterior h_det and parameters
hs = fit.result.posterior.h_det.stack(samples=('chain', 'draw'))
A = fit.result.posterior.A.stack(samples=('chain','draw'))
phi = fit.result.posterior.phi.stack(samples=('chain','draw'))
M = fit.result.posterior.M.stack(samples=('chain','draw'))
chi = fit.result.posterior.chi.stack(samples=('chain','draw'))
import sys,importlib
sys.path.append("/work/yifan.wang/ringdown/GW150914/pyring/gitlab-ringdown-systematics/utils/")
from pyring_par import *
import data
import wheel
import plot
test_par = input_par
test_par.update({'f-min-bp':20,
'f-max-bp':8191,
'sampling-rate': 16384,
'analysis-duration': 0.2,
'analysis-duration-n': int(0.2*16384),
'likelihood-method':'direct-inversion',
'split-inner-products':0,
'debug':0})
model = pyRing.KerrModel(modes=input_par['kerr-modes'],**test_par)
Trigtime in H1: 1126259462.423000 Reading data... Using GWPY to download data. Fetched 1 URLs from www.gw-openscience.org for [1126257414 .. 1126261510)) Reading data... [Done] Loaded channel GWOSC starting at 1126257414.0 length 4096.0s. Bandpassing the raw strain between [20, 8191] Hz. Computing the one-sided PSD with the Welch method and the standard ACF for comparison. No ACF was passed. Estimating ACF. Plancherel theorem E(f)/E(t) (expected value: 1) = 1.0661345786416638 Reading data... Using GWPY to download data. Fetched 1 URLs from www.gw-openscience.org for [1126257414 .. 1126261510)) Reading data... [Done] Loaded channel GWOSC starting at 1126257414.0 length 4096.0s. Bandpassing the raw strain between [20, 8191] Hz. Computing the one-sided PSD with the Welch method and the standard ACF for comparison. No ACF was passed. Estimating ACF. Plancherel theorem E(f)/E(t) (expected value: 1) = 1.0639064116263919 Using equatorial sky coordinates. ----Prior section - General---- I'll be running with the following prior bounds: Noise evidence: 297883.5826673859 Running the Kerr model with modes: [(2, 2, 2, 0), (2, 2, 2, 1)] ----Prior section - Kerr model---- I'll be running with the following prior bounds: Mf : [20.0, 200.0] af : [0.0, 0.99] ----Prior section - Kerr Model - Testing-GR parameters---- I'll be running with the following prior bounds: No TGR parameter was considered. ----Prior section - Kerr Model - Amplitudes and phases---- I'll be running with the following prior bounds: A2220 : [0.0, 50.0] A2221 : [0.0, 50.0] phi2220 : [0.0, 6.283185307179586] phi2221 : [0.0, 6.283185307179586] ----Prior section - Kerr Model - Fixed parameters---- I'll be running with the following fixed parameters: ra : 1.95 dec : -1.27 psi : 0.82 t : 0.0 cosiota : -1.0 phi : 0.0
result = wheel.wheel(model)
i = np.argmax(mfsnr)
prefactor = np.sqrt(16*np.pi/5)
pyring_par = {'Mf': M[i].values,
'af': chi[i].values,
'A2220': A[0][i].values/1e-21*prefactor,
'A2221': A[1][i].values/1e-21*prefactor,
'phi2220': -phi[0][i].values,
'phi2221': -phi[1][i].values}
pr_data = {}
pr_time = {}
for d in model.detectors.keys():
pr_data[d],pr_time[d] = data.local_load_data(d,**test_par)
Reading data... Using GWPY to download data. Fetched 1 URLs from www.gw-openscience.org for [1126257414 .. 1126261510)) Reading data... [Done] Loaded channel GWOSC starting at 1126257414.0 length 4096.0s. Bandpassing the raw strain between [20, 8191] Hz. Reading data... Using GWPY to download data. Fetched 1 URLs from www.gw-openscience.org for [1126257414 .. 1126261510)) Reading data... [Done] Loaded channel GWOSC starting at 1126257414.0 length 4096.0s. Bandpassing the raw strain between [20, 8191] Hz.
pr_optsnr_net,pr_mfsnr_net = wheel.compute_multiple_snr(
model,pr_time,pr_data,M,chi,A,phi,network=True,acf_from_psd=False)
100%|████████████████████████████████████████████████████████████████████████████████████| 4000/4000 [06:15<00:00, 10.64it/s]
wheel.plotsnr(optsnr,pr_optsnr_net,snr='Optimal SNR')
wheel.plotsnr(mfsnr,pr_mfsnr_net,snr='Matched-Filter SNR')
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
sns.set_context('notebook')
sns.set_palette('colorblind')
df_pyring = plot.read_pos('/work/yifan.wang/more_ringdown/runs/pyring/0p1-16384/')
mIMR = 67.4
chiIMR = 0.67
df = pd.DataFrame({
r'$M / M_\odot$': fit.result.posterior.M.values.flatten(),
r'$\chi$': fit.result.posterior.chi.values.flatten()
})
pg = sns.PairGrid(df, diag_sharey=False)
pg.map_diag(sns.kdeplot);
pg.map_upper(sns.scatterplot);
pg.map_lower(ringdown.kdeplot_2d_clevels, levels=[0.9, 0.5, 0.1])
pg.axes[0,0].axvline(mIMR, color='k')
pg.axes[0,1].axvline(chiIMR, color='k')
pg.axes[0,1].axhline(mIMR, color='k')
pg.axes[1,0].axhline(chiIMR, color='k')
pg.axes[1,0].axvline(mIMR, color='k')
pg.axes[1,1].axvline(chiIMR, color='k')
<matplotlib.lines.Line2D at 0x7fa77fd7d6a0>
sns.displot(fit.result.posterior.A[:,:,1].values.flatten()/1e-21*prefactor, kde=True, stat='density')
xlabel(r'$A_1$')
Text(0.5, 7.9599999999999795, '$A_1$')
sns.displot(df_pyring['A2221'],kde=True,stat='density')
xlabel(r'$A_1$')
Text(0.5, 7.9599999999999795, '$A_1$')
rd_a221 = fit.result.posterior.A[:,:,1].values.flatten()
median(rd_a221)/std(rd_a221)
2.3021402
median(df_pyring['A2221'])/std(df_pyring['A2221'])
2.1633854394424157