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

add comparison for duration, fix labeling issues

parent 9fba7cd6
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
1-compare-ringdown-pyring-snr.ipynb
+ 1181
− 1127
View file @ 3c29c440
This diff is collapsed.
utils/wheel.py
+ 73
− 35
View file @ 3c29c440
...@@ -73,7 +73,7 @@ def local_loglikelihood(model, ...@@ -73,7 +73,7 @@ def local_loglikelihood(model,
class wheel(): class wheel():
def __init__(self,model,fit=None): def __init__(self,model,fit=None,acf_from_psd=None):
'''A wheel to reproduce the PyRing loglikelihood and SNR computation '''A wheel to reproduce the PyRing loglikelihood and SNR computation
Parameter: Parameter:
...@@ -82,8 +82,6 @@ class wheel(): ...@@ -82,8 +82,6 @@ class wheel():
The KerrModel for pyring The KerrModel for pyring
''' '''
self.model = model
self.data = {} self.data = {}
self.cinv = {} self.cinv = {}
self.acf = {} self.acf = {}
...@@ -102,32 +100,64 @@ class wheel(): ...@@ -102,32 +100,64 @@ class wheel():
'phi2220': -0.9737754582321143, 'phi2220': -0.9737754582321143,
'phi2221': 1.652119726595113} 'phi2221': 1.652119726595113}
_ = model.log_likelihood(pyring_par) _ = model.log_likelihood(pyring_par)
waveform_model = model.get_waveform(pyring_par)
t_start = model.fixed_params['t'] t_start = model.fixed_params['t']
duration_n = model.duration_n duration_n = model.duration_n
for d in model.detectors.keys(): for d in model.detectors.keys():
dt = model.time_delay['{0}_'.format(self.model.ref_det)+d] dt = model.time_delay['{0}_'.format(model.ref_det)+d]
tref = LIGOTimeGPS(t_start+dt+self.model.tevent) tref = LIGOTimeGPS(t_start+dt+model.tevent)
# crop data # crop data
time_array_raw = self.model.detectors[d].time - (self.model.tevent+dt) time_array_raw = model.detectors[d].time - (model.tevent+dt)
time_array = time_array_raw[time_array_raw >= t_start][:duration_n] time_array = time_array_raw[time_array_raw >= t_start][:duration_n]
self.data[d] = self.model.detectors[d].time_series[time_array_raw >= t_start][:duration_n] self.data[d] = model.detectors[d].time_series[time_array_raw >= t_start][:duration_n]
self.cinv[d] = self.model.detectors[d].inverse_covariance self.cinv[d] = model.detectors[d].inverse_covariance
# read in the acf # read in the acf
self.acf[d] = np.loadtxt( if acf_from_psd:
'./GW150914_PROD1_Kerr_221_0M/Noise/ACF_TD_cropped_'+str(d)+'_1126257414_4096_4.0_2048_0.2.txt') self.acf[d] = np.loadtxt('/work/yifan.wang/ringdown/GW150914/pyring/compare-pyring-ringdown-pycbc/GW150914_PROD1_Kerr_221_0M/Noise/ACF_TD_cropped_'+str(d)+'_1126257414_4096_4.0_2048_0.2.txt')
psd = np.loadtxt( psd = np.loadtxt('/work/yifan.wang/ringdown/GW150914/pyring/compare-pyring-ringdown-pycbc/GW150914_PROD1_Kerr_221_0M/Noise/PSD_'+str(d)+'_1126257414_4096_4.0_2048.txt')
'./GW150914_PROD1_Kerr_221_0M/Noise/PSD_'+str(d)+'_1126257414_4096_4.0_2048.txt') acf_psd = 0.5*np.fft.irfft(psd[:,1]) * self.model.srate
acf_psd = 0.5*np.fft.irfft(psd[:,1]) * self.model.srate c = toeplitz(acf_psd[:model.duration_n])
c = toeplitz(acf_psd[:model.duration_n]) self.cinv_acf_from_psd[d] = inv(c)
self.cinv_acf_from_psd[d] = inv(c)
if fit is not None: if fit is not None:
acf = fit.acfs[d].values[:model.duration_n] acf = fit.acfs[d].values[:model.duration_n]
c = toeplitz(acf_psd[:model.duration_n]) c = toeplitz(acf[:model.duration_n])
self.cinv_acf_from_ringdown[d] = inv(c) self.cinv_acf_from_ringdown[d] = inv(c)
self.model = model
def get_data(self,detector_time,detector_rawdata):
'''Get the date corredsponding to waveform after correcting the labeling issue
Parameters:
-----------
detector_time: dict
By default, 4s duration
detector_rawdata: dict
By default, 4s data
Return:
-----------
cropped_time: dict
cropped time duration corredspoinding to waveform
cropeed_data: dict
cropped data corresponding to waveform
'''
cropped_time = {}
cropped_data = {}
t_start = self.model.fixed_params['t']
duration_n = self.model.duration_n
for d in self.model.detectors.keys():
dt = self.model.time_delay['{0}_'.format(self.model.ref_det)+d]
tref = LIGOTimeGPS(t_start+dt+self.model.tevent)
# crop data
lcrop = detector_time[d] >= self.model.tevent + dt + t_start
cropped_time[d] = detector_time[d][lcrop][:duration_n]
cropped_data[d] = detector_rawdata[d][lcrop][:duration_n]
return cropped_time,cropped_data
def get_hstrain(self,pyring_par,detector_time): def get_hstrain(self,pyring_par,detector_time):
'''Compute the waveform given parameters '''Compute the waveform given parameters
...@@ -140,6 +170,7 @@ class wheel(): ...@@ -140,6 +170,7 @@ class wheel():
''' '''
h = {} h = {}
waveform_time = {} waveform_time = {}
#To activate model.time_delay we need to run loglikelihood first #To activate model.time_delay we need to run loglikelihood first
_ = self.model.log_likelihood(pyring_par) _ = self.model.log_likelihood(pyring_par)
waveform_model = self.model.get_waveform(pyring_par) waveform_model = self.model.get_waveform(pyring_par)
...@@ -153,13 +184,14 @@ class wheel(): ...@@ -153,13 +184,14 @@ class wheel():
for d in self.model.detectors.keys(): for d in self.model.detectors.keys():
dt = self.model.time_delay['{0}_'.format(self.model.ref_det)+d] dt = self.model.time_delay['{0}_'.format(self.model.ref_det)+d]
tref = LIGOTimeGPS(t_start+dt+self.model.tevent) tref = LIGOTimeGPS(t_start+dt+self.model.tevent)
# crop data # crop data
#time_array_raw = self.model.detectors[d].time - (self.model.tevent+dt) #time_array_raw = self.model.detectors[d].time - (self.model.tevent+dt)
time_array_raw = detector_time[d] - (self.model.tevent+dt) time_array_raw = detector_time[d] - (self.model.tevent+dt)
time_array = time_array_raw[time_array_raw >= t_start][:duration_n] time_array = time_array_raw[time_array_raw >= t_start][:duration_n]
wf_model = waveform_model.waveform(time_array) wf_model = waveform_model.waveform(time_array)
hs, hvx, hvy, hp, hc = wf_model[0], wf_model[1], wf_model[2], wf_model[3], wf_model[4] hs, hvx, hvy, hp, hc = wf_model[0], wf_model[1], wf_model[2], wf_model[3], wf_model[4]
h[d] = project(hs, hvx, hvy, hp, hc, self.model.detectors[d].lal_detector, ra, dec, psi, tref) h[d] = project(hs, hvx, hvy, hp, hc, self.model.detectors[d].lal_detector, ra, dec, psi, tref)
waveform_time[d] = time_array + self.model.tevent + dt waveform_time[d] = time_array + self.model.tevent + dt
return waveform_time,h return waveform_time,h
...@@ -171,7 +203,7 @@ class wheel(): ...@@ -171,7 +203,7 @@ class wheel():
loglikelihood += loglikelihood_core(residuals, self.cinv[d],self.model.detectors[d].log_normalisation) loglikelihood += loglikelihood_core(residuals, self.cinv[d],self.model.detectors[d].log_normalisation)
return loglikelihood return loglikelihood
def optsnr(self,pyring_par,network=True,acf_from_psd=False,acf_from_ringdown=False): def optsnr(self,pyring_par,rawtime,network=True,acf_from_psd=False,acf_from_ringdown=False):
''' '''
Optimal SNR Optimal SNR
...@@ -182,13 +214,13 @@ class wheel(): ...@@ -182,13 +214,13 @@ class wheel():
If true, use the L from Ringdown, where C=LL^T (cholesky decomposition) If true, use the L from Ringdown, where C=LL^T (cholesky decomposition)
----------- -----------
''' '''
h = self.get_hstrain(pyring_par) _, h = self.get_hstrain(pyring_par,rawtime)
snr = {} snr = {}
for d in self.model.detectors.keys(): for d in self.model.detectors.keys():
if acf_from_psd == True: if acf_from_psd == True:
cinv = self.cinv_acf_from_psd[d] cinv = self.cinv_acf_from_psd[d]
elif acf_from_ringdown: elif acf_from_ringdown == True:
cinv = self.cinv_acf_from_ringdown[d] cinv = self.cinv_acf_from_ringdown[d]
else: else:
cinv = self.cinv[d] cinv = self.cinv[d]
...@@ -204,24 +236,30 @@ class wheel(): ...@@ -204,24 +236,30 @@ class wheel():
snr[d] = np.sqrt(snr[d]) snr[d] = np.sqrt(snr[d])
return snr return snr
def mfsnr(self,pyring_par,network=True,acf_from_psd=False,acf_from_ringdown=False): def mfsnr(self,pyring_par,rawtime,rawdata,network=True,acf_from_psd=False,acf_from_ringdown=False):
''' '''Matched-filter SNR
Matched-filter SNR
Parameters: Parameters:
----------- -----------
pyring_par: dict
Source parameters
rawtime: dict
The time duration for one chunk of signal
rawdata: dict
The data for one chunk of signal
''' '''
h = self.get_hstrain(pyring_par) _, h = self.get_hstrain(pyring_par,rawtime)
_, data = self.get_data(rawtime,rawdata)
snr = {} snr = {}
for d in self.model.detectors.keys(): for d in self.model.detectors.keys():
if acf_from_psd == True: if acf_from_psd:
cinv = self.cinv_acf_from_psd[d] cinv = self.cinv_acf_from_psd[d]
elif acf_from_ringdown: elif acf_from_ringdown:
cinv = self.cinv_acf_from_ringdown[d] cinv = self.cinv_acf_from_ringdown[d]
else: else:
cinv = self.cinv[d] cinv = self.cinv[d]
snr[d] = np.dot(self.data[d],np.dot(cinv,h[d]))**2 \ snr[d] = np.dot(data[d],np.dot(cinv,h[d]))**2 \
/np.dot(h[d],np.dot(cinv,h[d])) /np.dot(h[d],np.dot(cinv,h[d]))
if network: if network:
...@@ -289,7 +327,7 @@ class wheel(): ...@@ -289,7 +327,7 @@ class wheel():
/np.dot(self.hstrain[d],solve_toeplitz(self.acf[d],self.hstrain[d])) /np.dot(self.hstrain[d],solve_toeplitz(self.acf[d],self.hstrain[d]))
return np.sqrt(snr) return np.sqrt(snr)
def compute_pyring_snr(model,M,chi,A,phi, def compute_multiple_snr(model,pr_time,pr_data,M,chi,A,phi,
fit=None,network=True,acf_from_psd=False,acf_from_ringdown=False): fit=None,network=True,acf_from_psd=False,acf_from_ringdown=False):
''' '''
Loop compute the optimal SNR and matched-filter SNR given a PyRing Model Loop compute the optimal SNR and matched-filter SNR given a PyRing Model
...@@ -299,7 +337,7 @@ def compute_pyring_snr(model,M,chi,A,phi, ...@@ -299,7 +337,7 @@ def compute_pyring_snr(model,M,chi,A,phi,
------- -------
Optimal SNR, Matched-Filtering SNR: np.array() Optimal SNR, Matched-Filtering SNR: np.array()
''' '''
#Initialization #Initialization
if network: if network:
optsnr = [] optsnr = []
...@@ -321,12 +359,12 @@ def compute_pyring_snr(model,M,chi,A,phi, ...@@ -321,12 +359,12 @@ def compute_pyring_snr(model,M,chi,A,phi,
'phi2221': -phi[1][i].values} 'phi2221': -phi[1][i].values}
if network: if network:
optsnr.append(result.optsnr(pyring_par,network,acf_from_psd,acf_from_ringdown)) optsnr.append(result.optsnr(pyring_par,pr_time,network,acf_from_psd,acf_from_ringdown))
mfsnr.append(result.mfsnr(pyring_par,network,acf_from_psd,acf_from_ringdown)) mfsnr.append(result.mfsnr(pyring_par,pr_time,pr_data,network,acf_from_psd,acf_from_ringdown))
else: else:
for d in model.detectors.keys(): for d in model.detectors.keys():
optsnr[d].append(result.optsnr(pyring_par,network,acf_from_psd,acf_from_ringdown)[d]) optsnr[d].append(result.optsnr(pyring_par,pr_data,network,acf_from_psd,acf_from_ringdown)[d])
mfsnr[d].append(result.mfsnr(pyring_par,network,acf_from_psd,acf_from_ringdown)[d]) mfsnr[d].append(result.mfsnr(pyring_par,pr_time,pr_data,network,acf_from_psd,acf_from_ringdown)[d])
return optsnr,mfsnr return optsnr,mfsnr
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint